text:
# In this script, we will design a PID controller by optimization. The system model has uncertain parameters, and we pay a price not only for poor performance of the closed-loop system, but also for a high variance in the performance. In addition to this, we place a constraint on the 90:th percentile of the maximum of the sensitivity function. This way, we will get a doubly robust controller as a result :) To avoid excessive amplification of measurement noise, we penalize noise amplification above a certain frequency.
# We start by defining the system and some initial controller parameters
using MonteCarloMeasurements, Optim, ControlSystems, Plots
using MonteCarloMeasurements: ∓
unsafe_comparisons(true)
p = 1 + 0.1*Particles(100)
ζ = 0.3 + 0.05*Particles(100)
ω = 1 + 0.05*Particles(100)
const P = tf([p*ω], [1, 2ζ*ω, ω^2]) |> ss
const w = exp10.(LinRange(-2,3,100))
params = log.([1,0.1,0.1])
const Msc = 1.2 # Constraint on Ms

# We now define the cost function, which includes the constraint on the maximum sensitivity function

function systems(params::AbstractVector{T}) where T
    kp,ki,kd = exp.(params)
    C = convert(StateSpace{Continuous, T}, pid(kp=kp,ki=ki,kd=kd)*tf(1, [0.05, 1])^2, balance=false)
    G = feedback(P*C) # Closed-loop system
    S = 1/(1 + P*C)   # Sensitivity function
    CS = C*S          # Noise amplification
    local Gd
    try
        Gd = c2d(G,0.1) # Discretize the system. This might fail for some parameters, so we catch these cases and return a high value
    catch
        return T(10000)
    end
    y,t,_  = step(Gd,15) .|> vec # This is the time-domain simulation
    C, G, S, CS, y, t
end

function cost(params::AbstractVector{T}) where T
    C, G, S, CS, y, t = systems(params)
    Ms = maximum(bode(S, w)[1]) # Maximum of the sensitivity function
    q  = pquantile(Ms, 0.9)
    performance = mean(abs, 1 .- y)   # This is our performance measure
    robustness = (q > Msc ? 10000(q-Msc) : zero(T)) # This is our robustness constraint
    variance = pstd(performance)     # This is the price we pay for high variance in the performance
    noise = pmean(sum(bode(CS, w[end-30:end])[1]))
    100pmean(performance) + robustness  + 10variance + 0.002noise
end

# We are now ready to test the cost function.
@time cost(params)
#
res = Optim.optimize(cost, params, NelderMead(), Optim.Options(iterations=1000, show_trace=true, show_every=20));
println("Final cost: ", res.minimum)
# We can now perform the same computations as above to visualize the found controller
fig = plot(layout=2)
for params = (params, res.minimizer)
    C, G, S, CS, y, t = systems(params)
    mag   = bode(S, w)[1][:]
    plot!(t,y[:], title="Time response", subplot=1, legend=false)
    plot!(w, mag, title="Sensitivity function", xscale=:log10, yscale=:log10, subplot=2, legend=false)
end
hline!([Msc], l=(:black, :dash), subplot=2)
display(fig)

# Other things that could potentially be relevant is adding a probabilistic constraint on the time-domain output, such as the probability of having the step response go above 1.5 must be < 0.05 etc.


text:
#!/usr/bin/env julia

"""
Methods for finding and identifying prime numbers.
"""
module Primes

include("Sequences.jl")

# Currently computed prime number terms (in sorted order)
_prime_sequence = Int128[2]

"""
Precomputes and stores at least the first `n` prime numbers.
"""
function _compute_primes(n::Integer)
    global _prime_sequence

    # check if first n primes have already been computed
    if n < length(_prime_sequence)
        return
    end

    if n ≤ 25
        _compute_primes_up_to(97)
    else
        # set estimate and increment values according to PNT
        approx_gap = _estimate_prime_gap(n)
        estimate = trunc(BigInt, n * approx_gap)
        increment = trunc(BigInt, approx_gap)

        # get primes up to estimate, then step forward until n are found
        _compute_primes_up_to(estimate)
        while length(_prime_sequence) < n
            estimate += increment
            _compute_primes_up_to(estimate)
        end
    end
end

"""
Precomputes and stores the prime numbers up to `n`.
"""
function _compute_primes_up_to(n::Integer)
    global _prime_sequence

    # check if primes up to n have already been computed
    prime_max = last(_prime_sequence)
    if prime_max ≥ n
        return
    end

    # prepare sieve for numbers prime_max + 1 to n
    sieve_size = n - prime_max
    sieve = repeat([true], sieve_size)

    # sift out composite numbers using previously computed primes
    prime_count = length(_prime_sequence)
    for p in _prime_sequence
        p_start = Sequences.next_multiple(p, prime_max + 1)
        for j in p_start:p:n
            sieve[j - prime_max] = false
        end
    end

    # sift out composite numbers using previously computed primes
    for i in 1:sieve_size
        if sieve[i]
            p = i + prime_max
            push!(_prime_sequence, p)
            for j in (p^2 - prime_max):p:sieve_size
                sieve[j] = false
            end
        end
    end
end

"""
Returns an estimate for the average gap between the first `n` primes.

This formula is based on a result derived from the Prime Number Theorem:
https://en.wikipedia.org/wiki/Prime_number_theorem
"""
function _estimate_prime_gap(n::Integer)
    log_n = log(n)
    loglog_n = log(log_n)
    log_n2 = log(n^2)
    loglog_n2 = log(log_n^2)
    return (
        log_n + loglog_n - 1 + (loglog_n - 2)/log_n
        - (loglog_n2 - 6loglog_n + 11)/(2log_n2) + exp(1)/log_n2
    )
end

"""
Resets the currently cached list of prime numbers.
"""
function _reset_prime_cache()
    global _prime_sequence
    _prime_sequence = Int128[2]
    return
end

"""
    prime(n::Integer)

Finds the `n`th prime number.

## Preconditions
- `n > 0`

## Examples
```jldoctest
julia> prime(1)
2

julia> prime(6)
13
```
"""
function prime(n::Integer)
    _compute_primes(n)
    global _prime_sequence
    return _prime_sequence[n]
end

"""
    prime_factorization(n::Integer)

Computes the prime factorization of the natural number `n`.

Result is a 1D array of base-exponent tuple pairs, containing each prime factor
and its power in the prime factorization of n.

## Preconditions
- `n ≥ 2`

## Examples
```jldoctest
julia> prime_factorization(UInt8(2))
1-element Array{Tuple{UInt8,UInt8},1}:
 (0x02, 0x01)

julia> prime_factorization(Int64(12))
2-element Array{Tuple{Int64,Int64},1}:
 (2, 2)
 (3, 1)

julia> prime_factorization(BigInt(600))
3-element Array{Tuple{BigInt,BigInt},1}:
 (2, 3)
 (3, 1)
 (5, 2)
```
"""
function prime_factorization(n::Integer)
    n_type = typeof(n)

    # populate array with (factor, power) pairs
    factor = 2 * oneunit(n_type)
    factorization = Array{Tuple{n_type, n_type}, 1}()
    while factor ≤ floor(n_type, sqrt(n))
        # compute power of factor in factorization
        power = zero(n_type)
        (d, r) = divrem(n, factor)
        while r == 0
            n = d
            power += 1
            (d, r) = divrem(n, factor)
        end

        # add (factor, power) pair to factorization
        if power > zero(n_type)
            push!(factorization, (factor, power))
        end

        factor += 1
    end

    # append remaining prime factor if necessary
    if n > 1
        push!(factorization, (n, 1))
    end

    return factorization
end

end

ext: jl

max_stars_repo_name: curtislb/ProjectEuler

max_issues_repo_name: curtislb/ProjectEuler

max_forks_repo_name: curtislb/ProjectEuler


text:
@testset "Filter tests" begin
  using Test
  using CLIMA
  import GaussQuadrature
  using MPI
  using LinearAlgebra

  let
    # Values computed with:
    #   https://github.com/tcew/nodal-dg/blob/master/Codes1.1/Codes1D/Filter1D.m
    W = [0x3fe98f3cd0d725e8   0x3fddfd863c6c9a44   0xbfe111110d0fd334   0x3fddbe357bce0b5c   0xbfc970267f929618
         0x3fb608a150f6f927   0x3fe99528b1a1cd8d   0x3fcd41d41f8bae45   0xbfc987d5fabab8d5   0x3fb5da1cd858af87
         0xbfb333332eb1cd92   0x3fc666666826f178   0x3fe999999798faaa   0x3fc666666826f176   0xbfb333332eb1cd94
         0x3fb5da1cd858af84   0xbfc987d5fabab8d4   0x3fcd41d41f8bae46   0x3fe99528b1a1cd8e   0x3fb608a150f6f924
         0xbfc970267f929618   0x3fddbe357bce0b5c   0xbfe111110d0fd333   0x3fddfd863c6c9a44   0x3fe98f3cd0d725e8]
    W = reinterpret.(Float64, W)

    N = size(W, 1) - 1

    topology = CLIMA.Mesh.Topologies.BrickTopology(MPI.COMM_SELF, -1.0:2.0:1.0);

    grid = CLIMA.Mesh.Grids.DiscontinuousSpectralElementGrid(topology;
                                                        polynomialorder = N,
                                                        FloatType = Float64,
                                                        DeviceArray = Array)

    filter = CLIMA.Mesh.Filters.ExponentialFilter(grid, 0, 32)
    @test filter.filter ≈ W
  end

  let
    # Values computed with:
    #   https://github.com/tcew/nodal-dg/blob/master/Codes1.1/Codes1D/Filter1D.m
    W = [0x3fd822e5f54ecb62   0x3fedd204a0f08ef8   0xbfc7d3aa58fd6968   0xbfbf74682ac4d276
         0x3fc7db36e726d8c1   0x3fe59d16feee478b   0x3fc6745bfbb91e20   0xbfa30fbb7a645448
         0xbfa30fbb7a645455   0x3fc6745bfbb91e26   0x3fe59d16feee478a   0x3fc7db36e726d8c4
         0xbfbf74682ac4d280   0xbfc7d3aa58fd6962   0x3fedd204a0f08ef7   0x3fd822e5f54ecb62]
    W = reinterpret.(Float64, W)

    N = size(W, 1) - 1

    topology = CLIMA.Mesh.Topologies.BrickTopology(MPI.COMM_SELF, -1.0:2.0:1.0);
    grid = CLIMA.Mesh.Grids.DiscontinuousSpectralElementGrid(topology;
                                                        polynomialorder = N,
                                                        FloatType = Float64,
                                                        DeviceArray = Array)

    filter = CLIMA.Mesh.Filters.ExponentialFilter(grid, 1, 4)
    @test filter.filter ≈ W
  end

  let
    T = Float64
    N = 5
    Nc = 4

    topology = CLIMA.Mesh.Topologies.BrickTopology(MPI.COMM_SELF, -1.0:2.0:1.0);
    grid = CLIMA.Mesh.Grids.DiscontinuousSpectralElementGrid(topology;
                                                        polynomialorder = N,
                                                        FloatType = T,
                                                        DeviceArray = Array)

    ξ = CLIMA.Mesh.Grids.referencepoints(grid)
    a, b = GaussQuadrature.legendre_coefs(T, N)
    V = GaussQuadrature.orthonormal_poly(ξ, a, b)

    Σ = ones(T, N+1)
    Σ[(Nc:N).+1] .= 0

    W = V*Diagonal(Σ)/V

    filter = CLIMA.Mesh.Filters.CutoffFilter(grid, Nc)
    @test filter.filter ≈ W
  end
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
module CompileBot

export BotConfig, snoop_bot, snoop_bench

using Core: MethodInstance, CodeInfo
using YAML

if VERSION <=  v"1.1"
    isnothing(x) = x === nothing
end
if VERSION <=  v"1.2"
    Base.print(io, ::Nothing) = Base.print(io, "")
    Base.print(::Base.GenericIOBuffer{Array{UInt8,1}}, ::Nothing) = Base.print(io, "")
    Base.string(::Nothing) = ""
end

################################################################
const UStrings = Union{AbstractString,Regex,AbstractChar}

