text:
using Test
using OceanTurbulenceParameterEstimation
using Oceananigans
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ColumnEnsembleSize
using Oceananigans.TurbulenceClosures: ConvectiveAdjustmentVerticalDiffusivity
using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity, MixingLength, SurfaceTKEFlux
using Suppressor: @suppress

using OceanTurbulenceParameterEstimation.TurbulenceClosureParameters: closure_with_parameters, update_closure_ensemble_member!

const CAVD = ConvectiveAdjustmentVerticalDiffusivity
const CATKE = CATKEVerticalDiffusivity

@testset "TurbulenceClosureParameters tests" begin
    @info "Testing closure_with_parameters on ConvectiveAdjustmentVerticalDiffusivity"
    closures = [CAVD(background_κz=1.0) CAVD(background_κz=1.1)
                CAVD(background_κz=1.2) CAVD(background_κz=1.3)
                closure_with_parameters(CAVD(), (background_κz=1.4,)) closure_with_parameters(CAVD(), (background_κz=1.5,))]

    @test closures[3, 1].background_κz == 1.4
    @test closures[3, 2].background_κz == 1.5

    closure_tuple = closure_with_parameters((CAVD(), CAVD()), (background_κz = 14.0,))

    for closure in closure_tuple
        @test closure.background_κz == 14.0
    end

    @info "Testing closure_with_parameters on CATKEVerticalDiffusivity"
    catke_closure = @suppress CATKE()
    new_Cᴷuʳ = 1e3
    new_catke_closure = closure_with_parameters(catke_closure, (; Cᴷuʳ = new_Cᴷuʳ))
    @test new_catke_closure.mixing_length.Cᴷuʳ == new_Cᴷuʳ

    @info "Testing update_closure_ensemble_member! on ConvectiveAdjustmentVerticalDiffusivity"
    new_background_κz, new_convective_νz = 100.0, 10.0
    ensemble_member_no = 2
    update_closure_ensemble_member!(closures, ensemble_member_no, (background_κz = new_background_κz, convective_νz = new_convective_νz))

    for j in 1:2
        @test closures[ensemble_member_no, j].background_κz == new_background_κz
        @test closures[ensemble_member_no, j].convective_νz == new_convective_νz
    end

    @info "Testing update_closure_ensemble_member! on CATKEVerticalDiffusivity"
    closures = @suppress [CATKE() CATKE()
                          CATKE() CATKE()
                          CATKE() CATKE()]

    new_CᴷRiʷ = 100.0
    ensemble_member_no = 2
    update_closure_ensemble_member!(closures, ensemble_member_no, (; CᴷRiʷ = new_CᴷRiʷ))

    for j in 1:2
        @test closures[ensemble_member_no, j].mixing_length.CᴷRiʷ == new_CᴷRiʷ
    end
end


text:
"""
Transfer ITensor networks into networks with type tensortype
"""
function abstract_network!(tensortype, network::Vector{ITensor}, dict::Dict; kwargs...)
  for tensor in network
    if !haskey(dict, tensor)
      dict[tensor] = tensortype(tensor; kwargs...)
    end
  end
  return Vector{tensortype}([dict[t] for t in network])
end

function abstract_network!(tensortype, tree::SubNetwork, dict::Dict; kwargs...)
  for node in tree.inputs
    if !haskey(dict, node)
      if node isa ITensor
        dict[node] = tensortype(node; kwargs...)
      else
        dict[node] = abstract_network!(tensortype, node, dict; kwargs...)
      end
    end
  end
  return SubNetwork([dict[n] for n in tree.inputs])
end

function abstract_network(tensortype, networks::Vector{Vector{ITensor}}, vars; kwargs...)
  dict = Dict{ITensor,tensortype}()
  if length(vars) != 0
    vars = abstract_network!(tensortype, collect(vars), dict; kwargs...)
  end
  networks = Vector{Vector{tensortype}}([
    abstract_network!(tensortype, n, dict; kwargs...) for n in networks
  ])
  return networks, Tuple(vars)
end

function ChainRulesCore.rrule(
  ::typeof(abstract_network), tensortype, networks::Vector{Vector{ITensor}}, vars; kwargs...
)
  function pullback(v)
    d_networks, dvars = v[1], v[2]
    dvars_itensor = map(x -> ITensor(x), dvars)
    return (NoTangent(), NoTangent(), NoTangent(), dvars_itensor)
  end
  return abstract_network(tensortype, networks, vars; kwargs...), pullback
end

function abstract_network(tensortype, trees::Vector{SubNetwork}, vars; kwargs...)
  dict = Dict()
  if length(vars) != 0
    vars = abstract_network!(tensortype, collect(vars), dict; kwargs...)
  end
  trees = Vector{SubNetwork}([
    abstract_network!(tensortype, t, dict; kwargs...) for t in trees
  ])
  return trees, Tuple(vars)
end

function ChainRulesCore.rrule(
  ::typeof(abstract_network), tensortype, trees::Vector{SubNetwork}, vars; kwargs...
)
  function pullback(v)
    d_trees, dvars = v[1], v[2]
    dvars_itensor = map(x -> ITensor(x), dvars)
    return (NoTangent(), NoTangent(), NoTangent(), dvars_itensor)
  end
  return abstract_network(tensortype, trees, vars; kwargs...), pullback
end

ext: jl


text:
using Test, Random, LinearAlgebra, SparseArrays

using Yao
using QuAlgorithmZoo

@testset "RotBasis" begin
    rt = RotBasis(0.5, 0.4)
    crt = chain(rt)

    dispatch!(crt, [2., 3.])
    @test nparameters(crt) == 2

    for (t1, t2, t3) in zip(parameters(rt), (2, 3), parameters(crt))
        @test t1 == t2 == t3
    end

    # check consistency
    rb = roll(1, RotBasis(0.1, 0.3))#rot_basis(1)
    angles = randpolar(1)
    # prepair a state in the angles direction.
    psi = angles |> polar2u |> register

    # rotate to the same direction for measurements.
    dispatch!(rb, vec(angles))
    @test state(apply!(psi, rb)) ≈ [1, 0]

    @test nparameters(rot_basis(3)) == 6
    dispatch!(rb, :zero)
    @test parameters(rb)[1] == 0
    dispatch!(rb, :random)
    @test parameters(rb)[1] != 0
end

@testset "polar and u" begin
    polar = randpolar(10)
    @test size(polar) == (2, 10)
    @test polar |> polar2u |> u2polar ≈ polar
end


text:
using Test
using Statistics
using CompositeConvex
using OptimPackNextGen

