text:
module FESpaceFactoriesTests

using Test

using Gridap.Geometry
using Gridap.Arrays
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.FESpaces
using Gridap.CellData

domain = (0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)
order = 3

V = FESpace(model=model,reffe=QUAD4,conformity=GradConformity())
v1 = FEFunction(V,rand(num_free_dofs(V)))
@test isa(V,UnconstrainedFESpace)

V = FESpace(model=model,reffe=QUAD4,conformity=L2Conformity())
v2 = FEFunction(V,rand(num_free_dofs(V)))
@test isa(V,UnconstrainedFESpace)

V = FESpace(model=model,reffe=QUAD4,conformity=CDConformity((CONT,DISC)))
v3 = FEFunction(V,rand(num_free_dofs(V)))
@test isa(V,UnconstrainedFESpace)

#using Gridap.Visualization
#writevtk(Triangulation(model),"results",cellfields=["v1"=>v1,"v2"=>v2,"v3"=>v3])

V = FESpace(
 model=model,
 valuetype=Float64,
 reffe=:Lagrangian,
 order=order,
 conformity=:L2)

@test isa(V,UnconstrainedFESpace)

V = FESpace(
 triangulation=get_triangulation(model),
 valuetype=Float64,
 reffe=:Lagrangian,
 order=order,
 conformity=:L2)

@test isa(V,UnconstrainedFESpace)

V = FESpace(
 model=model,
 valuetype=Float64,
 reffe=:SLagrangian,
 order=order,
 conformity=:H1)

@test isa(V,UnconstrainedFESpace)

V = FESpace(
 model=model,
 valuetype=Float64,
 reffe=:QLagrangian,
 conformity=:H1,
 order=order,
 dirichlet_tags="boundary")

@test isa(V,UnconstrainedFESpace)

D = num_point_dims(model)

V = FESpace(
 model=model,
 valuetype=VectorValue{D,Float64},
 reffe=:Lagrangian,
 order=order,
 conformity=:H1,
 dirichlet_tags="boundary",
 dirichlet_masks=(true,false))

@test isa(V,UnconstrainedFESpace)

V = FESpace(
 model=model,
 valuetype=VectorValue{D,Float64},
 reffe=:SLagrangian,
 order=order,
 conformity=:H1,
 dirichlet_tags=[1,2],
 dirichlet_masks=[(true,false),(false,true)])

@test isa(V,UnconstrainedFESpace)

V = FESpace(
 model=model,
 valuetype=Float64,
 reffe=:PLagrangian,
 order=order,
 conformity=:L2)

#using Gridap.Visualization
#uh = FEFunction(V,rand(num_free_dofs(V)))
#writevtk(get_triangulation(model),"trian",nsubcells=20,cellfields=["uh"=>uh])

@test isa(V,UnconstrainedFESpace)

V = FESpace(
 model=model,
 valuetype=Float64,
 reffe=:PLagrangian,
 order=order,
 conformity=:L2,
 constraint=:zeromean)

@test isa(V,ZeroMeanFESpace)

trian = reindex(Triangulation(model),1:4)
quad = CellQuadrature(trian,order)

V = FESpace(
 model=model,
 valuetype=Float64,
 reffe=:PLagrangian,
 order=order,
 conformity=:L2,
 constraint=:zeromean,
 zeromean_trian = trian,
 zeromean_quad = quad)

@test abs(V.vol - 4/9) < 1.0e-9

V = FESpace(
 model=model,
 valuetype=Float64,
 reffe=:PLagrangian,
 order=order,
 conformity=:L2,
 constraint=:zeromean,
 zeromean_trian = trian)

@test abs(V.vol - 4/9) < 1.0e-9

@testset "FESpace argchecks" begin
    model = CartesianDiscreteModel((0, 2), 10)
    trian = Triangulation(model)
    TestFESpace(model=model, reffe=:Lagrangian, order=1, valuetype=Float64, dirichlet_tags="boundary")
    FESpace(model=model, reffe=:Lagrangian, order=1, valuetype=Float64, dirichlet_tags="boundary")

    # dirichlet_tags spelling error!
    @test_throws ArgumentError TestFESpace(model=model, reffe=:Lagrangian, order=1, valuetype=Float64,
        dirichlet_tag="boundary"
    )
    @test_throws ArgumentError FESpace(model=model, reffe=:Lagrangian, order=1, valuetype=Float64,
        dirichlet_tag="boundary"
    )
end

end # module


text:
"""
    solver_name(model::MOI.ModelLike, config::TestConfig)

Test that the [`MOI.SolverName`](@ref) attribute is implemented for `model`.
"""
function solver_name(model::MOI.ModelLike, config::TestConfig)
    if config.solve
        @test MOI.get(model, MOI.SolverName()) isa AbstractString
    end
end
unittests["solver_name"] = solver_name

"""
    raw_status_string(model::MOI.ModelLike, config::TestConfig)

Test that the [`MOI.RawStatusString`](@ref) attribute is implemented for
`model`.
"""
function raw_status_string(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, MOI.SingleVariable(x), MOI.GreaterThan(0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(),
            MOI.SingleVariable(x))
    test_model_solution(model, config, objective_value = 0.0,
                        variable_primal = [(x, 0.0)])
    if config.solve
        @test MOI.get(model, MOI.RawStatusString()) isa AbstractString
    end
end
unittests["raw_status_string"] = raw_status_string

"""
    solve_time(model::MOI.ModelLike, config::TestConfig)

Test that the [`MOI.SolveTime`](@ref) attribute is implemented for `model`.
"""
function solve_time(model::MOI.ModelLike, config::TestConfig)
    MOI.empty!(model)
    @test MOI.is_empty(model)
    x = MOI.add_variable(model)
    MOI.add_constraint(model, MOI.SingleVariable(x), MOI.LessThan(0.0))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(),
            MOI.SingleVariable(x))
    test_model_solution(model, config, objective_value = 0.0,
                        variable_primal = [(x, 0.0)])
    if config.solve
        time = MOI.get(model, MOI.SolveTime())
        @test time ≥ 0.0
    end
end
unittests["solve_time"] = solve_time

ext: jl


text:
# Run tests on nep_gallery (not tested elsewhere)

using NonlinearEigenproblemsTest
using NonlinearEigenproblems
using LinearAlgebra
using SparseArrays
using Test