"""
    BotConfig(package_name::AbstractString ; exclusions, os, else_os, version, else_version, package_path, precompiles_rootpath, subst, tmin)

Construct a SnoopCompile bot configuration. `package_name` is the name of the package. This object is passed to [`snoop_bot`](@ref)
and [`snoop_bench`](@ref).

You may supply the following optional **keyword** arguments:

- `exclusions` : A vector of of Strings (or RegExp) to exclude some functions from being precompiled

- `os`: A vector of of Strings (or RegExp) to support with precompile statements.

Example: `os = ["windows", "linux"]`

- `else_os`: If you want to use a specific operating system's precompile file as the default,
  set `else_os` to the name of that os.
  Not passing this argument skips precompilation on any operating system other than those
  explicitly listed in `os`.

Example: `else_os = "linux"`

- `version`: A vector of of Julia versions used to generate precompile signatures.

Example: `version = [v"1.1",v"1.5.0", "nightly"]`

It is assumed that the generated precompile signatures are valid for patch versions of Julia (e.g. givingv"1.5.0" supports v"1.4.0" to v"1.4.9").

- `else_version`: the Julia version used to generate the default signatures for other `version`s.
Not passing this argument skips precompilation on any Julia version other than those explicitly listed in `version`.

Example: `else_version =v"1.5.0"`

- `yml_path`: instead of directly passing `os` and `version` to BotConfig, you can pass `yml_path` which should be the GitHub actions YAML path or file name.
It assumes that the job name is `SnoopCompile`.

Example: `yaml_path = "SnoopCompile.yml"`

- `package_path`: path to the main `.jl` file of the package (similar to `pathof`).
  Default path is `pathof_noload(package_name)`.

- `precompiles_rootpath`: the path where precompile files are stored.
  Default path is "\$(dirname(dirname(package_path)))/deps/SnoopCompile/precompile".

- `subst` : A vector of pairs of Strings (or RegExp) to replace a packages precompile statements
  with another's package like `["ImageTest" => "Images"]`.

- `tmin`: Methods that take less time than `tmin` to be inferred will not be added to the
  precompile statements. Defaults to 0.01. Set it to 0.0 if you want to include all the sentences.

- `check_eval`: By default, the bot discards the precompile statements that cannot be `eval`ed.

In rare cases (when snooping is very time consuming), you may want to do this manually by using the printed errors to add the problematic functions to `exclusions` and then set `check_eval=false` for the future runs.

# Example
```julia
botconfig1 = BotConfig(
  "Zygote";                            # package name (the one this configuration lives in)
  os = ["linux", "windows", "macos"],  # operating systems for which to precompile
  version = [v"1.5.0", v"1.3.1"],      # supported Julia versions
  exclusions = ["SqEuclidean"],         # exclude functions (by name) that would be problematic if precompiled
)

botconfig2 = BotConfig(
  "Zygote";                            # package name (the one this configuration lives in)
  yml_path = "SnoopCompile.yml"        # parse `os` and `version` from `SnoopCompile.yml`
  exclusions = ["SqEuclidean"],         # exclude functions (by name) that would be problematic if precompiled
)

# A full example:
BotConfig("MatLang", exclusions = ["badfunction"], os = ["linux", "windows", "macos"], else_os = "linux", version = ["1.5.0", "1.2", "1.0.5"], else_version = "1.5.0" )

# Different examples for other possibilities:
BotConfig("MatLang")

BotConfig("MatLang", exclusions = ["badfunction"])

BotConfig("MatLang", os = ["linux", "windows"])

BotConfig("MatLang", os = ["windows", "linux"], else_os = "linux")

BotConfig("MatLang", version = [v"1.1",v"1.5.0"])

BotConfig("MatLang", version = [v"1.1",v"1.5.0"], else_version =v"1.5.0")

BotConfig("MatLang", os = ["linux", "windows"], version = [v"1.1",v"1.5.0"])
```
"""
struct BotConfig
    package_name::AbstractString
    exclusions::Vector{UStrings}
    check_eval::Bool
    os::Union{Vector{String}, Nothing}
    else_os::Union{String, Nothing}
    version::Union{Vector{VersionNumber}, Nothing}
    else_version::Union{VersionNumber, Nothing}
    package_path::AbstractString
    precompiles_rootpath::AbstractString
    subst::Vector{Pair{UStrings, UStrings}}
    tmin::AbstractFloat
end

function BotConfig(
    package_name::AbstractString;
    exclusions::AbstractVector = String[],
    check_eval::Bool = true,
    os::Union{Vector{String}, Nothing} = nothing,
    else_os::Union{String, Nothing} = nothing,
    version::Union{Vector{<:Union{VersionNumber,String}}, Nothing} = nothing,
    else_version::Union{VersionNumber, String, Nothing} = nothing,
    package_path::AbstractString = pathof_noload(package_name),
    precompiles_rootpath::AbstractString = "$(dirname(dirname(package_path)))/deps/SnoopCompile/precompile",
    subst::AbstractVector = Vector{Pair{UStrings, UStrings}}(),
    tmin::AbstractFloat = 0.01,
    yml_path::Union{String, Nothing} = nothing,
    blacklist = nothing       # deprecated keyword
    )

    if blacklist !== nothing
        Base.depwarn("`blacklist` is deprecated, please use `exclusions` to pass a list of excluded names", :BotConfig)
        append!(exclusions, blacklist)
    end

    package_root_path = dirname(dirname(package_path))

    # Error for workflow file
    yml_path_error = "SnoopCompile.yml"
    try
        yml_path_error = searchdirsboth([ pwd(), package_root_path, "$package_root_path/.github/workflows/"], yml_path_error)
    catch
        # file not found
    finally
        if isfile(yml_path_error) && !occursin("CompileBot.postprocess()", Base.read(yml_path_error, String))
            error("""
            CompileBot is now in a separate repository, and the API is changed because of that. Call `using CompileBot` directly in your snoop scripts and update your workflow based on this guide: [Configure the bot to run with a GitHub Action file]( https://aminya.github.io/CompileBot.jl/dev/#Configure-the-bot-to-run-with-a-GitHub-Action-file-1)
            """)
        end
    end

    # Parse os and version from the yaml file
    if !isnothing(yml_path)
        yml_path = searchdirsboth([ pwd(), package_root_path, "$package_root_path/.github/workflows/"], yml_path)
        if !isfile(yml_path)
            error("$yml_path not found")
        end

        # TODO This can be an option
        workflow_job = "SnoopCompile"

        yml = YAML.load_file(yml_path)
        matrix = yml["jobs"][workflow_job]["strategy"]["matrix"]
        if haskey(matrix, "os")
            os = matrix["os"]
        end
        if haskey(matrix, "version")
            version = matrix["version"]
        end
    end

    if !isnothing(version)
        version = JuliaVersionNumber.(version)
    end
    if !isnothing(else_version)
        else_version = JuliaVersionNumber(else_version)
    end

    return BotConfig(package_name, exclusions, check_eval, os, else_os, version, else_version, GoodPath(package_path), GoodPath(precompiles_rootpath), subst, tmin)
end

include("botutils.jl")
include("precompile_include.jl")
include("precompile_activation.jl")
include("snoop_bot.jl")
include("snoop_bench.jl")
include("postprocess.jl")

# Using GitHubActions logger in CI
using Logging: global_logger
using GitHubActions: GitHubActionsLogger
function __init__()
    get(ENV, "GITHUB_ACTIONS", "false") == "true" && global_logger(GitHubActionsLogger())
end

end # module


text:
using Distance
using NearestNeighbors
using StatsBase

function benchmark{T <: NearestNeighborTree, V}(k, d, n, r, warmup, ::Type{V}, ::Type{T})
	E = eltype(V)
	X = rand(E, d, n)
	tree = T(X, Euclidean())

	skip = warmup
	times = Float64[]
	print("pass: ")
	for i = 1 : r
		numReps = 50
	    while true
	    	time = time_ns()
	    	for ii = 1 : numReps
	  			nearest(tree, rand(E, 3), k)
	  		end
	    	time = time_ns() - time;

	    	t = float(time)/1.0e+09

	    	if t < 1.0
	    		numReps *= 5
	    		continue
	    	else
	    		if skip > 0
	    			skip -= 1
	    			print("$(i)* ")
	    		else
	    			push!(times, numReps/t)
	    			print("$(i) ")
	    		end
	    		break
	    	end
		end
	end

	# Statistics
	println("\nImpl: ", T,
			", Key type: ", V,
			", nn = $(k), searches/sec = ",int(mean(times)), " ± ", int(1.96*sem(times)), " [CI:95%]")
	return times
end

d = 3
n = 1000
r = 51
w = 1

k = 1
times = benchmark(k, d, n, 3, w, Float32, NaiveNeighborTree)
times = benchmark(k, d, n, r, w, Float64, NaiveNeighborTree)
times = benchmark(k, d, n, r, w, Float32, KDTree)
times = benchmark(k, d, n, r, w, Float64, KDTree)

k = 5
times = benchmark(k, d, n, r, w, Float32, NaiveNeighborTree)
times = benchmark(k, d, n, r, w, Float64, NaiveNeighborTree)
times = benchmark(k, d, n, r, w, Float32, KDTree)
times = benchmark(k, d, n, r, w, Float64, KDTree)

# Results [CI:95%]
# java kd-tree, nn = 1, searches/sec = 242236 ± 10386
# NaiveNeighborTree, Float32, nn = 1, searches/sec = 3600 ± 14
# NaiveNeighborTree, Float64, nn = 1, searches/sec = 3451 ± 12
# KDTree, Float32, nn = 1, searches/sec = 291513 ± 1474
# KDTree, Float64, nn = 1, searches/sec = 301366 ± 778

# java kd-tree, nn = 5, searches/sec =  102855 ± 10829
# NaiveNeighborTree, Float32, nn = 5, searches/sec = 2818 ± 25
# NaiveNeighborTree, Float64, nn = 5, searches/sec = 2757 ± 26
# KDTree, Float32, nn = 5, searches/sec = 120687 ± 432
# KDTree, Float64, nn = 5, searches/sec = 126815 ± 312


text:
##
"""
IDEAL OBSERVER MODEL: Probability of switch compared to no switch at poke k after consecutive failures
k = consecutive failure
C = estimated state (L or R, turned in High and Low as in being already in the rewarded side)
Pr = probability

Qk💠 = Pr(Ck+1 ≠ Ck|t[k0:k],C[k0:k]) / Pr(Ck+1 = Ck|t[k0:k],C[k0:k]) {P(switch)/P(no switch)}

i = trial
t = confidence, depends on
    yi =  p of correct discrimination(ABSENT IN THE Flipping)  equal to Prew; they used the subjective psychometric function
    λi =  p of transition at trial i given the hazard of Psw

P of not reward (in absence of switch) = λk0 +  Σ[i = k0:k](λi(Π[j = k0:i-1])) / Π[i = k0:k](1 - yi)(1 - λi)

CONFIDENCE-BASED MODEL
χΣ = trial by trial evidence for a switch, increment with a magnitude dipendent on the expected accuracy (confidence)
    it has a gaussian distribution at each trial with μχ and σχ  such that CBM was equal to IOM for 1B-Er. σχ grows linearly
    with the number of error trial (failures) (σ²χ for N failure = N σχ of first error)
Θ = threshold

Making CBM close to IOM
    Qk🌸 = P(χΣ > = Θ) \ P(χΣ < Θ) = Integral(Θ,inf)(N(μχ,σχ)dχΣ) / 1 - Integral(Θ,inf)(N(μχ,σχ)dχΣ) = Qk💠

    REARRANGED IN
    Integral(-inf, Θ) (N(μχ,σχ)dχΣ) = 1 - (Qk💠 / 1 - Qk💠)

α = perseverance factor

being μχΣ fixed by the IOM and its relationship with σχΣ the model can be reduced to 2 parameter σχΣ and α Ψ = {σχΣ, α}
probability of switch become

    p(Xy/n = SW|ψ) Integral(θ,inf) (N(α(μχ),σχ)dχΣ)
"""
##
using QuadGK
using Distributions
using Optim
using FileIO
using DataFrames
using StatsBase
using StatsFuns
using StructArrays
using Plots
import Flipping.nextcount
using LinearAlgebra
##
quadgk(f, a,b,c...; rtol=sqrt(eps), atol=0, maxevals=10^7, order=7, norm=norm)
fun(x) = exp(x)
quadgk(fun,-1.0,1.0)

plot(pdf.(Normal(), -1.0:0.2:1.0))

rand(Normal(1,0.5))
#########

function χ(p, leave::Bool, fail)
    failure_index = Int64(fail) + 1
    μ,σ = p
    if σ < 0
        return 10000000
    else
        store = Vector{Float64}(undef,failure_index)
        for i in 1:length(failure_index)
            store[i] = rand(Normal(μ,i*σ))
        end
        accumulation = Base.accumulate(+,store)
        result = accumulation[end]
        return - log(exp(-result)+1) - !leave * result
    end
end

function smart_data_evidence(p, leaves, failure_indices)
    soa = StructArray((leaves, failure_indices))
    cm = countmap(soa)
    smart_data_evidence(p, cm)
end

function smart_data_evidence(p, cm)
    sum(n*χ(p, value...) for (value, n) in cm)
end
###############
function failure_idx(count::T, rewarded) where {T <: Number}
    if rewarded
        count = 0.0
    else
        count + 1.0
    end
end
###############

filename = "/home/beatriz/mainen.flipping.5ht@gmail.com/Flipping/Datasets/Stimulations/DRN_Opto_Flipping/pokesDRN_Opto_Flipping.csv"
data = FileIO.load(filename) |> DataFrame
data[!,:LastPoke] .= false
data[!,:FailuresIdx] .= 0.0
by(data,[:Session,:Streak]) do dd
    dd.LastPoke[end] = true
    dd[:,:FailuresIdx] = Base.accumulate(failure_idx, occursin.("true",dd.Reward);init=0.0)
end
checkpoint = by(data,:Session) do d
    d.PokeIn[1] == d.PokeIn[2]
end
any(checkpoint[:,2])
println(describe(data))
filtered = filter(t -> t[:Protocol].!="45/15", data)
#################
any(isnan.(filtered[:,:LastPoke]))
any(isnan.(filtered[:,:FailuresIdx]))
#################
coll =by(filtered,[:Protocol,:Wall,:MouseID]) do dd
    leaves = dd[:,:LastPoke]
    failures_indices = dd[:,:FailuresIdx]
    ϵ = 10^(-2)
    res = optimize(p -> -smart_data_evidence(p, leaves, failures_indices) + ϵ*norm(p)^2, [1.0,1.0])
    mins = Optim.minimizer(res)
    DataFrame(Mu_G = mins[1], Sigma_N = mins[2])
    #DataFrame(Inverse_Temp = mins[1], Integration = mins[2])