@testset "Robust L1 regression test" begin

    # Generate problem data
    m,n = 1000,100
    num_noise = 10
    A = randn(m,n)
    x = randn(n)
    b = A*x
    e = zeros(m)
    e[rand(1:m,num_noise)]  = 1000*randn(num_noise)
    b_noisy = b + e

    # Choose norm for residuals
    residual_norm = pNorm{1}()

    # Construct operators for the problem
    c = (y,x)->affine_model!(y,A,x,b)
    ∇c_adjoint = (y,x,r)->mul!(y,A',r)
    gauge_projection = (y,x,t)->compute_norm_projection!(y,x,t,residual_norm)
    dual_gauge = x->compute_norm(x,dual_norm(residual_norm))
    params = ProblemParams(c,∇c_adjoint,gauge_projection,dual_gauge)
    info = CompositeConvex.MinComposite.IntermediateVars(zeros(m),zeros(m),0.0)

    # Function evaluated at the noiseless solution should be zero
    g = similar(x)
    f = CompositeConvex.MinComposite.value_function_objective!(x, g, params, info)
    @test f < 1e-10
    @test norm(g) < 1e-10

    f = CompositeConvex.MinComposite.value_function_objective!(zeros(n), g, params, info)
    @test norm(g - (-A'*b)) < 1e-10

    # Gradient test
    info = CompositeConvex.MinComposite.IntermediateVars(zeros(m),zeros(m),1.0)
    x0 = randn(n)
    obj! = (x,g)->CompositeConvex.MinComposite.value_function_objective!(x,g,params,info)
    f0 = obj!(x0,g)
    dx = 1e-2*randn(n)
    df = dot(g,dx)
    h = 10 .^ (-6.0:0.0)
    e1 = zeros(length(h))
    for (i,hi) in enumerate(h)
        fi = obj!(x0+hi*dx,g)
        e1[i] = abs(fi - f0 - hi*df)
    end
    taylor_exponent = median(diff(log10.(e1)))
    @test abs(taylor_exponent - 2) < 0.1

    # Recover the solution from noisy data, knowing the true τ value
    c = (y,x)->affine_model!(y,A,x,b_noisy)
    params = ProblemParams(c,∇c_adjoint,gauge_projection,dual_gauge)
    info = CompositeConvex.MinComposite.IntermediateVars(zeros(m),zeros(m),norm(e,1))

    obj! = (x,g)->CompositeConvex.MinComposite.value_function_objective!(x,g,params,info)
    x1 = vmlmb(obj!, zeros(n), mem=10, ftol=(0,1e-3), maxiter=20, maxeval=100)
    @test norm(x1-x) < 1e-3*norm(x)

    # Ensure the value function gives the correct value + gradient at the true τ value
    subproblem_opts = SubproblemOptParams{Float64}()
    f,df = CompositeConvex.MinComposite.value_function!(x, params, info, subproblem_opts)
    @test f < 1e-8
    @test df < 1e-8

    # Validate derivative of the value function vs the finite different gradient
    info.tau = norm(e,1)/2
    f,df = CompositeConvex.MinComposite.value_function!(zeros(n), params, info, subproblem_opts)
    h = 1e-6*info.tau
    info.tau += h
    f1,df1 = CompositeConvex.MinComposite.value_function!(zeros(n), params, info, subproblem_opts)
    @test abs((f1-f)/h - df) < 1e-2*abs(df)

    # Minimum working example - test interface
    residual_norm = pNorm{1}()
    gauge_projection = (y,x,t)->compute_norm_projection!(y,x,t,residual_norm)
    dual_gauge = x->compute_norm(x,dual_norm(residual_norm))
    subproblem_opts = SubproblemOptParams{Float64}()
    x1 = min_composite(A,b_noisy,gauge_projection,subproblem_opts,dual_gauge,verbose=false)
    @test norm(x1 - x) < 1e-2*norm(x)
end


text:
using Distributions
using Turing
using Turing: gradient_r, invlink, link, getval, realpart
using ForwardDiff
using ForwardDiff: Dual
using Test

# Define model
@model ad_test() = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0,sqrt(s))
  1.5 ~ Normal(m, sqrt(s))
  2.0  ~ Normal(m, sqrt(s))
  return s, m
end
# Turing.TURING[:modelex]
# Call Turing's AD
# The result out is the gradient information on R
ad_test_f = ad_test()
vi = ad_test_f(Turing.VarInfo(), nothing)
svn = collect(Iterators.filter(vn -> vn.sym == :s, keys(vi)))[1]
mvn = collect(Iterators.filter(vn -> vn.sym == :m, keys(vi)))[1]
_s = realpart(getval(vi, svn)[1])
_m = realpart(getval(vi, mvn)[1])

x = map(_->Float64(_), vi[nothing])
∇E = gradient_r(x, vi, ad_test_f)
# println(vi.vns)
# println(∇E)
grad_Turing = sort(∇E)

dist_s = InverseGamma(2,3)

# Hand-written logp
function logp(x::Vector)
  s = x[2]
  # s = invlink(dist_s, s)
  m = x[1]
  lik_dist = Normal(m, sqrt.(s))
  lp = logpdf(dist_s, s, false) + logpdf(Normal(0,sqrt.(s)), m, false)
  lp += logpdf(lik_dist, 1.5) + logpdf(lik_dist, 2.0)
  lp
end

# Call ForwardDiff's AD
g = x -> ForwardDiff.gradient(logp, x);
# _s = link(dist_s, _s)
_x = [_m, _s]
grad_FWAD = sort(-g(_x))

# Compare result
println(grad_Turing)
println(grad_FWAD)
@test grad_Turing ≈ grad_FWAD atol=1e-9


text:
using Turing
using Turing: Sampler
using Test

alg = SGLD(1000, 0.25)
sampler = Sampler(alg)
@test isa(alg, SGLD)
@test isa(sampler, Sampler{<:Turing.SGLD})

alg = SGLD(200, 0.25, :m)
sampler = Sampler(alg)
@test isa(alg, SGLD)
@test isa(sampler, Sampler{<:Turing.SGLD})

alg = SGLD(1000, 0.25, :s)
sampler = Sampler(alg)
@test isa(alg, SGLD)
@test isa(sampler, Sampler{<:Turing.SGLD})


text:
##=================####==============##
## TEST SIMULATION ABM NETWORK-EPIDEMICS
##=================####==============##

##=================####==============##

using StatsBase     # For sample function
using SparseArrays  # For Adyacency Matrix
using Distributions # For probability distributions
using DelimitedFiles # For reading/writing files
using LightGraphs

using Plots         # plotting
using LaTeXStrings  # latex labels in plots

##=================####==============##

include("NetABM.jl") # Include file with functions and definitions

##=================####==============##
# PATH TO SAVE IMGS (MODIFY WITH YOUR OWN!)

img_out_path = joinpath(homedir(), "Code", "Collaborations", "networkEpidemiology_MX", "julia_ABM", "figures")

##=================####==============##
## SYSTEM'S PARAMETERS

params = Params()

params.num_agents    = 2^16 # Number of agents (32768)
params.p_infected_t0 = 0.01 # Initial fraction of infected agents
params.attack_rate   = 0.25 # Virus attack rate
params.repetition    = 1 # Repetition ID for ensemble statistics

params.μ_cop_agents  = 0.15 # Mean of cooperation probability
params.σ_cop_agents  = 0.1 # Std deviation of cooperation probability
params.μ_recovT_agents = 7 # Mean recovery time

# Link probability (Erdos-Renyi)
# params.p_link = ((1+5exp10(-2))*log(params.num_agents))/params.num_agents

##=================####==============##
# INITIALIZE ADJACENCY MATRIX (LFR Model)

the_root = joinpath(repo_path, "src", "binary_networks")

println("BUILDING NETWORK...")

the_net, true_com = get_network(
    the_root,
    "julia_ABM",
    N    = params.num_agents,
    t2   = 0.25,
    minc = 100,
    maxc = 1000);

g, Nodes = lectura_uw(the_net)

adj_mat = adjacency_matrix(g)
rows = rowvals(adj_mat)

println("NETWORK DONE")

##=================####==============##
# INITIALIZE MEETS, COOPERATION AND RECOVERY TIMES DISTRIBUTION

rt_μ = params.μ_recovT_agents # MEAN RECOVERY TIME
pc_μ = params.μ_cop_agents # MEAN COOPERATION PROBABILITY
pc_σ = params.σ_cop_agents # STD_DEV COOPERATION PROBABILITY

# meets_dist = Exponential()
meets_dist  = Geometric()
coop_dist   = truncated(Normal(pc_μ, pc_σ), 0.0, 1.0)
recovt_dist = truncated(Poisson(rt_μ), 7, 35)

##=================####==============##

# DEFINE AGENTS POPULATION (DEFAULT VALUES)
agents = [Agent(id) for id in 1:params.num_agents]

# MODIFY AGENT'S CHARACTERISTICS ACCORDING TO ITS DEMOGRAPHICS
initialize_demographics!(agents, params, coop_dist, meets_dist, recovt_dist)

# ASSGINS CONTACTS TO agent FROM adj_mat
for ag in agents
    assign_contacts!(ag, agents, adj_mat, rows)
end

##=================####==============##

time_lockdown = 15

tau = []
num_S = []
num_I = []
num_R = []

# INITIALIZE TIME AND POPULATIONS
params.now_t = 0
populations = get_populations(agents, params)

while populations[2] > 0.0 # WHILE THERE ARE INFECTED AGENTS

    # RECORD DATA
    push!(tau, params.now_t)
    push!(num_S, populations[1])
    push!(num_I, populations[2])
    push!(num_R, populations[3])

    # EACH AGENTS DECIDES TO STAY AT HOME/ NOT TO STAY AT HOME
    if params.now_t >= time_lockdown
        set_coop_agents!(agents, params)
    end

    # FINDS NEXT STATE AND UPDATES ALL AGENTS
    update_all_agents!(agents, params)

    # GET CURRENT POPULATIONS AND UPDATE TIME
    global populations = get_populations(agents, params);
    params.now_t += 1
end

println("END")

##=================####==============##
# POPULATIONS GRAPH
#
plt = plot(
        xlabel = L"t",
        ylabel = L"N(s)/N",
        guidefontsize = 16,
        tickfontsize = 10
    )

plot!(plt, tau, num_S, color = :blue,  lab = "S", lw = 1.25)
plot!(plt, tau, num_I, color = :red,   lab = "I", lw = 1.25)
plot!(plt, tau, num_R, color = :green, lab = "R", lw = 1.25)
# plot!(plt, tau, 1.0 .- num_R, color = :green, lab = "R", lw = 1.25)

vline!(plt, [time_lockdown], ls = :dash, color = :red, lab = "")

display(plt)

##=================####==============##
savefig(plt, joinpath(img_out_path, "SIR_N_$(params.num_agents)_AR_$(params.attack_rate).png"))
##=================####==============##
# Distribución de Números de Encuentros

plt = plot(
        xlabel = L"N_m",
        ylabel = L"N (N_m)",
        title = "Encounters per agent per time step distribution",
        guidefontsize = 16,
        tickfontsize = 10
    )

vals = [ag.num_meets for ag in agents]
histogram!(plt, vals, lab = false, bar_width = 0.85)

xticks!(plt, 1:maximum(vals))
display(plt)

##=================####==============##
savefig(plt, joinpath(img_out_path, "meetsDist_N_$(params.num_agents).png"))
##=================####==============##
# Distribución de Tiempo de Recuperación

plt = plot(
        xlabel = L"R_t",
        ylabel = L"N(R_t)",
        title = "Agent's recovery times distribution",
        guidefontsize = 16,
        tickfontsize = 10
    )

vals = [ag.recovery_t for ag in agents]
histogram!(plt, vals, lab = false)

# histogram!(plt, vals, normalized = true, lab = false)

# xticks!(plt, unique(vals))
display(plt)

##=================####==============##
savefig(plt, joinpath(img_out_path, "recovtDist_N_$(params.num_agents).png"))
##=================####==============##
# Distribución de Probabilidad de Coperación

plt = plot(
        xlabel = L"p_C",
        ylabel = L"N(p_C)",
        title = "Agent's cooperation probability distribution",
        guidefontsize = 16,
        tickfontsize = 10
    )

vals = [ag.p_cop for ag in agents]
histogram!(plt, vals, lab = false)

display(plt)

##=================####==============##
savefig(plt, joinpath(img_out_path, "cooProb_N_$(params.num_agents).png"))
##=================####==============##
# Tiempo de Infección

plt = plot(
        xlabel = L"t_I",
        ylabel = L"p(t_I)",
        title = "Time of infection distribution",
        guidefontsize = 16,
        tickfontsize = 10
    )

histogram!(plt, [ag.infection_t for ag in agents], normalized = false, lab = false)
vline!(plt, [time_lockdown], color = :red, ls = :dash, lw = 1.25, lab = "lockdown")
display(plt)

##=================####==============##
savefig(plt, joinpath(img_out_path, "tiempoI_N_$(params.num_agents)_AR_$(params.attack_rate).png"))
##=================####==============##
# Distribución de Grado

plt = plot(
        xlabel = L"\kappa_i",
        ylabel = L"p\;(\kappa_i)",
        title = "Network's degree distribution ",
        guidefontsize = 16,
        tickfontsize = 10
    )

histogram!(plt, [ag.degree_t for ag in agents], normalized = true, lab = false)
display(plt)

##=================####==============##
savefig(plt, joinpath(img_out_path, "degDist_N_$(params.num_agents)_ER.png"))
##=================####==============##


text:
export RawBLAS

module RawBLAS
using Base.LinAlg # force built-in BLAS library initialization
using ..Mocha.blasfunc
using Compat

for (gemm, elty) in ((:dgemm_, Float64), (:sgemm_, Float32))
  @eval begin
    function gemm!(transA::Char, transB::Char, M::Int, N::Int, K::Int, alpha::$elty,
        A, lda, B, ldb, beta::$elty, C, ldc)
      ccall(($(blasfunc(gemm)), Base.libblas_name), Void,
          (Ptr{UInt8}, Ptr{UInt8}, Ptr{Base.LinAlg.BlasInt}, Ptr{Base.LinAlg.BlasInt},
          Ptr{Base.LinAlg.BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{Base.LinAlg.BlasInt},
          Ptr{$elty}, Ptr{Base.LinAlg.BlasInt}, Ptr{$elty}, Ptr{$elty},
          Ptr{Base.LinAlg.BlasInt}),
          &transA, &transB, &M, &N, &K, &alpha, A, &lda, B, &ldb, &beta, C, &ldc)
    end

    function gemm!(transA::Char, transB::Char, M::Int, N::Int, K::Int, alpha::$elty,
        A, B, beta::$elty, C)
      lda = (transA == 'N') ? M : K;
      ldb = (transB == 'N') ? K : N;
      ldc = M;
      gemm!(transA, transB, M, N, K, alpha, A, lda, B, ldb, beta, C, ldc)
    end
  end # eval
end

for (gemv, elty) in ((:dgemv_, Float64), (:sgemv_, Float32))
  @eval begin
    function gemv!(trans::Char, M::Int, N::Int, alpha::$elty, A, lda, x, incx, beta::$elty, y, incy)
      ccall(($(blasfunc(gemv)), Base.libblas_name), Void,
          (Ptr{UInt8}, Ptr{Base.LinAlg.BlasInt}, Ptr{Base.LinAlg.BlasInt}, Ptr{$elty},
          Ptr{$elty}, Ptr{Base.LinAlg.BlasInt}, Ptr{$elty}, Ptr{Base.LinAlg.BlasInt},
          Ptr{$elty}, Ptr{$elty}, Ptr{Base.LinAlg.BlasInt}),
          &trans, &M, &N, &alpha, A, &lda, x, &incx, &beta, y, &incy)
    end
    function gemv!(trans::Char, M::Int, N::Int, alpha::$elty, A, x, beta::$elty, y)
      lda = M
      incx = 1
      incy = 1
      gemv!(trans, M, N, alpha, A, lda, x, incx, beta, y, incy)
    end
  end
end

end # module RawBLAS


text:
"""
    KnotVector{T<:Number}

Vector to define knots used for [`BSplineCurve`](@ref) and [`BSplineSurface`](@ref).
"""
struct KnotVector{T<:Number}
    knots::Array{T,1}

    function KnotVector{T}(a::Array{T,1}) where {T<:Number}
        temp = a[1]
        for i in 2:lastindex(a) # ensure knot vector is non-decreasing
            @assert a[i] >= temp
            temp = a[i]
        end
        return new(copy(a))
    end
end
export KnotVector

Base.copy(knots::KnotVector{S}) where {S} = KnotVector{S}(knots.knots)

urange(knots::KnotVector) = (knots.knots[1], knots.knots[end])
numknots(knots::KnotVector) = size(knots.knots)[1]

struct InvalidParameterValue <: Exception
    invalid_u::Any
    validrange::Any
end

Base.showerror(io::IO, e::InvalidParameterValue) = print(io, "invalid parameter value: $(e.invalid_u) allowable parameter range: $(e.validrange)")

function findspan(knots::KnotVector, curveorder, u)
    pmin = curveorder + 1
    nknots = numknots(knots)
    pmax = nknots - (curveorder + 1)
    if knots.knots[pmax + 1] == u
        return pmax
    end

    for i in pmin:pmax # if u == ui+1 then the span goes from i to i+1.
        if u <= knots.knots[i + 1]
            return i
        end
    end
    # this might be too extreme. round off error could cause parameter value to go a little above max value. Might need to loosen this later.
    throw(InvalidParameterValue(u, (knots.knots[1], knots.knots[end]))) # return the bad u value and the legal range of u values
end

function basisfunctions(knots::KnotVector{T}, u::T, curveorder) where {T<:Number}
    # computes values of p+1 basis functions where p is the curve order. knotsegment is the index of the knot vector segment, u is the value of the curve parameter
    knotsegment = findspan(knots, curveorder, u)
    U = knots.knots
    tempsize = curveorder + 1
    left = zeros(T, tempsize) # later rewrite this code to have a struct that stores these temps so they aren't being allocated on every curve evaluation
    right = zeros(T, tempsize)
    N = zeros(T, tempsize)

    N[1] = one(T)

    for j in 1:curveorder
        left[j + 1] = u - U[knotsegment + 1 - j]
        right[j + 1] = U[knotsegment + j] - u
        saved = zero(T)
        for r in 0:(j - 1)
            leftindex = j - r + 1
            temp = N[r + 1] / (right[r + 2] + left[leftindex])
            N[r + 1] = saved + right[r + 2] * temp
            saved = left[leftindex] * temp
        end
        N[j + 1] = saved
    end
    return N
end


text:
# This file was generated, do not modify it. # hide
rpt = report(mach_pipe)
keys(rpt.evo_tree_classifier)

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using OhMyQSIM
using Compat.Test
using TestSetExtensions

@testset "All the tests" begin
    @includetests ARGS
end


text:
function InitialPleiades(T)

    q=[parse(BigFloat,"3.0"), parse(BigFloat,"3.0"),parse(BigFloat,"-1.0"),
       parse(BigFloat,"-3.0"),parse(BigFloat,"2.0"),parse(BigFloat,"-2.0"),
       parse(BigFloat,"2.0"), parse(BigFloat,"3.0"),parse(BigFloat,"-3.0"),
       parse(BigFloat,"2.0"),parse(BigFloat,"0.0"),parse(BigFloat,"0.0"),
       parse(BigFloat,"-4.0"),parse(BigFloat,"4.0")]

    v=[parse(BigFloat,"0.0"), parse(BigFloat,"0.0"),parse(BigFloat,"0.0"),
       parse(BigFloat,"0.0"),parse(BigFloat,"0.0"),parse(BigFloat,"1.75"),
       parse(BigFloat,"-1.5"), parse(BigFloat,"0.0"),parse(BigFloat,"0.0"),
       parse(BigFloat,"0.0"),parse(BigFloat,"-1.25"),parse(BigFloat,"1.0"),
       parse(BigFloat,"0.0"),parse(BigFloat,"0.0")]

    u0 = Array{T}(undef,1)

    u0=convert.(T,[q;v])

    return u0

end


text:
module TableViewAdvanced

using ..TableView, WebIO, Observables, JSExpr

export  Table,
        table_scope, table_scope_id, table_data, table_options, table_prepare_options, table_user_data, set_table_user_data,
        JS_loader_node,
        tva_example

const EVENTS_DEBUG = false

"""
    mutable struct Table

The main `Table` structure representing the data and scope of the table.
"""
mutable struct Table
    _table_data                 # the data source - currently a DataFrame only
    _scope                      # the table WebIO scope
    _options                    # the options for th ag-grid
    _prepare_options            # options for the Table object
    _user_data                  # aditional data that will pass between the different events and call back function - can be anything the user choose,
                                # such as the front end window in the Blink case or the JS context that is needed to run JS functions.

    function Table(table_data)
        return new(table_data, nothing, nothing, nothing)
    end
end

table_scope(t::Table) = t._scope
table_scope_id(t::Table) = t._scope.id
table_data(t::Table) = t._table_data
table_options(t::Table) = t._options
table_prepare_options(t::Table) = t._prepare_options

table_user_data(t::Table) = t._user_data
set_table_user_data(t::Table, data) = t._user_data = data

function Base.show(io::IO, obj::Table)
        # compact = get(io, :compact, false)

        println(io, obj._scope !== nothing ? "Scope is Initialized" : "Scope is nothing")
        Base.show(io, obj._table_data)
end

"""
    function prepare_table_setup_events!(table::Table)

Go over all the user's supplied events and register them.
This function is internal and should not be called directly.
"""
function prepare_table_setup_events!(table::Table)
    options = table_options(table)
    # scope = table_scope(table)

    prepare_options = options[:PrepareTableOptions]
    if (haskey(prepare_options, :events))
        for (event, func) in prepare_options[:events]
            prepare_table_setup_event!(string(event), table)
        end
    end
end

"""
    prepare_table_mutator!(options::Dict{Symbol, Any}, scope::WebIO.Scope)

Tegister a specific event by creating an observable if necessary and creating a message with certain parameters taken from the JS side of the grid.
The basic message contains the following parameters:
    row         - row number, starting at 1 for Julia
    col         - column name - the field property of the grid column (TODO: need to cdecide what to pass for calculated rows)
    oldValue    - this property exist in the cellValueChange event
    newValue    - this property exist in the cellValueChange event
    rowData     - a Dict object containing all the fields and their values in the specific row
    node        - additional information about the state of the grid row  - currently, just the indication if the row is selected.
                  this indication is important for the rowSelected event, which is triggered twice, for selecting and unselecting the previous row
This function is internal and should not be called directly.
"""
function prepare_table_setup_event!(event, table::Table)
    options = table_options(table)
    scope = table_scope(table)

    prepare_options = options[:PrepareTableOptions]

    # table_data = table_data(table)

    func = get(prepare_options[:events], Symbol(event), nothing)
    if (func !== nothing)
        if (EVENTS_DEBUG)
            @info "Registering Grid Event [$event]"
        end
        obs = scope[event]
        on(obs) do msg
            func(msg=msg, table=table)
        end

        js_func = @js function (ev)
            if ($EVENTS_DEBUG)
                console.group("Grid Event [" + $event + "]");
                console.log("Event Data: ")
                console.log(ev)
            end
            @var x = Dict()
            if ev.rowIndex !== undefined
                x["row"] = ev.rowIndex + 1
            end
            if ev.column !== undefined && ev.column.colDef !== undefined
                x["col"] = ev.column.colDef.field
            end
            if ev.oldValue !== undefined
                x["oldValue"] = ev.oldValue
            end
            if ev.newValue !== undefined
                x["newValue"] = ev.newValue
            end
            if ev.data !== undefined
                x["row_data"] = ev.data
            end
            if ev.node !== undefined
                x["node"] = Dict(
                    "selected" => ev.node.selected,
                )
            end
            $obs[] = x
            if ($EVENTS_DEBUG)
                console.groupEnd();
            end

            return
        end
        options[Symbol("on$event")] = js_func
    end
end

"""
    prepare_table_mutator!(options::Dict{Symbol, Any}, scope::WebIO.Scope)

This function is being called by the original `TableView.showtable` function. it allows us to change the grid options and
table scope before creating the JS grid object, This function is where we do most of our modifications over the original TableView process.
This function is internal and should not be called directly.
"""
function prepare_table_mutator!(options::Dict{Symbol, Any}, scope::WebIO.Scope)
    #@info keys(options)

    # try to include java script from a local file
    import!(scope, joinpath(dirname(@__FILE__), "TableViewAdvanced.js"))
    # pushfirst!(scope.imports, Asset(joinpath(dirname(@__FILE__), "TableViewAdvanced.js")))

    col_defs = options[:columnDefs]
    prepare_options = options[:PrepareTableOptions]

    table = prepare_options[:Table]
    table._scope = scope
    table._options = options

    # table_data = table_data(table)
    columns_options = get(prepare_options, :columnDefs, Dict())

    options[:rowSelection] = "single"
    options[:suppressRowDeselection] = false

    options[:defaultColDef] = get(prepare_options, :defaultColDef, Dict())

    # disable sorting on all rows and remove the row id
    for c in col_defs
        # skip the internal row field
        if (c[:field] == "__row__")
            continue
        end

        # c[:sortable] = false
        # remove certain default definition that the original TableView added per field
        for key in [:sortable, :resizable, :headerTooltip, :editable]
            if (haskey(c, key))
                delete!(c, key)
            end
        end

        # hide the ID field - the GUID field that identify each row
        # if (c[:field] == "ID")
        #     c[:hide] = true
        #     @show c
        # end

        column_options = get(columns_options, Symbol(c[:field]), Dict())

        for (col_opt_key, col_opt_value) in column_options
            c[col_opt_key] = col_opt_value
        end

        # c[:headerName] = get(column_options, :Title, c[:headerName])
        # c[:headerTooltip] = get(column_options, :TitleTooltip, c[:headerTooltip])
        # c[:editable] = get(column_options, :Editable, true)
        # c[:resizable] = get(column_options, :Resizable, c[:resizable])
        # c[:sortable] = get(column_options, :Sortable, c[:sortable])

        editor = get(column_options, :CellEditor, nothing)
        renderer = get(column_options, :CellRenderer, nothing)
        if (editor == "ComboBox")
            c[:cellEditor] = "agSelectCellEditor"
            c[:cellEditorParams] = Dict( :values => get(column_options, :CellEditorParams, []))
        elseif (editor !== nothing)
            c[:cellEditor] = editor
        end

        if (renderer !== nothing)
            c[:cellRenderer] = renderer
        end

    end

    # register the user's grid events
    prepare_table_setup_events!(table)
    # remove the prepare table information from the options Dict - was just used to pass the information to the mutator
    delete!(options, :PrepareTableOptions)

end

"""
    prepare_table(table_data, table_options::Dict{Symbol, Any})

Return a `Table` object that contains a reference to the table data and a scope for displaying the provided `table`.
This function is internal and should not be called directly.
"""
function prepare_table(table_data, table_options::Dict{Symbol, Any})
    # create a copy of the options and add the dataframe to pass to the mutator
    table_options = copy(table_options)
    # table_options[:TableData] = table_data

    # scope will be updated in the mutator function
    table = Table(table_data)
    table._prepare_options = table_options
    table_options[:Table] = table

    # these are the options we are sending to the original table view package in order to handle the data our way in the option mutator
    options = Dict{Symbol, Any}(
        :PrepareTableOptions => table_options
    )

    TableView.showtable(
        table_data,
        options = options,
        option_mutator! = prepare_table_mutator!,
    )

    return table
end

"""
    Table(data, table_options::Dict{Symbol, Any})

Return a `Table` object that contains a reference to the table data and a scope for displaying the provided `table`.
"""
function Table(data, table_options::Dict{Symbol, Any})
    table = prepare_table(data, table_options)
    return table
end

"""
    JS_loader_node(filenames...)

Return a `WebIO.Node` object that contains all the supplied files as children nodes. Recognized extentions files such
as '.js' and '.css' are contains in a :script and :style nodes respectivly.
"""
function JS_loader_node(filenames...)
    @info filenames
    nodes = []
    for filename in filenames
        ext = lowercase(splitext(filename)[2])
        @info ext
        data = read(filename)
        data = String(UInt8.(data))   # conver to a string from vector of bytes
        if (ext == ".js")
            push!(nodes, Node(:script, data))
        elseif (ext == ".css")
            push!(nodes, Node(:style, data))
        else
            push!(nodes, Node(:dome, data))
        end
    end

    return Node(:dom, nodes...)
end

function tva_example(example_name)
    if (splitext(example_name)[2] == "")
        example_name = "$(example_name).jl"
    end
    filename = abspath(joinpath(dirname(Base.find_package("TableView")), "..", "examples", example_name))
    return filename;
end

end # module TableViewAdvanced


text: norm(A) / sqrt(length(A))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
    parsegtf!(gencode::DataFrame)

Parse `gencode` by extracting `gene_id`, `gene_name`, `gene_type`, `transcript_id`,
`transcript_support_level` information from the `info` column.
"""
function parsegtf!(gencode::DataFrame)
    if !("info" in names(gencode))
        error("Cannot find :info column in the annotation.")
    end
    for col in [:gene_id, :gene_name, :gene_type]
        gencode[:, col] =
            [getindex(m.captures, 1) for m in match.(Regex("$(col) \"(.*?)\";"), gencode.info)]
    end
    for col in [:transcript_id, :transcript_support_level]
        storage = Vector(undef, size(gencode, 1))
        for (n, m) in enumerate(match.(Regex("$(col) \"(.*?)\";"), gencode.info))
            if isnothing(m) && col == :transcript_id
                storage[n] = gencode.gene_id[n]
            elseif isnothing(m) && col == :transcript_support_level
                storage[n] = missing
            else
                storage[n] = getindex(m.captures, 1)
            end
        end
        gencode[:, col] = storage
    end
    gencode.gene_id = [getindex(i, 1) for i in split.(gencode.gene_id, ".")]
    gencode.transcript_id = [getindex(i, 1) for i in split.(gencode.transcript_id, ".")]
    gencode.seqnames = replace.(gencode.seqnames, "chr" => "")
    return
end

"""
    findgene(gene::AbstractString, gencode::DataFrame)

Find chromosome, gene start, and gene stop sites for the `gene` of interest.
"""
function findgene(gene::AbstractString, gencode::DataFrame)
    ind = findfirst(isequal(gene), gencode.gene_name)
    if isnothing(ind)
        error("Cannot find $(gene) in the annotation.")
    else
        return gencode.seqnames[ind], gencode.start[ind], gencode[ind, :end]
    end
end


text:
export Agent

using Flux
using BSON
using JLD
using Setfield

"""
    Agent(;kwargs...)

One of the most commonly used [`AbstractAgent`](@ref).

Generally speaking, it does nothing but update the trajectory and policy appropriately in different stages.

# Keywords & Fields

- `policy`::[`AbstractPolicy`](@ref): the policy to use
- `trajectory`::[`AbstractTrajectory`](@ref): used to store transitions between an agent and an environment
- `role=:DEFAULT_PLAYER`: used to distinguish different agents
"""
Base.@kwdef mutable struct Agent{P<:AbstractPolicy,T<:AbstractTrajectory,R} <: AbstractAgent
    policy::P
    trajectory::T
    role::R = :DEFAULT_PLAYER
    is_training::Bool = true
end

# avoid polluting trajectory
(agent::Agent)(obs) = agent.policy(obs)

Flux.functor(x::Agent) = (policy = x.policy,), y -> @set x.policy = y.policy

function save(dir::String, agent::Agent; is_save_trajectory = true)
    mkpath(dir)
    @info "saving agent to $dir ..."

    t = @elapsed begin
        save(joinpath(dir, "policy.bson"), agent.policy)
        if is_save_trajectory
            JLD.save(joinpath(dir, "trajectory.jld"), "trajectory", agent.trajectory)
        else
            @warn "trajectory is skipped since you set `is_save_trajectory` to false"
        end
        BSON.bson(
            joinpath(dir, "agent_meta.bson"),
            Dict(
                :role => agent.role,
                :is_training => agent.is_training,
                :policy_type => typeof(agent.policy),
            ),
        )
    end

    @info "finished saving agent in $t seconds"
end

function load(dir::String, ::Type{<:Agent})
    @info "loading agent from $dir"
    BSON.@load joinpath(dir, "agent_meta.bson") role is_training policy_type
    policy = load(joinpath(dir, "policy.bson"), policy_type)
    JLD.@load joinpath(dir, "trajectory.jld") trajectory
    Agent(policy, trajectory, role, is_training)
end

get_role(agent::Agent) = agent.role

function Flux.testmode!(agent::Agent, mode = true)
    agent.is_training = !mode
    testmode!(agent.policy, mode)
end

(agent::Agent)(stage::AbstractStage, obs) =
    agent.is_training ? agent(Training(stage), obs) : agent(Testing(stage), obs)

(agent::Agent)(::Testing, obs) = nothing
(agent::Agent)(::Testing{PreActStage}, obs) = agent.policy(obs)

#####
# EpisodicCompactSARTSATrajectory
#####
function (agent::Agent{<:AbstractPolicy,<:EpisodicCompactSARTSATrajectory})(
    ::Training{PreEpisodeStage},
    obs,
)
    empty!(agent.trajectory)
    nothing
end

function (agent::Agent{<:AbstractPolicy,<:EpisodicCompactSARTSATrajectory})(
    ::Training{PreActStage},
    obs,
)
    action = agent.policy(obs)
    push!(agent.trajectory; state = get_state(obs), action = action)
    update!(agent.policy, agent.trajectory)
    action
end

function (agent::Agent{<:AbstractPolicy,<:EpisodicCompactSARTSATrajectory})(
    ::Training{PostActStage},
    obs,
)
    push!(agent.trajectory; reward = get_reward(obs), terminal = get_terminal(obs))
    nothing
end

function (agent::Agent{<:AbstractPolicy,<:EpisodicCompactSARTSATrajectory})(
    ::Training{PostEpisodeStage},
    obs,
)
    action = agent.policy(obs)
    push!(agent.trajectory; state = get_state(obs), action = action)
    update!(agent.policy, agent.trajectory)
    action
end

#####
# Union{CircularCompactSARTSATrajectory, CircularCompactPSARTSATrajectory}
#####

function (
    agent::Agent{
        <:AbstractPolicy,
        <:Union{CircularCompactSARTSATrajectory,CircularCompactPSARTSATrajectory},
    }
)(
    ::Training{PreEpisodeStage},
    obs,
)
    if length(agent.trajectory) > 0
        pop!(agent.trajectory, :state, :action)
    end
    nothing
end

function (
    agent::Agent{
        <:AbstractPolicy,
        <:Union{CircularCompactSARTSATrajectory,CircularCompactPSARTSATrajectory},
    }
)(
    ::Training{PreActStage},
    obs,
)
    action = agent.policy(obs)
    push!(agent.trajectory; state = get_state(obs), action = action)
    update!(agent.policy, agent.trajectory)
    action
end

function (
    agent::Agent{
        <:AbstractPolicy,
        <:Union{CircularCompactSARTSATrajectory,CircularCompactPSARTSATrajectory},
    }
)(
    ::Training{PostActStage},
    obs,
)
    push!(agent.trajectory; reward = get_reward(obs), terminal = get_terminal(obs))
    nothing
end

function (
    agent::Agent{
        <:AbstractPolicy,
        <:Union{CircularCompactSARTSATrajectory,CircularCompactPSARTSATrajectory},
    }
)(
    ::Training{PostEpisodeStage},
    obs,
)
    action = agent.policy(obs)
    push!(agent.trajectory; state = get_state(obs), action = action)
    update!(agent.policy, agent.trajectory)
    action
end

#####
# VectorialCompactSARTSATrajectory
#####

function (agent::Agent{<:AbstractPolicy,<:VectorialCompactSARTSATrajectory})(
    ::Training{PreEpisodeStage},
    obs,
)
    if length(agent.trajectory) > 0
        pop!(agent.trajectory, :state, :action)
    end
    nothing
end

function (agent::Agent{<:AbstractPolicy,<:VectorialCompactSARTSATrajectory})(
    ::Training{PreActStage},
    obs,
)
    action = agent.policy(obs)
    push!(agent.trajectory; state = get_state(obs), action = action)
    update!(agent.policy, agent.trajectory)
    action
end

function (agent::Agent{<:AbstractPolicy,<:VectorialCompactSARTSATrajectory})(
    ::Training{PostActStage},
    obs,
)
    push!(agent.trajectory; reward = get_reward(obs), terminal = get_terminal(obs))
    nothing
end

function (agent::Agent{<:AbstractPolicy,<:VectorialCompactSARTSATrajectory})(
    ::Training{PostEpisodeStage},
    obs,
)
    action = agent.policy(obs)
    push!(agent.trajectory; state = get_state(obs), action = action)
    update!(agent.policy, agent.trajectory)
    action
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "COINGLPKBuilder"
version = v"1.10.5"

# Collection of sources required to build COINGLPKBuilder
sources = [
    "https://github.com/coin-or-tools/ThirdParty-Glpk/archive/releases/1.10.5.tar.gz" =>
    "85292f5f1816a5717b384a1f010d1edf75cd03249511b709642b9cab78974819",

]

# Bash recipe for building across all platforms
script = raw"""
set -e
cd $WORKSPACE/srcdir
cd ThirdParty-Glpk-releases-1.10.5/

# Get Glpk from mirror. gnu.org can be unreliable. Note: this will fail if wget or fetch are not present
cp get.Glpk get.Glpk.orig
cat > get.Glpk.patch <<'END'
--- get.Glpk.orig 2018-08-23 15:09:22.682356366 -0400
+++ get.Glpk2018-08-23 15:09:51.562429801 -0400
@@ -25,7 +25,7 @@
 rm -f glpk*.tar.gz

 echo "Downloading the source code from ftp.gnu.org..."
-$wgetcmd ftp://ftp.gnu.org/gnu/glpk/glpk-${glpk_ver}.tar.gz
+$wgetcmd http://ftpmirror.gnu.org/gnu/glpk/glpk-${glpk_ver}.tar.gz

 echo "Uncompressing the tarball..."
 gunzip -f glpk-${glpk_ver}.tar.gz
END
patch -l get.Glpk.orig get.Glpk.patch -o get.Glpk

./get.Glpk
update_configure_scripts
mkdir build
cd build/
../configure --prefix=$prefix --with-pic --disable-pkg-config --host=${target} --enable-shared --disable-static --enable-dependency-linking lt_cv_deplibs_check_method=pass_all
make -j${nproc}
make install

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libcoinglpk", :libcoinglpk)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: juan-pablo-vielma/COINGLPKBuilder

max_issues_repo_name: juan-pablo-vielma/COINGLPKBuilder

max_forks_repo_name: juan-pablo-vielma/COINGLPKBuilder


text:
using Test

#--- parameters
para = ( ϵ  = 0.1, μ = 0.1)  # physical parameters
paraX= ( N  = 6, L  = 3)   # mesh with 6 collocation points: x=[-3.0, -2.0, -1.0, 0.0, 1.0, 2.0]
paraT= ( T  = 2.5, dt = 1) # timegrid with three instants: t=[0.0, 1.0, 2.0]
param = merge(para,paraX)  # used to construct models
parap = merge(paraX,paraT) # used to construct problems

#--- initial data
init     = Init(x->exp.(-x.^2),x-> x )

#--- models
models=[]
push!(models,Boussinesq(param;
                    a=-1,b=3,
                    dealias=2,ktol=1e-12,
                    label="Boussinesq",
                    verbose=false) )

push!(models,WhithamBoussinesq(param;
				α = 1.5, a = -1, b = 3,
				dealias = 2,
				ktol	= 1e-12,
				label 	= "Whitham-Boussinesq",
				verbose	=false) )

push!(models,WhithamBoussinesq(param;
				Boussinesq=true,
				α = 1.5, a = -1, b = 3,
				dealias = 2,
				ktol	= 1e-12,
				label 	= "(Whitham-)Boussinesq",
				verbose	=false) )

push!(models,DeepQuadratic_fast( param;
                    dealias=true, verbose=false, label="fast deep quadratic") )

push!(models,DeepQuadratic( param;
                    dealias=false, verbose=false, label="naive deep quadratic") )

push!(models,Matsuno_fast( param;
                    dealias=true, verbose=false, label="fast Matsuno") )

push!(models,Matsuno( param;
                    dealias=false, verbose=false, label="naive Matsuno") )

push!(models,IsobeKakinuma(param;
            dealias = 2,
            ktol	= 1e-10,
            iterate = true,
            gtol	= 1e-14,
            precond = false,
            restart	= 1,
            maxiter	= 2,
            label	= "Isobe-Kakinuma with GMRES",
            verbose	= false) )

push!(models, modifiedMatsuno(param;
			ν=2,ktol=01e-10,dealias=2,
			verbose=false,
			label="modified Matsuno") )

push!(models,IsobeKakinuma(param::NamedTuple;
			dealias = 0,
			ktol	= 0,
			iterate = false,
			gtol	= 0,
			precond = false,
			restart	= nothing,
			maxiter	= nothing,
			label	= "Isobe-Kakinuma with LU",
			verbose	= false) )

push!(models,NonHydrostatic(param;
            dealias = 2,
            ktol	= 1e-10,
            iterate = true,
            gtol	= 1e-14,
            precond = false,
            restart	= 1,
            maxiter	= 2,
            label	= "non-hydrostatic with GMRES",
            verbose	= false) )

push!(models,NonHydrostatic(param::NamedTuple;
			dealias = 0,
			ktol	= 0,
			iterate = false,
			gtol	= 0,
			precond = false,
			restart	= nothing,
			maxiter	= nothing,
			label	= "non-hydrostatic with LU",
			verbose	= false) )

push!(models,SquareRootDepth(param;
            dealias = 2,
            ktol	= 1e-10,
            iterate = true,
            gtol	= 1e-14,
            precond = false,
            restart	= 1,
            maxiter	= 2,
            label	= "square-root depth with GMRES",
            verbose	= false) )

push!(models,SquareRootDepth(param::NamedTuple;
			dealias = 0,
			ktol	= 0,
			iterate = false,
			gtol	= 0,
			precond = false,
			restart	= nothing,
			maxiter	= nothing,
			label	= "square-root depth with LU",
			verbose	= false) )

push!(models,SerreGreenNaghdi(param;
            dealias = 2,
            ktol	= 1e-10,
            iterate = true,
            gtol	= 1e-14,
            precond = false,
            restart	= 1,
            maxiter	= 2,
            label	= "Green-Naghdi with GMRES",
            verbose	= false) )

push!(models,SerreGreenNaghdi(param::NamedTuple;
			dealias = 0,
			ktol	= 0,
			iterate = false,
			gtol	= 0,
			precond = false,
			restart	= nothing,
			maxiter	= nothing,
			label	= "Green-Naghdi with LU",
			verbose	= false) )

WhithamGreenNaghdi

push!(models,WhithamGreenNaghdi(param;
            dealias = 2,
            ktol	= 1e-10,
            iterate = true,
            gtol	= 1e-14,
            precond = false,
            restart	= 1,
            maxiter	= 2,
            label	= "Whitham-Green-Naghdi with GMRES",
            verbose	= false) )

push!(models,WhithamGreenNaghdi(param::NamedTuple;
			dealias = 0,
			ktol	= 0,
			iterate = false,
			gtol	= 0,
			precond = false,
			restart	= nothing,
			maxiter	= nothing,
			label	= "Whitham-Green-Naghdi with LU",
			verbose	= false) )

push!(models,WhithamGreenNaghdi(param;
			SGN=true,
            dealias = 2,
            ktol	= 1e-10,
            iterate = true,
            gtol	= 1e-14,
            precond = false,
            restart	= 1,
            maxiter	= 2,
            label	= "(Whitham-)Green-Naghdi with GMRES",
            verbose	= false) )

push!(models,WhithamGreenNaghdi(param::NamedTuple;
			SGN=true,
			dealias = 0,
			ktol	= 0,
			iterate = false,
			gtol	= 0,
			precond = false,
			restart	= nothing,
			maxiter	= nothing,
			label	= "(Whitham-)Green-Naghdi with LU",
			verbose	= false) )

push!(models, WWn(param;
			ν		= 2,
			n		= 1,
			δ		= 0.1,
			m		= -12,
			ktol	= 1e-12,
			dealias	= 2,
			label	= "WW1",
			verbose	= false) )

push!(models, WWn(param;
			ν		= 1/2,
			n		= 7,
			δ		= 0.1,
			m		= -Inf,
			ktol	= 1e-12,
			dealias	= 2,
			label	= "WW2",
			verbose	= false) )

push!(models, WWn(param;
			ν		= 0,
			n		= 3,
			δ		= 0.1,
			m		= (-1,2),
			ktol	= 1e-12,
			dealias	= 2,
			label	= "WW3",
			verbose	= false) )

push!(models, WaterWaves(param;
					ν	    = 1,
					IL	    = true,
					method  = 1,
					tol	    = 1e-8,
					maxiter = 5,
					dealias	= 2,
					ktol	= 1e-10,
					label	= "water waves 1",
					verbose	= false) )

push!(models, WaterWaves(param;
					ν	    = 2,
					IL	    = false,
					method  = 2,
					tol	    = 1e-8,
					maxiter = 5,
					dealias	= 1,
					ktol	= 1e-10,
					label	= "water waves 2",
					verbose	= false) )

push!(models, WaterWaves(param;
					ν	    = 1/2,
					IL	    = false,
					method  = 3,
					tol	    = 1e-8,
					maxiter = 5,
					dealias	= 0,
					ktol	= 1e-10,
					label	= "water waves 3",
					verbose	= false) )

#--- tests
for model in models
    # build the initial-value problem
    problem = Problem( model, init, parap; verbose = false )
    # solve the initial-value problem
    solve!(problem ; verbose=false)
    # check everything went well
    @testset "Initial value problem with the model $(model.label)" begin
        @test length(problem.data.U)==length(Times(paraT).tc)
        @test size(problem.data.U[end])==(length(Mesh(paraX).x),2)
        @test !any(isnan,problem.data.U[end][1])
        @test !any(isnan,problem.data.U[end][2])
    end
end


text:
module Proposals
using Distributions: Distribution, logpdf
using Random: AbstractRNG
using ..Traits, ..RNG, ..Var, ..Tagging, ..Space
using ..Util

abstract type ProposalAlg end

export propose

"""
  `propose(rng::AbstractRNG, x, Ω = defΩ())`

Propose `ω::Ω` such that `x(ω)` is well defined with corresnponding proposal probability/density

# Returns
- `ret = (ω = ω_, ℓ = ℓ_)::NamedTuple` where
  - `ret.ω`:
  - `ret.ℓ`:
"""
function propose end

@inline taglogpdf(ω, logpdf_ = 0.0, seen = Set{idtype(ω)}()) =
  tag(ω, (logpdf = (ℓ = Box(logpdf_), seen = seen),))

function propose(rng::AbstractRNG, x, Ω = defΩ())
  ω = tagrng(Ω(), rng)
  ω = taglogpdf(ω)
  ret = x(ω)
  (ω = ω, ℓ = ω.tags.logpdf.ℓ.val)
  # FIXME, should we return ω with rng tagged? Feel like i answers this with rejsample
end

function Var.posthook(::trait(LogPdf), ret, f::Distribution, id, ω)
  #ω.tags.logpdf.val += logpdf(f, ret)
  # FIXME: scope doesn't exist.
  # What should it be?
  if id ∉ ω.tags.logpdf.seen
    ω.tags.logpdf.ℓ.val += logpdf(f, ret)
    push!(ω.tags.logpdf.seen, id)
  end
  nothing
end

## Probs
# id or without id
# avoid duplicats, need seen
# generalize to anything with a logpdf

end


text:
module LatexSVG

@doc read(joinpath(@__DIR__, "..", "README.md"), String) LatexSVG

using Printf: @sprintf
using EzXML: parsexml, root
using Preferences

export
    preamble, reset_preamble!, add_preamble!,
    PDFLaTeX, XeLaTeX, LuaLaTeX, texengine, texengine!,
    LaTeXSVG, latexsvg, savesvg,
    @Lsvg_str,
    config, config!

include("latexdocument.jl")

include("latexengines/interface.jl")
include("latexengines/pdflatex.jl")
include("latexengines/xelatex.jl")
include("latexengines/lualatex.jl")

include("svgtype.jl")
include("svghtmldisplay.jl")
include("svgrender.jl")
include("svgmacros.jl")

include("configurations.jl")
include("initialize.jl")

function __init__()
    _initialize_latex()
    return nothing
end

end # module


text:
function loc2cs3D(loc1::Union{Int64,Array{Int64}},loc2::Union{Int64,Array{Int64}},loc3::Union{Int64,Array{Int64}},n::Array{Int64,1})
@inbounds cs = loc1 .+ (loc2.-1)*n[1] .+ (loc3.-1)*n[1]*n[2];
return cs;
end

# function getDippingMatrix(Mesh::RegularMesh)
# S = kron(speye(Mesh.n[3]),spones(1,Mesh.n[1]*Mesh.n[2]));
# return S;
# end
export generateSamplingMatrix
function generateSamplingMatrix(Mesh::RegularMesh,theta_phi_rad::Array{Float64,2},samplingBinning::Int64 = 1,calcSR = 1)
X = getCellCenteredGrid(Mesh);

mid_domain = [(Mesh.domain[1]+Mesh.domain[2])/2,(Mesh.domain[3]+Mesh.domain[4])/2,(Mesh.domain[5]+Mesh.domain[6])/2];
half_length_domain = [(Mesh.domain[2]-Mesh.domain[1])/2,(Mesh.domain[4]-Mesh.domain[3])/2,(Mesh.domain[6]-Mesh.domain[5])/2];
X[:,1] .= X[:,1] .- mid_domain[1];
X[:,2] .= X[:,2] .- mid_domain[2];
X[:,3] .= X[:,3] .- mid_domain[3];

num_rotations = size(theta_phi_rad,1);
offset = 0;
n =  Mesh.n;
trace_length = Mesh.n[1]*Mesh.n[2];
volume = prod(Mesh.n);

globalRowIdx = zeros(Int32,(volume*num_rotations+5));
globalColIdx = zeros(Int32,(volume*num_rotations+5));

R = 0;
offset = 0;
# println("Begin generating matrix");
# tic()
for k=1:num_rotations
	Rk = inv(getRotate3D(theta_phi_rad[k,1],theta_phi_rad[k,2]));
	Xk = X*Rk';
	# Here we move the domain to be [0,L] to compute the coordinates.
	Xk[:,1] .= (Xk[:,1] .+ half_length_domain[1])*(1.0 ./Mesh.h[1]) .+ 0.5;
	Xk[:,2] .= (Xk[:,2] .+ half_length_domain[2])*(1.0 ./Mesh.h[2]) .+ 0.5;
	Xk[:,3] .= (Xk[:,3] .+ half_length_domain[3])*(1.0 ./Mesh.h[3]) .+ 0.5;
	Xk = round.(Int64,Xk);
	II = findall((Xk[:,1].>=1) .& (Xk[:,1].<=Mesh.n[1]) .& (Xk[:,2].>=1) .& (Xk[:,2].<=Mesh.n[2]) .& (Xk[:,3].>=1) .& (Xk[:,3].<=Mesh.n[3]));
	Xk = Xk[II,:];

	ColIdx_k = loc2cs3D(Xk[:,1],Xk[:,2],Xk[:,3],Mesh.n);
	# Rk = sparse(II,ColIdx_k,ones(size(II)),volume,volume);
	# SRk = S*Rk;

	### IF doing SR
	if calcSR == 1
		II .-= 1;
		II = div.(II,n[1]*n[2]*samplingBinning).+1; ## here we "bin"
		globalRowIdx[(offset+1):(offset + length(II))] = II .+ (k-1)*div(n[3],samplingBinning);
	else ### IF doing R
		globalRowIdx[(offset+1):(offset + length(II))] = II .+ (k-1)*volume;
	end
	globalColIdx[(offset+1):(offset + length(II))] = ColIdx_k;
	offset = offset + length(II);
end
# println("Done collecting idxs for matrix");
X = 0;
II = globalRowIdx .> 0;
globalRowIdx = globalRowIdx[II];
globalColIdx = globalColIdx[II];
II = 0;
V = ones(Float32,size(globalColIdx));
V[:] .= prod(Mesh.h);
## We now define (SR) transposed
if calcSR == 1
	# println("Defining the sparse matrix");
	RT = sparse(globalColIdx,globalRowIdx,V,volume,div(num_rotations*n[3],samplingBinning));
	# println("Done defining the sparse matrix");
else
	RT = sparse(globalColIdx,globalRowIdx,V,volume,num_rotations*volume);
end
# println("Done generating matrix");
return RT;
end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
using LinearAlgebra
import AtomicLevels: angular_integral

@testset "jj2ℓsj" begin
    @testset "Transform matrix" begin
        R = jj2ℓsj(ros"1[s] 2[s-p] k[s-d]"...)
        # Since it is a rotation matrix, its inverse is simply the
        # transpose.
        @test norm(inv(Matrix(R)) - R') < 10eps()
    end

    @testset "Transform individual spin-orbitals" begin
        for o in [o"5s", o"5p", o"5d"]
            for so in spin_orbitals(o)
                ℓ = so.orb.ℓ
                mℓ,ms = so.m
                mj = mℓ + ms
                pure = abs(mj) == ℓ + half(1)
                linear_combination = jj2ℓsj(so)
                @test length(linear_combination) == (pure ? 1 : 2)
                @test norm(last.(linear_combination)) ≈ 1
                @test all(isequal(mj), map(o -> first(o.m), first.(linear_combination)))
            end
        end
    end

    @testset "Angular integrals" begin
        @test isone(angular_integral(SpinOrbital(o"ks", (0,half(1))), SpinOrbital(o"1s", (0,half(1)))))
        @test iszero(angular_integral(SpinOrbital(o"ks", (0,half(1))), SpinOrbital(o"2p", (0,half(1)))))

        @test isone(angular_integral(SpinOrbital(ro"ks", half(1)), SpinOrbital(ro"1s", half(1))))
        @test iszero(angular_integral(SpinOrbital(ro"ks", half(1)), SpinOrbital(ro"2p-", half(1))))

        @test isone(angular_integral(SpinOrbital(o"ks", (0,half(1))), SpinOrbital(ro"1s", half(1))))
        @test isone(angular_integral(SpinOrbital(ro"1s", half(1)), SpinOrbital(o"ks", (0,half(1)))))

        @test angular_integral(SpinOrbital(o"2p", (0,half(1))), SpinOrbital(ro"2p-", half(1))) == clebschgordan(1, 0, half(1), half(1), half(1), half(1))
    end

    @testset "#orbitals = $(length(ros))" for (ros,nb) in ((rsos"l[p]", 4), (rsos"k[s-d] l[d]", 18), (filter(o -> o.m[1]==1//2, rsos"l[d]"), 1))
        os, bs = jj2ℓsj(ros)
        @test length(os) == length(ros)
        for (ro,o) in zip(ros, os)
            @test sum(o.m) == ro.m[1]
            @test o.orb.n == ro.orb.n
            @test o.orb.ℓ == ro.orb.ℓ
        end
        @test length(bs) == nb
        for (i,j,b) in bs
            if size(b,1) == 1
                @test i == j
            else
                @test i ≠ j
            end
            @test b == b'
            @test b^2 ≈ I
        end
    end
end


text:
# implements sparsity-exploiting versions of the main functions that
# take advantage of neighborhoods in the data

# computes a sparse approximation to the conditional neighbor matrix based on
# neighborhood graph "neighbors" computed with a fast nearsest neighbor algorithm (see util)
function conditional_neighbor(X::AbstractMatrix, σ²::AbstractVector{<:Real},
        neighbors::AbstractVecOfVec{Int}, kernel = gaussian, distance = euclidean)
    x = [c for c in eachcol(X)]
    n = length(x)
    K = spzeros(n, n)
    for i in 1:n
        xn, sn = @views x[neighbors[i]], σ²[neighbors[i]]
        @. K[neighbors[i], i] = kernel(distance(xn, (x[i],)) / sqrt(σ²[i]))
    end
    sumK = sum(K, dims = 1) # this is now a sparse operation!
    # @. K ./= sumK # this yields a dense array!
    for (i, k) in enumerate(eachcol(K))
        k ./= sumK[i]
    end
    return K
end

function optimize_variances(X::AbstractMatrix, u::Real, neighbors::AbstractVecOfVec{Int}, kernel = gaussian, distance = euclidean;
                            tol::Real = 1e-6, max_iter::Int = 32)
    n = size(X, 2)
    σ² = fill(5one(eltype(X)), n) # if this is too small initially, Newton has problems
    @threads for i in 1:n
        σ²[i] = optimize_variances(X, u, i, σ²[i], kernel, distance, tol = tol, max_iter = max_iter)
    end
    any(isnan, σ²) && throw(DomainError("Warning: variance optimization yielded NaN: ", σ²))
    return σ²
end

# optimizing variances using sparse operations
# neighbors_i are the indices of the neighbors of x[i]
function optimize_variances(X::AbstractMatrix, u::Real, i::Int, σ²::Real,
                            neighbors_i::AbstractVector{Int},
                            kernel = gaussian, distance = euclidean;
                            tol::Real = 1e-6, max_iter::Int = 32)
    x = [c for c in eachcol(X)]
    n = length(x)
    xn = @view x[neighbors_i]
    function objective(log_σ²::Real) # allocating version, is faster for high d, equal for low
        local σ² = exp(log_σ²)
        p = @. kernel(distance((x[i],), xn)/σ²) # only evaluate kernel function on neighbors
        p ./= sum(p)
        p = SparseVector(n, neighbors_i, p) # create sparse vectors
        entropy(p) - log2(u)
    end
    log_σ² = newtons_method(objective, log(σ²), tol = tol, max_iter = max_iter)
    return exp(log_σ²)
end


text:
export Neighbor
export NeighborSet

"""
    Neighbor{T, F}(idx::T, dist::F)
Simple struct for storing a _neighbor_ index
along with its distance to some query element.
"""
struct Neighbor{T, F}
    idx::T
    dist::F
end

"""
    NeighborSet(n::Neighbor)
An ordered List of ['Neighbor'](@ref), sorted by
the `dist` field.
Elements can be added by calling

    insert!(ns::NeighborsSet, n::Neighbor)

and retrieved with

    nearest(ns::NeighborSet) → Neighbor
    nearest(ns::NeighborSet, k) → k nearest Neighbors
    furthest(ns::NeighborSet) → Neighbor
    pop_nearest!(ns::NeighborSet) → Neighbor
    pop_furthest!(ns::NeighborSet) → Neighbor
"""
struct NeighborSet{T <: Integer, F <: Real}
    neighbor::Vector{Neighbor{T,F}}
end

NeighborSet{T,F}() where {T,F} = NeighborSet{T,F}(Neighbor{T,F}[])
NeighborSet(idx::Integer, dist::Real) = NeighborSet([Neighbor(idx,dist)])
NeighborSet(n::Neighbor) = NeighborSet([n])

nearest(ns::NeighborSet) = ns.neighbor[1]
furthest(ns::NeighborSet) = ns.neighbor[end]
Base.getindex(ns::NeighborSet, i) = ns.neighbor[i]

pop_nearest!(ns::NeighborSet) = popfirst!(ns.neighbor)
pop_furthest!(ns::NeighborSet) = pop!(ns.neighbor)
Base.length(ns::NeighborSet) = length(ns.neighbor)
Base.iterate(ns::NeighborSet) = iterate(ns.neighbor)
Base.iterate(ns::NeighborSet, i) = iterate(ns.neighbor,i)

function Base.insert!(ns::NeighborSet, n::Neighbor)
    idx = searchsortedfirst(ns.neighbor, n, by=x->x.dist)
    insert!(ns.neighbor, idx, n)
    return nothing
end

nearest(ns::NeighborSet, k) = ns.neighbor[1:min(k, length(ns))]


text:
importall ReinforcementLearning

type GymRandomAgent <: AbsAgent end

function play(agent::GymRandomAgent, state::GymState, env::GymEnv; learn=false)
	return sample(env)
end

#env = GymEnv("Pong-v0")
env = GymEnv("MountainCarContinuous-v0")
agent = GymRandomAgent()
threshold = 1000
for i=1:2
	state = getInitialState(env)
	totalRewards = 0.0
	numOfStates = 1.0

	while !isTerminal(state, env) && numOfStates < threshold
		render_env(env)
		action = play(agent, state, env; learn=false)
		state, reward = transfer(env, state, action)
		totalRewards += reward
		numOfStates += 1
	end

	println("Episode i: $i, rewards: $(totalRewards)")
end

max_stars_repo_name: ozanarkancan/ReinforcementLearning

max_issues_repo_name: ozanarkancan/ReinforcementLearning

max_forks_repo_name: ozanarkancan/ReinforcementLearning


text:
export AbstractSampleArray, nframes, nchannels, data, rate, rateHz, domain
export names, names!
export AbstractSpectrumArray, nfreqs, data_no0, domain_no0, slice, zerophase, zerophase!, delay, delay!

import Base: BroadcastStyle, IndexStyle, getindex, setindex!, show, similar, to_index, to_indices, ==
import Base: cat, hcat, vcat, names

# ----- AbstractSampleArray ---------------------

FrameIndex = Union{Colon,AbstractRange,BitArray{1},Vector{Bool},Vector{Int},ClosedInterval{<:Time},ClosedInterval{<:Frequency}}
SingleChannelIndex = Union{Int,Symbol}
MultipleChannelIndex = Union{Vector{Symbol},ClosedInterval{Symbol},Colon,AbstractRange,Vector{Bool},BitArray{1},Vector{Int}}
ChannelIndex = Union{SingleChannelIndex, MultipleChannelIndex}

abstract type AbstractSampleArray{T} <: AbstractMatrix{T} end;

nframes(a::AbstractSampleArray) = size(data(a), 1)
nchannels(a::AbstractSampleArray) = size(data(a), 2)
data(a::AbstractSampleArray) = a.data
rate(a::AbstractSampleArray) = a.rate

names(a::AbstractSampleArray)::Vector{Symbol} = copy(a.names)

# TODO maybe use OrderedDict for names?
function names!(a::AbstractSampleArray, names::Vector{Symbol})
    _check_channel_names(names)
    if length(names) != nchannels(a)
        throw(DimensionMismatch("the number of names given ($(length(names))) does not match the number of channels ($(nchannels(a)))!"))
    end
    a.names .= names
end

function names!(a::AbstractSampleArray, names_::Union{Symbol, Vector{Symbol}}, I::ChannelIndex)::Vector{Symbol}
    cnames = names(a)
    cnames[tochannelidx(a, I)] .= names_
    names!(a, cnames)
end

_default_channel_names(nchannels::Int) = [Symbol(i) for i in 1:nchannels]

function _check_channel_names(names::Vector{Symbol})
    length(names) == length(unique(names)) || throw(ArgumentError("channel names not unique: $(names)!"))
end

function _unique_channel_names(X::AbstractSampleArray...)
    allnames = vcat(map(names, X)...)
    ncnt = counter(Symbol)
    newnames = similar(allnames)
    for (i, name) in enumerate(allnames)
        cnt = ncnt[name]
        if cnt == 0
            newnames[i] = name
            inc!(ncnt, name)
        else
            j = 1
            while true
                newname = Symbol("$(String(name))_$(cnt+j)")
                if ncnt[newname] == 0
                    newnames[i] = newname
                    inc!(ncnt, newname)
                    break
                end
                j += 1
            end
        end
    end
    newnames
end

rateHz(a) = rate(a) * Hz

function domain end

function nameindex(a::AbstractSampleArray, name::Symbol)
   idx = findfirst(isequal(name), a.names)
   isnothing(idx) && throw(ArgumentError("channel :$(name) does not exist!"))
   idx
end

Base.size(a::AbstractSampleArray, dim...) = size(data(a), dim...)

Base.IndexStyle(::Type{T}) where {T <: AbstractSampleArray} = Base.IndexLinear()

@inline Base.:(==)(a::T, b::T) where {T <: AbstractSampleArray} = (rate(a) == rate(b)) && (names(a) == names(b)) && (data(a) == data(b))

@inline toframeidx(::AbstractSampleArray{T}, I) where T = I
@inline toframeidx(a::AbstractSampleArray{T}, i::Int) where T = toframeidx(a, i:i)

@inline tochannelidx(::AbstractSampleArray{T}, I) where T  = I
@inline tochannelidx(a::AbstractSampleArray{T}, i::Int) where T  = tochannelidx(a, i:i)
@inline tochannelidx(a::AbstractSampleArray{T}, name::Symbol) where T  = tochannelidx(a, nameindex(a, name))
@inline tochannelidx(a::AbstractSampleArray{T}, names::Vector{Symbol}) where T  = [nameindex(a, name) for name in names]
@inline tochannelidx(a::AbstractSampleArray{T}, names::ClosedInterval{Symbol}) where T  = nameindex(a, minimum(names)):nameindex(a, maximum(names))

# define Base.to_indices to allow for new type of indexing for Base.IndexLinear()
function Base.to_indices(a::AbstractSampleArray{T}, I::Tuple{FrameIndex}) where T
    # automatically append channel index, so all channels are selected
    to_indices(data(a), (toframeidx(a, I[1]), tochannelidx(a, :)))
end

function Base.to_indices(a::AbstractSampleArray{T}, I::Tuple{FrameIndex, ChannelIndex}) where T
    to_indices(data(a), (toframeidx(a, I[1]), tochannelidx(a, I[2])))
end

# needed for linear indices (Base.IndexLinear())
Base.getindex(a::AbstractSampleArray{T}, i::Int) where {T} = data(a)[i]::T
Base.setindex!(a::AbstractSampleArray{T}, v, i::Int) where T = setindex!(data(a), v, i)::Matrix{T}

Base.BroadcastStyle(::Type{<:AbstractSampleArray}) = Broadcast.ArrayStyle{AbstractSampleArray}()

function _iscompatible(saslist::Vector{<:AbstractSampleArray})
    doms = unique(map(domain, saslist))
    if length(doms) > 1
        throw(DimensionMismatch("can't broadcast different domains: $(doms)!"))
    end
    rates = unique(map(rate, saslist))
    if length(rates) > 1
        throw(ArgumentError("can't broadcast different rates: $(rates)!"))
    end
end

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{T}}, ::Type{ElType}) where {T <: AbstractSampleArray,ElType}
    collectsas(bc::Broadcast.Broadcasted, rest...) = [collectsas(bc.args...)..., collectsas(rest...)...]
    collectsas(a::AbstractSampleArray{T}, rest...) where T = [a, collectsas(rest...)...]
    collectsas(x, rest...) = collectsas(rest...)
    collectsas() = []

    saslist = collectsas(bc)
    _iscompatible(saslist)

    dims = length.(axes(bc))
    imax = argmax(map(nchannels, saslist)) # take the one with most channels to get all channel names
    a = saslist[imax]
    sim = similar(a, ElType, dims)
    for o in saslist
        if names(o) != names(a)[1:nchannels(o)]
            names!(sim, _default_channel_names(nchannels(a)))
            break
        end
    end
    return sim
end

function Base.hcat(X::AbstractSampleArray...)
    length(unique(rate.(X))) == 1 || throw(ArgumentError("hcat: non-unique rates!"))
    length(unique(nframes.(X))) == 1 || throw(ArgumentError("hcat: non-unique number of frames!"))
    newnames = _unique_channel_names(X...)
    data_ = hcat(map(data, X)...)
    return eltype(X)(data_, rate(X[1]), newnames) # eltype gives common supertype
end

function Base.vcat(X::AbstractSampleArray...)
    length(unique(rate.(X))) == 1 || throw(ArgumentError("vcat: non-unique rates!"))
    namelists = names.(X)
    length(unique(namelists)) == 1 || throw(ArgumentError("vcat: non-unique channel names!"))
    data_ = vcat(map(data, X)...)
    return eltype(X)(data_, rate(X[1]), namelists[1])
end

function Base.cat(X::AbstractSampleArray...; dims)
    throw(ErrorException("cat: not implemented for $(eltype(X))"))
end

# ----- AbstractSpectrumArray ---------------------------
abstract type AbstractSpectrumArray{T} <: AbstractSampleArray{T} end

# removed DC
function data_no0 end
function domain_no0 end

@inline toframeidx(::AbstractSpectrumArray{T}, ::R) where {T, R <: ClosedInterval} = throw(ArgumentError("only Frequency intervals allowed!"))
@inline toframeidx(X::AbstractSpectrumArray{T}, ti::R) where {T, R <: ClosedInterval{<:Frequency}} = toindex(X, toHz(minimum(ti))):toindex(X, toHz(maximum(ti)))

@inline Base.:(==)(a::T, b::T) where {T <: AbstractSpectrumArray} = (rate(a) == rate(b)) && (names(a) == names(b)) && (domain(a) == domain(b)) && (data(a) == data(b))

function slice(X::AbstractSpectrumArray, f::Frequency)
    error("missing TEST! is this used anywhere?")
    data(X)[toindex(X, f), :]
end

function zerophase!(X::AbstractSpectrumArray{<:Complex})
    data(X) .= zerophase.(data(X))
    X
end

zerophase(a) = zerophase!(deepcopy(a))

function delay(X::AbstractSpectrumArray{T}, Δt::Time) where T
    T.(MagPhase.(abs.(X), angle.(X) .- 2π .* domain(X) .* tos(Δt)))
end

function delay!(X::AbstractSpectrumArray, shift)
    X .= delay(X, shift)
end


text:
"""
    newton(f,dfdx,x₁[;maxiter,ftol,xtol])

Use Newton's method to find a root of `f` starting from `x₁`, where
`dfdx` is the derivative of `f`. Returns a vector of root estimates.

The optional keyword parameters set the maximum number of iterations
and the stopping tolerance for values of `f` and changes in `x`.
"""
function newton(f,dfdx,x₁;maxiter=40,ftol=100*eps(),xtol=100*eps())
    x = [float(x₁)]
    y = f(x₁)
    Δx = Inf   # for initial pass below
    k = 1

    while (abs(Δx) > xtol) && (abs(y) > ftol)
        dydx = dfdx(x[k])
        Δx = -y/dydx            # Newton step
        push!(x,x[k]+Δx)        # append new estimate

        k += 1
        y = f(x[k])
        if k==maxiter
            @warn "Maximum number of iterations reached."
            break   # exit loop
        end
    end
    return x
end

"""
    secant(f,x₁,x₂[;maxiter,ftol,xtol])

Use the secant method to find a root of `f` starting from `x₁` and
`x₂`. Returns a vector of root estimates.

The optional keyword parameters set the maximum number of iterations
and the stopping tolerance for values of `f` and changes in `x`.
"""
function secant(f,x₁,x₂;maxiter=40,ftol=100*eps(),xtol=100*eps())
    x = [float(x₁),float(x₂)]
    y₁ = f(x₁)
    Δx,y₂ = Inf,Inf   # for initial pass in the loop below
    k = 2

    while (abs(Δx) > xtol) && (abs(y₂) > ftol)
        y₂ = f(x[k])
        Δx = -y₂ * (x[k]-x[k-1]) / (y₂-y₁)   # secant step
        push!(x,x[k]+Δx)        # append new estimate

        k += 1
        y₁ = y₂    # current f-value becomes the old one next time

        if k==maxiter
            @warn "Maximum number of iterations reached."
            break   # exit loop
        end
    end
    return x
end

"""
    newtonsys(f,jac,x₁[;maxiter,ftol,xtol])

Use Newton's method to find a root of a system of equations,
starting from `x₁`. The functions `f` and `jac should return the
residual vector and the Jacobian matrix, respectively. Returns the
history of root estimates as a vector of vectors.

The optional keyword parameters set the maximum number of iterations
and the stopping tolerance for values of `f` and changes in `x`.

"""
function newtonsys(f,jac,x₁;maxiter=40,ftol=1000*eps(),xtol=1000*eps())
    x = [float(x₁)]
    y,J = f(x₁),jac(x₁)
    Δx = Inf   # for initial pass below
    k = 1

    while (norm(Δx) > xtol) && (norm(y) > ftol)
        Δx = -(J\y)             # Newton step
        push!(x,x[k] + Δx)    # append to history
        k += 1
        y,J = f(x[k]),jac(x[k])

        if k==maxiter
            @warn "Maximum number of iterations reached."
            break
        end
    end
    return x
end

"""
    fdjac(f,x₀[,y₀])

Compute a finite-difference approximation of the Jacobian matrix for
`f` at `x₀`, where `y₀`=`f(x₀)` may be given.
"""
function fdjac(f,x₀,y₀=f(x₀))
    δ = sqrt(eps())*max(norm(x₀),1)   # FD step size
    m,n = length(y₀),length(x₀)
    if n==1
        J = (f(x₀+δ) - y₀) / δ
    else
        J = zeros(m,n)
        x = copy(x₀)
        for j in 1:n
            x[j] += δ
            J[:,j] = (f(x) - y₀) / δ
            x[j] -= δ
        end
    end
    return J
end

"""
    levenberg(f,x₁[;maxiter,ftol,xtol])

Use Levenberg's quasi-Newton iteration to find a root of the system
`f`, starting from `x₁`, with `tol` as the stopping tolerance in
both step size and residual norm. Returns the history of root estimates
as a vector of vectors.

The optional keyword parameters set the maximum number of iterations
and the stopping tolerance for values of `f` and changes in `x`.

"""
function levenberg(f,x₁;maxiter=40,ftol=1e-12,xtol=1e-12)
    x = [float(x₁)]
    yₖ = f(x₁)
    k = 1;  s = Inf;
    A = fdjac(f,x[k],yₖ)   # start with FD Jacobian
    jac_is_new = true

    λ = 10;
    while (norm(s) > xtol) && (norm(yₖ) > ftol)
        # Compute the proposed step.
        B = A'*A + λ*I
        z = A'*yₖ
        s = -(B\z)

        x̂ = x[k] + s
        ŷ = f(x̂)

        # Do we accept the result?
        if norm(ŷ) < norm(yₖ)    # accept
            λ = λ/10   # get closer to Newton
            # Broyden update of the Jacobian.
            A += (ŷ-yₖ-A*s)*(s'/(s'*s))
            jac_is_new = false

            push!(x,x̂)
            yₖ = ŷ
            k += 1
        else                       # don't accept
            # Get closer to gradient descent.
            λ = 4λ
            # Re-initialize the Jacobian if it's out of date.
            if !jac_is_new
                A = fdjac(f,x[k],yₖ)
                jac_is_new = true
            end
        end

        if k==maxiter
            @warn "Maximum number of iterations reached."
            break
        end

    end
    return x
end


text:
###############################################################################
## Adrian Lewis' OWOS problem
#    min_(x,y) = 2x² + y² + |x²-y|

struct quadmaxquadAL <: CompositeProblem
    regularizer::maxquadAL
end

quadmaxquadAL() = quadmaxquadAL(maxquadAL())

problem_dimension(::quadmaxquadAL) = 2

## f
# 0th order
function f(::quadmaxquadAL, xy)
    x, y = xy
    return 2*x^2 + y^2
end

# 1st order
function ∇f!(::quadmaxquadAL, res, xy)
    x, y = xy
    return res .= [4x, 2y]
end

# 2nd order
function ∇²f_h!(::quadmaxquadAL, res, x, h)
    res[1] = 4*h[1]
    res[2] = 2*h[2]
    return res
end

# conditioning
# get_gradlips(::quadmaxquadAL)
# get_μ_cvx(::quadmaxquadAL) = (svdvals(pb.A)[end])^2

function firstorder_optimality_tangnorm(pb::quadmaxquadAL, x, M::PlaneParabola, ∇f_x)
    return norm(project(M, x, ∇f_x)), -1
end

function firstorder_optimality_tangnorm(::quadmaxquadAL, x, M::Euclidean, ∇f_x)
    ∇g_x = x[1]^2-x[2] > 0 ? [2*x[1], -1] : [-2*x[1], 1]
    return norm(∇f_x + ∇g_x), -1
end


text:
abstract type AbstractNexResource end

struct Font <: AbstractNexResource
    id::UInt8
end
Font(id::Number) = Font(UInt8(id))

struct Picture <: AbstractNexResource
    id::UInt8
end
Picture(id::Number) = Picture(UInt8(id))
# Int32(pic::Picture) = Int32(pic.id)


text:
function CircleShape(body::Body, radius::Real, offset::Vect)
    CircleShape(ccall(dlsym(libchipmunk, :cpCircleShapeNew), Ptr{Void}, (Ptr{Void}, Cdouble, Vect,), body.ptr, radius, offset))
end

function get_offset(shape::CircleShape)
    ccall(dlsym(libchipmunk, :cpCircleShapeGetOffset), Vect, (Ptr{Void},), shape.ptr)
end

function get_radius(shape::CircleShape)
    ccall(dlsym(libchipmunk, :cpCircleShapeGetRadius), Cdouble, (Ptr{Void},), shape.ptr)
end

function SegmentShape(body::Body, a::Vect, b::Vect, radius::Real)
    SegmentShape(ccall(dlsym(libchipmunk, :cpSegmentShapeNew), Ptr{Void}, (Ptr{Void}, Vect, Vect, Cdouble,), body.ptr, a, b, radius))
end

function set_neighbors(shape::SegmentShape, prev::Vect, next::Vect)
    ccall(dlsym(libchipmunk, :cpSegmentShapeSetNeighbors), Void, (Ptr{Void}, Vect, Vect,), shape.ptr, prev, next)
end

function get_a(shape::SegmentShape)
    ccall(dlsym(libchipmunk, :cpSegmentShapeGetA), Vect, (Ptr{Void},), shape.ptr)
end

function get_b(shape::SegmentShape)
    ccall(dlsym(libchipmunk, :cpSegmentShapeGetB), Vect, (Ptr{Void},), shape.ptr)
end

function get_normal(shape::SegmentShape)
    ccall(dlsym(libchipmunk, :cpSegmentShapeGetNormal), Vect, (Ptr{Void},), shape.ptr)
end

function get_radius(shape::SegmentShape)
    ccall(dlsym(libchipmunk, :cpSegmentShapeGetRadius), Cdouble, (Ptr{Void},), shape.ptr)
end

function free(shape::Shape)
    ccall(dlsym(libchipmunk, :cpShapeFree), Void, (Ptr{Void},), shape.ptr)
end

function get_space(shape::Shape)
    Space(ccall(dlsym(libchipmunk, :cpShapeGetSpace), Ptr{Void}, (Ptr{Void},), shape.ptr))
end

function get_body(shape::Shape)
    Body(ccall(dlsym(libchipmunk, :cpShapeGetBody), Ptr{Void}, (Ptr{Void},), shape.ptr))
end

function set_body(shape::Shape, body::Body)
    ccall(dlsym(libchipmunk, :cpShapeSetBody), Void, (Ptr{Void}, Ptr{Void},), shape.ptr, body.ptr)
end

function get_mass(shape::Shape)
    ccall(dlsym(libchipmunk, :cpShapeGetMass), Cdouble, (Ptr{Void},), shape.ptr)
end

function set_mass(shape::Shape, mass::Real)
    ccall(dlsym(libchipmunk, :cpShapeSetMass), Void, (Ptr{Void}, Cdouble,), shape.ptr, mass)
end

function get_density(shape::Shape)
    ccall(dlsym(libchipmunk, :cpShapeGetDensity), Cdouble, (Ptr{Void},), shape.ptr)
end

function set_density(shape::Shape, density::Real)
    ccall(dlsym(libchipmunk, :cpShapeSetDensity), Void, (Ptr{Void}, Cdouble,), shape.ptr, density)
end

function get_moment(shape::Shape)
    ccall(dlsym(libchipmunk, :cpShapeGetMoment), Cdouble, (Ptr{Void},), shape.ptr)
end

function get_area(shape::Shape)
    ccall(dlsym(libchipmunk, :cpShapeGetArea), Cdouble, (Ptr{Void},), shape.ptr)
end

function get_center_of_gravity(shape::Shape)
    ccall(dlsym(libchipmunk, :cpShapeGetCenterOfGravity), Vect, (Ptr{Void},), shape.ptr)
end

function get_sensor(shape::Shape)
    Bool(ccall(libchipmunk, :cpShapeGetSensor), Uint8, (Ptr{Void},), shape.ptr)
end

function set_sensor(shape::Shape, sensor::Bool)
    ccall(dlsym(libchipmunk, :cpShapeSetSensor), Void, (Ptr{Void}, Int32,), shape.ptr, sensor)
end

function get_elasticity(shape::Shape)
    ccall(dlsym(libchipmunk, :cpShapeGetElasticity), Cdouble, (Ptr{Void}, Int32,), shape.ptr)
end

function set_elasticity(shape::Shape, elasticity::Real)
    ccall(dlsym(libchipmunk, :cpShapeSetElasticity), Void, (Ptr{Void}, Cdouble,), shape.ptr, elasticity)
end

function get_friction(shape::Shape)
    ccall(dlsym(libchipmunk, :cpShapeGetFriction), Cdouble, (Ptr{Void},), shape.ptr)
end

function set_friction(shape::Shape, friction::Real)
    ccall(dlsym(libchipmunk, :cpShapeSetFriction), Void, (Ptr{Void}, Cdouble,), shape.ptr, friction)
end

function get_surface_velocity(shape::Shape)
    ccall(dlsym(libchipmunk, :cpShapeGetSurfaceVelocity), Vect, (Ptr{Void},), shape.ptr)
end

function set_surface_velocity(shape::Shape, velocity::Vect)
    ccall(dlsym(libchipmunk, :cpShapeSetSurfaceVelocity), Void, (Ptr{Void}, Vect,), shape.ptr, velocity)
end

function get_userdata(shape::Shape)
    ccall(dlsym(libchipmunk, :cpShapeGetUserData), Ptr{Void}, (Ptr{Void},), shape.ptr)
end

function set_userdata(shape::Shape, userdata::Ptr{Void})
    ccall(dlsym(libchipmunk, :cpShapeSetUserData), Void, (Ptr{Void}, Ptr{Void},), shape.ptr, userdata)
end

function get_collision_type(shape::Shape)
    ccall(dlsym(libchipmunk, :cpShapeGetCollisionType), Uint32, (Ptr{Void},), shape.ptr)
end

function set_collision_type(shape::Shape, ctype::Integer)
    ccall(dlsym(libchipmunk, :cpShapeSetCollisionType), Void, (Ptr{Void}, Uint32,), shape.ptr, ctype)
end

export free, get_space, get_body, set_body, get_mass, set_mass, get_density, set_density, get_moment, get_area,
get_center_of_gravity, get_sensor, set_sensor, get_elasticity, set_elasticity, get_friction, set_friction,
CircleShape, SegmentShape, get_surface_velocity, set_surface_velocity, get_userdata, set_userdata,
get_collision_type, set_collision_type, get_offset, get_radius, get_a, get_b, get_normal


text:
##################################################
# Computation of the local velocity deflections
##################################################

### Todo
# change vmax by wmax
# for coherence with notes and paper

"""
    _vmax(r,vr,vt,sinth,costh,cosph)

Computes the maximum allowed velocity in Rosenbluth integrals, for a test star at position `r` and velocity `(vr,vt)`.
At fixed integration variable (over field) `(θ,φ)` such that `sin(θ)=sinth`, `cos(θ)=costh` and `cos(φ)=cosph`.
Returns false if no such vmax exists (i.e. non-positive discriminant).

# Remarks
- Integrands in Rosenbluth integrals vanishes for v > _vmax.

# Arguments
- `r::Float64`: Radial position of the test star.
- `vr::Float64`: Radial velocity of the test star.
- `vt::Float64`: Tangential velocity of the test star.
- `sinth::Float64`: Value `sin(θ)=sinth` corresponding to `θ` integration variable.
- `costh::Float64`: Value `cos(θ)=costh` corresponding to `θ` integration variable.
- `cosph::Float64`: Value `cos(φ)=cosph` corresponding to `φ` integration variable.
"""
function _vmax(r::Float64, vr::Float64, vt::Float64,
    sinth::Float64,costh::Float64,cosph::Float64)

    vSq = vr^2 + vt^2
    deltav = (vr*costh+vt*sinth*cosph)^2 - (2*psi(r)+vSq)

    if (deltav > 0.0)
        return (vr*costh+vt*sinth*cosph) + sqrt(deltav)
    else
        return false # in order to test whether the discriminant is positive or not
    end

end

"""
    _EaLa(r,var2,vat2)

Computes the energy and angular momentum used in Rosenbluth integrands.

# Arguments
- `r::Float64`: Radial position of the test star.
- `var2::Float64`: Square radial velocity in Rosenbluth integrand.
- `vat2::Float64`: Square tangential velocity in Rosenbluth integrand.
"""
function _EaLa(r::Float64, var2::Float64, vat2::Float64)
    Ea = psi(r) + (1/2)*(var2+vat2)
    La = r*sqrt(vat2)
    return Ea, La
end

"""
    _fa(r,vr,vt,vp,snth,csth,snph,csph)

Distribution function used in Rosenbluth integrand

# Arguments
- `r::Float64`: Radial position of the test star.
- `vr::Float64`: Radial velocity of the test star.
- `vt::Float64`: Tangential velocity of the test star.
- `vp::Float64`: Velocity integration variable.
- `snth::Float64`: Value `sin(θ)=snth` corresponding to `θ` integration variable.
- `csth::Float64`: Value `cos(θ)=csth` corresponding to `θ` integration variable.
- `snph::Float64`: Value `sin(φ)=snph` corresponding to `φ` integration variable.
- `csph::Float64`: Value `cos(φ)=csph` corresponding to `φ` integration variable.
"""
function _fa(r::Float64, vr::Float64, vt::Float64, vp::Float64,
                snth::Float64,csth::Float64,snph::Float64,csph::Float64)

    pref_r = vr-vp*csth
    pref_t = vt-vp*snth*csph
    pref_t2 = vp*snth*snph

    var2 = pref_r^2
    vat2 = pref_t^2+pref_t2^2
    Ea, La = _EaLa(r,var2,vat2)
    DF = _F(Ea,La)

    return DF
end

"""
    RosenbluthPotentials(r,vr,vt,[nbK=nbK_default])

Computation of the Rosenbluth potential and its derivatives.
Returns a tuple `(dh/dvr, dh/dvt, dg/dvt, d2g/dvr2, d2g/dvrdvt, d2g/dvt2)`.

# Arguments
- `r::Float64`: Radial position of the test star.
- `vr::Float64`: Radial velocity of the test star.
- `vt::Float64`: Tangential velocity of the test star.
- `nbK::Int64`: Sampling number of Rosenbluth integrals.
"""
function RosenbluthPotentials(r::Float64, vr::Float64, vt::Float64, nbK::Int64=nbK_default)
    #X = (v',x,phi)

    sumhr = 0
    sumht = 0
    sumgt = 0
    sumh = 0
    sumgrrInt = 0
    sumgrt = 0
    sumgttInt = 0
    for ith=1:nbK
        thtemp = (2*ith-1)/(2*nbK)
        sinthtemp,costhtemp = sincos(PI*thtemp)
        sumhr_ph = 0
        sumht_ph = 0
        sumgt_ph = 0
        sumh_ph = 0
        sumgrrInt_ph = 0
        sumgrt_ph = 0
        sumgttInt_ph = 0

        for iph=1:nbK
            phtemp = (2*iph-1)/(2*nbK)
            sinphtemp, cosphtemp = sincos(2*PI*phtemp)
            vmax = _vmax(r,vr,vt,sinthtemp,costhtemp,cosphtemp)

            if (vmax != false)
                # if discriminant > 0 (i.e there exists an acceptable w-range where E'<0)
                # this is always the case for bound stars

                if (vmax > 0)
                    # if acceptable range has some part where w>0
                    # always the case for bounds stars
                    # happens for some unbound stars (in that case, one can optimize this code as vmin>0)

                    sumhr_v = 0
                    sumht_v = 0
                    sumgt_v = 0
                    sumh_v = 0
                    sumgrrInt_v = 0
                    sumgrt_v = 0
                    sumgttInt_v = 0

                    for iv=1:nbK
                        vtemp = (2*iv-1)/(2*nbK)
                        floc = _fa(r,vr,vt,vmax*vtemp,sinthtemp,costhtemp,sinphtemp,cosphtemp)
                        sumhr_v += -floc
                        sumht_v += -floc
                        sumgt_v += vtemp^2*floc

                        # compute h and the other integrals separately
                        sumh_v += vtemp*floc
                        sumgrrInt_v += -vtemp*floc
                        sumgrt_v += -vtemp*floc
                        sumgttInt_v += -vtemp*floc
                    end
                    sumhr_v *= vmax
                    sumht_v *= vmax
                    sumgt_v *= vmax^3
                    sumh_v *= vmax^2
                    sumgrrInt_v *= vmax^2
                    sumgrt_v *= vmax^2
                    sumgttInt_v *= vmax^2

                    sumhr_ph += sumhr_v
                    sumht_ph += sumht_v*cosphtemp
                    sumgt_ph += sumgt_v*cosphtemp
                    sumh_ph += sumh_v
                    sumgrrInt_ph += sumgrrInt_v
                    sumgrt_ph += sumgrt_v*cosphtemp
                    sumgttInt_ph += sumgttInt_v*cosphtemp^2
                end

                # when (vmax <= 0), there is no positive w-range, hence w-integral yields 0
                # sometimes happens for unbound stars
                # in that case, add nothing

            end

            # if discriminant < 0 (i.e for all w, E'>=0)
            # sometimes happens for unbound stars
            # hence w-integral yields 0
            # in that case, add nothing

        end
        sumhr += sinthtemp*costhtemp*sumhr_ph
        sumht += sinthtemp^2*sumht_ph
        sumgt += sinthtemp^2*sumgt_ph
        sumh += sinthtemp*sumh_ph
        sumgrrInt += sinthtemp*costhtemp^2*sumgrrInt_ph
        sumgrt += sinthtemp^2*costhtemp*sumgrt_ph
        sumgttInt += sinthtemp^3*sumgttInt_ph
    end
    sumhr *= PI*2*PI/(nbK^3)
    sumht *= PI*2*PI/(nbK^3)
    sumgt *= PI*2*PI/(nbK^3)
    sumh *= PI*2*PI/(nbK^3)
    sumgrrInt *= PI*2*PI/(nbK^3)
    sumgrt *= PI*2*PI/(nbK^3)
    sumgttInt *= PI*2*PI/(nbK^3)

    return sumhr, sumht, sumgt, sumh+sumgrrInt, sumgrt, sumh+sumgttInt
end

"""
    localVelChange(r,vr,vt,m_field,[nbK=nbK_default,m_test=m_field])

Computation of the local velocity deflections along and perpendicular to the trajectory of the test star.
Returns a tuple `(dvPar, dvPar2, dvPerp2)`.

# Arguments
- `r::Float64`: Radial position of the test star.
- `vr::Float64`: Radial velocity of the test star.
- `vt::Float64`: Tangential velocity of the test star.
- `m_field::Float64`: Mass of field star.
- `nbK::Int64`: Sampling number of Rosenbluth integrals.
- `m_test::Float64`: Mass of test star.
"""
function localVelChange(r::Float64, vr::Float64, vt::Float64,
                        m_field::Float64, nbK::Int64=nbK_default,
                        m_test::Float64=m_field)

    dhdvr, dhdvt, dgdvt, d2gdvr2, d2gdvrdvt, d2gdvt2 = RosenbluthPotentials(r,vr,vt,nbK)

    cst       = 4.0*PI*_G^2*logCoulomb
    v         = sqrt(vr^2+vt^2)

    if (vt>0)
        vr_v   = vr/v
        vt_v   = vt/v

        dvPar   = cst*(m_field+m_test) *(vr_v*dhdvr+vt_v*dhdvt)
        dvPar2  = cst* m_field         *(vr_v^2*d2gdvr2+(2*vr_v*vt_v)*d2gdvrdvt
                                        +vt_v^2*d2gdvt2)
        dvPerp2 = cst* m_field         *(vt_v^2*d2gdvr2-(2*vr_v*vt_v)*d2gdvrdvt
                                        +vr_v^2*d2gdvt2+(1/vt)*dgdvt)
    else # vt = 0

        dvPar   = cst*(m_field+m_test) *dhdvr
        dvPar2  = cst* m_field         *d2gdvr2
        dvPerp2 = cst* m_field         *2*d2gdvt2

    end

    return dvPar, dvPar2, dvPerp2
end

"""
    localOrbitChange(r,E,L,m_field,[nbK=nbK_default,m_test=m_field])

Computation of the local energy and angular momentum deflections.
Returns a tuple `(dE, dE2, dL, dEdL, dL2)`.

# Arguments
- `r::Float64`: Radial position of the test star.
- `E::Float64`: Energy of the test star.
- `L::Float64`: Angular momentum of the test star.
- `m_field::Float64`: Mass of field star.
- `nbK::Int64`: Sampling number of Rosenbluth integrals.
- `m_test::Float64`: Mass of test star.
"""
function localOrbitChange(r::Float64, E::Float64, L::Float64,
                           m_field::Float64, nbK::Int64=nbK_default,
                           m_test::Float64=m_field)

    vr = radialVelocity(r,E,L)
    vt = tangentVelocity(r,E,L)
    v = sqrt(vr^2+vt^2)

    dvPar, dvPar2, dvPerp2 = localVelChange(r,vr,vt,m_field,nbK,m_test)

    dE   = dvPar*(v)   + dvPar2*(1/2)     + dvPerp2*(1/2)
    dE2  =               dvPar2*(v^2)
    dL   = dvPar*(L/v)                    + dvPerp2*(r^2/(4.0*L))
    dL2  =               dvPar2*(L^2/v^2) + dvPerp2*(1/2 * (r^2 - L^2/v^2))
    dEdL =               dvPar2*(L)

    return dE, dE2, dL, dEdL, dL2
end

max_stars_repo_name: KerwannTEP/CAT

max_issues_repo_name: KerwannTEP/CAT

max_forks_repo_name: KerwannTEP/CAT


text:
module BOMoD

    using Kronecker
    using StatsBase
    using Statistics
    using StringDistances
    using Distances
    using DataFrames
    using Random
    using LinearAlgebra
    using Optim
    using Stheno
    using Distributions

    import Base: +, *, getindex,length, eltype,isequal,isless,push!,in,summary,show,==
    import Stheno:ew,pw, BaseKernel,Kernel
    import LinearAlgebra: norm, eigvals!, Diagonal,dot
    import Optim: minimizer
    import Kronecker: KroneckerPower
    import StringDistances: QGramDistance

    # include construction of space
    include(joinpath("design", "Mod.jl"))
    include(joinpath("design", "Combinatroics.jl"))
    include(joinpath("design", "Constraints.jl"))
    include(joinpath("design", "Space.jl"))
    include(joinpath("design", "Construct.jl"))
    include(joinpath("design", "Designspace.jl"))
    include(joinpath("design", "Sample.jl"))
    #BO
    include(joinpath("BO", "util_graph.jl"))
    include(joinpath("BO", "Kernels.jl"))
    include(joinpath("BO", "GP_model.jl"))
    include(joinpath("BO", "batch_sample.jl"))

    export groupmod, getspace, Mod, constructdesign, OrderedConstraint
    export UnOrderedConstraint, OrderedConstruct,  UnorderedConstruct
    export GPModel, GPpredict, fit_gp, predict_gp, ts_sampler, ei_sampler
    export pi_sampler, gpucb_sampler, QGramKernel, EditDistancesKernel
    export PRandomKernel, DiffusionKernel, sample, CosRule, LevRule, Linear

end

max_stars_repo_name: SamDier/BOMoD

max_issues_repo_name: SamDier/BOMoD

max_forks_repo_name: SamDier/BOMoD


text:
"""
    WhittedIntegrator

The standard (and currently only) SurfaceIntegrator. Traces rays from surface point to light
sources and adds together their contributions. Recursion to ray depths greater than one is not
currently implemented yet.

"""
struct WhittedIntegrator <: SurfaceIntegrator
    "maxdepth, the maximum ray depth to trace recursively."
    maxdepth::Int64
end

"""
    intensity(integrator, scene, isect, ray, sample)

Return intensity given a WhittedIntegrator, a scene, intersection point, ray and camera
sample.

# Arguments
- `integrator::WhittedIntegrator`: a WhittedIntegrator instance.
- `scene::Scene`: the Scene to trace illumination rays in.
- `isect::Intersection`: an Intersection, giving the local geometry at ray intersection.
- `ray::Ray`, the incident Ray (not used by this function but part of SurfaceIntegrator API)
- `sample::Sample`, a Sample (not used by this function but part of SurfaceIntegrator API)

"""
function intensity(intgr::WhittedIntegrator, scene::Scene, isect::Intersection,
                   ray::Ray, sample::Sample)

    # add contribution of each light source
    L = nolight
    for light in scene.lights
        L = L .+ inner_int(light, isect.geometry, isect.material, -ray.direction, scene)
    end
    return L

    # TODO: recursive for depth
#    if ray.depth < maxDepth
        # trace rays for transmission and reflection
#    end
end

const RAY_EPS = 1e-5

"""
    inner_int(light::LightSource, dg::DifferentialGeometry, mat::BxDF, w1::Vector3, scene::Scene)

The inner loop of the `intensity` function, separated for better optimization.

"""
function inner_int(light::LightSource, dg::DifferentialGeometry, mat::BxDF, w1::Vector3, scene::Scene)
    light_spectrum, pdf, light_ray = sample_L(light, dg.p + RAY_EPS * dg.n)
    w0 = light_ray.direction
    if dot(w0, dg.n) < 0.0 || intersectP(light_ray, scene)
        return nolight
    else
        brdf_value = evaluate(mat, dg.T, w1, w0)
        return (brdf_value .* light_spectrum) .* dot(w0, dg.n)
    end
end

"""
    PressureIntegrator

A Whitted integrator that computes radiation pressure force and torque.

"""
struct PressureIntegrator <: SurfaceIntegrator
    maxdepth::Int64
    force::Array{Vector3,1}
    torque::Array{Vector3,1}
    counts::Array{Int64,1}
end

"""

"""
function compute_pressure(intgr::PressureIntegrator, light::LightSource, isect::Intersection,
                   ray::Ray)
    world_to_local = isect.geometry.T
    local_to_world = inv(world_to_local)
    incident_direction = world_to_local(-ray.direction)
    f = compute_pressure(isect.material, incident_direction, light.intensity)
    f = local_to_world(f)
    t = cross(f, isect.geometry.p)

    return f, t
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
# This file was generated, do not modify it. # hide
foo(x) = sum(abs.(x)) / length(x)
d = describe(boston, :mean, :median, foo => :foo)
first(d, 3)


text:
module TestStorage

using LibCloud.Storage
using Base.Test
using PyCall

const TEST_BKT_NAME = "julia-libcloud-test"

function s3_like(storage, do_download=true)
    println("Regions:\n", list_regions(storage))

    println("Creating test container ($TEST_BKT_NAME)...")
	cont = create_container(storage, TEST_BKT_NAME)
    @test isa(cont, Container)
    @test cont.name == TEST_BKT_NAME

    println("Creating temp file...")
    fname, io = mktemp()
    for x in 1:10
        println(io, "hello world")
    end
    close(io)

    println("Uploading test files...")
    for idx in 1:5
        meta_data = Dict{Any,Any}("idx" => idx)
        extra = Dict{Any,Any}("meta_data" => meta_data)
        println("    - ", meta_data["idx"])
        oname = "storage_test_$idx"
        obj = upload_object(cont, fname, oname, extra)
        @test isa(obj, Object)
        @test obj.name == oname
    end

    if do_download
        println("Downloading test files...")
        for idx in 1:5
            obj = get_object(cont, "storage_test_$idx")
            @test isa(obj, Object)
            println("    - ", obj.meta_data["idx"])
            @test download_object(cont, obj, fname, overwrite_existing=true)
        end
    end

    println("Listing files...")
    for container in iterate_containers(storage)
        @test isa(container, Container)
        println(container)
        try
            objcount = 0
            for obj in iterate_container_objects(storage, container)
                @test isa(obj, Object)
                println("    - ", obj)
                objcount += 1
                if objcount > 20
                    break
                end
            end
        end
    end

    println("Deleting test files...")
    for idx in 1:5
        obj = get_object(cont, "storage_test_$idx")
        @test isa(obj, Object)
        println("    - ", obj.meta_data["idx"])
        @test delete_object(cont, obj)
    end

    println("Deleting test container ($TEST_BKT_NAME)...")
	@test delete_container(storage, cont)

    println("Deleting temp file ($fname)...")
    rm(fname)

    println("Done.")
end

function s3(key, secret)
    storage = StorageDriver(StorageProvider.S3, key, secret)
    s3_like(storage)
end

function gce(service_account_id, key, project)
	storage = StorageDriver(StorageProvider.GOOGLE_STORAGE, service_account_id, key, project=project)
    s3_like(storage)
end

function dummy()
	storage = StorageDriver(StorageProvider.DUMMY, "key", "secret")
    s3_like(storage, false)
end

end # module TestStorage