@testset "Gallery" begin
    @testset "Gallery extra loading" begin
        # Test if it is possible to load some extra gallery
        @test (using GalleryWaveguide) == nothing
    end

    @testset "Basic random examples" begin
        @info "Testing dep0"
        nep = nep_gallery("dep0")
        nep = nep_gallery("dep0", 15)
        A = compute_Mder(nep, 3.0)
        @test isa(nep, DEP) # Not so sophisticated, but improves code coverage
        @test_throws MethodError nep_gallery("dep0", 15, 8)
        @test_throws MethodError nep_gallery("dep0", 15, t=8)

        @info "Testing dep0_sparse"
        nep = nep_gallery("dep0_sparse")
        nep = nep_gallery("dep0_sparse", 20)
        nep = nep_gallery("dep0_sparse", 20, 0.25)
        A = compute_Mder(nep, 3.0)
        @test isa(nep, DEP) # Not so sophisticated, but improves code coverage
        @test_throws MethodError nep_gallery("dep0_sparse", 20, 0.25, 8)
        @test_throws MethodError nep_gallery("dep0_sparse", 20, 0.25, n=8)

        @info "Testing dep0_tridiag"
        nep = nep_gallery("dep0_tridiag")
        nep = nep_gallery("dep0_tridiag", 15)
        A = compute_Mder(nep, 3.0)
        @test isa(nep, DEP) # Not so sophisticated, but improves code coverage
        @test_throws MethodError nep_gallery("dep0", 15, 8)
        @test_throws MethodError nep_gallery("dep0", 15, t=8)

        @info "pep0"
        nep = nep_gallery("pep0")
        nep = nep_gallery("pep0", 15)
        A = compute_Mder(nep, 3.0)
        @test isa(nep, PEP) # Not so sophisticated, but improves code coverage
        @test_throws MethodError nep_gallery("pep0", 15, 8)
        @test_throws MethodError nep_gallery("pep0", 15, t=8)

        @info "pep0_sym"
        nep = nep_gallery("pep0_sym")
        nep = nep_gallery("pep0_sym", 15)
        A = compute_Mder(nep, 3.0)
        @test isa(nep, PEP) # Not so sophisticated, but improves code coverage
        @test_throws MethodError nep_gallery("pep0_sym", 15, 8)
        @test_throws MethodError nep_gallery("pep0_sym", 15, t=8)

        @info "Testing pep0_sparse"
        nep=nep_gallery("pep0_sparse")
        nep=nep_gallery("pep0_sparse", 15)
        nep=nep_gallery("pep0_sparse", 15, 0.12)
        A=compute_Mder(nep,3.0)
        @test issparse(A) # Not so sophisticated, but improves code coverage
        @test_throws MethodError nep_gallery("pep0_sparse", 15, 0.12, 8)
        @test_throws MethodError nep_gallery("pep0_sparse", 15, 0.12, t=8)

        @info "Testing pdde_stability"
        pep=nep_gallery("nlevp_native_pdde_stability",3);
        λ=-0.5064627366803098 - 0.00019049029618525343im; # Reference solution
        @test minimum(svdvals!(Matrix(compute_Mder(pep,λ)))) < eps()*100;

    end

    @info "Testing dep_symm_double"
    nep = nep_gallery("dep_symm_double")
    nep = nep_gallery("dep_symm_double", 15)
    A = compute_Mder(nep, 3.0)
    @test isa(nep, DEP) # Not so sophisticated, but improves code coverage
    @test_throws MethodError nep_gallery("dep_symm_double", 15, 8)
    @test_throws MethodError nep_gallery("dep_symm_double", 15, t=8)

    @info "Testing dep_double"
    nep = nep_gallery("dep_double")
    A = compute_Mder(nep, 3.0)
    @test isa(nep, DEP) # Not so sophisticated, but improves code coverage
    @test_throws MethodError nep_gallery("dep_double", 15)
    @test_throws MethodError nep_gallery("dep_double", t=15)

    @info "real_quadratic"
    nep = nep_gallery("real_quadratic")
    A = compute_Mder(nep, 3.0)
    @test isa(nep, PEP) # Not so sophisticated, but improves code coverage
    @test_throws MethodError nep_gallery("real_quadratic", 15)
    @test_throws MethodError nep_gallery("real_quadratic", t=15)

    for i = [0, 1]
        nep_name = "qdep" * string(i)
        @info nep_name
        nep = nep_gallery(nep_name)
        A = compute_Mder(nep, 3.0)
        @test isa(nep, SPMF_NEP) # Not so sophisticated, but improves code coverage
        @test_throws MethodError nep_gallery(nep_name, 15)
        @test_throws MethodError nep_gallery(nep_name, t=15)
    end

    @info "Testing sine"
    nep=nep_gallery("sine")
    tol=1e-10
    λ0=λ=-1.4;
    v0=normalize!(compute_Mder(nep,λ0)\ones(size(nep,1)))
    λ,v=quasinewton(Float64,nep,λ=-1.4566,v=v0,tol=tol,
                    errmeasure=ResidualErrmeasure(nep),
                    logger=displaylevel,armijo_factor=0.5,armijo_max=3)
    normalize!(v)
    @test norm(compute_Mlincomb(nep,λ,v))<sqrt(tol)
    @test isa(nep,SPMFSumNEP)
    @test_throws MethodError nep_gallery("sine", 15)
    @test_throws MethodError nep_gallery("sine", t=15)

    @info "Testing dep_symm_double"
    nep=nep_gallery("dep_symm_double")
    A=compute_Mder(nep,3.0)
    @test norm(A-A')==0

    @info "Testing qep_fixed_eig"
    nep=nep_gallery("qep_fixed_eig")
    nep=nep_gallery("qep_fixed_eig",3)
    nep=nep_gallery("qep_fixed_eig",3,1:6)
    A=compute_Mder(nep,3.0)
    λstar=5.0 # this problem is constructed such that this is an eigenvalue
    @test minimum(svdvals(compute_Mder(nep,λstar)))<100*eps()
    @test_throws MethodError nep_gallery("qep_fixed_eig",3,1:6, 8)
    @test_throws MethodError nep_gallery("qep_fixed_eig",3,1:6, t=8)

    @info "Testing neuron0"
    nep=nep_gallery("neuron0")
    A=compute_Mder(nep,3.0)
    @test_throws MethodError nep_gallery("neuron0", 8)
    @test_throws MethodError nep_gallery("neuron0", t=8)

    @info "Testing schrodinger_movebc"
    n=5; λ=-3.0;
    nep=nep_gallery("schrodinger_movebc",n)
    A=compute_Mder(nep,λ)
    h=1/(n-1);
    @test A[1,1]== -2/(h^2)-λ-1 # Check (1,1)-point in discretization

    @info "Testing beam"
    nep=nep_gallery("beam")
    nep=nep_gallery("beam", 15)
    A=compute_Mder(nep,3.0)
    @test_throws MethodError nep_gallery("beam", 15, 8)
    @test_throws MethodError nep_gallery("beam", 15, t=8)

    @info "Testing hadeler"
    nep=nep_gallery("nlevp_native_hadeler")
    A=compute_Mder(nep,3.0);
    @test isreal(A)

    @info "Testing loaded_string"
    # Compare with hard-coded value
    λ=3.0+4im;
    Zref=[9.600000000000000 - 0.533333333333333im -5.100000000000000 - 0.133333333333333im  0.000000000000000 + 0.000000000000000im  0.000000000000000 + 0.000000000000000im  0.000000000000000 + 0.000000000000000im
 -5.100000000000000 - 0.133333333333333im  9.600000000000000 - 0.533333333333333im -5.100000000000000 - 0.133333333333333im  0.000000000000000 + 0.000000000000000im  0.000000000000000 + 0.000000000000000im
  0.000000000000000 + 0.000000000000000im -5.100000000000000 - 0.133333333333333im  9.600000000000000 - 0.533333333333333im -5.100000000000000 - 0.133333333333333im  0.000000000000000 + 0.000000000000000im
  0.000000000000000 + 0.000000000000000im  0.000000000000000 + 0.000000000000000im -5.100000000000000 - 0.133333333333333im  9.600000000000000 - 0.533333333333333im -5.100000000000000 - 0.133333333333333im
       0.000000000000000 + 0.000000000000000im  0.000000000000000 + 0.000000000000000im  0.000000000000000 + 0.000000000000000im -5.100000000000000 - 0.133333333333333im  9.137811900191938 - 0.880870121561100im ];
    nep=nep_gallery("nlevp_native_loaded_string",5,4,5);
    Z=compute_Mder(nep,λ);
    @test norm(Z-Zref)<1e-10

    @info "testing bem_fichera"
    nep=nep_gallery("bem_fichera",1);

    λ_ref=8.790558462139456 - 0.010815457827738698im
    M=compute_Mder(nep,λ_ref);
    @test rank(M)<size(M,1)
    @onlybench @testset "bem_fichera + IAR" begin
       v=ones(size(nep,1));
       (λ,vv)=iar(nep,σ=9,v=v,logger=displaylevel,neigs=1,maxit=50,tol=1e-6) # normally takes 30 iterations
       @test norm(compute_Mlincomb(nep,λ[1],vv[:,1]))<sqrt(eps())*100;
    end
    @info "non-existing example"
    @test_throws ErrorException nep_gallery("non-existing example")

end


text:
# wraps the solver functions
function ODEjl_wrapper(tr::TestRun)
    tc  = tr.tc
    so = tr.solver
    ###
    # 0) Wrap tc.fn!, tc.jac!, tc.mass! if necessary
    fn = nonmod_fn(tc)
    if hasjacobian(tc)
        jac = nonmod_jac(tc)
    end
    if hasmass(tc)
        error("ODE.jl solvers do not support a mass matrix")
    end

    ###
    # 1) Make call signature
    if isadaptive(so)
        args = (fn, tc.ic, tc.tspan)
        # TODO: also use h0
        kwargs  = ((:reltol, tr.reltol), (:abstol, tr.abstol), (:points, :specified))
        if hasjacobian(tc) && isimplicit(so)
            kwargs  = ((:jacobian, jac), kwargs...)
        end
    else
        args = (fn, tc.ic, tr.tsteps)
        if hasjacobian(tc) && isimplicit(tr.solver)
            kwargs  = ((:jacobian, jac),)
        else
            kwargs  = ()
        end
    end

    ###
    # 2) Call solver, if it does not succeed throw an error (if that
    # is not done anyway)
    (t,y) = so.solverfn(args...; kwargs...)

    ###
    # 3) Transform output to conform to standard:
    # tend -- end time reached
    # yend -- solution at tend
    # stats -- statistics, if available: (steps_total,steps_accepted, fn_evals, jac_evals, linear_solves)
    #                         otherwise  (-1, -1, -1, -1, -1)
    tend = t[end]
    yend = y[end]
    stats = (-1, -1, -1, -1, -1)
    return tend, yend, stats
end

import ODE
##############################################################################
#List of all ODE.jl solvers avaible for testing
##############################################################################
## Non-stiff fixed step solvers
nonstiff_fixedstep= [
           ODE.ode1,
#          ODE.ode2_midpoint,# add back in after 8 color issues is resolved
           ODE.ode2_heun,
           ODE.ode4,
           ODE.ode4ms,
           ODE.ode5ms
#          ODE.ode113 Adam Bashforth under construction
           ]
## Non-stiff adaptive step solvers
nonstiff_adaptive=[
#          ODE.ode21, # this fails on Travis with 0.4?! TODO revert once fixed.
           ODE.ode23,
           ODE.ode45,
           ODE.ode45_dp,
           ODE.ode45_fe,
           ODE.ode78,
#          ODE.ode_ab_adaptive #Adaptive Adam Bashforth under construction
           ]
## Stiff fixed-step solvers
stiff_fixedstep=[
           ODE.ode4s_s,
           ODE.ode4s_kr
           ]
## Stiff adaptive solvers
stiff_adaptive = [
           ODE.ode23s
 #          ODE.odeRadauIIA #RADAU methods under construction
           ]

ode_only = 0 # dae index
pkg = "ODE.jl"
#    ode23s = Solver{:im}(ODE.ode23s, stiff)

ODEsolvers = Dict{Any,Solver}()
sl = 1 # to make it global so it works with eval
# adaptive non-stiff solvers
for fn in nonstiff_adaptive
    sl = Solver{:ex}(fn, ODE, ODEjl_wrapper, nonstiff, adaptive, ode_only, explicit_eq)
    ODEsolvers[fn] = sl
end

# fixed step non-stiff solvers
for fn in nonstiff_fixedstep
    sl = Solver{:ex}(fn, ODE, ODEjl_wrapper, nonstiff, nonadaptive, ode_only, explicit_eq)
    ODEsolvers[fn] = sl
end

# adaptive stiff solvers
for fn in stiff_adaptive
    sl = Solver{:im}(fn, ODE, ODEjl_wrapper, stiff, adaptive, ode_only, explicit_eq)
    ODEsolvers[fn] = sl
end

# fixed step stiff solvers
for fn in stiff_fixedstep
    sl = Solver{:im}(fn, ODE, ODEjl_wrapper, stiff, nonadaptive, ode_only, explicit_eq)
    ODEsolvers[fn] = sl
end
merge!(allsolvers, ODEsolvers)


text:
using Soss

get_distname(x::Symbol) = Symbol(:_, x, :_dist)

"""
    withmeasures(m::Model) -> Model

julia> m = @model begin
    σ ~ HalfNormal()
    y ~ For(10) do j
        Normal(0,σ)
    end
end;

julia> m_dists = Soss.withmeasures(m)
@model begin
        _σ_dist = HalfNormal()
        σ ~ _σ_dist
        _y_dist = For(10) do j
                Normal(0, σ)
            end
        y ~ _y_dist
    end

julia> ydist = rand(m_dists())._y_dist
For{GeneralizedGenerated.Closure{function = (σ, M, j;) -> begin
    M.Normal(0, σ)
end,Tuple{Float64,Module}},Tuple{Int64},Normal{Float64},Float64}(GeneralizedGenerated.Closure{function = (σ, M, j;) -> begin
    M.Normal(0, σ)
end,Tuple{Float64,Module}}((0.031328640120683524, Main)), (10,))

julia> rand(ydist)
10-element Array{Float64,1}:
  0.03454891487870426
  0.008832782323408313
 -0.007395186925623771
 -0.030669004243492004
 -0.01728630026691135
  0.011892877715064682
  0.025576319363013512
 -0.029323425779917773
 -0.020502677724193594
  0.04612690097957398
"""
function withmeasures(m::Model)
    theModule = getmodule(m)
    m_init = Model(theModule, m.args, NamedTuple(), NamedTuple(), nothing)

    function proc(st::Sample)
        distname = get_distname(st.x)
        assgn = Model(theModule, Assign(distname, st.rhs))
        sampl = Model(theModule, Sample(st.x, distname))
        return merge(assgn, sampl)
    end

    proc(st::Arg) = nothing
    proc(st) = Model(theModule, st)

    # Rewrite the statements of the model one by one.
    m_new = foldl(statements(m); init=m_init) do m0,st
        merge(m0, proc(st))
    end
    return m_new
end

function withmeasures(d::ConditionalModel)
    withmeasures(Model(d))(argvals(d)) | observations(d)
end

# TODO: Finish this
# function predict_measure(rng::AbstractRNG, d::ConditionalModel, post::AbstractVector{<:NamedTuple{N}}) where {N}


text:
using Test
using LinearAlgebra
using NonNegLeastSquares
using SparseArrays

#test specific
using Random
using PyCall
const pyopt = pyimport_conda("scipy.optimize", "scipy")

function test_algorithm(fh, ε=1e-5)
	# Solve A*x = b for x, subject to x >=0
	A = [ 0.53879488  0.65816267
	      0.12873446  0.98669198
	      0.24555042  0.00598804
	      0.80491791  0.32793762 ]

	b = [0.888,  0.562,  0.255,  0.077]
	x = [0.15512102, 0.69328985] # approx solution from scipy
	@test norm(fh(A,b)-x) < ε

	## A second test case
	A2 = [ -0.24  -0.82   1.35   0.36   0.35
	       -0.53  -0.20  -0.76   0.98  -0.54
	        0.22   1.25  -1.60  -1.37  -1.94
	       -0.51  -0.56  -0.08   0.96   0.46
	        0.48  -2.25   0.38   0.06  -1.29 ]
	b2 = [-1.6,  0.19,  0.17,  0.31, -1.27]
	x2 = [2.2010416, 1.19009924, 0.0, 1.55001345, 0.0]
	@test norm(fh(A2,b2)-x2) < ε

	## Test a bunch of random cases
	for i = 1:100
		m,n = rand(1:10),rand(1:10)
		A3 = randn(m,n)
		b3 = randn(m)
		x3,resid = pyopt.nnls(A3,b3)
		if resid > ε
	        @test norm(fh(A3,b3)-x3) < ε
	    else
	    	@test norm(A3*fh(A3,b3)-b3) < ε
	    end
	end
end

nnls(A,b) = nonneg_lsq(A, b; alg=:nnls)
nnls_gram(A,b) = nonneg_lsq(A'*A, A'*b; alg=:nnls, gram=true)
fnnls(A,b) = nonneg_lsq(A, b; alg=:fnnls)
fnnls_gram(A,b) = nonneg_lsq(A'*A, A'*b; alg=:fnnls, gram=true)
pivot(A,b) = nonneg_lsq(A, b; alg=:pivot)
pivot_comb(A,b) = nonneg_lsq(A, b; alg=:pivot, variant=:comb)
pivot_cache(A,b) = nonneg_lsq(A, b; alg=:pivot, variant=:cache)

algs = [nnls, nnls_gram, fnnls, fnnls_gram, pivot, pivot_comb, pivot_cache]
errs = fill(1e-5, length(algs))

for (f, ε) in zip(algs, errs)
 	print("testing ")
 	@show f
 	test_algorithm(f, ε)
 	println("done")
end

@testset "NNLS" begin include("nnls_test.jl") end
@testset "FNNLS" begin include("fnnls_test.jl") end
@testset "Pivot" begin include("pivot_test.jl") end
@testset "Sparse" begin include("sparse_test.jl") end
@testset "CoordinateDescent" begin include("cd_test.jl") end


text:
# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
# licensing details.
"""
    connected_components!(label, g)

Fill `label` with the `id` of the connected component in the undirected graph
`g` to which it belongs. Return a vector representing the component assigned
to each vertex. The component value is the smallest vertex ID in the component.

### Performance
This algorithm is linear in the number of edges of the graph.
"""
function connected_components!(label::AbstractVector, g::AbstractGraph{T}) where T
    nvg = nv(g)

    for u in vertices(g)
        label[u] != zero(T) && continue
        label[u] = u
        Q = Vector{T}()
        push!(Q, u)
        while !isempty(Q)
            src = popfirst!(Q)
            for vertex in all_neighbors(g, src)
                if label[vertex] == zero(T)
                    push!(Q, vertex)
                    label[vertex] = u
                end
            end
        end
    end
    return label
end

"""
    components_dict(labels)

Convert an array of labels to a map of component id to vertices, and return
a map with each key corresponding to a given component id
and each value containing the vertices associated with that component.
"""
function components_dict(labels::Vector{T}) where T <: Integer
    d = Dict{T,Vector{T}}()
    for (v, l) in enumerate(labels)
        vec = get(d, l, Vector{T}())
        push!(vec, v)
        d[l] = vec
    end
    return d
end

"""
    components(labels)

Given a vector of component labels, return a vector of vectors representing the vertices associated
with a given component id.
"""
function components(labels::Vector{T}) where T <: Integer
    d = Dict{T,T}()
    c = Vector{Vector{T}}()
    i = one(T)
    for (v, l) in enumerate(labels)
        index = get!(d, l, i)
        if length(c) >= index
            push!(c[index], v)
        else
            push!(c, [v])
            i += 1
        end
    end
    return c, d
end

"""
    connected_components(g)

Return the [connected components](https://en.wikipedia.org/wiki/Connectivity_(graph_theory))
of an undirected graph `g` as a vector of components, with each element a vector of vertices
belonging to the component.

For directed graphs, see [`strongly_connected_components`](@ref) and
[`weakly_connected_components`](@ref).

# Examples
```jldoctest
julia> g = SimpleGraph([0 1 0; 1 0 1; 0 1 0]);

julia> connected_components(g)
1-element Array{Array{Int64,1},1}:
 [1, 2, 3]

julia> g = SimpleGraph([0 1 0 0 0; 1 0 1 0 0; 0 1 0 0 0; 0 0 0 0 1; 0 0 0 1 0]);

julia> connected_components(g)
2-element Array{Array{Int64,1},1}:
 [1, 2, 3]
 [4, 5]
```
"""
function connected_components(g::AbstractGraph{T}) where T
    label = zeros(T, nv(g))
    connected_components!(label, g)
    c, d = components(label)
    return c
end

"""
    is_connected(g)

Return `true` if graph `g` is connected. For directed graphs, return `true`
if graph `g` is weakly connected.

# Examples
```jldoctest
julia> g = SimpleGraph([0 1 0; 1 0 1; 0 1 0]);

julia> is_connected(g)
true

julia> g = SimpleGraph([0 1 0 0 0; 1 0 1 0 0; 0 1 0 0 0; 0 0 0 0 1; 0 0 0 1 0]);

julia> is_connected(g)
false

julia> g = SimpleDiGraph([0 1 0; 0 0 1; 1 0 0]);

julia> is_connected(g)
true
```
"""
function is_connected(g::AbstractGraph)
    mult = is_directed(g) ? 2 : 1
    return mult * ne(g) + 1 >= nv(g) && length(connected_components(g)) == 1
end

"""
    weakly_connected_components(g)

Return the weakly connected components of the graph `g`. This
is equivalent to the connected components of the undirected equivalent of `g`.
For undirected graphs this is equivalent to the [`connected_components`](@ref) of `g`.

# Examples
```jldoctest
julia> g = SimpleDiGraph([0 1 0; 1 0 1; 0 0 0]);

julia> weakly_connected_components(g)
1-element Array{Array{Int64,1},1}:
 [1, 2, 3]
```
"""
weakly_connected_components(g) = connected_components(g)

"""
    is_weakly_connected(g)

Return `true` if the graph `g` is weakly connected. If `g` is undirected,
this function is equivalent to [`is_connected(g)`](@ref).

# Examples
```jldoctest
julia> g = SimpleDiGraph([0 1 0; 0 0 1; 1 0 0]);

julia> is_weakly_connected(g)
true

julia> g = SimpleDiGraph([0 1 0; 1 0 1; 0 0 0]);

julia> is_connected(g)
true

julia> is_strongly_connected(g)
false

julia> is_weakly_connected(g)
true
```
"""
is_weakly_connected(g) = is_connected(g)

"""
    strongly_connected_components(g)

Compute the strongly connected components of a directed graph `g`.

Return an array of arrays, each of which is the entire connected component.

### Implementation Notes
The order of the components is not part of the API contract.

# Examples
```jldoctest
julia> g = SimpleDiGraph([0 1 0; 1 0 1; 0 0 0]);

julia> strongly_connected_components(g)
2-element Array{Array{Int64,1},1}:
 [3]
 [1, 2]
```
"""
function strongly_connected_components end
# see https://github.com/mauro3/SimpleTraits.jl/issues/47#issuecomment-327880153 for syntax
@traitfn function strongly_connected_components(g::AG::IsDirected) where {T, AG <: AbstractGraph{T}}
    zero_t = zero(T)
    one_t = one(T)
    nvg = nv(g)
    count = one_t
    index = zeros(T, nvg)         # first time in which vertex is discovered
    stack = Vector{T}()           # stores vertices which have been discovered and not yet assigned to any component
    onstack = zeros(Bool, nvg)    # false if a vertex is waiting in the stack to receive a component assignment
    lowlink = zeros(T, nvg)       # lowest index vertex that it can reach through back edge (index array not vertex id number)
    parents = zeros(T, nvg)       # parent of every vertex in dfs
    components = Vector{Vector{T}}()    # maintains a list of scc (order is not guaranteed in API)
    sizehint!(stack, nvg)
    sizehint!(components, nvg)

    for s in vertices(g)
        if index[s] == zero_t
            index[s] = count
            lowlink[s] = count
            onstack[s] = true
            parents[s] = s
            push!(stack, s)
            count = count + one_t

            # start dfs from 's'
            dfs_stack = Vector{T}([s])
            while !isempty(dfs_stack)
                v = dfs_stack[end] #end is the most recently added item
                u = zero_t
                for n in outneighbors(g, v)
                    if index[n] == zero_t
                        # unvisited neighbor found
                        u = n
                        break
                        #GOTO A push u onto DFS stack and continue DFS
                    elseif onstack[n]
                        # we have already seen n, but can update the lowlink of v
                        # which has the effect of possibly keeping v on the stack until n is ready to pop.
                        # update lowest index 'v' can reach through out neighbors
                        lowlink[v] = min(lowlink[v], index[n])
                    end
                end
                if u == zero_t
                    # All out neighbors already visited or no out neighbors
                    # we have fully explored the DFS tree from v.
                    # time to start popping.
                    popped = pop!(dfs_stack)
                    lowlink[parents[popped]] = min(lowlink[parents[popped]], lowlink[popped])
                    if index[v] == lowlink[v]
                        # found a cycle in a completed dfs tree.
                        component = Vector{T}()
                        sizehint!(component, length(stack))
                        while !isempty(stack) #break when popped == v
                            # drain stack until we see v.
                            # everything on the stack until we see v is in the SCC rooted at v.
                            popped = pop!(stack)
                            push!(component, popped)
                            onstack[popped] = false
                            # popped has been assigned a component, so we will never see it again.
                            if popped == v
                                # we have drained the stack of an entire component.
                                break
                            end
                        end
                        push!(components, reverse(component))
                    end
                else #LABEL A
                    # add unvisited neighbor to dfs
                    index[u] = count
                    lowlink[u] = count
                    onstack[u] = true
                    parents[u] = v
                    count = count + one_t
                    push!(stack, u)
                    push!(dfs_stack, u)
                    # next iteration of while loop will expand the DFS tree from u.
                end
            end
        end
    end

    return components
end

"""
    is_strongly_connected(g)

Return `true` if directed graph `g` is strongly connected.

# Examples
```jldoctest
julia> g = SimpleDiGraph([0 1 0; 0 0 1; 1 0 0]);

julia> is_strongly_connected(g)
true
```
"""
function is_strongly_connected end
@traitfn is_strongly_connected(g::::IsDirected) = length(strongly_connected_components(g)) == 1

"""
    period(g)

Return the (common) period for all vertices in a strongly connected directed graph.
Will throw an error if the graph is not strongly connected.

# Examples
```jldoctest
julia> g = SimpleDiGraph([0 1 0; 0 0 1; 1 0 0]);

julia> period(g)
3
```
"""
function period end
# see https://github.com/mauro3/SimpleTraits.jl/issues/47#issuecomment-327880153 for syntax
@traitfn function period(g::AG::IsDirected) where {T, AG <: AbstractGraph{T}}
    !is_strongly_connected(g) && throw(ArgumentError("Graph must be strongly connected"))

    # First check if there's a self loop
    has_self_loops(g) && return 1

    g_bfs_tree  = bfs_tree(g, 1)
    levels      = gdistances(g_bfs_tree, 1)
    tree_diff   = difference(g, g_bfs_tree)
    edge_values = Vector{T}()

    divisor = 0
    for e in edges(tree_diff)
        @inbounds value = levels[src(e)] - levels[dst(e)] + 1
        divisor = gcd(divisor, value)
        isequal(divisor, 1) && return 1
    end

    return divisor
end

"""
    condensation(g[, scc])

Return the condensation graph of the strongly connected components `scc`
in the directed graph `g`. If `scc` is missing, generate the strongly
connected components first.

# Examples
```jldoctest
julia> g = SimpleDiGraph([0 1 0 0 0; 0 0 1 0 0; 1 0 0 1 0; 0 0 0 0 1; 0 0 0 1 0])
{5, 6} directed simple Int64 graph

julia> strongly_connected_components(g)
2-element Array{Array{Int64,1},1}:
 [4, 5]
 [1, 2, 3]

julia> foreach(println, edges(condensation(g)))
Edge 2 => 1
```
"""
function condensation end
@traitfn function condensation(g::::IsDirected, scc::Vector{Vector{T}}) where T <: Integer
    h = DiGraph{T}(length(scc))

    component = Vector{T}(undef, nv(g))

    for (i, s) in enumerate(scc)
        @inbounds component[s] .= i
    end

    @inbounds for e in edges(g)
        s, d = component[src(e)], component[dst(e)]
        if (s != d)
            add_edge!(h, s, d)
        end
    end
    return h
end
@traitfn condensation(g::::IsDirected) = condensation(g, strongly_connected_components(g))

"""
    attracting_components(g)

Return a vector of vectors of integers representing lists of attracting
components in the directed graph `g`.

The attracting components are a subset of the strongly
connected components in which the components do not have any leaving edges.

# Examples
```jldoctest
julia> g = SimpleDiGraph([0 1 0 0 0; 0 0 1 0 0; 1 0 0 1 0; 0 0 0 0 1; 0 0 0 1 0])
{5, 6} directed simple Int64 graph

julia> strongly_connected_components(g)
2-element Array{Array{Int64,1},1}:
 [4, 5]
 [1, 2, 3]

julia> attracting_components(g)
1-element Array{Array{Int64,1},1}:
 [4, 5]
```
"""
function attracting_components end
# see https://github.com/mauro3/SimpleTraits.jl/issues/47#issuecomment-327880153 for syntax
@traitfn function attracting_components(g::AG::IsDirected) where {T, AG <: AbstractGraph{T}}
    scc  = strongly_connected_components(g)
    cond = condensation(g, scc)

    attracting = Vector{T}()

    for v in vertices(cond)
        if outdegree(cond, v) == 0
            push!(attracting, v)
        end
    end
    return scc[attracting]
end

"""
    neighborhood(g, v, d, distmx=weights(g))

Return a vector of each vertex in `g` at a geodesic distance less than or equal to `d`, where distances
may be specified by `distmx`.

### Optional Arguments
- `dir=:out`: If `g` is directed, this argument specifies the edge direction
with respect to `v` of the edges to be considered. Possible values: `:in` or `:out`.

# Examples
```jldoctest
julia> g = SimpleDiGraph([0 1 0 0 0; 0 0 1 0 0; 1 0 0 1 0; 0 0 0 0 1; 0 0 0 1 0]);

julia> neighborhood(g, 1, 2)
3-element Array{Int64,1}:
 1
 2
 3

julia> neighborhood(g, 1, 3)
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> neighborhood(g, 1, 3, [0 1 0 0 0; 0 0 1 0 0; 1 0 0 0.25 0; 0 0 0 0 0.25; 0 0 0 0.25 0])
5-element Array{Int64,1}:
 1
 2
 3
 4
 5
```
"""
neighborhood(g::AbstractGraph{T}, v::Integer, d, distmx::AbstractMatrix{U}=weights(g); dir=:out) where T <: Integer where U <: Real =
    first.(neighborhood_dists(g, v, d, distmx; dir=dir))

"""
    neighborhood_dists(g, v, d, distmx=weights(g))

Return a a vector of tuples representing each vertex which is at a geodesic distance less than or equal to `d`, along with
its distance from `v`. Non-negative distances may be specified by `distmx`.

### Optional Arguments
- `dir=:out`: If `g` is directed, this argument specifies the edge direction
with respect to `v` of the edges to be considered. Possible values: `:in` or `:out`.

# Examples
```jldoctest
julia> g = SimpleDiGraph([0 1 0 0 0; 0 0 1 0 0; 1 0 0 1 0; 0 0 0 0 1; 0 0 0 1 0]);

julia> neighborhood_dists(g, 1, 3)
4-element Array{Tuple{Int64,Int64},1}:
 (1, 0)
 (2, 1)
 (3, 2)
 (4, 3)

julia> neighborhood_dists(g, 1, 3, [0 1 0 0 0; 0 0 1 0 0; 1 0 0 0.25 0; 0 0 0 0 0.25; 0 0 0 0.25 0])
5-element Array{Tuple{Int64,Float64},1}:
 (1, 0.0)
 (2, 1.0)
 (3, 2.0)
 (4, 2.25)
 (5, 2.5)

julia> neighborhood_dists(g, 4, 3)
2-element Array{Tuple{Int64,Int64},1}:
 (4, 0)
 (5, 1)

julia> neighborhood_dists(g, 4, 3, dir=:in)
5-element Array{Tuple{Int64,Int64},1}:
 (4, 0)
 (3, 1)
 (5, 1)
 (2, 2)
 (1, 3)
```
"""
neighborhood_dists(g::AbstractGraph{T}, v::Integer, d, distmx::AbstractMatrix{U}=weights(g); dir=:out) where T <: Integer where U <: Real =
    (dir == :out) ? _neighborhood(g, v, d, distmx, outneighbors) : _neighborhood(g, v, d, distmx, inneighbors)

function _neighborhood(g::AbstractGraph{T}, v::Integer, d::Real, distmx::AbstractMatrix{U}, neighborfn::Function) where T <: Integer where U <: Real
    Q = Vector{Tuple{T,U}}()
    d < zero(U) && return Q
    push!(Q, (v,zero(U),) )
    seen = fill(false,nv(g)); seen[v] = true #Bool Vector benchmarks faster than BitArray
    for (src,currdist) in Q
        currdist >= d && continue
        for dst in neighborfn(g,src)
            if !seen[dst]
                seen[dst]=true
                if currdist+distmx[src,dst] <= d
                    push!(Q, (dst , currdist+distmx[src,dst],))
                end
            end
        end
    end
    return Q
end

"""
    isgraphical(degs)

Return true if the degree sequence `degs` is graphical.
A sequence of integers is called graphical, if there exists a graph where the degrees of its vertices form that same sequence.

### Performance
Time complexity: ``\\mathcal{O}(|degs|*\\log(|degs|))``.

### Implementation Notes
According to Erdös-Gallai theorem, a degree sequence ``\\{d_1, ...,d_n\\}`` (sorted in descending order) is graphic iff the sum of vertex degrees is even and the sequence obeys the property -
```math
\\sum_{i=1}^{r} d_i \\leq r(r-1) + \\sum_{i=r+1}^n min(r,d_i)
```
for each integer r <= n-1
"""
function isgraphical(degs::Vector{<:Integer})
    iseven(sum(degs)) || return false
    sorted_degs = sort(degs, rev = true)
    n = length(sorted_degs)
    cur_sum = zero(UInt64)
    mindeg = Vector{UInt64}(undef, n)
    @inbounds for i = 1:n
        mindeg[i] = min(i, sorted_degs[i])
    end
    cum_min = sum(mindeg)
    @inbounds for r = 1:(n - 1)
        cur_sum += sorted_degs[r]
        cum_min -= mindeg[r]
        cond = cur_sum <= (r * (r - 1) + cum_min)
        cond || return false
    end
    return true
end


text:
using BenchmarkTools
using Random
using AugmentedGaussianProcesses

const SUITE = BenchmarkGroup()

# SUITE["Kernel"] = BenchmarkGroup()
# include("kernel.jl")
SUITE["Models"] = BenchmarkGroup()
include("models.jl")


text:
let
    # fitting noisy data to an exponential model
    # TODO: Change to `.-x` when 0.5 support is dropped
    @. model(x, p) = p[1] * exp(-x * p[2])
    model_inplace(F, x, p) = (@. F = p[1] * exp(-x * p[2]))

    # some example data
    Random.seed!(12345)
    xdata = range(0, stop=10, length=500000)
    ydata = model(xdata, [1.0, 2.0]) + 0.01*randn(length(xdata))
    p0 = [0.5, 0.5]

    # if your model is differentiable, it can be faster and/or more accurate
    # to supply your own jacobian instead of using the finite difference
    function jacobian_model(x,p)
        J = Array{Float64}(undef, length(x), length(p))
        @. J[:,1] = exp(-x*p[2])     #dmodel/dp[1]
        @. @views J[:,2] = -x*p[1]*J[:,1]
        J
    end

    function jacobian_model_inplace(J::Array{Float64,2},x,p)
        @. J[:,1] = exp(-x*p[2])     #dmodel/dp[1]
        @. @views J[:,2] = -x*p[1]*J[:,1]
    end

    f(p) = model(xdata, p) - ydata
    g(p) = jacobian_model(xdata, p)
    df = OnceDifferentiable(f, g, p0, similar(ydata); inplace = false);
    evalf(x) = NLSolversBase.value!!(df, x)
    evalg(x) = NLSolversBase.jacobian!!(df, x)
    r = evalf(p0);
    j = evalg(p0);

    f_inplace = (F,p) -> (model_inplace(F, xdata, p); @. F = F - ydata)
    g_inplace = (G,p) -> jacobian_model_inplace(G, xdata, p)
    df_inplace = OnceDifferentiable(f_inplace, g_inplace, p0, similar(ydata); inplace = true);
    evalf_inplace(x) = NLSolversBase.value!!(df_inplace,x)
    evalg_inplace(x) = NLSolversBase.jacobian!!(df_inplace,x)
    r_inplace = evalf_inplace(p0);
    j_inplace = evalg_inplace(p0);

    @test r == r_inplace
    @test j == j_inplace

    println("--------------\nPerformance of non-inplace")
    println("\t Evaluation function")

    stop = 8 #8 because the tests afterwards will call the eval function 8 or 9 times, so it makes it easy to compare
    step = 1

    @time for i=range(0,stop=stop,step=step)
        evalf(p0);
    end

    println("\t Jacobian function")
    @time for i=range(0,stop=stop,step=step)
        evalg(p0);
    end

    println("--------------\nPerformance of inplace")
    println("\t Evaluation function")
    @time for i=range(0,stop=stop,step=step)
        evalf_inplace(p0);
    end

    println("\t Jacobian function")
    @time for i=range(0,stop=stop,step=step)
        evalg_inplace(p0);
    end

    curve_fit(model, xdata, ydata, p0; maxIter=100); #warmup
    curve_fit(model_inplace, xdata, ydata, p0; inplace = true, maxIter=100);

    #explicit jac
    curve_fit(model, jacobian_model, xdata, ydata, p0; maxIter=100);
    curve_fit(model_inplace, jacobian_model_inplace, xdata, ydata, p0; inplace = true, maxIter=100);

    println("--------------\nPerformance of curve_fit")

    println("\t Non-inplace")
    fit = @time curve_fit(model, xdata, ydata, p0; maxIter=100)
    @test fit.converged

    println("\t Inplace")
    fit_inplace = @time curve_fit(model_inplace, xdata, ydata, p0; inplace = true, maxIter=100)
    @test fit_inplace.converged

    @test fit_inplace.param == fit.param

    println("\t Non-inplace with jacobian")
    fit_jac = @time curve_fit(model, jacobian_model, xdata, ydata, p0; maxIter=100)
    @test fit_jac.converged

    println("\t Inplace with jacobian")
    fit_inplace_jac = @time curve_fit(model_inplace, jacobian_model_inplace, xdata, ydata, p0; inplace = true, maxIter=100)
    @test fit_inplace_jac.converged

    @test fit_jac.param == fit_inplace_jac.param

    # some example data
    yvars = 1e-6*rand(length(xdata))
    ydata = model(xdata, [1.0, 2.0]) + sqrt.(yvars) .* randn(length(xdata))

    println("--------------\nPerformance of curve_fit with weights")

    curve_fit(model, xdata, ydata, 1 ./ yvars, [0.5, 0.5]);
    curve_fit(model_inplace, xdata, ydata, 1 ./ yvars, [0.5, 0.5]; inplace = true, maxIter=100);

    curve_fit(model, jacobian_model, xdata, ydata, 1 ./ yvars, [0.5, 0.5]);
    curve_fit(model_inplace, jacobian_model_inplace, xdata, ydata, 1 ./ yvars, [0.5, 0.5]; inplace = true, maxIter=100);

    println("\t Non-inplace with weights")
    fit_wt = @time curve_fit(model, jacobian_model, xdata, ydata, 1 ./ yvars, [0.5, 0.5]; maxIter=100)
    @test fit_wt.converged

    println("\t Inplace with weights")
    fit_inplace_wt = @time curve_fit(model_inplace, xdata, ydata, 1 ./ yvars, [0.5, 0.5]; inplace = true, maxIter=100)
    @test fit_inplace_wt.converged

    @test maximum(abs.(fit_wt.param - fit_inplace_wt.param)) < 1e-15

    println("\t Non-inplace with jacobian with weights")
    fit_wt_jac = @time curve_fit(model, jacobian_model, xdata, ydata, 1 ./ yvars, [0.5, 0.5]; maxIter=100)
    @test fit_wt_jac.converged

    println("\t Inplace with jacobian with weights")
    fit_inplace_wt_jac = @time curve_fit(model_inplace, jacobian_model_inplace, xdata, ydata, 1 ./ yvars, [0.5, 0.5]; inplace = true, maxIter=100)
    @test fit_inplace_wt_jac.converged

    @test maximum(abs.(fit_wt_jac.param - fit_inplace_wt_jac.param)) < 1e-15

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, Pkg

name = "SpeexDSP"
version = v"1.2.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://ftp.osuosl.org/pub/xiph/releases/speex/speexdsp-1.2.0.tar.gz", "682042fc6f9bee6294ec453f470dadc26c6ff29b9c9e9ad2ffc1f4312fd64771")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd speexdsp-1.2.0
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make
make install
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

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

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

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using Bio.Align
using Glob
using GZip