end
println(describe(coll))
f = plot(; legend = :bottomright)
param_plot(coll,:Wall,[:Mu_G,:Sigma_N])
###############
function param_plot(df::AbstractDataFrame,split_var::Symbol,params::Union{Tuple,AbstractVector};groups = :MouseID)
    z = []
    sort!(df, [groups,split_var])
    dof = union(df[:,split_var])
    idxs_x = df[:,split_var].== dof[1]
    idxs_y = df[:,split_var].== dof[end]
    for par in params
        o = scatter(df[idxs_x,par],df[idxs_y,par],
        group = df[idxs_x,groups],
        legend = false,
        xlabel = dof[1],
        ylabel = dof[end],
        title = string(par),
        markersize = 1)
        Plots.abline!(1,0)
        push!(z,o)
    end
    plot(z...)
end


text:
using SeisProcessing
using Test

d = SeisLinearEvents();
snr_i =1.2
dn = SeisAddNoise(d,snr_i);

snr = MeasureSNR(d,dn)

# test that snr = snr_i
quality_factor = abs(snr-snr_i)
@test quality_factor < 1


text:
export ExponentialLinearQuadratic

import Distributions: pdf, logpdf, ContinuousUnivariateDistribution

struct ExponentialLinearQuadratic{ A <: AbstractApproximationMethod, T <: Real } <: ContinuousUnivariateDistribution
    approximation :: A
    a :: T
    b :: T
    c :: T
    d :: T
end

ExponentialLinearQuadratic(approximation, a::Real, b::Real, c::Real, d::Real)             = ExponentialLinearQuadratic(approximation, promote(a, b, c, d)...)
ExponentialLinearQuadratic(approximation, a::Integer, b::Integer, c::Integer, d::Integer) = ExponentialLinearQuadratic(approximation, float(a), float(b), float(c), float(d))

Distributions.pdf(dist::ExponentialLinearQuadratic, x::Real)    = exp(logpdf(dist, x))
Distributions.logpdf(dist::ExponentialLinearQuadratic, x::Real) = -0.5 * (dist.a * x + dist.b * exp(dist.c * x + dist.d * x ^ 2 / 2.0))
Distributions.mean(dist::ExponentialLinearQuadratic)            = approximate_meancov(dist.approximation, (z) -> pdf(dist, z) * exp(0.5 * z^2), NormalMeanVariance(0.0,1.0))[1]
Distributions.var(dist::ExponentialLinearQuadratic)             = approximate_meancov(dist.approximation, (z) -> pdf(dist, z) * exp(0.5 * z^2), NormalMeanVariance(0.0,1.0))[2]
Distributions.cov(dist::ExponentialLinearQuadratic)             = var(dist)

precision(dist::ExponentialLinearQuadratic)    = inv(var(dist))
weightedmean(dist::ExponentialLinearQuadratic) = precision(dist)*mean(dist)

prod_analytical_rule(::Type{ <: UnivariateNormalDistributionsFamily }, ::Type{ <: ExponentialLinearQuadratic }) = ProdAnalyticalRuleAvailable()

function prod(::ProdAnalytical, left::UnivariateNormalDistributionsFamily, right::ExponentialLinearQuadratic)
    mean, variance = approximate_meancov(right.approximation, (z) -> pdf(right, z), left)
    return NormalMeanVariance(mean, variance)
end


text:
@testset "Hostlist" begin
  slurm_hostlist = "linux[1-3,6]"
  hostlist = Hostlist(slurm_hostlist)
  @test all(collect(hostlist) .== ["linux1", "linux2", "linux3", "linux6"])
end


text:
# Toy problem to test the solver

using luteos

P = P1() # Polynomial order of solution

println("Generate mesh")
@time mesh   = Mesh3D( "cube", P, N = 13)#9)
println("Generate master")
@time master = Master3D( P )

# source function
function funcS( p::Array{Float64} )
  return -3*pi^2 * sin.(pi * (p[:,1]-0.5)) .* sin.(pi * (p[:,2]-0.5)) .* sin.(pi * (p[:,3]-0.5))
end

function funcB( p::Array{Float64} )
  return fill(0.0, size(p,1), 1)
end

function funcB2( p::Array{Float64} )
  return cos.( π * (p[:,1] - 0.5) ).^2 .* cos.( π * (p[:,3] - 0.5) ).^2
end

function funcU( p::Array{Float64} )
  return sin.(pi * (p[:,1]-0.5)) .* sin.(pi * (p[:,2]-0.5)) .* sin.(pi * (p[:,3]-0.5))
end

println("Setup problem")
@time prob = Problem( "Example 4 - Poisson", 1.0, [0.0,0.0,0.0], funcS, [1,1,1,2,1,1], true, [funcB, funcB, funcB, funcB, funcB, funcB] )

println("Solve problem")
@time (uhath, uh, qh, uhathTri ) = hdgSolve( master, mesh, prob )

# write solution
println("Write solution")
@time writeTecplot( "blaCD3D.dat", prob, mesh, uh, qh )


text:
"""
    AlgoDeforNonlinearModule

Module for algorithms used in nonlinear deformation models.
"""
module AlgoDeforNonlinearModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtools.AlgoBaseModule: dcheck!
import Arpack: eigs
import SparseArrays: spzeros
import LinearAlgebra: mul!
my_A_mul_B!(C, A, B) = mul!(C, A, B)
import FinEtools.FieldModule: AbstractField, ndofs, setebc!, numberdofs!, applyebc!, scattersysvec!, wipe!, copyto!, incrscattersysvec!
import FinEtools.NodalFieldModule: NodalField, nnodes
import FinEtools.FEMMBaseModule: associategeometry!, distribloads, fieldfromintegpoints, elemfieldfromintegpoints
import ..FEMMDeforNonlinearBaseModule: stiffness, geostiffness, nzebcloads, restoringforce
import FinEtoolsDeforLinear.DeforModelRedModule: stresscomponentmap
import FinEtools.ForceIntensityModule: ForceIntensity, settime!
import FinEtools.MeshModificationModule: meshboundary
import FinEtools.MeshExportModule.VTK: vtkexportmesh
import LinearAlgebra: eigen, qr, dot, cholesky, sum
# using Debugger

"""
    AlgoDeforNonlinearModule.nonlinearstatics(modeldata::FDataDict)

Algorithm for static nonlinear deformation (stress) analysis.

The algorithm chooses steps from the array of load multipliers: the step
takes it precisely from the preceding step to the next step in one go.

# Argument
`modeldata` = dictionary with values for keys

* `"fens"`  = finite element node set
* `"regions"`  = array of region dictionaries
* `"essential_bcs"` = array of essential boundary condition dictionaries
* `"traction_bcs"` = array of traction boundary condition dictionaries
* `"temperature_change"` = dictionary of data for temperature change

For each region (connected piece of the domain made of a particular
material), mandatory, the  region dictionary  contains values for keys:
* `"femm"` = finite element model machine (mandatory);

For essential boundary conditions (optional) each dictionary
would hold
  + `"displacement"` = when this key is not present, the assumption is
        that the displacement is fixed at zero (0). Otherwise, this needs to
        be set to a function with signature `f(x, lambda)`, where `x` is the
        location of the node in the reference configuration, and `lambda` is the
        load factor. In other words, whenever this quantity is supplied, it is
        implied that the displacement depends on the load factor.
  + `"component"` = which component is prescribed  (1, 2, 3)?
  + `"node_list"` = list of nodes on the boundary to which the condition
        applies (mandatory)

For traction boundary conditions (optional) each dictionary
would hold
  + `"femm"` = finite element model machine (mandatory);
  + `"traction_vector"` = traction vector, a force-intensity
    (`ForceIntensity`) object.

Control parameters
The following attributes  may be supplied:
* `"load_multipliers"` = For what load multipliers should the solution be
          calculated? Array of monotonically increasing numbers.
* `"line_search"` = Should we use line search? Boolean.  Default = true.
* `"maxdu_tol"` = Tolerance on the magnitude  of the largest incremental
          displacement component.
* `"maxbal_tol"` = Tolerance on the magnitude of the largest out-of-balance
          force component.
* `"iteration_observer"` = observer function to be called
      after each iteration.  Default is to do nothing.
      The observer function has a signature
                `iteration_observer(lambda,iter,du,modeldata)`
      where `lambda` is the current load factor, `iter` is the iteration
      number, `du` is the nodal field of current displacement increments.
* `"increment_observer"` = observer function
      to be called after convergence is reached in each step (optional)
      The observer function has a signature
                `output(lambda, modeldata)`
      where `lambda` is the current load factor. Default is to do nothing.
      The increment observer can refer to the following key-value pairs in `modeldata`:
      - "un1": converged displacement in current step
      - "un": converged displacement in the last step
      - "t": current value of load factor
      - "dt": increment of the load factor

# Output
`modeldata` = the dictionary on input is augmented with the keys
* `"geom"` = the nodal field that is the geometry
* `"u"` = the nodal field that is the computed displacement
* `"reactions"` = computed reaction nodal field
* `"timing"` = dictionary with timing results
"""
function nonlinearstatics(modeldata::FDataDict)
    # Lists of recognized keys for the data dictionaries:
    modeldata_recognized_keys = ["fens", "regions", "essential_bcs",  "traction_bcs", "factorize"]
    essential_bcs_recognized_keys = ["displacement", "node_list", "component"]
    traction_bcs_recognized_keys = ["femm", "traction_vector"]
    regions_recognized_keys = ["femm", "body_load"]

    # Should we use line search?
    line_search = get(modeldata, "line_search", false);

    # For what load multipliers should dissolution be calculated?
    load_multipliers = get(modeldata, "load_multipliers", 1.0);

    # Is something to be done after each iteration?
    iteration_observer = get(modeldata, "iteration_observer", nothing);

    increment_observer = get(modeldata, "increment_observer", (lambda, modeldata) -> println("Load multiplier = $(lambda)"));

    # Tolerance on the magnitude  of the largest incremental displacement component
    maxdu_tol = get(modeldata, "maxdu_tol", 0.0);

    # Tolerance on the magnitude  of the out-of-balance force
    maxbal_tol = get(modeldata, "maxbal_tol", 0.0);

    # Maximum number of iterations allowed
    maxiter = get(modeldata, "maxiter", 4);

    # Extract the nodes
    fens=get(()->error("Must get fens!"), modeldata, "fens")

    # Construct the geometry field
    geom = NodalField(fens.xyz)

    # Construct the displacement field
    un1 = NodalField(zeros(nnodes(geom), ndofs(geom)))

    modeldata["timing"] = FDataDict()

    # Apply the essential boundary conditions on the displacement field
    essential_bcs = get(modeldata, "essential_bcs", nothing);
    if (essential_bcs != nothing)
        for j = 1:length(essential_bcs)
            ebc = essential_bcs[j]
            dcheck!(ebc, essential_bcs_recognized_keys)
            fenids = get(()->error("Must get node list!"), ebc, "node_list");
            displacement = get(ebc, "displacement", nothing);
            u_fixed = zeros(FFlt, length(fenids)); # default is  zero displacement
            if (displacement != nothing) # if it is nonzero, it is dependent on the load factor
                ebc["load_factor_dependent"] = true
            end
            component = get(ebc, "component", 0); # which component?
            setebc!(un1, fenids[:], true, component, u_fixed);
        end
        applyebc!(un1);
    end

    # Number the equations
    numberdofs!(un1)
    # Converged displacement at the beginning of the current step, i. e.  u_{n}
    un = deepcopy(un1);
    # Converged displacement one step back, i. e.  u_{n-1}
    unm1 = deepcopy(un);

    # Create the necessary FEMMs
    regions = get(()->error("Must get region list!"), modeldata, "regions")
    for i = 1:length(regions)
        region = regions[i]
        dcheck!(region, regions_recognized_keys)
        femm = region["femm"];
        # Give the  FEMM a chance  to precompute  geometry-related quantities
        femm = associategeometry!(femm, geom);
    end

    # Initially the stiffness matrix is empty.  We can check it  in the
    # treatment of nonzero essential boundary conditions and construct it if
    # needed.
    K = nothing;

    # Increment magnitudes
    load_increments = diff(unique(sort(vcat([0], vec(load_multipliers)))));

    # For all  load increments
    for incr in 1:length(load_increments) # Load-incrementation loop

        lambda = sum(load_increments[1:incr]);
        dlambda = load_increments[incr];

        # Initial value of the displacement  in the current step
        copyto!(un1, un)

        # Process load-factor-dependent essential boundary conditions:
        aany_nonzero_EBC = false;
        essential_bcs = get(modeldata, "essential_bcs", nothing);
        if (essential_bcs != nothing)
            for j = 1:length(essential_bcs)
                ebc = essential_bcs[j]
                if haskey(ebc, "load_factor_dependent")
                    fenids = get(()->error("Must get node list!"), ebc, "node_list");
                    displacement = ebc["displacement"]
                    component = get(ebc, "component", 0); # which component?
                    for k in 1:length(fenids)
                        u_fixed = displacement(geom.values[fenids[k],:], lambda);
                        setebc!(un1, [fenids[k]], true, component, u_fixed[1]);
                        aany_nonzero_EBC = aany_nonzero_EBC || (u_fixed[1] != 0.0);
                    end
                    applyebc!(un1);
                    numberdofs!(un1)
                end
            end
        end

        # Initialization
        applyebc!(un1); # Apply EBC
        numberdofs!(un1)
        du = deepcopy(un1); # Displacement increment
        applyebc!(du);

        # Initialize the load vector
        F = fill(0.0, un1.nfreedofs);
        FL = similar(F);

        # If any boundary conditions are inhomogeneous, calculate  the force
        # vector due to the displacement increment. Then update the guess of
        # the new displacement.
        if (aany_nonzero_EBC)
            du.fixed_values .= un1.fixed_values .- un.fixed_values
            for i = 1:length(regions)
                femm = deepcopy(regions[i]["femm"])
                F .+= nzebcloads(femm, geom, un, unm1, du, lambda, dlambda);
            end
            # Provided we got converged results  in the last step, we already
            # have a usable stiffness matrix.
            if K == nothing # we don't have a stiffness matrix
                K = spzeros(un1.nfreedofs, un1.nfreedofs);
                for i = 1:length(regions)
                    femm = deepcopy(regions[i]["femm"])
                    K = K + stiffness(femm, geom, un, unm1, lambda, dlambda);
                    K = K + geostiffness(femm, geom, un, unm1, lambda, dlambda);
                end
            end
            incrscattersysvec!(un1, K\F);
        end

        # Iteration loop
        iterationsuccessful = false
        iter = 1;
        while true #  Iteration loop

            # Initialize the load vector
            fill!(F, 0.0)

            # Construct the system stiffness matrix.
            K =  spzeros(un1.nfreedofs,un1.nfreedofs);
            for i = 1:length(regions)
                femm = regions[i]["femm"];
                K = K + stiffness(femm, geom, un1, un, lambda, dlambda);
                K = K + geostiffness(femm, geom, un1, un, lambda, dlambda);
            end

            # # Process the body load
            # if (isfield(model_data, 'body_load' ))
            #     for j=1:length(model_data.body_load)
            #         body_load =model_data.body_load{j};
            #         femm = femm_deformation (struct ('material',[],...
            #             'fes',body_load.fes,...
            #             'integration_rule',body_load.integration_rule));
            #         fi= force_intensity(struct('magn',body_load.force));
            #         F = F + distrib_loads(femm, sysvec_assembler, geom, un1, fi, 3);
            #     end
            #     clear body_load fi  femm
            # end

            # Process the traction boundary condition
            traction_bcs = get(modeldata, "traction_bcs", nothing);
            if (traction_bcs != nothing)
                for j=1:length(traction_bcs)
                    tractionbc = traction_bcs[j]
                    dcheck!(tractionbc, traction_bcs_recognized_keys)
                    femm = tractionbc["femm"]
                    fi = tractionbc["traction_vector"];
                    settime!(fi, lambda)
                    F .= F .+ distribloads(femm, geom, un1, fi, 2);
                end
            end

            # # Process the nodal force boundary condition
            # if (isfield(model_data.boundary_conditions, 'nodal_force' ))
            #     for j=1:length(model_data.boundary_conditions.nodal_force)
            #         nodal_force =model_data.boundary_conditions.nodal_force{j};
            #         femm = femm_deformation (struct ('material',[],...
            #             'fes',fe_set_P1(struct('conn',reshape(nodal_force.node_list,[],1))),...
            #             'integration_rule',point_rule));
            #         fi= force_intensity(struct('magn',nodal_force.force));
            #         F = F + distrib_loads(femm, sysvec_assembler, geom, un1, fi, 0);
            #     end
            #     clear nodal_force fi femm
            # end

            #  External loads vector
            copyto!(FL, F);

            # The restoring  force vector
            FR = fill(0.0, un1.nfreedofs);
            for i = 1:length(regions)
                femm = regions[i]["femm"];
                FR .+= restoringforce(femm, geom, un1, un, lambda, dlambda);
            end

            # Solve the system of linear algebraic equations
            F .= (FL .+ FR);
            dusol = K\F;

            #  Distribute the solution
            fill!(du.values, 0.0)
            scattersysvec!(du, dusol);

            # Do we need line search?
            eta = 1.0;
            if (line_search)
                un1p = deepcopy(un1)
                R0 = dot(F, dusol);
                for  linesrch in 1:1 # How many  searches should we do?
                    fill!(FR, 0.0);
                    un1p.values .= un1.values .+ eta.*du.values
                    for i = 1:length(regions)
                        femm = regions[i]["femm"];
                        FR .+= restoringforce(femm, geom, un1p, un, lambda, dlambda);
                    end
                    F .= FL .+ FR;
                    R1 = dot(F, dusol);
                    a = R0/R1;
                    if ( a < 0 )
                        eta = a/2 +sqrt((a/2)^2 -a);
                    else
                        eta = a/2;
                    end
                    eta = min(eta, 1.0);
                end
            end

            # Increment the displacements
            un1.values .= un1.values .+ eta*du.values;

            # Compute the increment data
            maxdu  = eta * maximum(abs.(du.values));
            maxbal = maximum(abs.(F));

            # Report  iteration results?
            if iteration_observer != nothing
                setindex!(modeldata, geom, "geom");
                setindex!(modeldata, un1, "un1");
                setindex!(modeldata, un, "un");
                setindex!(modeldata, lambda, "t");
                setindex!(modeldata, dlambda, "dt");
                setindex!(modeldata, maxdu, "maxdu");
                setindex!(modeldata, maxbal, "maxbal");
                iteration_observer(lambda, iter, du, modeldata);
            end

            # Converged?
            if (maxdu <= maxdu_tol)  || (maxbal <= maxbal_tol)
                iterationsuccessful = true
                break; # Breakout of the iteration loop
            end

            iter = iter+1;
            if iter > maxiter
                @error "Maximum number of iterations reached"
                break
            end

        end #  Iteration loop

        if !iterationsuccessful
            break
        end

        # Update the model data
        setindex!(modeldata, un1, "un1");
        setindex!(modeldata, un, "un");
        setindex!(modeldata, lambda, "t");
        setindex!(modeldata, dlambda, "dt");

        # Calculate the reactions
        # This field holds the current converged displacements;
        un1Allfree = deepcopy(un1);
        fill!(un1Allfree.is_fixed, false);# All nodes, all dofs are free.
        numberdofs!(un1Allfree);
        copyto!(un1Allfree.values, un1.values);
        unAllfree = deepcopy(un1Allfree);
        copyto!(unAllfree.values, un.values);
        # Since all the degrees of freedom are free, the restoring force method
        # can be used to compute the reactions. Note the negative sign: the
        # reactions are the opposite of the resisting forces of the
        # material.
        FR = fill(0.0, un1Allfree.nfreedofs);
        for i = 1:length(regions)
            femm = deepcopy(regions[i]["femm"])
            FR .-= restoringforce(femm, geom, un1Allfree, unAllfree, lambda, dlambda);
        end
        reactions = deepcopy(un1Allfree);
        scattersysvec!(reactions, FR);
        setindex!(modeldata, reactions, "reactions");

        # Converged
        # Update the FEMMs: save the material state vectors
        for i = 1:length(regions)
            femm = regions[i]["femm"];
            restoringforce(femm, geom, un1, un, lambda, dlambda, true);
        end

        # Report results
        increment_observer(lambda, modeldata);

        # Reset  the displacement field for the next load step
        copyto!(unm1, un);
        copyto!(un, un1);
        fill!(F, 0.0);

    end # Load-incrementation loop

    # Return some computed quantities: we stored most of them inside the
    # incrementation loop, now we add a reference to the final displacement
    setindex!(modeldata, modeldata["un1"], "u");

    return modeldata