function coveragemap(reader, chrom, region)
    cov = zeros(Int, length(region))

    aln = Bio.Align.Alignment("")

    for rec in intersect(reader, chrom, region)
        if !ismapped(rec)
            continue
        end

        aln = alignment(rec)
        for i in 1:seqlength(rec)
            j, op = seq2ref(aln, i)
            if ismatchop(op) && j in region
                cov[j - first(region) + 1] += 1
            end
        end
    end
    return cov
end

function batchOpenBamReader(bamFiles)
    map(bf -> open(BAMReader, bf, index=bf * ".bai"), bamFiles)
end

function batchCloseBamReader(bamReaders)
    map(x -> close(x), bamReaders)
end

# function getCoverage(bamReaders, chrom, region)
#     mat = zeros(length(bamReaders), length(region))
#     for i in 1:length(bamReaders)
#         mat[i,:] = coveragemap(bamReaders[i], chrom, region)
#     end
#     return mat
# end

function getNormalizedCoverageSum(bamReaders, chrom, region, normFac)
    arr = zeros(length(region))
    cov = zeros(Int, length(region))
    for i in 1:length(bamReaders)
        cov = coveragemap(bamReaders[i], chrom, region)
        for j in 1:length(cov)
            arr[j] +=  cov[j] / normFac[i]
        end
    end
    return arr