end

end


text:
list_test_str = "light red bags contain 1 bright white bag, 2 muted yellow bags.
dark orange bags contain 3 bright white bags, 4 muted yellow bags.
bright white bags contain 1 shiny gold bag.
muted yellow bags contain 2 shiny gold bags, 9 faded blue bags.
shiny gold bags contain 1 dark olive bag, 2 vibrant plum bags.
dark olive bags contain 3 faded blue bags, 4 dotted black bags.
vibrant plum bags contain 5 faded blue bags, 6 dotted black bags.
faded blue bags contain no other bags.
dotted black bags contain no other bags."

list_test = split(list_test_str, '\n')

list = readlines("day7_input.txt")

bag = "shiny gold"

#regxrule = r"( bags contain | bag\, | bags\, | bag\.| bags\.)"

regxrule = r"( bags contain | bags?, | bags?.)"

regxrule = r" bags contain | bags?[,.] ?"

regxrule = r" bags contain | bags?(, |\.)"

function get_outer_bag_rules(list)
    outer = Dict()
    for line in list
        parts = split(line, regxrule, keepempty=false)
            for part in parts[2:end]
                if part ≠ "no other"
                    num, bag = split(part, ' ', limit=2)
                    if bag in keys(outer)
                        push!(outer[bag], parts[1])
                    else
                        push!(outer, bag => Set([parts[1]]))
                    end
                end
            end
    end
    return outer
end

function get_outer_bags(bag, outer_bag_rules, outer_set = Set())
    if bag in keys(outer_bag_rules)
        outermost_bags = outer_bag_rules[bag]
        if outermost_bags ∉ outer_set
            union!(outer_set,outermost_bags)
            for outer_bag ∈ outermost_bags
                union!(outer_set, get_outer_bags(outer_bag, outer_bag_rules, outer_set))
            end
        end
    end
    return outer_set
end

get_number_outer_bags(bag, list) = length(get_outer_bags(bag, get_outer_bag_rules(list)))

println("Number of outer bags containing `$bag` bag in test list is $(get_number_outer_bags(bag, list_test))\n")
println("Number of outer bags containing `$bag` bag in problem list is $(get_number_outer_bags(bag, list))\n")

#=
```julia
julia> @btime get_number_outer_bags("shiny gold", list)
  5.744 ms (99019 allocations: 3.22 MiB)
185

julia> @btime get_outer_bag_rules(list);
  1.970 ms (17606 allocations: 1.25 MiB)

julia> outer_bag_rules = get_outer_bag_rules(list);

julia> @btime length(get_outer_bags(bag, outer_bag_rules));
  3.742 ms (81414 allocations: 1.97 MiB)
```
=#

# Outra versão com mais regex, mas só um pouco mais rápida.

# A versão regex_web foi testada em https://regexr.com/5hutn
regexrule_web = r"(.+)(?: bags contain )((\d+) ([^,\n]*) (bag,|bags,) )?((\d+) ([^,\n]*) (bag,|bags,) )?((\d+) ([^,\n]*) (bag,|bags,) )?((\d+) ([^,\n]*) (bag,|bags,))?((\d+|no) ([^,\n]*) ((bag\.|bags\.)))\n?"
regexrule = r"(.+)(?: bags contain )((\d+) ([^,]*) (bag,|bags,) )?((\d+) ([^,]*) (bag,|bags,) )?((\d+) ([^,]*) (bag,|bags,) )?((\d+) ([^,]*) (bag,|bags,))?((\d+|no) ([^,]*) ((bag\.|bags\.)))"
#regexrule = r"(.+)(?: bags contain )((\d+) ([^,]*) (bag,|bags,) )?{4}((\d+|no) ([^,]*) ((bag\.|bags\.)))"

function get_outer_bag_rules2(list)
    outer = Dict()
    for line in list
        m = match(regexrule, line)
            for j in 3:4:length(m.captures)
                if m[j] ∉ ("", "no") && !(m[j+1] === nothing)
                    num = m[j]
                    bag = m[j+1]
                    if bag in keys(outer)
                        push!(outer[bag], m[1])
                    else
                        push!(outer, bag => Set([m[1]]))
                    end
                end
            end
    end
    return outer
end

# part 2

list_test2_str = "shiny gold bags contain 2 dark red bags.
dark red bags contain 2 dark orange bags.
dark orange bags contain 2 dark yellow bags.
dark yellow bags contain 2 dark green bags.
dark green bags contain 2 dark blue bags.
dark blue bags contain 2 dark violet bags.
dark violet bags contain no other bags."

list_test2 = split(list_test2_str, '\n')

function get_inner_bag_rules(list)
    inner = Dict()
    for line in list
        parts = split(line, regxrule, keepempty=false)
            for part in parts[2:end]
                if part ≠ "no other"
                    num, bag = split(part, ' ', limit=2)
                    if parts[1] in keys(inner)
                        push!(inner[parts[1]], (bag, parse(Int, num)))
                    else
                        push!(inner, parts[1] => [(bag, parse(Int, num))])
                    end
                end
            end
    end
    return inner
end

function get_next_inner_bags(bag, inner_bag_rules)
    if bag in keys(inner_bag_rules)
        innermost_stuff = inner_bag_rules[bag]
        inner_sum = 1 # start with containing bag
        for (innermost_bag, innermost_num) in innermost_stuff
            # recursively count inner bags
            inner_sum += innermost_num * get_next_inner_bags(innermost_bag, inner_bag_rules)
        end
        return inner_sum
    else
        return 1
    end
end

get_number_inner_bags(bag, list) = get_next_inner_bags(bag, get_inner_bag_rules(list)) - 1 # discount starting bag

println("Number of bags contained in `$bag` bag in test list is $(get_number_inner_bags(bag, list_test))\n")
println("Number of bags contained in `$bag` bag in second test list is $(get_number_inner_bags(bag, list_test2))\n")
println("Number of bags contained in `$bag` bag in problem list is $(get_number_inner_bags(bag, list))\n")

#=
# In this case, most of the time is taken extracting the rules.
# We don't need to keep growing with the set of color bags, just the number
```julia
julia> @btime get_number_inner_bags("shiny gold", list)
  1.970 ms (17252 allocations: 1.10 MiB)
89084

julia> inner_bag_rules = get_inner_bag_rules(list);

julia> @btime get_inner_bag_rules(list);
  1.906 ms (16404 allocations: 1.07 MiB)

julia> @btime get_next_inner_bags(bag, inner_bag_rules) - 1;
  61.513 μs (848 allocations: 33.66 KiB)

```
=#