end

function getNormalizedCoverageBase(bamReaders, chrom, region, normFac)
    mat = zeros(length(bamReaders), length(region))
    for i in 1:length(bamReaders)
        mat[i,:] = coveragemap(bamReaders[i], chrom, region) / normFac[i]
    end
    return vec(sum(mat, 1))
end

function parsePathBmFiles(pathBamFiles)
    bamFiles = Array{Any,1}()
    normFac = Array{Float64,1}()

    open(pathBamFiles) do f
        readline(f)
        for line in eachline(f)
            arr = split(chomp(line), "\t")
            push!(bamFiles, arr[1])
            push!(normFac, parse(Float64, arr[2]))
        end
    end

    return bamFiles, normFac
end

function parsePathRegions(pathRegions)
    I_CHROM = 1
    I_LPOS = 2
    I_RPOS = 3
    I_whichIsExonAnchor = 5
    I_NAME = 4
    I_REF_LPOS = 8
    I_REF_RPOS = 9

    ids = Array{Any,1}()
    chroms = Array{Any,1}()
    regions = Array{UnitRange{Int64},1}()
    rsss = Array{Int64,1}()

    rss = 0

    open(pathRegions) do f
        readline(f)
        for line in eachline(f)
            arr = split(chomp(line), "\t")

            if arr[I_whichIsExonAnchor] == "1"
                rss = parse(Int,arr[I_RPOS])
            else
                rss = parse(Int,arr[I_LPOS]) + 1
            end

            region = (parse(Int,arr[I_REF_LPOS])+1):(parse(Int,arr[I_REF_RPOS]))
            id = arr[I_NAME]

            push!(ids, id)
            push!(regions, region)
            push!(rsss, rss)
            push!(chroms, arr[I_CHROM])
        end
    end

    return ids, chroms, regions, rsss
end

function writeResult(ofile, region, cov, rss)
    RSS = ""
    GZip.open(ofile, "w") do fw
        write(fw, "bp\tdepth\tRSS\n")
        for i in 1:length(region)
            if region[i] < rss
                RSS = "TRUE"
            else
                RSS = "FALSE"
            end

            x = @sprintf "%d\t%.5f\t%s\n" region[i] cov[i] RSS
            write(fw, x)
        end
    end
end

function getNormalizedCoverage(pathBamFiles, pathRegions, odirPrefix)

    bamFiles, normFac = parsePathBmFiles(pathBamFiles)
    println(bamFiles)
    ids, chroms, regions, rsss = parsePathRegions(pathRegions)

    bamReaders = batchOpenBamReader(bamFiles)

    for i in 1:length(ids)
        println((@sprintf "Processing %s" ids[i]))

        id = ids[i]
        chrom = chroms[i]
        region = regions[i]
        rss = rsss[i]

        cov = getNormalizedCoverageSum(bamReaders, chrom, region, normFac)

        ofile = @sprintf "%s/%s.cov.txt.gz" odirPrefix id

        writeResult(ofile, region, cov, rss)
    end

    batchCloseBamReader(bamReaders)
end


text:
using Test

include("pascals-triangle.jl")

const rows = [[1], [1, 1], [1, 2, 1], [1, 3, 3, 1], [1, 4, 6, 4, 1], [1, 5, 10, 10, 5, 1], [1, 6, 15, 20, 15, 6, 1], [1, 7, 21, 35, 35, 21, 7, 1], [1, 8, 28, 56, 70, 56, 28, 8, 1], [1, 9, 36, 84, 126, 126, 84, 36, 9, 1]]

@testset "$i row(s)" for i in eachindex(rows)
    @test triangle(i) == @view rows[1:i]
end

@testset "special cases" begin
    @test_throws DomainError triangle(-1)
    @test isempty(triangle(0))
end


text:
struct UniformWave{N,T<:AbstractFloat} <: Wave{N,T}
    c::T    # constant wave speed
    fmax::T # max supported frequency
    dx::T   # spatial step (c*dt/λ)
    dt::T   # time step (dx*λ/c)
    λ::T    # courant number c*dt/dx
end

function (::Type{UniformWave{N}})(::Type{T}=Float64; c=340, fmax=-1, lambda=-1, dt=-1, dx=-1) where {N, T <: AbstractFloat}
    λ = T(lambda == -1 ? √(1/N) : lambda)
    if fmax == -1 && dt == -1 && dx == -1
        dₜ = T(λ / (2*2e3))
        dₓ = T(c * dₜ / λ)
    elseif dt == -1 && dx == -1
        dₜ = T(λ / (2*fmax))
        dₓ = T(c * dₜ / λ)
    elseif dt != -1 && dx == -1
        dₜ = T(dt)
        dₓ = T(c * dₜ / λ)
    elseif dt == -1 && dx != -1
        dₓ = T(dx)
        dₜ = T(dₓ * λ / c)
    elseif dt != -1 && dx != -1
        dₜ = T(dt)
        dₓ = T(dx)
        λ = T(c * dₜ / dₓ)
    end
    fₘₐₓ = T(λ / (2*dₜ))
    fmax == -1 || fₘₐₓ ≈ fmax || warn("Overwriting specified fₘₐₓ = $(round(fmax/1000,2)) kHz with implied fₘₐₓ = $(round(fₘₐₓ/1000,2)) kHz")
    UniformWave{N,T}(c, fₘₐₓ, dₓ, dₜ, λ)
end


text:
"""
        data_path(region::String="GER_1")
- Returning AbstractString where the example data is stored for regions:
- `"GER_1"`: Germany 1 node
- `"GER_18"`: Germany 18 nodes
- `"CA_1"`: California 1 node
- `"CA_14"`: California 14 nodes
- `"TX_1"`: Texas 1 node
"""
function data_path(region::String="GER_1")
    # Check for existance of the region in data
    region in readdir(normpath(joinpath(@__DIR__,"..","..","data"))) || throw(error("The region $region is not found. The provided regions are: GER_1: Germany 1 node, GER_18: Germany 18 nodes, CA_1: California 1 node, CA_14: California 14 nodes, TX_1: Texas 1 node"))
    # Generate the data path based on application and region
    data_path=normpath(joinpath(@__DIR__,"..","..","data",region,"TS"))
    return data_path
end


text:
module Utils

    export @define, @big

    include("utils/macros.jl")

    export compensated_summation, compensated_summation!

    include("utils/summation.jl")

    export L2norm, l2norm

    include("utils/norms.jl")

    export relative_l2_error, relative_maximum_error

    include("utils/errors.jl")

end


text:
import Base: view

using Base: ScalarIndex, ViewIndex, Slice, @_inline_meta, @boundscheck,
            to_indices, compute_offset1, unsafe_length, _maybe_reshape_parent, index_ndims

struct Contiguous end
struct NonContiguous end

# Detect whether the view is contiguous or not
CuIndexStyle() = Contiguous()
CuIndexStyle(I...) = NonContiguous()
CuIndexStyle(i1::Colon, ::ScalarIndex...) = Contiguous()
CuIndexStyle(i1::AbstractUnitRange, ::ScalarIndex...) = Contiguous()
CuIndexStyle(i1::Colon, I...) = CuIndexStyle(I...)

cuviewlength() = ()
cuviewlength(::Real, I...) = (@_inline_meta; cuviewlength(I...)) # skip scalars
cuviewlength(i1::AbstractUnitRange, I...) = (@_inline_meta; (unsafe_length(i1), cuviewlength(I...)...))
cuviewlength(i1::AbstractUnitRange, ::ScalarIndex...) = (@_inline_meta; (unsafe_length(i1),))

view(A::CuArray, I::Vararg{Any,N}) where {N} = (@_inline_meta; _cuview(A, I, CuIndexStyle(I...)))

function _cuview(A, I, ::Contiguous)
    @_inline_meta
    J = to_indices(A, I)
    @boundscheck checkbounds(A, J...)
    _cuview(_maybe_reshape_parent(A, index_ndims(J...)), J, cuviewlength(J...))
end

# for contiguous views just return a new CuArray
_cuview(A::CuArray{T}, I::NTuple{N,ViewIndex}, dims::NTuple{M,Integer}) where {T,N,M} =
    CuArray{T,M}(A.buf, dims; offset=A.offset + compute_offset1(A, 1, I) * sizeof(T), own=A.own)

# fallback to SubArray when the view is not contiguous
_cuview(A, I, ::NonContiguous) where {N} = invoke(view, Tuple{AbstractArray, typeof(I).parameters...}, A, I...)


text:
export CoolProp_parameters, CoolProp_fluids;
"""
# CoolProp parameters table, to build run `CoolProp.buildparameters()`

$(isfile(abspath(@__FILE__, "..", "parameters.table")) ? readstring(abspath(@__FILE__, "..", "parameters.table")) : "")
"""
const CoolProp_parameters = "Type `?CoolProp_arameters` to get a list of all CoolProp parameters."
buildparameters() = begin
    logf = open("parameters.table", "w");
    println(logf, "Paramerer |Description |Unit |Comment ");
    println(logf, ":---------|:-----------|:----|:-------" );
    counter = 0;
    longunits = Set();
    for p in coolpropparameters
        longunit = get_parameter_information_string(p, "long") * " | " * get_parameter_information_string(p, "units");
        note = "";
        if (!in(longunit, longunits))
            push!(longunits, longunit);
        else
            note = " *Duplicated* "
        end
        for fluid in coolpropfluids
            try
                res = ("$(PropsSI(p, fluid))");
                note *= " **Constant Property** "
            break;
            catch err
            end
        end
        println(logf, "$p" * " | " * longunit * " | " * note);
    end
    close(logf);
end
"""
# CoolProp fluids table, to build run `CoolProp.buildfluids()`

$(isfile(abspath(@__FILE__, "..", "fluids.table")) ? readstring(abspath(@__FILE__, "..", "fluids.table")) : "")
"""
const CoolProp_fluids = "Type `?CoolProp_fluids` to get a list of all CoolProp fluids."
buildfluids() = begin
    logf = open("fluids.table", "w");
    println(logf, "ID |Name |Alias |CAS |Pure |Formula |BibTeX ");
    println(logf, ":--|:----|:-----|:---|:----|:-------|:------");
    id = 0;
    for fluid in coolpropfluids
        id+=1;
        print(logf, "$id | $fluid | $(get_fluid_param_string(fluid, "aliases"))");
        print(logf, " | $(get_fluid_param_string(fluid, "CAS"))");
        pure = get_fluid_param_string(fluid, "pure");
        print(logf, " | $pure");
        print(logf, " | $(get_fluid_param_string(fluid, "formula")) | ");
        for bi in ["BibTeX-CONDUCTIVITY", "BibTeX-EOS", "BibTeX-CP0", "BibTeX-SURFACE_TENSION","BibTeX-MELTING_LINE","BibTeX-VISCOSITY"]
            print(logf, " $bi:$(get_fluid_param_string(fluid, bi))");
        end
        print(logf, "\n");
     end
     close(logf);
end
# ---------------------------------
#       Information functions
# ---------------------------------

"""
    get_global_param_string(key::AbstractString)

Get a globally-defined string.

# Ref
ref CoolProp::get_global_param_string

# Arguments
* `key`: A string represents parameter name, could be one of $inputs_to_get_global_param_string
"""
function get_global_param_string(key::AbstractString)
    val = ccall( (:get_global_param_string, "CoolProp"), Clong, (Cstring, Ptr{UInt8}, Int), key, message_buffer::Array{UInt8, 1}, buffer_length)
    if val == 0
        error("CoolProp: ", get_global_param_string("errstring"))
    end
    return unsafe_string(convert(Ptr{UInt8}, pointer(message_buffer::Array{UInt8, 1})))
end

"""
    get_parameter_information_string(key::AbstractString, outtype::AbstractString)
    get_parameter_information_string(key::AbstractString)

Get information for a parameter.

# Arguments
* `key`: A string represents parameter name, to see full list check "Table of string inputs to PropsSI function": http://www.coolprop.org/coolprop/HighLevelAPI.html#parameter-table, or simply type `get_global_param_string("parameter_list")`
* `outtype="long"`: Output type, could be one of the `["IO", "short", "long", "units"]`, with a default value of "long"

# Example
```julia
julia> get_parameter_information_string("HMOLAR")
"Molar specific enthalpy"
julia> get_parameter_information_string("HMOLAR", "units")
"J/mol"
```

# Note
A tabular output for this function is available with `?CoolProp_parameters`
"""
function get_parameter_information_string(key::AbstractString, outtype::AbstractString)
    message_buffer[1:length(outtype)+1] = [Vector{UInt8}(outtype); 0x00]
    val = ccall( (:get_parameter_information_string, "CoolProp"), Clong, (Cstring, Ptr{UInt8}, Int), key, message_buffer::Array{UInt8, 1}, buffer_length)
    if val == 0
        error("CoolProp: ", get_global_param_string("errstring"))
    end
    return unsafe_string(convert(Ptr{UInt8}, pointer(message_buffer::Array{UInt8, 1})))
end

function get_parameter_information_string(key::AbstractString)
    return get_parameter_information_string(key, "long")
end

"""
    get_fluid_param_string(fluid::AbstractString, param::AbstractString)

Get a string for a value from a fluid (numerical values for the fluid can be obtained from Props1SI function).

# Arguments
* `fluid`: The name of the fluid that is part of CoolProp, for instance "n-Propane"
* `param`: A string, can be in one of the terms described in the following table

ParamName                    | Description
:----------------------------|:----------------------------------------
"aliases"                    | A comma separated list of aliases for the fluid
"CAS", "CAS_number"          | The CAS number
"ASHRAE34"                   | The ASHRAE standard 34 safety rating
"REFPROPName", "REFPROP_name"| The name of the fluid used in REFPROP
"Bibtex-XXX"                 | A BibTeX key, where XXX is one of the bibtex keys used in get_BibTeXKey
"pure"                       | "true" if the fluid is pure, "false" otherwise
"formula"                    | The chemical formula of the fluid in LaTeX form if available, "" otherwise

# Note
A tabular output for this function is available with `?CoolProp_fluids`
"""
function get_fluid_param_string(fluid::AbstractString, param::AbstractString)
    val = ccall( (:get_fluid_param_string, "CoolProp"), Clong, (Cstring, Cstring, Ptr{UInt8}, Int), fluid, param, message_buffer::Array{UInt8, 1}, buffer_length)
    if val == 0
        error("CoolProp: ", get_global_param_string("errstring"))
    end
    return unsafe_string(convert(Ptr{UInt8}, pointer(message_buffer::Array{UInt8, 1})))