text:
###
### The advantage of the diagonal versions is improved performance for operations
### like sqrt / exp / log of diagonal elements, which is O(1) -- asymptotically irrelevant.
### Therefore, I am not planning on bothering with implementing loops over kernels for
### diagonal matrices within the forseeable future.
### For now, for medium to large sizes, packed matrices are the way to go.
### The switch here simply decides which sort to return.
###


text:
"""
From Andrej Karpathy's implementation:

- the initial stem consists of a combination of token encoding and a positional encoding
- the meat of it is a uniform sequence of Transformer blocks
    - each Transformer is a sequential combination of a 1-hidden-layer MLP block and a self-attention block
    - all blocks feed into a central residual pathway similar to resnets
- the final decoder is a linear projection into a vanilla Softmax classifier
"""

"""
for the Julia implementation
    Use Transformer.jl's Transformer.
    For some reason PositionEmbedding was not learning, defined my own
    Didn't find a simple linear layer, defined my own (TODO look for one in FLUX)
"""

import Flux: Dense
using Transformers: Transformer, PositionEmbedding, Embed
using Transformers.Basic: Positionwise, PwFFN, Dropout
using Flux:params, @functor, LayerNorm, Chain
using NNlib: gelu
using Random
using Distributions: Normal, UnivariateDistribution
using BSON: @save, @load

smallNormalInit(out, in) = Float32.(rand(Normal(0.0, 0.02), out, in))
zeroBias(x) = Float32.(zeros(x))
Dense(in::Integer, out::Integer)= Dense(in,out, initW=smallNormalInit,initb=zeroBias)
mingptTransformer(hiddenSize=512, attentionHeads=4, pdrop=.1) = Transformer(hiddenSize,attentionHeads, 4*hiddenSize, future = false, act = gelu, pdrop = 0.1)
gptEmbed(d::UnivariateDistribution, size::Int, vocab_size::Int; scale = one(Float32)) = Embed(Float32(scale), Float32.(rand(d, size, vocab_size)))

# Trying to see if a linear layer works better than dense layer for the last layer in gpt
struct Linear{S<:AbstractArray}
  W::S
end

Linear(W) = Linear(W)

function Linear(in::Integer, out::Integer)
  return Linear(smallNormalInit(out, in))
end

@functor Linear

function (a::Linear)(x::AbstractArray{T}) where T
  a.W*x
end

struct GPTPositionEmbedding{F, W <: AbstractArray{F}}
	    embedding::W
	end

@functor GPTPositionEmbedding

function GPTPositionEmbedding(size::Int, max_len::Int = 1024)
    GPTPositionEmbedding( zeros(Float32,size, max_len))
end

function (pe::GPTPositionEmbedding{F})(x::AbstractArray{F}) where F
    pe.embedding[:, 1:size(x,2)]  # just return the needed position embedding
end

struct minGPT{ EM<:Embed,
		PM<:GPTPositionEmbedding,
		C<:Chain,DP<:Dropout,
		LN<:LayerNorm,
		PW<:Positionwise }
    em::EM
    pe::PM
	drop::DP
	ln::LN
    c::C
	pw::PW
end

@functor minGPT

function (g::minGPT)(x::A) where {T, N, A<:AbstractArray{T, N}}
	x=g.em(x)
	x=g.drop(x .+g.pe(x))
	x = g.ln(x)
	x=g.c(x)
	g.pw(x)
end

function minGPT(;numberAttnLayers=2, hiddenSize=512, attentionHeads=4, vocabSize=65, maxBlockSize=50, pdrop=.1)
    embed = gptEmbed(Normal(0.0, 0.02), hiddenSize, vocabSize)
	posembd = GPTPositionEmbedding(hiddenSize, maxBlockSize)
	ln = LayerNorm(hiddenSize)
	pw = Positionwise(Dense(hiddenSize, vocabSize))   #Linear(hiddenSize, vocabSize)) #Logits.
    c=Chain(
        [mingptTransformer(hiddenSize,attentionHeads, pdrop) for i in 1:numberAttnLayers]...,
    )

    minGPT(embed,posembd, Dropout(pdrop), ln, c, pw)
end

## ---

function noDecayList(t::Transformer, objs::IdDict=IdDict())
	#MultiheadAttention is not exported from Transformers.  Get the noDecayList
	#directly.
	for o in [t.mh.iqproj.b,t.mh.ikproj.b,t.mh.ivproj.b,t.mh.oproj.b]
							get!(objs, o, 1)
	end

	noDecayList(t.mhn, objs)
	noDecayList(t.pw, objs)
	noDecayList(t.pwn, objs)
	noDecayList(t.drop, objs)
	objs
end

function noDecayList(c::Chain, objs::IdDict=IdDict())
	for i in 1:length(c.layers)
		noDecayList(c.layers[i], objs)
	end
	objs
end

function noDecayList(pe::GPTPositionEmbedding, objs::IdDict=IdDict())
	objs
end

function noDecayList(em::Embed, objs::IdDict=IdDict())

	get!(objs, em.embedding, 1)
	objs
end

function noDecayList(ln::LayerNorm, objs::IdDict=IdDict())
	get!(objs, ln.diag.α, 1)
	get!(objs, ln.diag.β, 1)
	objs
end

function noDecayList(pw::PwFFN, objs::IdDict=IdDict())
	get!(objs, pw.din.b, 1)
	get!(objs, pw.dout.b, 1)
	objs
end

function noDecayList(drop::Dropout, objs::IdDict=IdDict())
	objs
end

function noDecayList(d::Dense, objs::IdDict=IdDict())
	get!(objs, d.b, 1)
	objs
end

function noDecayList(g::minGPT, objs::IdDict=IdDict())
	noDecayList(g.em, objs)
	noDecayList(g.pe, objs)
	noDecayList(g.drop, objs)
	noDecayList(g.ln, objs)
	noDecayList(g.c, objs)
	#noDecayList(g.pw)...,
	for o in g.pw.models   # positionwise layers
		noDecayList(o, objs)
	end
	objs

end

## ---
using Zygote
using Flux
using Transformers
function functionDummy( m)
	## Don't call the above function, just put here to show how to save and restore model
	weights = params(model);
	@save "mymodel.bson" weights

	#model needs same parameters
	newModel=minGPT(maxBlockSize=lengthOfSentences, vocabSize=length(vocab))
	 @load "mymodel.bson"  weights   ## this must be the same name as the one used in the save!!!!
	Flux.loadparams!(newModel, weights)
end

function functionDummy2()
	WWW## Don't call the above function, just put here to show how to save and restore model
	##assuming mytestModel is set  flux, zygote and transformers must be included
	@save "testmodel.bson" mytestmodel

	@load "testmodel.bson"  mytestmodel   ## this must be the same name as the one used in the save!!!!
	mytestmodel(vocab(collect("Where ")))
end

max_stars_repo_name: daryoush/minGPT

max_issues_repo_name: daryoush/minGPT

max_forks_repo_name: daryoush/minGPT


text:
using Documenter, FastaIO

makedocs(
    modules  = [FastaIO],
    format   = :html,
    sitename = "FastaIO.jl",
    pages    = Any[
        "Home" => "index.md",
       ]
    )

deploydocs(
    repo   = "github.com/carlobaldassi/FastaIO.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
    julia  = "0.7"
)


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %%
using LinearAlgebra
using ForwardDiff
using Symbolics
using StaticArrays
using BenchmarkTools

# %%
function f(x)
    @show x
    @show zero(x)
    @assert x > zero(x)
    sin(x)/x
end

# %%
ForwardDiff.derivative(f, 1.0)

# %%
ForwardDiff.Dual{ForwardDiff.Tag{typeof(f), Float64}}(1.0,1.0) > 0

# %%
module My

struct Foo{T}
    a::T
    b::T
end

end

# %%
methods(My.Foo)

# %%
methods(My.Foo{Int})

# %%
My.Foo(a::T, b::T) where T<:Real = My.Foo{T}(a, b)

# %%
methods(My.Foo)

# %%
@which My.Foo(1im, 2im)

# %%
@which My.Foo(1, 2)

# %%
G(x, A) = exp(-dot(x, A, x)/2)

# %%
struct Bar{F, G}
    f::F
    ∇f::G
end

# %%
function Bar(f)
    ∇f(x, param) = ForwardDiff.gradient(x -> f(x, param), x)
    Bar(f, ∇f)
end

# %%
bar = Bar(G)

# %%
A = SA[
    2 -1
    -1 2
]

# %%
bar.f(SA[1, 2], A)

# %%
bar.∇f(SA[1, 2], A)

# %%
@variables a[axes(A)...], x[axes(SA[1, 2])...]

# %%
aa = collect(a)

# %%
collect(aa) - aa

# %%
xx = collect(x)

# %%
_G = G(xx, aa) # |> expand |> simplify

# %%
_∇G = Symbolics.gradient(_G, xx)

# %%
_∇G / _G

# %%
∇G = build_function(_∇G, xx, vec(aa); expression=Val(false))[1]

# %%
∇G(SA[1, 2], A)

# %%
function Bar(f, x, param)
    x_sym = collect(x)
    param_sym = collect(param)
    _f = simplify(expand(f(x_sym, param_sym)))
    _∇f = Symbolics.gradient(_f, x_sym)
    ff = build_function(_f, vec(x_sym), vec(param_sym); expression=Val(false))
    ∇f = build_function(_∇f, vec(x_sym), vec(param_sym); expression=Val(false))[1]
    Bar(ff, ∇f)
end

# %%
@variables x[1:2] a[1:2, 1:2]
car = Bar(G, x, a)

# %%
x = SA[1, 2]
A = SA[2 -1; -1 2]
car.f(x, A)

# %%
car.∇f(x, A)

# %%
@btime G($x, $A)
@btime $bar.f($x, $A)
@btime $car.f($x, $A)

# %%
@btime $∇G($x, $A)
@btime $bar.∇f($x, $A)
@btime $car.∇f($x, $A)