end

"""
    F2K(tf::Real)

Convert from degrees Fahrenheit to Kelvin (useful primarily for testing).
"""
function F2K(tf::Real)
    return ccall( (:F2K, "CoolProp"), Cdouble, (Cdouble,), tf)
end

"""
    K2F(tk::Real)

Convert from Kelvin to degrees Fahrenheit (useful primarily for testing).
"""
function K2F(tk::Real)
    return ccall( (:K2F, "CoolProp"), Cdouble, (Cdouble,), tk)
end

"""
    get_param_index(param::AbstractString)

Get the index as a long for a parameter "T", "P", etc, for `abstractstate_keyed_output()` function.

# Arguments
* `param`: A string represents parameter name, to see full list check "Table of string inputs to PropsSI function": http://www.coolprop.org/coolprop/HighLevelAPI.html#parameter-table, or simply type `get_global_param_string("parameter_list")`
"""
function get_param_index(param::AbstractString)
    val = ccall( (:get_param_index, "CoolProp"), Clong, (Cstring,), param)
    if val == -1
        error("CoolProp: Unknown parameter: ", param)
    end
    return val
end

"""
    get_input_pair_index(pair::AbstractString)

Get the index for an input pair "PT_INPUTS", "HmolarQ_INPUTS", etc, for `abstractstate_update()` function.

# Arguments
* `pair::AbstractString`: The name of an input pair, described in the following table

Input Pair          |Description
:-------------------|:------------------------------------------------
QT_INPUTS           |Molar quality, Temperature in K
QSmolar_INPUTS      |Molar quality, Entropy in J/mol/K
QSmass_INPUTS       |Molar quality, Entropy in J/kg/K
HmolarQ_INPUTS      |Enthalpy in J/mol, Molar quality
HmassQ_INPUTS       |Enthalpy in J/kg, Molar quality
DmassQ_INPUTS       |Molar density kg/m^3, Molar quality
DmolarQ_INPUTS      |Molar density in mol/m^3, Molar quality
PQ_INPUTS           |Pressure in Pa, Molar quality
PT_INPUTS           |Pressure in Pa, Temperature in K
DmassT_INPUTS       |Mass density in kg/m^3, Temperature in K
DmolarT_INPUTS      |Molar density in mol/m^3, Temperature in K
HmassT_INPUTS       |Enthalpy in J/kg, Temperature in K
HmolarT_INPUTS      |Enthalpy in J/mol, Temperature in K
SmassT_INPUTS       |Entropy in J/kg/K, Temperature in K
SmolarT_INPUTS      |Entropy in J/mol/K, Temperature in K
TUmass_INPUTS       |Temperature in K, Internal energy in J/kg
TUmolar_INPUTS      |Temperature in K, Internal energy in J/mol
DmassP_INPUTS       |Mass density in kg/m^3, Pressure in Pa
DmolarP_INPUTS      |Molar density in mol/m^3, Pressure in Pa
HmassP_INPUTS       |Enthalpy in J/kg, Pressure in Pa
HmolarP_INPUTS      |Enthalpy in J/mol, Pressure in Pa
PSmass_INPUTS       |Pressure in Pa, Entropy in J/kg/K
PSmolar_INPUTS      |Pressure in Pa, Entropy in J/mol/K
PUmass_INPUTS       |Pressure in Pa, Internal energy in J/kg
PUmolar_INPUTS      |Pressure in Pa, Internal energy in J/mol
DmassHmass_INPUTS   |Mass density in kg/m^3, Enthalpy in J/kg
DmolarHmolar_INPUTS |Molar density in mol/m^3, Enthalpy in J/mol
DmassSmass_INPUTS   |Mass density in kg/m^3, Entropy in J/kg/K
DmolarSmolar_INPUTS |Molar density in mol/m^3, Entropy in J/mol/K
DmassUmass_INPUTS   |Mass density in kg/m^3, Internal energy in J/kg
DmolarUmolar_INPUTS |Molar density in mol/m^3, Internal energy in J/mol
HmassSmass_INPUTS   |Enthalpy in J/kg, Entropy in J/kg/K
HmolarSmolar_INPUTS |Enthalpy in J/mol, Entropy in J/mol/K
SmassUmass_INPUTS   |Entropy in J/kg/K, Internal energy in J/kg
SmolarUmolar_INPUTS |Entropy in J/mol/K, Internal energy in J/mol

# Example
```julia
julia> get_input_pair_index("PT_INPUTS")
9
```
"""
function get_input_pair_index(pair::AbstractString)
    val = ccall( (:get_input_pair_index, "CoolProp"), Clong, (Cstring,), pair)
    if val == -1
        error("CoolProp: Unknown input pair: ", pair)
    end
    return val
end
const coolpropparameters = map(Compat.String, split(get_global_param_string("parameter_list"),','));
const coolpropfluids = map(Compat.String, split(get_global_param_string("FluidsList"),','));


text:
# starting to add exports here

global WORKERPOOL = WorkerPool()

"""
    $SIGNATURES

For use with `multiproc`, nominal use is a worker pool of all processes available above and including 2..., but will return single process [1;] if only the first processes is available.
"""
function setWorkerPool!(pool::Vector{Int}=1 < nprocs() ? setdiff(procs(), [1;]) : [1;])
  global WORKERPOOL
  WORKERPOOL = WorkerPool(pool)
end

function getWorkerPool()
  global WORKERPOOL
  return WORKERPOOL
end

## =============================================================================================
# Iterate over variables in clique
## =============================================================================================

"""
    $SIGNATURES

Dev Notes
- part of refactoring fmcmc.
- function seems excessive
"""
function compileFMCMessages(fgl::AbstractDFG,
                            lbls::Vector{Symbol},
                            solveKey::Symbol,
                            logger=ConsoleLogger())
  #
  d = Dict{Symbol,TreeBelief}()
  for vsym in lbls
    vari = DFG.getVariable(fgl,vsym)
    d[vsym] = TreeBelief(vari, solveKey)
  end
  return d
end

function doFMCIteration(fgl::AbstractDFG,
                        vsym::Symbol,
                        solveKey::Symbol,
                        cliq::TreeClique,
                        fmsgs,
                        N::Int,
                        dbg::Bool,
                        needFreshMeasurements::Bool=true,
                        logger=ConsoleLogger()  )
  #

  vert = DFG.getVariable(fgl, vsym)
  if !getSolverData(vert, solveKey).ismargin
    # potprod = nothing
    densPts, ipc = predictbelief(fgl, vsym, :, needFreshMeasurements=needFreshMeasurements, N=N, dbg=dbg, logger=logger)

    if 0 < length(densPts)
      # TODO --  can we remove this duplicate getVert?
      # updvert = DFG.getVariable(fgl, vsym)
      setValKDE!(vert, densPts, true, ipc)
      # TODO perhaps more debugging inside `predictbelief`?
    end
  end
  nothing
end

"""
    $(SIGNATURES)

Iterate successive approximations of clique marginal beliefs by means
of the stipulated proposal convolutions and products of the functional objects
for tree clique `cliq`.
"""
function fmcmc!(fgl::AbstractDFG,
                cliq::TreeClique,
                fmsgs::Vector{<:LikelihoodMessage},
                lbls::Vector{Symbol},
                solveKey::Symbol,
                N::Int,
                MCMCIter::Int,
                dbg::Bool=false,
                logger=ConsoleLogger(),
                multithreaded::Bool=false  )
  #
  with_logger(logger) do
    @info "---------- successive fnc approx ------------$(getLabel(cliq))"
  end
  # repeat several iterations of functional Gibbs sampling for fixed point convergence
  if length(lbls) == 1
      MCMCIter=1
  end
  # mcmcdbg = Array{CliqGibbsMC,1}()

  # burn-in loop for outer Gibbs
  for iter in 1:MCMCIter
    # iterate through each of the variables, KL-divergence tolerence would be nice test here
    with_logger(logger) do
      @info "#$(iter)\t -- "
    end
    # dbgvals = !dbg ? nothing : CliqGibbsMC([], Symbol[])

    needFreshMeasurements = iter == 1 || getSolverParams(fgl).alwaysFreshMeasurements

    # outer Gibbs cycle
    for vsym in lbls
        doFMCIteration(fgl, vsym, solveKey, cliq, fmsgs, N, dbg, needFreshMeasurements, logger)
    end
    # !dbg ? nothing : push!(mcmcdbg, dbgvals)
  end

  # populate dictionary for return NBPMessage in multiple dispatch
  msgdict = compileFMCMessages(fgl, lbls, solveKey, logger)

  return msgdict
end

## =============================================================================================
# Up solve
## =============================================================================================

"""
    $(SIGNATURES)

Perform computations required for the upward message passing during belief propation on the Bayes (Junction) tree.
This function is usually called as via remote_call for multiprocess dispatch.

Notes
- `fg` factor graph,
- `tree` Bayes tree,
- `cliq` which cliq to perform the computation on,
- `parent` the parent clique to where the upward message will be sent,
- `childmsgs` is for any incoming messages from child cliques.

DevNotes
- FIXME total rewrite with AMP #41 and RoME #244 in mind
"""
function upGibbsCliqueDensity(dfg::AbstractDFG,
                              cliq::TreeClique,
                              solveKey::Symbol,
                              inmsgs,
                              N::Int=getSolverParams(dfg).N,
                              dbg::Bool=false,
                              iters::Int=3,
                              logger=ConsoleLogger()  ) # where {T, T2}
  #
  with_logger(logger) do
    @info "up w $(length(inmsgs)) msgs"
  end
  # TODO -- some weirdness with: d,. = d = ., nothing
  # mcmcdbg = Array{CliqGibbsMC,1}()
  d = Dict{Symbol,TreeBelief}()

  # priorprods = Vector{CliqGibbsMC}()

  cliqdata = getCliqueData(cliq)

  # use nested structure for more efficient Chapman-Kolmogorov solution approximation
  if false
    IDS = [cliqdata.frontalIDs;cliqdata.separatorIDs] #inp.cliq.attributes["frontalIDs"]
    d = fmcmc!(dfg, cliq, inmsgs, IDS, solveKey, N, iters, dbg, logger)
  else
    # NOTE -- previous mistake, must iterate over directsvarIDs also (or incorporate once at the right time)
    # NOTE -- double up on directs to allow inflation to take proper affect, see #1051
    d = fmcmc!(dfg, cliq, inmsgs, cliqdata.directFrtlMsgIDs, solveKey, N, 1, dbg, logger, true)

    if length(cliqdata.msgskipIDs) > 0
      dd = fmcmc!(dfg, cliq, inmsgs, cliqdata.msgskipIDs, solveKey, N, 1, dbg, logger, true)
      for md in dd d[md[1]] = md[2]; end
    end
    if length(cliqdata.itervarIDs) > 0
      ddd = fmcmc!(dfg, cliq, inmsgs, cliqdata.itervarIDs, solveKey, N, iters, dbg, logger, false)
      for md in ddd d[md[1]] = md[2]; end
    end
    if length(cliqdata.directPriorMsgIDs) > 0
      doids = setdiff(cliqdata.directPriorMsgIDs, cliqdata.msgskipIDs)
      dddd = fmcmc!(dfg, cliq, inmsgs, doids, solveKey, N, 1, dbg, logger, true)
      for md in dddd d[md[1]] = md[2]; end
    end
  end

  return d
end

## =============================================================================================
# Down solve
## =============================================================================================

## ============================================================================
# Initialization is slightly different and likely to be consolidated
## ============================================================================

"""
    $SIGNATURES

Cycle through var order and initialize variables as possible in `subfg::AbstractDFG`.
Return true if something was updated.

Notes:
- assumed `subfg` is a subgraph containing only the factors that can be used.
  - including the required up or down messages
- intended for both up and down initialization operations.

Dev Notes
- Should monitor updates based on the number of inferred & solvable dimensions
"""
function cycleInitByVarOrder!(subfg::AbstractDFG,
                              varorder::Vector{Symbol};
                              solveKey::Symbol=:default,
                              logger=ConsoleLogger()  )
  #
  with_logger(logger) do
    @info "cycleInitByVarOrder! -- varorder=$(varorder)"
  end
  retval = false
  count = 1
  while count > 0
    count = 0
    for vsym in varorder
      var = DFG.getVariable(subfg, vsym)
      isinit = isInitialized(var, solveKey)
      with_logger(logger) do
        @info "var.label=$(var.label) is initialized=$(isinit)"
      end
      doautoinit!(subfg, [var;], solveKey=solveKey, logger=logger)
      if isinit != isInitialized(var, solveKey)
        count += 1
        retval = true
      end
    end
  end
  with_logger(logger) do
    @info "cycleInitByVarOrder!, retval=$(retval)"
  end
  flush(logger.stream)
  return retval