# %%
dG(x, A) = (A + A')*x/2*G(x, A)
#dG2(x, A) = A*x*G(x, A)
@btime $dG2($x, $A)

# %%
dG(xx, aa)

# %%

max_stars_repo_path: 0020/ForwardDiff.jl vs. Symbolics.jl.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0020/ForwardDiff.jl vs. Symbolics.jl.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0020/ForwardDiff.jl vs. Symbolics.jl.jl

max_forks_repo_name: genkuroki/public


text:
#
# Dimensions & Factors IO
#
# Printing and stuff

####################################################################################################
#                                   Dimensions
Base.show(io::IO, d::Dimension) =
    print(io, name(d), ":  ", _support_repr(values(d)))

_support_repr(v::AbstractVector) = repr(v) * " (" * repr(length(v)) * ")"

_support_repr(v::Range) =
    repr(first(v)) * ":" * repr(step(v)) * ":" * repr(last(v)) * " (" * repr(length(v)) * ")"

_support_repr(v::UnitRange) =
    repr(first(v)) * ":" * repr(last(v)) * " (" * repr(length(v)) * ")"

_support_repr(v::Base.OneTo) = "1:" * repr(last(v))

Base.mimewritable(::MIME"text/html", d::Dimension) = false
Base.show(io::IO, a::MIME"text/html", d::Dimension) =
        print(io, replace(replace(repr(d), "\n", "<br>"), "  ", " "))

####################################################################################################
#                                   Factors
function Base.show(io::IO, ϕ::Factor)
    if length(ϕ) > 1
        print(io, length(ϕ), " instantiations:")
        for d in ϕ.dimensions
            println(io, "")
            print(io, "\t", d)
        end
    else
        print(io, "1 instantiation: ", ϕ.potential[1])
    end
end

Base.mimewritable(::MIME"text/html", ϕ::Factor) = false
Base.show(io::IO, a::MIME"text/html", ϕ::Factor) = show(io, a, DataFrame(ϕ))


text:
@testset "Regions" begin

@testset "Point" begin
    @test SI.Point((1, 2)) isa SI.Point{Int, 2}
    @test SI.Point((1.0, 2.0, 3.0)) isa SI.Point{Float64, 3}
    @test SI.Point{Float64, 3}((1, 2, 3)) isa SI.Point{Float64, 3}
    @test isequal(SI.empty(SI.Point{Float64, 2}).coord, (NaN, NaN))
    @test isequal(SI.empty(SI.Point{Int, 1}).coord, (0,))
    @test SI.isvalid(SI.Point((1, 2)))
    @test SI.isvalid(SI.Point((1.0, 2.0, 3.0)))
    @test !SI.isvalid(SI.Point((NaN, 1.0)))

    @test_throws MethodError SI.Point((1.0, 2, 3.0))
    @test_throws InexactError SI.Point{Int,2}((1.5, 3.0))
    @test_throws MethodError SI.Point{Int,3}((1, 3))

    @test SI.area(SI.Point((1,2))) === 0
    @test SI.area(SI.Point((1.0,2.0))) === 0.0
    @test SI.perimeter(SI.Point((1,2))) === 0
    @test SI.perimeter(SI.Point((1.0,2.0))) === 0.0
    @test SI.sqrdistance(SI.Point((1,)), SI.Point((3,))) === 4
    @test SI.sqrdistance(SI.Point((1, 0)), SI.Point((3, -3))) === 13
end

@testset "Rect" begin
    @test SI.Rect((1, 2), (2, 3)) isa SI.Rect{Int, 2}
    @test SI.Rect((1.0, 2.0, 3.0), (2, 3, 4)) isa SI.Rect{Float64, 3}
    @test SI.Rect{Float64, 3}((1, 2, 3), (2, 3, 4)) isa SI.Rect{Float64, 3}
    @test isequal(SI.empty(SI.Rect{Float64, 2}), SI.Rect((NaN, NaN), (NaN, NaN)))
    @test isequal(SI.empty(SI.Rect{Int, 1}), SI.Rect((typemax(Int),), (typemax(Int),)))
    @test SI.isvalid(SI.Rect((1, 2), (2, 3)))
    @test SI.isvalid(SI.Rect((1.0, 2.0), (2.0, 3.0)))
    @test SI.isvalid(SI.Rect((1.0, 2.0), (Inf, Inf)))
    @test SI.isvalid(SI.Rect((-Inf, -Inf), (2.0, 3.0)))
    @test !SI.isvalid(SI.Rect((NaN, 2.0), (2.0, 3.0)))
    @test !SI.isvalid(SI.Rect((1.0, 2.0), (2.0, NaN)))
    @test !SI.isvalid(SI.Rect((2.5, 2.0), (2.0, 3.0)))
    @test !SI.isvalid(SI.Rect((1.0, 2.0), (2.0, -3.0)))
    @test !SI.isvalid(SI.Rect((1.0, Inf), (2.0, 3.0)))

    @test_throws MethodError SI.Rect((1.0, 2, 3.0), (2, 3.0, 4))
    @test_throws MethodError SI.Rect((1.0, 3.0), (2, 3.0, 4))
    @test_throws InexactError SI.Rect{Int,2}((1.0, 3.0), (1.5, 2.0))

    @test SI.area(SI.Rect((1.0,), (2.0,))) === 1.0
    @test SI.area(SI.Rect((1,), (0,))) === -1
    @test SI.area(SI.Rect((1.0, 2.0), (2.5, 4.0))) === 3.0

    @test SI.perimeter(SI.Rect((1.0,), (2.0,))) === 1.0
    @test SI.perimeter(SI.Rect((1,), (0,))) === -1
    @test SI.perimeter(SI.Rect((1.0, 2.0), (2.5, 4.0))) === 7.0

    @test SI.center(SI.Rect((1.0,), (2.0,))) == SI.Point((1.5,))
    @test SI.center(SI.Rect((1.0, 2.0), (2.5, 4.0))) == SI.Point((1.75, 3.0))

    ai = SI.Rect((0, 0), (1, 1))
    bi = SI.Rect((1, 1), (2, 2))
    ci = SI.Rect((2, 2), (3, 3))
    @test SI.intersect(ai, bi) == SI.Rect((1,1), (1,1))
    @test SI.intersect(ai, ci) === SI.empty(SI.Rect{Int, 2})
    @test SI.combine(ai, bi) == SI.Rect((0,0), (2,2))
    @test SI.center(SI.Rect((2, 2), (3, 4))) == SI.Point((2.5, 3.0))
    @test SI.center(SI.Rect((1, 2), (3, 4))) == SI.Point((2.0, 3.0))

    a = SI.Rect{Float64, 2}((0, 0), (1, 1))
    b = SI.Rect{Float64, 2}((1, 1), (2, 2))
    c = SI.Rect{Float64, 2}((2, 2), (3, 3))
    @test SI.intersect(a, b) == SI.Rect((1.0, 1.0), (1.0, 1.0))
    @test SI.intersect(a, c) === SI.empty(typeof(a))
    @test SI.combine(a, b) == SI.Rect((0.0,0.0), (2.0,2.0))
    @test SI.combined_area(a, b) == 4
    @test SI.combined_area(a, c) == 9
    @test SI.overlap_area(a, b) == 0
    @test SI.overlap_area(a, c) == 0
    @test SI.overlap_area(SI.combine(a, b), SI.combine(b, c)) == 1

    @test SI.touches(a, SI.combine(a, b))
    @test SI.touches(SI.combine(a, b), a)
    @test SI.touches(b, SI.combine(a, b))
    @test !SI.touches(SI.combine(a, b), c)
    @test SI.touches(SI.combine(a, c), a)
    @test SI.touches(SI.combine(a, c), c)
    @test !SI.touches(SI.combine(a, c), b)
    @test SI.touches(SI.combine(a, b), SI.Rect((1.0, 0.0), (2.0, 1.0)))
    @test SI.touches(SI.combine(a, b), SI.Rect((-1.0, 0.0), (0.0, 1.0))) #outside
    @test SI.touches(SI.combine(a, b), SI.Rect((0.0, -1.0), (1.0, 0.0))) # outside
    @test SI.touches(SI.combine(a, b), SI.Rect((2.0, 0.0), (3.0, 1.0))) # outside
    @test SI.touches(SI.combine(a, b), SI.Rect((0.0, 2.0), (1.0, 3.0))) # outside
end

end


text:
using BinaryBuilder

name = "Sundials"
version = v"5.2.1" # <-- There is no version 5.2.1, but we need to change version to build for Julia v1.7

# Collection of sources required to build Sundials
sources = [
    GitSource("https://github.com/LLNL/sundials.git",
              "b16d3d3995668c9a13c9f4bee8b0113ff6a9cf6d"),
    DirectorySource("../bundled@5"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/sundials*

# Set up CFLAGS
if [[ "${target}" == *-mingw* ]]; then
    atomic_patch -p1 $WORKSPACE/srcdir/patches/Sundials_windows.patch
    # Work around https://github.com/LLNL/sundials/issues/29
    export CFLAGS="-DBUILD_SUNDIALS_LIBRARY"
    # See https://github.com/LLNL/sundials/issues/35
    atomic_patch -p1 ../patches/Sundials_lapackband.patch
    # When looking for KLU libraries, CMake searches only for import libraries,
    # this patch ensures we look also for shared libraries.
    atomic_patch -p1 ../patches/Sundials_findklu_suffixes.patch
elif [[ "${target}" == powerpc64le-* ]]; then
    export CFLAGS="-Wl,-rpath-link,/opt/${target}/${target}/lib64"
fi

# Set up LAPACK
LAPACK_LIBRARIES="-lgfortran"
if [[ ${nbits} == 64 ]]; then
    atomic_patch -p1 $WORKSPACE/srcdir/patches/Sundials_Fortran.patch
    LAPACK_LIBRARIES="${LAPACK_LIBRARIES} ${libdir}/libopenblas64_.${dlext}"
else
    LAPACK_LIBRARIES="${LAPACK_LIBRARIES} ${libdir}/libopenblas.${dlext}"
fi
if [[ "${target}" == i686-* ]] || [[ "${target}" == x86_64-* ]]; then
    LAPACK_LIBRARIES="${LAPACK_LIBRARIES} -lquadmath"
elif [[ "${target}" == powerpc64le-* ]]; then
    LAPACK_LIBRARIES="${LAPACK_LIBRARIES} -lgomp -ldl -lm -lpthread"
fi

# Build
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
    -DEXAMPLES_ENABLE_C=OFF \
    -DKLU_ENABLE=ON \
    -DKLU_INCLUDE_DIR="${includedir}" \
    -DKLU_LIBRARY_DIR="${libdir}" \
    -DLAPACK_ENABLE=ON \
    -DLAPACK_LIBRARIES:STRING="${LAPACK_LIBRARIES}" \
    ..
make -j${nproc}
make install

# Move libraries to ${libdir} on Windows
if [[ "${target}" == *-mingw* ]]; then
    mv ${prefix}/lib/libsundials_*.${dlext} "${libdir}"
fi
"""

# We attempt to build for all defined platforms
platforms = filter!(p -> arch(p) != "powerpc64le", supported_platforms(; experimental=true))
platforms = expand_gfortran_versions(platforms)

products = [
    LibraryProduct("libsundials_arkode", :libsundials_arkode),
    LibraryProduct("libsundials_cvode", :libsundials_cvode),
    LibraryProduct("libsundials_cvodes", :libsundials_cvodes),
    LibraryProduct("libsundials_ida", :libsundials_ida),
    LibraryProduct("libsundials_idas", :libsundials_idas),
    LibraryProduct("libsundials_kinsol", :libsundials_kinsol),
    LibraryProduct("libsundials_nvecmanyvector", :libsundials_nvecmanyvector),
    LibraryProduct("libsundials_nvecserial", :libsundials_nvecserial),
    LibraryProduct("libsundials_sunlinsolband", :libsundials_sunlinsolband),
    LibraryProduct("libsundials_sunlinsoldense", :libsundials_sunlinsoldense),
    LibraryProduct("libsundials_sunlinsolklu", :libsundials_sunlinsolklu),
    LibraryProduct("libsundials_sunlinsollapackband", :libsundials_sunlinsollapackband),
    LibraryProduct("libsundials_sunlinsollapackdense", :libsundials_sunlinsollapackdense),
    LibraryProduct("libsundials_sunlinsolpcg", :libsundials_sunlinsolpcg),
    LibraryProduct("libsundials_sunlinsolspbcgs", :libsundials_sunlinsolspbcgs),
    LibraryProduct("libsundials_sunlinsolspfgmr", :libsundials_sunlinsolspfgmr),
    LibraryProduct("libsundials_sunlinsolspgmr", :libsundials_sunlinsolspgmr),
    LibraryProduct("libsundials_sunlinsolsptfqmr", :libsundials_sunlinsolsptfqmr),
    LibraryProduct("libsundials_sunmatrixband", :libsundials_sunmatrixband),
    LibraryProduct("libsundials_sunmatrixdense", :libsundials_sunmatrixdense),
    LibraryProduct("libsundials_sunmatrixsparse", :libsundials_sunmatrixsparse),
    LibraryProduct("libsundials_sunnonlinsolfixedpoint", :libsundials_sunnonlinsolfixedpoint),
    LibraryProduct("libsundials_sunnonlinsolnewton", :libsundials_sunnonlinsolnewton),
]

dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("OpenBLAS_jll"),
    Dependency("SuiteSparse_jll"; compat="~5.10.1"),
]

# Build the tarballs.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"6", julia_compat="1.7")

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module lap

using TimerOutputs

export pflap, single_lap,
       sparse_pflap, single_lap_at_points

include("lap_algs/helpers.jl")
include("lap_algs/single.jl")
include("lap_algs/multi.jl")

end # module


text:
using Base.Cartesian
export chfield

############# ancillas ################
export @fieldview
"""
    @fieldview fname(x::TYPE) = x.fieldname
    @fieldview fname(x::TYPE) = x[i]
    ...

Create a function fieldview that can be accessed by a reversible program

```jldoctest; setup=:(using NiLangCore)
julia> struct GVar{T, GT}
           x::T
           g::GT
       end

julia> @fieldview xx(x::GVar) = x.x

julia> chfield(GVar(1.0, 0.0), xx, 2.0)
GVar{Float64, Float64}(2.0, 0.0)
```
"""
macro fieldview(ex)
    @smatch ex begin
        :($f($obj::$tp) = begin $line; $ex end) => begin
            xval = gensym("value")
            esc(Expr(:block,
                :(Base.@__doc__ $f($obj::$tp) = begin $line; $ex end),
                :($NiLangCore.chfield($obj::$tp, ::typeof($f), $xval) = $(Expr(:block, assign_ex(ex, xval, false), obj)))
            ))
        end
        _ => error("expect expression `f(obj::type) = obj.prop`, got $ex")
    end
end

chfield(a, b, c) = error("chfield($a, $b, $c) not defined!")
chfield(x, ::typeof(identity), xval) = xval
chfield(x::T, ::typeof(-), y::T) where T = -y
chfield(x::T, ::typeof(adjoint), y) where T = adjoint(y)

############ dataview patches ############
export tget, subarray

"""
    tget(i::Int)

Get the i-th entry of a tuple.
"""
tget(i::Int) = x::Tuple -> x[i]

"""
    subarray(ranges...)

Get a subarray, same as `view` in Base.
"""
subarray(args...) = x -> view(x, args...)


text:
module SpecialFuncs

using PyCall
using SymPy

import Base: gamma, polygamma, beta,
             airyai, airybi,
             besseli, besselj, besselk, bessely

for meth in (
             :jacobi,
             :gegenbauer,
             :chebyshevt, :chebyshevu,
             :legendre, :assoc_legendre,
             :hermite,
             :laguerre, :assoc_laguerre,
             :Ynm, :Ynm_c,
             :hankel1, :hankel2,
             :jn, :yn
             )
    meth_name = string(meth)
    @eval begin
        @doc """
            `$($meth_name)`: a SymPy function.
                The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
                """ ->
        ($meth)(xs...;kwargs...) = sympy_meth($meth_name, xs...; kwargs...)
    end
    eval(Expr(:export, meth))
end

## these have (parameter, x) signature. Use symbolic x to call sympy version, othewise
## should dispatch to julia version.
for fn in (:besselj, :bessely, :besseli, :besselk)
    meth = string(fn)
#    eval(Expr(:import, :Base, fn))
    @eval ($fn)(nu::Number, x::Sym; kwargs...) = sympy_meth($meth, nu, x; kwargs...)
    @eval ($fn)(nu::Number, a::Array{Sym}) = map(x ->$fn(nu, x), a)
end

for meth in (:fresnels, :fresnelc, :Ei, :Si, :Ci,
             :airyai, :airybi
           )
    meth_name = string(meth)

    @eval begin
        @doc """
`$($meth_name)`: a SymPy function.
The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
""" ->
        ($meth)(x::Sym;kwargs...) = sympy_meth($meth_name, x; kwargs...)
    end
    @eval ($meth)(a::Array{Sym}) = map($meth, a)
    eval(Expr(:export, meth))
end

## :gamma, :beta... # need import
beta(x::Sym, y::Sym) = sympy_meth(:beta, x, y)
gamma(x::Sym) = sympy_meth(:gamma, x)
polygamma(nu::Sym, x::Sym) = sympy_meth(:polygamma, nu, x)

## Hyper and friends don't really have symbolic use...
"""

Evaluates the generalized hypergeometric function:
[hyper](http://docs.sympy.org/dev/modules/mpmath/functions/hypergeometric.html#hyper)

"""
hyper{T<:Number, S<:Number}(as::Vector{T}, bs::Vector{S}, z::Number) = sympy_meth(:hyper, as, bs, z)
export hyper

"""

[docs](http://docs.sympy.org/dev/modules/mpmath/functions/hypergeometric.html#meijer-g-function)

"""
function meijerg{T<:Number, S<:Number}(a1s::Vector{T}, a2s::Vector{T}, b1s::Vector{S}, b2s::Vector{S}, z::Number, r=1;kwargs...)
    sympy_meth(:meijerg, a1s, a2s, b1s, b2s, z;  kwargs...)
end
export meijerg

end


text:
using Random

Random.seed!(123)

# smooth case

f1 = SqrNormL2()
f2 = Translate(SqrNormL2(2.5), randn(10))
f = Sum(f1, f2)

predicates_test(f)

@test ProximalOperators.is_quadratic(f) == true
@test ProximalOperators.is_strongly_convex(f) == true
@test ProximalOperators.is_set(f) == false

xtest = randn(10)

result = f1(xtest) + f2(xtest)
@test f(xtest) ≈ result

grad1, val1 = gradient(f1, xtest)
grad2, val2 = gradient(f2, xtest)

gradsum = randn(size(xtest))
valsum = gradient!(gradsum, f, xtest)
@test gradsum ≈ grad1 + grad2

# nonsmooth case

g1 = NormL2()
g2 = Translate(SqrNormL2(2.5), randn(10))
g = Sum(g1, g2)

predicates_test(g)

@test ProximalOperators.is_smooth(g) == false
@test ProximalOperators.is_strongly_convex(g) == true
@test ProximalOperators.is_set(g) == false

xtest = randn(10)

result = g1(xtest) + g2(xtest)
@test g(xtest) ≈ result

grad1, val1 = gradient(g1, xtest)
grad2, val2 = gradient(g2, xtest)

gradsum = randn(size(xtest))
valsum = gradient!(gradsum, g, xtest)
@test gradsum ≈ grad1 + grad2


text:
function construct( e::HomographyMatrix,
                   𝐊₁::AbstractArray,
                   𝐑₁::AbstractArray,
                   𝐭₁::AbstractArray,
                   𝐊₂::AbstractArray,
                   𝐑₂::AbstractArray,
                   𝐭₂::AbstractArray,
                   𝐧::AbstractArray,
                   d::Real)

    if size(𝐊₁) != (3,3) || size(𝐊₂) != (3,3) ||
       size(𝐑₁) != (3,3) || size(𝐑₂) != (3,3)
        throw(ArgumentError("Expect 3 x 3 calibration and rotation matrices."))
    end
    if length(𝐭₁) != 3 || length(𝐭₂) != 3
        throw(ArgumentError("Expect length-3 translation vectors."))
    end
    if length(𝐧) != 3
        throw(ArgumentError("Expect length-3 normal vector."))
    end
    # TODO Check that camera center does not lie on the plane.
    𝐈 = SMatrix{3,3}(1.0I)
    𝐇 = 𝐊₂*𝐑₂*(𝐈 - (𝐧'*𝐭₁ + d)^-1 * (𝐭₁ - 𝐭₂) * 𝐧')/𝐑₁/𝐊₁
    SMatrix{3,3,Float64,3*3}(𝐇)
end


text:
# -----------------------------------------------------

function get_source_destiny_weight{T}(mat::AbstractArray{T,2})
    nrow, ncol = size(mat)    # rows are sources and columns are destinies

    nosymmetric = !issymmetric(mat) # plots only triu for symmetric matrices
    nosparse = !issparse(mat) # doesn't plot zeros from a sparse matrix

    L = length(mat)

    source  = Array(Int, L)
    destiny = Array(Int, L)
    weights  = Array(T, L)

    idx = 1
    for i in 1:nrow, j in 1:ncol
        value = mat[i, j]
        if !isnan(value) && ( nosparse || value != zero(T) ) # TODO: deal with Nullable

            if i < j
                source[idx]  = i
                destiny[idx] = j
                weights[idx]  = value
                idx += 1
            elseif nosymmetric && (i > j)
                source[idx]  = i
                destiny[idx] = j
                weights[idx]  = value
                idx += 1
            end

        end
    end

    resize!(source, idx-1), resize!(destiny, idx-1), resize!(weights, idx-1)
end

function get_source_destiny_weight(source::AVec, destiny::AVec)
    if length(source) != length(destiny)
        throw(ArgumentError("Source and destiny must have the same length."))
    end
    source, destiny, Float64[ 1.0 for i in source ]
end

function get_source_destiny_weight(source::AVec, destiny::AVec, weights::AVec)
    if !(length(source) == length(destiny) == length(weights))
        throw(ArgumentError("Source, destiny and weights must have the same length."))
    end
    source, destiny, weights
end

# -----------------------------------------------------

function get_adjacency_matrix(source::AVec{Int}, destiny::AVec{Int}, weights::AVec)
    n = max(maximum(source), maximum(destiny))
    full(sparse(source, destiny, weights, n, n))
end

function make_symmetric(A::AMat)
    A = copy(A)
    for i=1:size(A,1), j=i+1:size(A,2)
        A[i,j] = A[j,i] = A[i,j]+A[j,i]
    end
    A
end

function compute_laplacian(adjmat::AMat, node_weights::AVec)
    n, m = size(adjmat)
    # @show size(adjmat), size(node_weights)
    @assert n == m == length(node_weights)

    # scale the edge values by the product of node_weights, so that "heavier" nodes also form
    # stronger connections
    adjmat = adjmat .* sqrt(node_weights * node_weights')

    # D is a diagonal matrix with the degrees (total weights for that node) on the diagonal
    deg = vec(sum(adjmat,1)) - diag(adjmat)
    D = diagm(deg)

    # Laplacian (L = D - adjmat)
    L = Float64[i == j ? deg[i] : -adjmat[i,j] for i=1:n,j=1:n]

    L, D
end

# -----------------------------------------------------
# -----------------------------------------------------

# see: http://www.research.att.com/export/sites/att_labs/groups/infovis/res/legacy_papers/DBLP-journals-camwa-Koren05.pdf
# also: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.3.2055&rep=rep1&type=pdf

# this recipe uses the technique of Spectral Graph Drawing, which is an
# under-appreciated method of graph layouts; easier, simpler, and faster
# than the more common spring-based methods.
function spectral_graph(adjmat::AMat; node_weights::AVec = ones(size(adjmat,1)), kw...)
    adjmat = make_symmetric(adjmat)
    L, D = compute_laplacian(adjmat, node_weights)

    # get the matrix of eigenvectors
    v = eig(L, D)[2]

    # x, y, and z are the 2nd through 4th eigenvectors of the solution to the
    # generalized eigenvalue problem Lv = λDv
    vec(v[2,:]), vec(v[3,:]), vec(v[4,:])
end

function spectral_graph(source::AVec{Int}, destiny::AVec{Int}, weights::AVec; kw...)
    spectral_graph(get_adjacency_matrix(source, destiny, weights); kw...)
end

# -----------------------------------------------------

# Axis-by-Axis Stress Minimization -- Yehuda Koren and David Harel
# See: http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.437.3177&rep=rep1&type=pdf

# # NOTES:
# #   - dᵢⱼ = the "graph-theoretical distance between nodes i and j"
# #         = Aᵢⱼ
# #   - kᵢⱼ = dᵢⱼ⁻²
# #   - b̃ᵢ = ∑ᵢ≠ⱼ ((x̃ⱼ ≤ x̃ᵢ ? 1 : -1) / dᵢⱼ)
# #   - need to solve for x each iteration: Lx = b̃

# # Solve for one axis at a time while holding the others constant.
# # dims is 2 (2D) or 3 (3D).  free_dims is a vector of the dimensions to update (for example if you fix y and solve for x)
# function by_axis_stress_graph(adjmat::AMat, node_weights::AVec = ones(size(adjmat,1));
#                               dims = 2, free_dims = 1:dims,
#                               x = rand(length(node_weights)),
#                               y = rand(length(node_weights)),
#                               z = rand(length(node_weights)))
#     adjmat = make_symmetric(adjmat)
#     L, D = compute_laplacian(adjmat, node_weights)

#     n = length(node_weights)
#     maxiter = 100 # TODO: something else

#     @assert dims == 2

#     @show adjmat L

#     for _ in 1:maxiter
#         x̃ = x
#         b̃ = Float64[sum(Float64[(i==j || adjmat[i,j] == 0) ? 0.0 : ((x̃[j] <= x̃[i] ? 1.0 : -1.0) / adjmat[i,j]) for j=1:n]) for i=1:n]
#         @show x̃ b̃
#         x = L \ b̃

#         xdiff = x - x̃
#         @show norm(xdiff)
#         if norm(xdiff) < 1e-4
#             info("converged. norm(xdiff) = $(norm(xdiff))")
#             break
#         end
#     end
#     @show x y
#     x, y, z
# end

norm_ij(X, i, j) = sqrt(sum(Float64[(v[i]-v[j])^2 for v in X]))
stress(X, dist, w, i, j) = w[i,j] * (norm_ij(X, i, j) - dist[i,j])^2
function stress(X, dist, w)
    tot = 0.0
    for i=1:size(X,1), j=1:i-1
        tot += stress(X, dist, w, i, j)
    end
    tot
end

# follows section 2.3 from http://link.springer.com/chapter/10.1007%2F978-3-540-31843-9_25#page-1
# Localized optimization, updates: x
function by_axis_local_stress_graph(adjmat::AMat;
                              node_weights::AVec = ones(size(adjmat,1)),
                              dim = 2, free_dims = 1:dim,
                              x = rand(length(node_weights)),
                              y = rand(length(node_weights)),
                              z = rand(length(node_weights)),
                              maxiter = 300,
                              kw...)
    adjmat = make_symmetric(adjmat)
    n = length(node_weights)

    # graph-theoretical distance between node i and j (i.e. shortest path distance)
    # TODO: calculate a real distance
    dist = estimate_distance(adjmat)
    # @show dist

    # also known as kᵢⱼ in "axis-by-axis stress minimization".  the -2 could also be 0 or -1?
    w = dist .^ -2

    # in each iteration, we update one dimension/node at a time, reducing the total stress with each update
    X = dim == 2 ? (x, y) : (x, y, z)
    laststress = stress(X, dist, w)
    for k in 1:maxiter
        for p in free_dims
            for i=1:n
                numer, denom = 0.0, 0.0
                for j=1:n
                    i==j && continue
                    numer += w[i,j] * (X[p][j] + dist[i,j] * (X[p][i] - X[p][j]) / norm_ij(X, i, j))
                    denom += w[i,j]
                end
                if denom != 0
                    X[p][i] = numer / denom
                end
            end
        end

        # check for convergence of the total stress
        thisstress = stress(X, dist, w)
        if abs(thisstress - laststress) / abs(laststress) < 1e-5
            # info("converged. numiter=$k last=$laststress this=$thisstress")
            break
        end
        laststress = thisstress
    end

    dim == 2 ? (X..., nothing) : X
end

function by_axis_local_stress_graph(source::AVec{Int}, destiny::AVec{Int}, weights::AVec; kw...)
    by_axis_local_stress_graph(get_adjacency_matrix(source, destiny, weights); kw...)
end

# -----------------------------------------------------

if Plots.is_installed("LightGraphs")
    @eval begin
        import LightGraphs

        # TODO: so much wasteful conversion... do better
        function estimate_distance(adjmat::AMat)
            source, destiny, weights = get_source_destiny_weight(sparse(adjmat))
            n = size(adjmat,1)
            g = LightGraphs.Graph(n)
            for (si,di,wi) in zip(source,destiny,weights)
                LightGraphs.add_edge!(g, si, di)
            end
            convert(Matrix{Float64}, hcat(map(i->LightGraphs.dijkstra_shortest_paths(g, i).dists, 1:n)...))
        end
    end
else
    @eval function estimate_distance(adjmat::AMat)
        warn("Install LightGraphs for the best layout calculations.")
        map(a -> a==0 ? 5.0 : a, adjmat)
    end
end

function tree_graph(adjmat::AMat; kw...)
    tree_graph(get_source_destiny_weight(adjmat)...; kw...)
end

function tree_graph(source::AVec{Int}, destiny::AVec{Int}, weights::AVec;
                    node_weights::AVec = ones(max(maximum(source), maximum(destiny))),
                    root::Symbol = :top,  # flow of tree: left, right, top, bottom
                    layers_scalar = 1.0,
                    layers = nothing,
                    dim = 2,
                    kw...)
    extrakw = KW(kw)
    n = length(node_weights)

    # @show root

    # TODO: compute layers, which get bigger as you go away from the root
    if layers == nothing
        # layers = rand(1:4, n)
        layers = compute_tree_layers2(source, destiny, n)
    end

    # reverse direction?
    if root in (:top, :right)
        layers = -layers
    end

    # add noise
    layers = layers + 0.6rand(size(layers)...)

    # TODO: normalize layers somehow so it's in line with distances
    layers .*= layers_scalar
    if dim == 2
        if root in (:top, :bottom)
            extrakw[:y] = layers
            extrakw[:free_dims] = [1]
        elseif root in (:left, :right)
            extrakw[:x] = layers
            extrakw[:free_dims] = [2]
        else
            error("unknown root: $root")
        end
    else
        error("3d not supported")
    end

    # now that we've fixed one dimension, let the stress algo solve for the other(s)
    by_axis_local_stress_graph(get_adjacency_matrix(source, destiny, weights);
                               node_weights=node_weights,
                               dim=dim,
                               extrakw...)
end

function adjlist_and_degrees(source, destiny, n)
    # build a list of children (adjacency list)
    alist = Vector{Int}[Int[] for i=1:n]
    indeg, outdeg = zeros(Int, n), zeros(Int, n)
    for (si,di) in zip(source, destiny)
        push!(alist[si], di)
        indeg[di] += 1
        outdeg[si] += 1
    end
    alist, indeg, outdeg
end

function compute_tree_layers(source, destiny, n)
    alist, indeg, outdeg = adjlist_and_degrees(source, destiny, n)

    # choose root to be the node with lots going out, but few coming in
    netdeg = outdeg - 50indeg
    idxs = sortperm(netdeg, rev=true)
    # rootidx = findmax(netdeg)
    # @show outdeg indeg netdeg idxs alist
    placed = Int[]

    layers = zeros(n)
    for i=1:n
        idx = shift!(idxs)

        # first, place this after its parents
        for j in placed
            if idx in alist[j]
                layers[idx] = max(layers[idx], layers[j] + 1)
            end
        end

        # next, shift its children lower
        for j in idxs
            if j in alist[idx]
                layers[j] = max(layers[j], layers[idx] + 1)
            end
        end

        push!(placed, idx)
    end
    layers
end

# an alternative algo to pick tree layers... generate a list of roots,
# and for each root, make a pass through the tree (without recurrency)
# and push the children below their parents
function compute_tree_layers2(source, destiny, n)
    alist, indeg, outdeg = adjlist_and_degrees(source, destiny, n)
    roots = filter(i->indeg[i]==0, 1:n)
    if isempty(roots)
        roots = [1]
    end

    layers = zeros(Int,n)
    for i in roots
        shift_children!(layers, alist, Int[], i)
    end

    # now that we've shifted children out, move parents closer to their closest children
    while true
        shifted = false
        for parent=1:n
            if !(isempty(alist[parent]))
                minidx = minimum(layers[child] for child in alist[parent])
                if layers[parent] < minidx - 1
                    shifted = true
                    layers[parent] = minidx - 1
                end
            end
        end
        shifted || break
    end

    layers
end

function shift_children!(layers, alist, placed, parent)
    for idx in alist[parent]
        if !(idx in placed) && layers[idx] <= layers[parent]
            layers[idx] = layers[parent] + 1
        end
    end
    for idx in alist[parent]
        if idx != parent && !(idx in placed)
            push!(placed, idx)
            shift_children!(layers, alist, placed, idx)
        end
    end
end

# -----------------------------------------------------

# TODO: maybe also implement Catmull-Rom Splines? http://www.mvps.org/directx/articles/catmull/

# -----------------------------------------------------

# we want to randomly pick a point to be the center control point of a bezier
# curve, which is both equidistant between the endpoints and normally distributed
# around the midpoint
function random_control_point(xi, xj, yi, yj, curvature_scalar)
    xmid = 0.5 * (xi+xj)
    ymid = 0.5 * (yi+yj)

    # get the angle of y relative to x
    theta = atan((yj-yi) / (xj-xi)) + 0.5pi

    # calc random shift relative to dist between x and y
    dist = sqrt((xj-xi)^2 + (yj-yi)^2)
    dist_from_mid = curvature_scalar * (rand()-0.5) * dist

    # now we have polar coords, we can compute the position, adding to the midpoint
    (xmid + dist_from_mid * cos(theta),
     ymid + dist_from_mid * sin(theta))
end

# -----------------------------------------------------

const _graph_funcs = KW(
    :spectral => spectral_graph,
    :stress => by_axis_local_stress_graph,
    :tree => tree_graph,
)

# a graphplot takes in either an (N x N) adjacency matrix
#   note: you may want to pass node weights to markersize or marker_z
# A graph has N nodes where adj_mat[i,j] is the strength of edge i --> j.  (adj_mat[i,j]==0 implies no edge)

# NOTE: this is for undirected graphs... adjmat should be symmetric and non-negative

@userplot GraphPlot

@recipe function f(g::GraphPlot;
                   dim = 2,
                   free_dims = 1:dim,
                   T = Float64,
                   curves = true,
                   curvature_scalar = 0.2,
                   root = :top,
                   node_weights = nothing,
                   names = [],
                   x = nothing,
                   y = nothing,
                   z = nothing,
                   func = spectral_graph,
                   shorten = 0.1
                  )
    @assert dim in (2, 3)
    _3d = dim == 3

    source, destiny, weights = get_source_destiny_weight(g.args...)
    if !(eltype(source) <: Integer)
        names = unique(sort(vcat(source,destiny)))
        source = Int[findfirst(names, si) for si in source]
        destiny = Int[findfirst(names, di) for di in destiny]
    end
    n = max(maximum(source), maximum(destiny))

    if node_weights == nothing
        node_weights = ones(n)
    end
    @assert length(node_weights) == n

    # do we want to compute coordinates?
    if (_3d && (x == nothing || y == nothing || z == nothing)) || (!_3d && (x == nothing || y == nothing))
        if isa(func, Symbol)
            func = _graph_funcs[func]
        end
        x, y, z = func(
            source, destiny, weights;
            node_weights = node_weights,
            dim = dim,
            free_dims = free_dims,
            root = root
        )
    end

    # create a series for the line segments
    if get(d, :linewidth, 1) > 0
        @series begin
            xseg, yseg, zseg = Segments(), Segments(), Segments()
            for (si, di, wi) in zip(source, destiny, weights)
                # add a line segment
                xsi, ysi, xdi, ydi = shorten_segment(x[si], y[si], x[di], y[di], shorten)
                # ysi, ydi = shorten_segment(y[si], y[di])
                if curves
                    xpt, ypt = random_control_point(xsi, xdi,
                                                    ysi, ydi,
                                                    curvature_scalar)
                    push!(xseg, xsi, xpt, xdi)
                    push!(yseg, ysi, ypt, ydi)
                    _3d && push!(zseg, z[si], z[si], z[di])
                else
                    push!(xseg, xsi, xdi)
                    push!(yseg, ysi, ydi)
                    _3d && push!(zseg, z[si], z[di])
                end
            end

            # generate a list of colors, one per segment
            grad = get(d, :linecolor, nothing)
            if isa(grad, ColorGradient)
                line_z := weights
            end

            seriestype := (curves ? :curves : (_3d ? :path3d : :path))
            series_annotations := []
            linewidth --> 1
            markershape := :none
            markercolor := :black
            primary := false
            _3d ? (xseg.pts, yseg.pts, zseg.pts) : (xseg.pts, yseg.pts)
        end
    end

    xlims --> extrema_plus_buffer(x)
    ylims --> extrema_plus_buffer(y)
    if _3d
        zlims --> extrema_plus_buffer(z)
    end

    seriestype := (_3d ? :scatter3d : :scatter)
    linewidth := 0
    linealpha := 0
    foreground_color_border --> nothing
    grid --> false
    legend --> false
    ticks --> nothing
    series_annotations --> map(string,names)
    markersize --> 10 + 100node_weights / sum(node_weights)
    _3d ? (x, y, z) : (x, y)
end

function extrema_plus_buffer(v)
    vmin,vmax = extrema(v)
    vdiff = vmax-vmin
    buffer = vdiff * 0.1
    vmin - buffer, vmax + buffer
end

function shorten_segment(x1, y1, x2, y2, shorten)
    xshort = shorten * (x2-x1)
    yshort = shorten * (y2-y1)
    x1+xshort, y1+yshort, x2-xshort, y2-yshort
end

# ---------------------------------------------------------------------------
# Arc/Chord Diagrams

function arcvertices{T}(source::AVec{T}, destiny::AVec{T})
    values = unique(vcat(source, destiny))
    [(i, i) for i in values ]
end

function arcvertices{T<:Union{Char,Symbol,AbstractString}}(source::AVec{T},
                                                                   destiny::AVec{T})
    lab2x = Dict{T,Int}()
    n = 1
    for element in vcat(source, destiny)
        if !haskey(lab2x, element)
            lab2x[element] = n
            n += 1
        end
    end
    lab2x
end

@userplot ArcDiagram

@recipe function f(h::ArcDiagram)

    source, destiny, weights = get_source_destiny_weight(h.args...)

    vertices = arcvertices(source, destiny)

    # Box setup
    legend --> false
    aspect_ratio --> :equal
    grid --> false
    foreground_color_axis --> nothing
    foreground_color_border --> nothing
    ticks --> nothing

    usegradient = length(unique(weights)) != 1

    if usegradient
        colorgradient = ColorGradient(get(d,:linecolor,cgrad()))
        wmin,wmax = extrema(weights)
    end

    for (i, j, value) in zip(source,destiny,weights)
        @series begin

            xi = vertices[i]
            xj = vertices[j]

            if usegradient
                linecolor --> colorgradient[(value-wmin)/(wmax-wmin)]
            end

            legend --> false
            label := ""
            primary := false

            r  = (xj - xi) / 2
            x₀ = (xi + xj) / 2
            θ = linspace(0,π,30)
            x₀ + r * cos(θ), r * sin(θ)
        end
    end

    @series begin
        if eltype(keys(vertices)) <: Union{Char, AbstractString, Symbol}
            series_annotations --> collect(keys(vertices))
            markersize --> 0
            y = -0.1
        else
            y =  0.0
        end
        seriestype --> :scatter
        collect(values(vertices)) , Float64[ y for i in vertices ]
    end
end

# =================================================
# Arc and chord diagrams

# ---------------------------------------------------------------------------
# Chord diagram

function arcshape(θ1, θ2)
    Plots.shape_coords(Shape(vcat(
        Plots.partialcircle(θ1, θ2, 15, 1.1),
        reverse(Plots.partialcircle(θ1, θ2, 15, 0.9))
    )))
end

# """
# `chorddiagram(source, destiny, weights[, grad, zcolor, group])`

# Plots a chord diagram, form `source` to `destiny`,
# using `weights` to determine the edge colors using `grad`.
# `zcolor` or `group` can be used to determine the node colors.
# """

@userplot ChordDiagram

@recipe function f(h::ChordDiagram)
    source, destiny, weights = get_source_destiny_weight(h.args...)

    xlims := (-1.2,1.2)
    ylims := (-1.2,1.2)
    legend := false
    grid := false
    xticks := nothing
    yticks := nothing

    nodemin, nodemax = extrema(vcat(source, destiny))
    weightmin, weightmax = extrema(weights)

    A  = 1.5π # Filled space
    B  = 0.5π # White space (empirical)

    Δα = A / nodemax
    Δβ = B / nodemax
    δ = Δα  + Δβ

    grad = cgrad(get(d,:linecolor,:inferno), get(d, :linealpha, nothing))

    for i in 1:length(source)
        # TODO: this could be a shape with varying width
        @series begin
            seriestype := :curves
            x := [cos((source[i ]-1)*δ + 0.5Δα), 0.0, cos((destiny[i]-1)*δ + 0.5Δα)]
            y := [sin((source[i ]-1)*δ + 0.5Δα), 0.0, sin((destiny[i]-1)*δ + 0.5Δα)]
            linecolor := grad[(weights[i] - weightmin) / (weightmax - weightmin)]
            primary := false
            ()
        end
    end

    for n in 0:(nodemax-1)
        sx, sy = arcshape(n*δ, n*δ + Δα)
        @series begin
            seriestype := :shape
            x := sx
            y := sy
            ()
        end
    end
end


text:
#test passed

@testset "loreg_research.jl" begin

    data_surv = DataFrame(CSV.File("C:\\Users\\Alexandrine\\cso_raw.csv"))
    data_trudeau = DataFrame(CSV.File("C:\\Users\\Alexandrine\\Montreal-Trudeau_15min_precipitations.csv"))
    data = innerjoin(data_surv, data_trudeau, on=:Date)
    dropmissing!(data, :Duration)
    addColOverflow!(data, "Duration", "Surverse", 3)

    topredict = "Surverse"

    train, test = train_test_year(data, 2020)

    result = CR_logreg(topredict, [:d15min, :d45min, :d1h, :d24h], train, test)
    varnames = result[1, :Model]
    F1 = result[1, :F₁]
    variables = stepF1_logreg(topredict, [:d15min, :d45min, :d1h, :d24h], [], train, test, false)

    #The F1 score is well derived
    @test F1 == logreg("Surverse", varnames, train, test)

    #The df is well sorted
    @test F1 == maximum(result.F₁)

    #The higher F1 score of the CR_logreg corresponds to the higher found stepwise
    @test F1 == logreg("Surverse", variables, train, test)

end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
function increment(x)
    return x + 1
end

max_stars_repo_name: kallooa-bhr/dockerfiles

max_issues_repo_name: kallooa-bhr/dockerfiles

max_forks_repo_name: kallooa-bhr/dockerfiles