end

#


text:
using LinearAlgebra, Plots, JLD2
using Flux, Flux.Zygote, Optim, DiffEqFlux
using KitBase, ProgressMeter
import KitML

maxwell_boltzmann(x) = x * log(x) - x

function kinetic_entropy(α, m, weights)
    B = KitBase.maxwell_boltzmann_dual_prime.(α' * m)[:]
    return sum(maxwell_boltzmann.(B) .* weights)
end

cd(@__DIR__)
include("math.jl")
model = KitML.load_model("best_model.h5"; mode = :tf)

begin
    # space
    x0 = 0.0
    x1 = 1.0
    nx = 100
    dx = (x1 - x0) / nx

    pspace = KitBase.PSpace1D(x0, x1, nx)

    # time
    tEnd = 1.0
    cfl = 0.5
    dt = cfl * dx

    # quadrature
    quadratureorder = 5
    points, triangulation = KitBase.octa_quadrature(quadratureorder)
    weights = KitBase.quadrature_weights(points, triangulation)
    nq = size(points, 1)

    # moments
    L = 1
    ne = (L + 1)^2
    #m = KitBase.eval_spherharmonic(points, L)
    m = ComputeSphericalBasisAnalytical(L, 3, points)

    # collision
    σs = ones(Float32, nx)
    σa = zeros(Float32, nx)# * 0.1
    #σa = exp.(pspace.x[1:nx])

    σt = σs + σa
    σq = zeros(Float32, nx)
    #σq = 1e-10 .* exp.(-100 * pspace.x[1:nx])

    global t = 0.0
    flux = zeros(Float32, ne, nx + 1)
end

pdf = zeros(Float32, axes(weights))
begin
    α = zeros(Float32, ne, nx)
    phi = ones(Float32, ne, nx) .* 1e-6
    phi[1, 1] = 0.5
    #=
    for i = 1:nx
        #=for j = 1:nq
            if points[j, 1] > 0.
                pdf[j] = (1.0 + sin(2π * pspace.x[i])) ./ 4π
            end
        end
        phi[:, i] .= m * pdf=#

        #phi[1, i] = 1.0 #- 0.05 * (pspace.x[i] - x0)
        #phi[1, i] = 1.0 + 0.2 * sin(2π * pspace.x[i])
        #phi[2, i] = 0.2
    end=#
end

#αT = zeros(Float32, nx, ne)
#phiT = zeros(Float32, nx, ne)

#global X = zeros(Float32, 1, ne)
#global Y = zeros(Float32, 1, 1)

#res = KitBase.optimize_closure(X[1, :], m, weights, phi[:, nx÷2], KitBase.maxwell_boltzmann_dual)
#X[1, :] .= res.minimizer
#Y[1, 1] = kinetic_entropy(X[1, :], m, weights)

@showprogress for iter = 1:20#20
    #phi_old .= phi
#=
    # regularization
    @inbounds for i = 1:ne
        phiT[:, i] .= phi[i, :, :][:]
    end
    αT .= KitML.neural_closure(model, phiT)
    @inbounds for i = 1:ne
        α[i, :, :] .= reshape(αT[:, i], nx, ny)
    end
    @inbounds Threads.@threads for j = 1:ny
        for i = 1:nx
            phi_temp[:, i, j] .= KitBase.realizable_reconstruct(α[:, i, j], m, weights, KitBase.maxwell_boltzmann_dual_prime)
        end
    end

    @inbounds for j = 1:ny
        for i = 1:nx
            #@show norm(phi_temp[:, i, j] - phi_old[:, i, j], 2)

            if norm(phi_temp[:, i, j] - phi_old[:, i, j], 2) > 1e-1
                res = KitBase.optimize_closure(α[:, i, j], m, weights, phi[:, i, j], KitBase.maxwell_boltzmann_dual)
                α[:, i, j] .= res.minimizer
                phi[:, i, j] .= KitBase.realizable_reconstruct(res.minimizer, m, weights, KitBase.maxwell_boltzmann_dual_prime)

                if phi[1, i, j] > 0.01
                    X = vcat(X, permutedims(α[:, i, j]))
                    Y = vcat(Y, kinetic_entropy(α[:, i, j], m, weights))
                end
            else
                phi[:, i, j] .= phi_temp[:, i, j]
            end
        end
    end
=#

    res = KitBase.optimize_closure(α[:, 1], m, weights, phi[:, 1], KitBase.maxwell_boltzmann_dual)
    α[:, 1] .= res.minimizer
    #res = KitBase.optimize_closure(α[:, nx], m, weights, phi[:, nx], KitBase.maxwell_boltzmann_dual)
    #α[:, nx] .= res.minimizer

    @inbounds for i = 2:nx
        res = KitBase.optimize_closure(α[:, i], m, weights, phi[:, i], KitBase.maxwell_boltzmann_dual; optimizer=Newton())
        α[:, i] .= res.minimizer
        phi[:, i] .= KitBase.realizable_reconstruct(res.minimizer, m, weights, KitBase.maxwell_boltzmann_dual_prime)
    end

    # flux
    fη = zeros(nq)
    @inbounds for i = 2:nx
        KitBase.flux_kfvs!(fη, KitBase.maxwell_boltzmann_dual.(α[:, i-1]' * m)[:], KitBase.maxwell_boltzmann_dual.(α[:, i]' * m)[:], points[:, 1], dt)

        for k in axes(flux, 1)
            flux[k, i] = sum(m[k, :] .* weights .* fη)
        end
    end
    #=
    # periodic
    KitBase.flux_kfvs!(fη, KitBase.maxwell_boltzmann_dual.(α[:, nx]' * m)[:], KitBase.maxwell_boltzmann_dual.(α[:, 1]' * m)[:], points[:, 1], dt)
    for k in axes(flux, 1)
        flux[k, 1] = sum(m[k, :] .* weights .* fη)
        flux[k, nx+1] = flux[k, 1]
    end
    =#

    # update
    @inbounds for i = 2:nx-1
        for q = 1:1
            phi[q, i] =
                phi[q, i] +
                (flux[q, i] - flux[q, i+1]) / dx +
                (σs[i] * phi[q, i] - σt[i] * phi[q, i]) * dt +
                σq[i] * dt
        end

        for q = 2:ne
            phi[q, i] =
                phi[q, i] +
                (flux[q, i] - flux[q, i+1]) / dx +
                (-σt[i] * phi[q, i]) * dt
        end
    end

    phi[:, nx] .=  phi[:, nx-1]

    global t += dt

    #=
    if iter%9 == 0
        model.fit(X, Y, epochs=2)

        X = zeros(Float32, 1, ne)
        Y = zeros(Float32, 1, 1)
        res = KitBase.optimize_closure(X[1, :], m, weights, phi[:, nx÷2, ny÷2], KitBase.maxwell_boltzmann_dual)
        X[1, :] .= res.minimizer
        Y[1, 1] = kinetic_entropy(X[1, :], m, weights)
    end=#
end

scatter(pspace.x[1:nx], phi[1, :])
scatter(pspace.x[1:nx], α[1, :])


text:
mutable struct ZoomTransform{T <: AbstractFloat}
    # An optional (occasionally used) translation.
    position::Vector2D{T}

    # The zoom factor generally incremented in small steps.
    # For example, 0.1
    scale::Vector2D{T}

    # The focal point where zooming occurs
    zoom_at::Vector2D{T}

    # A "running" accumulating transform
    acc_transform::AffineTransform

    # A transform that includes position translation.
    transform::AffineTransform

    function ZoomTransform{T}() where {T <: AbstractFloat}
        o = new()

        o.position = Vector2D{Float64}()
        o.scale = Vector2D{Float64}(1.0, 1.0)
        o.zoom_at = Vector2D{Float64}()

        o.acc_transform = AffineTransform{Float64}()

        o.transform = AffineTransform{Float64}()
        o
    end
end

function get_transform(zoom::ZoomTransform)
    update!(zoom)
    zoom.transform
end

function update!(zoom::ZoomTransform)
    # Accumulate zoom transformations.
    # acc_transform is an intermediate accumulative matrix used for tracking the current zoom target.
    translate!(zoom.acc_transform, zoom.zoom_at.x, zoom.zoom_at.y)
    scale!(zoom.acc_transform, zoom.scale.x, zoom.scale.y)
    translate!(zoom.acc_transform, -zoom.zoom_at.x, -zoom.zoom_at.y)

    # We reset Scale because acc_transform is accumulative and has "captured" the information.
    set!(zoom.scale, 1.0, 1.0)

    # We want to leave acc_transform solely responsible for zooming.
    # "transform" is the final matrix.
    set!(zoom.transform, zoom.acc_transform)

    # Tack on translation. Note: we don't append it, but concat it into a separate matrix.
    translate!(zoom.transform, zoom.position.x, zoom.position.y)
end

# Use this if you want to manually set the positional value.
# You would typically use translateBy() instead.
function set_position!(zoom::ZoomTransform, x::T, y::T) where {T <: AbstractFloat}
    set!(zoom.position, x, y)
end

# A relative zoom.
# [delta] is relative to the current scale/zoom.
function zoom_by!(zoom::ZoomTransform, delta_x::T, delta_y::T) where {T <: AbstractFloat}
    set!(zoom.scale, zoom.scale.x + delta_x, zoom.scale.y + delta_y)
end

function translate_by!(zoom::ZoomTransform, delta_x::T, delta_y::T) where {T <: AbstractFloat}
    set!(zoom.position, zoom.position.x + delta_x, zoom.position.y + delta_y)
end

function set_scale!(zoom::ZoomTransform, scale::T) where {T <: AbstractFloat}
    update!(zoom)

    # We use dimensional analysis to set the scale. Remember we can't
    # just set the scale absolutely because acc_transform is an accumulating matrix.
    # We have to take its current value and compute a new value based
    # on the passed in value.

    # Also, I can use acc_transform.a because I don't allow rotations for zooms,
    # so the diagonal components correctly represent the matrix's current scale.
    # And because I only perform uniform scaling I can safely use just the "a" element.
    scale_factor = scale / zoom.acc_transform.a

    set!(zoom.scaleX, scale_factor, scale_factor)
end

function set_zoom_at!(zoom::ZoomTransform, x::T, y::T) where {T <: AbstractFloat}
    set!(zoom.zoom_at, x, y)
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
@marginalrule typeof(dot)(:in1_in2) (m_out::NormalDistributionsFamily, m_in1::PointMass, m_in2::NormalDistributionsFamily, meta::AbstractCorrection) = begin

    # Forward message towards `in2` edge
    mf_in2 = @call_rule typeof(dot)(:in2, Marginalisation) (m_out = m_out, m_in1 = m_in1, meta = meta)
    q_in2  = prod(ProdAnalytical(), m_in2, mf_in2)

    return (in1 = m_in1, in2 = q_in2)
end

@marginalrule typeof(dot)(:in1_in2) (m_out::NormalDistributionsFamily, m_in1::NormalDistributionsFamily, m_in2::PointMass, meta::AbstractCorrection) = begin
    symmetric = @call_marginalrule typeof(dot)(:in1_in2) (m_out = m_out, m_in1 = m_in2, m_in2 = m_in1, meta = meta)
    return (in1 = symmetric[:in2], in2 = symmetric[:in1])
end


text:
using Optim
using StatsBase
include("../dataset.jl")
include("../likelihoods.jl")

percent = ARGS[1]
split_prop = parse(Float64, percent)/100
# fname must be a SNAP dataset, preprocessed according to README.md
fname = ARGS[2]

# Compute sufficient statistics used by the likelihood functions
degs, ts, _, _, _, _ = trainTestSplitSnapData("$fname", split_prop)
Tend = sum(degs)
dmap = countmap(degs)
ds = collect(keys(dmap))
dcounts = collect(values(dmap))
K = sum(dcounts)

# Fit a PYP to this dataset
result = optimize(params -> -pyp_llikelihood(params, ds, dcounts, ts, K, Tend),
                  (storage, params) -> neg_grad_pyp_llikelihood!(storage, params, ds, dcounts, ts, K, Tend),
                  [0., .5], LBFGS())
tau = exp(result.minimizer[1])/(1+exp(result.minimizer[1]))
theta = result.minimizer[2]
println("Tau ", tau, " Theta ", theta, )
println("Optimized ll ", -result.minimum)
println("Expected powerlaw ", 1+tau)

# Fit a Geom(beta)Beta Neutral-to-the-left model to this dataset
result = optimize(a -> -ntl_llikelihood(a, ds, dcounts, ts, K, Tend),
                  (storage, params) -> neg_grad_ntl_llikelihood!(storage, params, ds, dcounts, ts, K, Tend),
                  [0.], LBFGS())
alpha = 1 - exp(result.minimizer[1])
println("NTL alpha ", alpha)
ll_ntl = result.minimum
result = optimize(g -> -geom_llikelihood(g, K, Tend), [.5], LBFGS())
g = exp(result.minimizer[1])/(1+exp(result.minimizer[1]))
println("g ", g)
ll = ll_ntl + result.minimum
println("Optimized ll ", -ll)
println("Expected powerlaw ", 1+(1/g - alpha)/(1/g - 1))


text:
"""
Holidays for TARGET Eurozone (Trans-European Automated Real-time Gross Settlement Express Transfer System)
"""
struct TARGET <: HolidayCalendar end
const TARGET2 = TARGET
const EuroZone = TARGET

function isholiday(::TARGET, dt::Dates.Date)
    yy = Dates.year(dt)
    mm = Dates.month(dt)
    dd = Dates.day(dt)

    # Bisection
    if mm < 8
        # Fixed Holidays
        if (
            # New Year's Day
            ((mm == 1) && (dd == 1))
            ||
            # Labour Day
            ((mm == 5) && (dd == 1) && (yy >= 2000))
            )
            return true
        end

        if yy >= 2000

            dt_rata::Int = Dates.days(dt)
            e_rata::Int = easter_rata(Dates.Year(yy))

            if (
                # Good Friday
                ( dt_rata == ( e_rata -  2 ))
                ||
                # Easter Monday
                ( dt_rata == ( e_rata + 1))
                )
                return true
            end
        end
    else
        # mm >= 8
        if (
            # Christmas
            ((mm == 12) && (dd == 25))
            ||
            # Day of Goodwill
            ((mm == 12) && (dd == 26) && (yy >= 2000))
            ||
            # End of year
            ((mm == 12) && (dd == 31) && (yy == 1998 || yy == 1999 || yy == 2001))
            )
            return true
        end
    end

    return false
end


text: struct a end


text:
include("util.jl")

using Distributions

function clean_comment(s::String)
    replace(s, r"(#\s+)"=>"")
end

function parse_node_label(s::String)
    return parse(Int64, s) + 1
end

function normalize_edgelist!(mat::Adj_Matrix)
    rows, cols = size(mat)
    for i in 1:rows
        total = 0
        @simd for j in 1:cols
            if mat[i, j] != 0.0
                total += mat[i, j]
            end
        end

        @simd for j in 1:cols
            if mat[i, j] != 0.0
                mat[i, j] /= total
            end
        end
    end
end

function scale_mat!(f::PSYCHOLOGICAL_FACTOR, mat::Adj_Matrix)
    rows, cols = size(mat)
    for i = 1:rows
        @simd for j = 1:cols
            mat[i, j] *= f[i]
        end
    end
end

function scale_mat!(f::PSYCHOLOGICAL_FACTOR, mat::Array{WEIGHT, 2})
    rows, cols = size(mat)
    for i = 1:rows
        @simd for j = 1:cols
            mat[i, j] *= f[i]
        end
    end
end

function read_edgelist(s::String,
        delim::Union{String, Regex})::Tuple{Adj_Matrix, Adj_Matrix, Int64}
    I_pos = Array{Int64, 1}()
    J_pos = Array{Int64, 1}()
    K_pos = Array{WEIGHT, 1}()

    I_neg = Array{Int64, 1}()
    J_neg = Array{Int64, 1}()
    K_neg = Array{WEIGHT, 1}()

    n = 0
    open(s) do file
        for (index, line) in enumerate(readlines(file))
            if startswith(line, "#") && (index == 1)
                comment = clean_comment(line)
                n = parse(Int64, comment)
            elseif !startswith(line, "#")
                contents = map(x -> String(x), split(line, delim))
                #(u, v, w) = (0, 0, WEIGHT(1.0))
                u = parse_node_label(contents[1])
                v = parse_node_label(contents[2])
                w = WEIGHT(1.0)
                if length(contents) > 2
                    w = parse(WEIGHT, contents[3])
                    if w > 0.0
                        push!(I_pos, u)
                        push!(J_pos, v)
                        push!(I_pos, v)
                        push!(J_pos, u)
                        push!(K_pos, w)
                        push!(K_pos, w)
                    elseif w < 0.0
                        push!(I_neg, u)
                        push!(J_neg, v)
                        push!(I_neg, v)
                        push!(J_neg, u)
                        push!(K_neg, w)
                        push!(K_neg, w)
                    end
                end
            end
        end
    end
    n_temp = max(maximum(I_pos), maximum(J_pos))
    n = max(n, n_temp)
    n_temp = max(maximum(I_neg), maximum(J_neg))
    n = max(n, n_temp)

    pos = sparse(I_pos, J_pos, K_pos, n, n)
    neg = sparse(I_neg, J_neg, K_neg, n, n)

    return (pos, neg, n)
end

function schatten_norm(A, p)
    inner = tr(A' * A)
    inner = inner ^ (p / 2)
    outer = inner ^ 1/p
end

function construct_graph(A)
    (n, m) = size(A)
    g = SimpleGraph(n)
    for i = 1:n, j = 1:n
        if A[i, j] != 0
            add_edge!(g, i, j)
        end
    end
    return g
end

max_stars_repo_name: InzamamRahaman/SignedDegroot

max_issues_repo_name: InzamamRahaman/SignedDegroot

max_forks_repo_name: InzamamRahaman/SignedDegroot


text:
function encode_decode (x)
    e, n = encode(x)
    d, m = decode(uint8(e))
    @assert n == m
    d, m
end

x, nbytes = encode_decode(uint8(255))
@test x == uint8(255)

x, nbytes = encode_decode(uint32(256))
@test x == uint32(256)

x, nbytes = encode_decode(float(1.0))
@test x == float(1.0)

x, nbytes = encode_decode(big(256256256256))
@test x == big(256256256256)

x, nbytes = encode_decode(:symbol)
@test x == :symbol

x, nbytes = encode_decode(tuple())
@test x == tuple()

x, nbytes = encode_decode((:a, uint8(255), float(1.0)))
@test x == (:a, uint8(255), float(1.0))

x, nbytes = encode_decode("a_string")
@test x == "a_string"

#x, nbytes = encode_decode([:a, :b, :c])
#@test x == [:a, :b, :c]

x, nbytes = encode_decode(uint8([1, 2, 3, 4]))
@test x == uint8([1, 2, 3, 4])

max_stars_repo_name: esl/FsBert

max_issues_repo_name: esl/FsBert

max_forks_repo_name: esl/FsBert


text:
module TestParallelReduce
include("preamble.jl")

struct Recorder
    records
    lock
end

Recorder() = Recorder([], Threads.SpinLock())

function (recorder::Recorder)(x)
    try
        lock(recorder.lock)
        push!(recorder.records, x)
        return x
    finally
        unlock(recorder.lock)
    end
end

@testset "early termination (runtime)" begin
    recorder = Recorder()
    records = recorder.records
    @test reduce(right, Map(recorder) |> ReduceIf(==(1)), 1:32;
                 basesize=1) === 1
    if VERSION >= v"1.3-alpha"
        if Threads.nthreads() == 1
            @test length(records) == 1
        else
            # Not sure if this safe, but let's check if _some_
            # tasks are canceled.
            @test length(records) < sum(2 .^ (0:5))
            # (though it seems `length(records)` is always 1?)
        end
    end
end

@testset "early termination (grid)" begin
    @testset for needle in 1:20, len in 1:20
        @test reduce(right, ReduceIf(x -> x >= needle), 1:len; basesize=1) ==
            min(needle, len)
    end
end

@testset "early termination (shuffle)" begin
    @testset for seed in rand(UInt32, 10), p in Any[50, 95, 99, Inf]
        rng = MersenneTwister(seed)
        basesize = 16
        xs = collect(enumerate(rand(rng, 1:100, 100 * basesize)))
        xf = ReduceIf(x -> x[2] >= p)
        @test reduce(right, xf, xs; basesize=basesize) == foldl(right, xf, xs)
    end
end

@testset "withprogress" begin
    xf = Map() do x
        x
    end
    @test reduce(+, xf, withprogress(1:100; interval=0); basesize=1) == 5050

    xf2 = ScanEmit(0) do u, x
        y = u + x
        y, y
    end
    @test reduce(right, xf2, withprogress(1:100; interval=0); basesize=1) == 5050
end

end  # module


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -17.39516631142554 <= q <= 12.050938789121952)

                     add_NL_constraint(m, :(softplus(0.1892165831624366 + -0.6797303891528506*softplus(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])) + -0.5563706493225857*softplus(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])) + -0.5881166064998418*softplus(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])) + -0.5324749908347566*softplus(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))) + softplus(-0.7087227741366906 + -0.3339546766337098*softplus(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])) + -0.8987724666577948*softplus(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])) + 0.3907505378456695*softplus(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])) + -0.7983638213649877*softplus(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))) + softplus(0.16171979073463483 + 0.31171990946083294*softplus(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])) + 0.2757921127232579*softplus(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])) + -0.732639957443523*softplus(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])) + -0.947011017876227*softplus(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))) + softplus(-0.17479526548317725 + -0.7121643152938684*softplus(0.23382152176968107 + 0.970568112692761*$(x[1]) + 0.9147961292868287*$(x[2]) + 0.25099364821370385*$(x[3])) + -0.20573796119512977*softplus(0.1824369939040409 + -0.636682724457545*$(x[1]) + 0.684450108697257*$(x[2]) + -0.3896097552087272*$(x[3])) + 0.920351043144791*softplus(-0.7017137443246066 + -0.579761469455264*$(x[1]) + -0.38090529589576505*$(x[2]) + -0.30474275998045197*$(x[3])) + -0.8770524767395065*softplus(0.495407339331571 + -0.11045663155169017*$(x[1]) + -0.19556503402268088*$(x[2]) + 0.8751710578699647*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using DataFrames, CSV
using Distributions
using RDatasets

df = DataFrame(
    name = ["todd", "jane", "yiting", "emilio", "betty"],
    attribute = randn(5),
    awesomeness = rand(1:5, 5),
    wealth = rand(0:0.1:5_000_000, 5),
    owns_hat = rand(Bool, 5),
)

CSV.write("session-1/example.csv", df)


text:
@testset "regression" begin
    Random.seed!(14)
    beta = [10; 20; 30]
    n = 1000
    X = rand(n, 3)
    y = X*beta + randn(n)

    model = regression(y, X)

    ss = statespace(model)

    beta_OLS = X\y

    # The result should be the same of OLS
    @test ss.filter.att[end, :] ≈ beta_OLS rtol = 1e-4
    @test ss.model.H[1] ≈ 1.0 rtol = 1e-2
end


text:
using DataTypesBasic
using Test

@test Try(4) == Identity(4)
@test Try(ErrorException("hi")) == Const(ErrorException("hi"))

@test istry(Identity(4))
@test istry(Const(ErrorException("hi")))
@test istry(@Try error("an error"))
@test istry(@Try :value)
@test istry(@TryCatch ErrorException error("an error"))
@test istry(@TryCatch ErrorException :value)
@test !istry(Some(4))

@test issuccess(@Try 5)
@test !issuccess(@Try error("hi"))
@test_throws MethodError issuccess(Const("other"))

@test isfailure(@Try error("hi"))
@test !isfailure(@Try 5)
@test_throws MethodError isfailure(Const("other"))

@test Thrown(ErrorException("one"), []) == Thrown(ErrorException("one"), [])
@test repr(Thrown(ErrorException("some"), [])) == """Thrown(ErrorException("some"))"""

@test (@Try error("hi")) isa Const{<:Thrown}
@test (@TryCatch ErrorException error("hi")) isa Const{<:Thrown}
@test_throws ArgumentError (@TryCatch ErrorException throw(ArgumentError("bad")))

me = MultipleExceptions(ErrorException("hi"))
@test me.exceptions == (ErrorException("hi"),)
@test MultipleExceptions(me, ErrorException("ho")).exceptions == (ErrorException("hi"), ErrorException("ho"))

e1 = Thrown(ErrorException("one"), [])
e2 = AssertionError("error")
@test merge(e1, e2) == MultipleExceptions(e1, e2)
@test merge(e2, e1) == MultipleExceptions(e2, e1)
@test merge(e1, e1) == MultipleExceptions(e1, e1)
m12 = MultipleExceptions(e1, e2)
@test merge(e2, m12) == MultipleExceptions(e2, e1, e2)
@test merge(m12, e2) == MultipleExceptions(e1, e2, e2)
@test merge(m12, m12) == MultipleExceptions(e1, e2, e1, e2)

@test eltype(Try{Int}) == Int
@test eltype(Identity{Int}) == Int
@test eltype(Thrown) == Any
@test eltype(Const{Thrown}) == Any
@test eltype(Try) == Any