text:
# Tests for Russell DEA Models
@testset "RussellDEAModel" begin

    # ------------------
    # Input oriented
    # ------------------
    X = [2 2; 1 4; 4 1; 4 3; 5 5; 6 1; 2 5; 1.6 8]
    Y = [1; 1; 1; 1; 1; 1; 1; 1]

    # Input oriented CRS
    dearussellio = dearussell(X, Y, orient = :Input, rts = :CRS)

    @test typeof(dearussellio) == RussellDEAModel

    @test nobs(dearussellio) == 8
    @test ninputs(dearussellio) == 2
    @test noutputs(dearussellio) == 1
    @test efficiency(dearussellio, :X) ≈ [1 1; 1 1; 1 1; 0.5 2/3; 0.4 0.4; 2/3 1; 0.5 0.8; 0.625 0.5]
    @test_throws ArgumentError efficiency(dearussellio, :Y) # No output efficiency in input oriented model
    @test efficiency(dearussellio) ≈ [1; 1; 1; 0.5833333333333334; 0.4; 0.8333333333333334; 0.65; 0.5625]
    @test convert(Matrix, peers(dearussellio)) ==
            [ 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
            1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
            1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0]
    @test slacks(dearussellio, :X) ≈ zeros(8,2)
    @test slacks(dearussellio, :Y) ≈ zeros(8,1)

    @test efficiency(dearussell(targets(dearussellio, :X), targets(dearussellio, :Y), orient = :Input, rts = :CRS)) ≈ ones(8,1)

    # Test default is :Input with :CRS
    efficiency(dearussell(X, Y)) == efficiency(dearussellio)

    # Test no slacks
    dearussellionoslack = dearussell(X, Y, slack = false)
    @test efficiency(dearussellionoslack) == efficiency(dearussellio)
    @test isempty(slacks(dearussellionoslack, :X)) == 1
    @test isempty(slacks(dearussellionoslack, :Y)) == 1

    # Input oriented VRS
    dearusselliovrs  = dearussell(X, Y, orient = :Input, rts = :VRS)
    @test efficiency(dearusselliovrs) ≈ [1; 1; 1; 0.5833333333333334; 0.4; 0.8333333333333334; 0.65; 0.5625]

    ## Test if one-by-one DEA using evaluation and reference sets match initial results
    dearussell_io_ref_eff = zeros(size(X, 1))
    dearussell_iovrs_ref_eff = zeros(size(X, 1))

    Xref = X[:,:]
    Yref = Y[:,:]

    for i = 1:size(X, 1)
        Xeval = X[i:i,:]
        Xeval = Xeval[:,:]
        Yeval = Y[i:i,:]
        Yeval = Yeval[:,:]

        dearussell_io_ref_eff[i] = efficiency(dearussell(Xeval, Yeval, orient = :Input, rts = :CRS, Xref = Xref, Yref = Yref))[1]
        dearussell_iovrs_ref_eff[i] = efficiency(dearussell(Xeval, Yeval, orient = :Input, rts = :VRS, Xref = Xref, Yref = Yref))[1]
    end

    @test dearussell_io_ref_eff == efficiency(dearussellio)
    @test dearussell_iovrs_ref_eff == efficiency(dearusselliovrs)

    # Add new DMU with slack
    dearussellioslack = dearussell([X; 2 2], [Y; 2], orient = :Input, rts = :VRS)
    @test slacks(dearussellioslack, :Y) ≈ [1; 0; 0; 1; 1; 0; 0; 0; 0]

    # Print
    show(IOBuffer(), dearussellio)
    show(IOBuffer(), dearussellionoslack)

    # ------------------
    # Output oriented
    # ------------------
    X = [1; 1; 1; 1; 1; 1; 1; 1]
    Y = [7 7; 4 8; 8 4; 3 5; 3 3; 8 2; 6 4; 1.5 5]

    # Output oriented CRS
    dearusselloo = dearussell(X, Y, orient = :Output, rts = :CRS)

    @test nobs(dearusselloo) == 8
    @test ninputs(dearusselloo) == 1
    @test noutputs(dearusselloo) == 2
    @test efficiency(dearusselloo, :Y) ≈ [1 1; 1 1; 1 1; 7/3 1.4; 7/3 7/3; 1 2; 7/6 7/4; 46/9 1]
    @test_throws ArgumentError efficiency(dearusselloo, :X) # No input efficiency in input oriented model
    @test efficiency(dearusselloo) ≈ [1.0; 1.0; 1.0; 1.8666666666666665; 2.333333333333333; 1.5; 1.4583333333333335; 3.0555555555555554]
    @test convert(Matrix, peers(dearusselloo)) ==
            [ 1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  1.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
            1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
            1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  0.0  0.0  0.0  0.0
            1.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
            1/3  0.0  2/3  0.0  0.0  0.0  0.0  0.0]
    @test slacks(dearusselloo, :X) ≈ zeros(8,1)
    @test slacks(dearusselloo, :Y) ≈ zeros(8,2) atol = 1e-14

    @test efficiency(dearussell(targets(dearusselloo, :X), targets(dearusselloo, :Y), orient = :Output, rts = :CRS)) ≈ ones(8,1)

    # Test no slacks
    dearusselloonoslack = dearussell(X, Y, orient = :Output, slack = false)
    @test efficiency(dearusselloonoslack) == efficiency(dearusselloo)
    @test isempty(slacks(dearusselloonoslack, :X)) == 1
    @test isempty(slacks(dearusselloonoslack, :Y)) == 1

    # Output oriented VRS
    dearusselloovrs  = dearussell(X, Y, orient = :Output, rts = :VRS)
    @test efficiency(dearusselloovrs) ≈ [1.0; 1.0; 1.0; 1.8666666666666665; 2.333333333333333; 1.5; 1.4583333333333335; 3.0555555555555554]

    ## Test if one-by-one DEA using evaluation and reference sets match initial results
    dearussell_oo_ref_eff = zeros(size(X, 1))
    dearussell_oovrs_ref_eff = zeros(size(X, 1))

    Xref = X[:,:]
    Yref = Y[:,:]

    for i = 1:size(X, 1)
        Xeval = X[i:i,:]
        Xeval = Xeval[:,:]
        Yeval = Y[i:i,:]
        Yeval = Yeval[:,:]

        dearussell_oo_ref_eff[i] = efficiency(dearussell(Xeval, Yeval, orient = :Output, rts = :CRS, Xref = Xref, Yref = Yref))[1]
        dearussell_oovrs_ref_eff[i] = efficiency(dearussell(Xeval, Yeval, orient = :Output, rts = :VRS, Xref = Xref, Yref = Yref))[1]
    end

    @test dearussell_oo_ref_eff == efficiency(dearusselloo)
    @test dearussell_oovrs_ref_eff == efficiency(dearusselloovrs)

    # Add new DMU with slack
    dearussellooslack = dearussell([X; 0.5], [Y; 7 7], orient = :Output, rts = :VRS)
    @test slacks(dearussellooslack, :X) ≈ [0.5; 0; 0; 0.5; 0.5; 0; 0.5; 1/6; 0]

    # Print
    show(IOBuffer(), dearusselloo)
    show(IOBuffer(), dearusselloonoslack)

    # ------------------
    # Graph
    # ------------------
    X = [2; 4; 8; 12; 6; 14; 14; 9.412]
    Y = [1; 5; 8; 9; 3; 7; 9; 2.353]

    # Graph CRS
    dearussellgr = dearussell(X, Y, orient = :Graph, rts = :CRS)

    @test nobs(dearussellgr) == 8
    @test ninputs(dearussellgr) == 1
    @test noutputs(dearussellgr) == 1
    @test efficiency(dearussellgr, :X) ≈ [ 0.632455;
        1.0;
        0.894427;
        0.774597;
        0.632456;
        0.632456;
        0.717137;
        0.447214] atol = 1e-5
    @test efficiency(dearussellgr, :Y) ≈ [ 1.581139;
        1.0;
        1.118034;
        1.290994;
        1.581139;
        1.581139;
        1.394433;
        2.236068] atol = 1e-5
    @test efficiency(dearussellgr) ≈ [ 0.632455;
        1.0;
        0.894427;
        0.774597;
        0.632455;
        0.632455;
        0.717137;
        0.447214] atol = 1e-5

    @test isempty(slacks(dearussellgr, :X)) == 1
    @test isempty(slacks(dearussellgr, :Y)) == 1

    @test efficiency(dearussell(targets(dearussellgr, :X), targets(dearussellgr, :Y), orient = :Graph, rts = :CRS)) ≈ ones(8,1) atol = 1e-6

    # Graph VRS
    dearussellgrvrs  = dearussell(X, Y, orient = :Graph, rts = :VRS)

    @test efficiency(dearussellgrvrs, :X) ≈ [  1.0;
        1.0;
        1.0;
        1.0;
        2/3;
        0.571429;
        0.857143;
        0.424989] atol = 1e-5
    @test efficiency(dearussellgrvrs, :Y) ≈ [ 1.0;
        1.0;
        1.0;
        1.0;
        5/3;
        1.142857;
        1.0;
        2.124947] atol = 1e-5
    @test efficiency(dearussellgrvrs) ≈ [ 1.0;
        1.0;
        1.0;
        1.0;
        0.633333;
        0.723214;
        0.928571;
        0.447795] atol = 1e-5

    ## Test if one-by-one DEA using evaluation and reference sets match initial results
    dearussell_gr_ref_eff = zeros(size(X, 1))
    dearussell_grvrs_ref_eff = zeros(size(X, 1))

    Xref = X[:,:]
    Yref = Y[:,:]

    for i = 1:size(X, 1)
        Xeval = X[i:i,:]
        Xeval = Xeval[:,:]
        Yeval = Y[i:i,:]
        Yeval = Yeval[:,:]

        dearussell_gr_ref_eff[i] = efficiency(dearussell(Xeval, Yeval, orient = :Graph, rts = :CRS, Xref = Xref, Yref = Yref))[1]
        dearussell_grvrs_ref_eff[i] = efficiency(dearussell(Xeval, Yeval, orient = :Graph, rts = :VRS, Xref = Xref, Yref = Yref))[1]
    end

    @test dearussell_gr_ref_eff == efficiency(dearussellgr)
    @test dearussell_grvrs_ref_eff == efficiency(dearussellgrvrs)

    # Print
    show(IOBuffer(), dearussellgr)

    # ------------------
    # Test errors
    # ------------------
    @test_throws DimensionMismatch dearussell([1; 2 ; 3], [4 ; 5]) #  Different number of observations
    @test_throws DimensionMismatch dearussell([1; 2], [4 ; 5], Xref = [1; 2; 3; 4]) # Different number of observations in reference sets
    @test_throws DimensionMismatch dearussell([1 1; 2 2], [4 4; 5 5], Xref = [1 1 1; 2 2 2]) # Different number of inputs
    @test_throws DimensionMismatch dearussell([1 1; 2 2], [4 4; 5 5], Yref = [4 4 4; 5 5 5]) # Different number of inputs
    @test_throws ArgumentError dearussell([1; 2; 3], [4; 5; 6], orient = :Error) # Invalid orientation
    @test_throws ArgumentError dearussell([1; 2; 3], [4; 5; 6], rts = :Error) # Invalid returns to scale

    @test_throws ArgumentError efficiency(dearussellio, :Error) # Invalid efficiency type

end


text:
using Random
using Statistics
using StatsBase
using Plots
import Base.+
import Base.-

const λ  = 0.002 #in mm => 150GHz
const θᵦ = (π / 180.) * 0.17   #0.17 deg for STRIP W-band
const ω₀ = λ / ( π * θᵦ ) # θᵦ is the beam opening angle in radians
const z_atm = 3E4 #depends on the observation site
const Tᵧ = 290 # ground temperature in kelvin
const χ₁_0 = 1 # in m^{-2}
const χ₂_0 = 1 # in m^{-2}
const z_0  = 2E3 # the quote where water vapor vanish
const L₀   = 100 # in m, the turbulance correlation length

struct Point{T}
	x::T
 	y::T
	z::T
end

function(p::Point)()
	x = p.x * cos(p.y) * cos(p.z)
	y = p.x * cos(p.y) * sin(p.z)
	z = p.x * sin(p.y)
	return Point(x, y, z)
end

Base.:+(a::Point{Float64}, b::Point{Float64}) = Point{Float64}(a.x+b.x , a.y+b.y, a.z+b.z)
Base.:-(a::Point{Float64}, b::Point{Float64}) = Point{Float64}(a.x-b.x , a.y-b.y, a.z-b.z)
Base.:*(a::Point{Float64}, b::Point{Float64}) = Point{Float64}(a.x*b.x , a.y*b.y, a.z*b.z)
Base.:/(a::Point{Float64}, b::Float64) = Point{Float64}(a.x/b , a.y/b, a.z/b)

r = 0:20000

punti = [Point{Float64}(i, π / 2 , 0.)()   for i in r]
z = [i.z for i in punti]
ww = [ ω₀ * ( 1 + ( (λ * z) / (π * ω₀^2)  )^2)^0.5  for z in r]

plot( ww[2:length(r)] ./ ww[2]  , z[2:length(r)], yaxis=:log, xlims=(0,2), ylims = (1,2E4))

ext: jl

max_stars_repo_name: gruffoc/AFO

max_issues_repo_name: gruffoc/AFO

max_forks_repo_name: gruffoc/AFO


text:
using Test
using Base.Iterators: take
using Random
using LinearAlgebra

using MeasureBase
using MeasureBase: logpdf

using Aqua
Aqua.test_all(MeasureBase; ambiguities=false, unbound_args=false)

d = ∫exp(x -> -x^2, Lebesgue(ℝ))

# if VERSION ≥ v"1.6"
#     @eval using JETTest

#     @eval begin
#         @test_nodispatch density(Lebesgue(ℝ), 0.3)

#         @test_nodispatch density(Dirac(0), 0.3)
#         @test_nodispatch density(Dirac(0), 0)

#         @test_nodispatch density(d, 3)

#         @test_nodispatch basemeasure(d)

#         @test_nodispatch logdensity(For(3) do j Dirac(j) end, [1,2,3])
#     end
# end

# function draw2(μ)
#     x = rand(μ)
#     y = rand(μ)
#     while x == y
#         y = rand(μ)
#     end
#     return (x,y)
# end

function test_testvalue(μ)
    logdensity(μ, testvalue(μ)) isa AbstractFloat
end

test_measures = [
    # Chain(x -> Normal(μ=x), Normal(μ=0.0))
    For(3) do j Dirac(j) end
    For(2,3) do i,j Dirac(i) + Dirac(j) end
    Lebesgue(ℝ) ^ 3
    Lebesgue(ℝ) ^ (2,3)
    3 * Lebesgue(ℝ)
    Dirac(π)
    Lebesgue(ℝ)
    Dirac(0.0) + Lebesgue(ℝ)
    SpikeMixture(Lebesgue(ℝ), 2)
    d ⊙ d
]

testbroken_measures = [
    # InverseGamma(2) # Not defined yet
    # MvNormal(I(3)) # Entirely broken for now
    CountingMeasure(Float64)
    Likelihood

    TrivialMeasure()
]

@testset "testvalue" begin
    for μ in test_measures
        @info "testing $μ"
        @test test_testvalue(μ)
    end

    for μ in testbroken_measures
        @info "testing $μ"
        @test_broken test_testvalue(μ)
    end

end

# @testset "Kernel" begin
#     κ = MeasureBase.kernel(MeasureBase.Dirac, identity)
#     @test rand(κ(1.1)) == 1.1
# end

@testset "SpikeMixture" begin
    @test rand(SpikeMixture(Dirac(0), 0.5)) == 0
    @test rand(SpikeMixture(Dirac(1), 1.0)) == 1
    w = 1/3
    m = SpikeMixture(d, w)
    bm = basemeasure(m)
    @test (bm.s*bm.w)*bm.m == 1.0*basemeasure(d)
    @test density(m, 1.0)*(bm.s*bm.w) ≈ w*density(d,1.0)
    @test density(m, 0)*(bm.s*(1-bm.w)) ≈ (1-w)
end

@testset "Dirac" begin
    @test rand(Dirac(0.2)) == 0.2
    @test logdensity(Dirac(0.3), 0.3) == 0.0
    @test logdensity(Dirac(0.3), 0.4) == -Inf
end

@testset "AffineTransform" begin
    f = AffineTransform((μ=3,σ=2))
    @test f(inv(f)(1)) == 1
    @test inv(f)(f(1)) == 1

    f = AffineTransform((μ=3,ω=2))
    @test f(inv(f)(1)) == 1
    @test inv(f)(f(1)) == 1

    f = AffineTransform((σ=2,))
    @test f(inv(f)(1)) == 1
    @test inv(f)(f(1)) == 1

    f = AffineTransform((ω=2,))
    @test f(inv(f)(1)) == 1
    @test inv(f)(f(1)) == 1

    f = AffineTransform((μ=3,))
    @test f(inv(f)(1)) == 1
    @test inv(f)(f(1)) == 1
end

@testset "Affine" begin
    unif = ∫(x -> 0<x<1, Lebesgue(ℝ))
    f1 = AffineTransform((μ=3.,σ=2.))
    f2 = AffineTransform((μ=3.,ω=2.))
    f3 = AffineTransform((μ=3.,))
    f4 = AffineTransform((σ=2.,))
    f5 = AffineTransform((ω=2.,))

    for f in [f1,f2,f3,f4,f5]
        par = getfield(f, :par)
        @test Affine(par)(unif) == Affine(f, unif)
        @test density(Affine(f, Affine(inv(f), unif)), 0.5) == 1
    end

    d = ∫exp(x -> -x^2, Lebesgue(ℝ))

    μ = randn(3)
    σ = LowerTriangular(randn(3,3))
    ω = inv(σ)

    x = randn(3)

    @test logdensity(Affine((μ=μ,σ=σ), d^3), x) ≈  logdensity(Affine((μ=μ,ω=ω), d^3), x)
    @test logdensity(Affine((σ=σ,), d^3), x) ≈  logdensity(Affine((ω=ω,), d^3), x)
    @test logdensity(Affine((μ=μ,), d^3), x) ≈  logdensity(d^3, x-μ)

    d = ∫exp(x -> -x^2, Lebesgue(ℝ))
    a = Affine((σ=[1 0]',), d^1)
    x = randn(2)
    y = randn(1)
    @test logpdf(a, x) ≈ logpdf(d, inv(a.f)(x)[1])
    @test logpdf(a, a.f(y)) ≈ logpdf(d^1, y)

    b = Affine((ω=[1 0]'',), d^1)
    @test logpdf(b, x) ≈ logpdf(d, inv(b.f)(x)[1])
    @test logpdf(b, b.f(y)) ≈ logpdf(d^1, y)
end

@testset "For" begin
    FORDISTS = [
        For(1:10) do j Dirac(j) end
        For(4,3) do i,j Dirac(i) ⊗ Dirac(j) end
        For(1:4, 1:4) do i,j Dirac(i) ⊗ Dirac(j) end
        For(eachrow(rand(4,2))) do x Dirac(x[1]) ⊗ Dirac(x[2]) end
        For(rand(4), rand(4)) do i,j Dirac(i) ⊗ Dirac(j) end
    ]

    for d in FORDISTS
        @info "testing $d"
        @test logdensity(d, rand(d)) isa Float64
    end
end

@testset "Half" begin
    Normal() = ∫exp(x -> -0.5x^2, Lebesgue(ℝ))
    @half Normal
    @test logdensity(HalfNormal(), Lebesgue(ℝ), -0.2) == -Inf
    @test logdensity(HalfNormal(), 0.2) == logdensity(Normal(), 0.2)
    @test density(HalfNormal(), Lebesgue(ℝ), 0.2) ≈ 2 * density(Normal(), Lebesgue(ℝ), 0.2)
end

# import MeasureBase.:⋅
# function ⋅(μ::Normal, kernel)
#     m = kernel(μ)
#     Normal(μ = m.μ.μ, σ = sqrt(m.μ.σ^2 + m.σ^2))
# end
# struct AffineMap{S,T}
#     B::S
#     β::T
# end
# (a::AffineMap)(x) = a.B*x + a.β
# (a::AffineMap)(p::Normal) = Normal(μ = a.B*mean(p) + a.β, σ = sqrt(a.B*p.σ^2*a.B'))

# @testset "Likelihood" begin
#     dps = [
#         (Normal()                             ,    2.0  )
#         # (Pushforward(as((μ=asℝ,)), Normal()^1), (μ=2.0,))
#     ]

#     ℓs = [
#         Likelihood(Normal{(:μ,)},              3.0)
#         Likelihood(kernel(Normal, x -> (μ=x, σ=2.0)), 3.0)
#     ]

#     for (d,p) in dps
#         for ℓ in ℓs
#             @test logdensity(d ⊙ ℓ, p) == logdensity(d, p) + logdensity(ℓ, p)
#         end
#     end
# end

# @testset "ProductMeasure" begin
#     d = For(1:10) do j Poisson(exp(j)) end
#     x = Vector{Int16}(undef, 10)
#     @test rand!(d,x) isa Vector
#     @test rand(d) isa Vector

#     @testset "Indexed by Generator" begin
#         d = For((j^2 for j in 1:10)) do i Poisson(i) end
#         x = Vector{Int16}(undef, 10)
#         @test rand!(d,x) isa Vector
#         @test_broken rand(d) isa Base.Generator
#     end

#     @testset "Indexed by multiple Ints" begin
#         d = For(2,3) do μ,σ Normal(μ,σ) end
#         x = Matrix{Float16}(undef, 2, 3)
#         @test rand!(d, x) isa Matrix
#         @test_broken rand(d) isa Matrix{Float16}
#     end
# end

@testset "Show methods" begin
    @testset "PowerMeasure" begin
        @test repr(Lebesgue(ℝ) ^ 5) == "Lebesgue(ℝ) ^ 5"
        @test repr(Lebesgue(ℝ) ^ (3, 2)) == "Lebesgue(ℝ) ^ (3, 2)"
    end
end

# @testset "Density measures and Radon-Nikodym" begin
#     x = randn()
#     let d = ∫(𝒹(Cauchy(), Normal()), Normal())
#         @test logdensity(d, x) ≈ logdensity(Cauchy(), x)
#     end

#     let f = 𝒹(∫(x -> x^2, Normal()), Normal())
#         @test f(x) ≈ x^2
#     end

#     let d = ∫exp(log𝒹(Cauchy(), Normal()), Normal())
#         @test logdensity(d, x) ≈ logdensity(Cauchy(), x)
#     end

#     let f = log𝒹(∫exp(x -> x^2, Normal()), Normal())
#         @test f(x) ≈ x^2
#     end
# end


text:
# Adapted from https://github.com/julia-vscode/julia-vscode/blob/1da85dc37b54723c1bc4b313d2c11a974c7656b4/scripts/tasks/task_compileenv.jl instead
import Pkg, Libdl, PackageCompiler

# crude CLI toggle. Either:
# julia ./deps/create_sysimage.jl   # defaults
# julia ./deps/create_sysimage.jl --only_package --no_precompile_tests   #change to only use the one package and not precompile tests, etc.
# julia ./deps/create_sysimage.jl --replace_default # replace the default sysimage.  Ony intended for batch environments with vscode/etc.

all_packages = !("--only_package" in ARGS)
precompile_tests = !("--no_precompile_tests" in ARGS)
replace_default = !("--replace_default" in ARGS)

project = dirname(@__DIR__)  # Must be one level deep below the project file
sysimage_path = joinpath(project, "JuliaSysimage.$(Libdl.dlext)")  # consistent naming with vscode
project_filename = isfile(joinpath(project, "JuliaProject.toml")) ? joinpath(project, "JuliaProject.toml") : joinpath(project, "Project.toml")
project_object = Pkg.API.read_project(project_filename)
packages = all_packages ? Symbol.(collect(keys(project_object.deps))) : Symbol(project_object.name)
precompile_execution_file = precompile_tests ? joinpath(project, "test/runtests.jl") : String[]

precompile_tests && @info "Using unit tests for precompile execution file."
@info "Precompiling packages: $packages: "
@info "Now building a custom sysimage for the environment '$project'."
if replace_default
    @info "Now replacing the default sysimage with one for the environment '$project'"
    PackageCompiler.create_sysimage(packages; replace_default, project, precompile_execution_file)
else
    @info "Now building a custom sysimage for the environment '$project' to $sysimage_path"
    PackageCompiler.create_sysimage(packages; sysimage_path, project, precompile_execution_file)
end


text:
# p11.jl - Chebyshev differentation of a smooth function

xx = -1:.01:1; uu = @. exp(xx)*sin(5*xx); clf();
for N = [10 20]
    (D,x) = cheb(N); u = @. exp(x)*sin(5*x);
    axes([.15,.66-.4*(N==20),.31,.28]);
    plot(x,u,".",markersize=6); grid(true);
    plot(xx,uu);
    title("u(x),  N=$N")
    error = D*u - @. exp(x)*(sin(5*x)+5*cos(5*x));
    axes([.55,.66-.4*(N==20),.31,.28]);
    plot(x,error,".-",markersize=6); grid(true);
    title("error in u'(x),  N=$N");
end


text:
module Float8s

    import Base: (-),(==),(<),(<=),isless,bitstring,
                isnan,iszero,one,zero,abs,isfinite,
                floatmin,floatmax,typemin,typemax,
                Float16,Float32,Float64,
                UInt8,Int8,Int16,Int32,Int64,
                (+), (-), (*), (/), (\), (^),
                sin,cos,tan,asin,acos,atan,sinh,cosh,tanh,asinh,acosh,
                atanh,exp,exp2,exp10,expm1,log,log2,log10,sqrt,cbrt,log1p,
                atan,hypot,round,show,nextfloat,prevfloat,
                promote_rule, sign, signbit

    export Float8, Float8_4, NaN8, Inf8, NaN8_4, Inf8_4

    include("float8.jl")
    include("float8_to_float32.jl")
    include("float32_to_float8.jl")
    include("float32_to_float8_old.jl")

end


text:
import COBRA

export AVCobraModel

"""
Minimal metabolic network exhibiting a switch, in COBRA's format.
The network objective is ATP production.
    Vg, Vl: uptake bounds of glucose and lactate
    R: maximum respiratory rate
    L: large bound, used in place of ∞ (default: 1000)
"""
function AVCobraModel(Vg::Real, Vl::Real = 0; R = 0.9, L = 1000)
    @assert 0 ≤ Vg < Inf && 0 ≤ Vl < Inf
    @assert 0 ≤ R < Inf
    @assert 0 < L < Inf

    rxns = ["vg", "vl", "vo", "vatp"]
    mets = ["pyr", "atp"]

    S = [2 1 -1 0; 2 0 18 -1]
    b = zeros(2)
    c = [0, 0, 0, 1]
    lb = [-L, -L, 0, -L]
    ub = [Vg, Vl, R, L]
    osense = 1
    csense = ['E', 'E']

    COBRA.LPproblem(S, b, c, lb, ub, osense, csense, rxns, mets)
end


text:
using Primes

const L = 2*10^4
acnt = 0

println("Amicable pairs not greater than ", L)

for i in 2:L
    !isprime(i) || continue
    j = divisorsum(i)
    j < i && divisorsum(j) == i || continue
    acnt += 1
    println(@sprintf("%4d", acnt), " => ", j, ", ", i)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
include("../src/advection.jl")
include("../src/spline.jl")
include("../src/bspline.jl")
include("../src/bsplinelu.jl")
include("../src/bsplinefft.jl")
include("../src/lagrange.jl")

function aff_graph(ind, fp)
    println("ind=$ind begin")
    minval = 0
    maxval = 0
    for (i, val) in enumerate(fp)
        #        println("$i\t$val")
        minval = min(val, minval)
        maxval = max(val, maxval)
    end
    println("minval=$minval maxval=$maxval")
    println("ind=$ind end")
end
function test_dirac(bsp::InterpolationType, order, len, nb, modval)
    #        println("lag=$lag")
    n = len
    # whennot circuar coef != 1 set the function unperiodic
    fp = zeros(BigFloat, n)
    fp[div(n, 2)] = big"1.0"
    fi = zeros(BigFloat, n)
    value = big"0.30529810681113334445566767713091"
    normax = 0.0
    aff_graph(1, fp)
    for i = 1:nb
        fi .= fp
        interpolate!(missing, fp, fi, value, bsp)
        if i % modval == 0
            aff_graph(i, fp)
        end
    end
end
function test_dirac_lagrange(order, len, nb, modval)
    lag = Lagrange(BigFloat, order; iscirc = true)
    test_dirac(lag, order, len, nb, modval)
end
function test_dirac_splu(order, len, nb, modval)
    bsp = B_SplineLU(order, len, BigFloat; iscirc = true)
    test_dirac(bsp, order, len, nb, modval)
end
function test_dirac_spfft(order, len, nb, modval)
    bsp = B_SplineFFT(order, len, BigFloat)
    test_dirac(bsp, order, len, nb, modval)
end

test_dirac_lagrange(31, 256, 500, 50)
test_dirac_splu(31, 256, 500, 50)
test_dirac_spfft(31, 256, 500, 50)


text:
module Scanivalve

using Sockets
using AbstractDAQs

export AbstractDAQ, DSA3217
export close, numchans, listscan, listzero, listgain, listoffset, listdelta
export scan!, stopscan, isreading, samplesread, scanconfig, daqparam
export clearbuffer!
export daqread, daqstart, daqacquire, daqconfig, daqconfigdev
export daqzero, daqstop, daqaddinput, numchannels, daqchannels
export DAQConfig, iparameters, fparameters, sparameters, devname
export daqdevip, daqdevmodel, daqdevserialnum, daqdevtag

"Base type for Scanivalve pressure scanners"
abstract type AbstractScanivalve <: AbstractPressureScanner end

include("dsa3200packets.jl")
include("dsa3217.jl")

end


text:
# István Maros. Computational Techniques of the Simplex Method. Kluwer Academic Publishers, Boston, 2003. [CTSM]
# Robert J. Vanderbei. Linear Programming: Foundations and Extensions. Springer, second edition, 2001. [LPFE]
# Suhl, Suhl. A fast LU update for linear programming. [LULP]

#=
	List of algorithms

	Linear Programming
		I. Solve, Revised, Canonical Simplex.
			1. revised simplex, basis inversion, no reuse
			2. revised simplex, basis factorization, no reuse
			3. revised simplex, basis factorization, reuse, no re-factorization
			4. revised simplex, basis factorization, reuse, re-factorization
			5. revised simplex, basis factorization, reuse, re-factorization, degen (CTSM p230), sing
			6. 5. with 'textbook' Phase-I
			7. 5. with Maros' Phase-I (p203)

		II. Solve, Revised, Sparse, Canonical Simplex.
			1. revised sparse simplex

		III. Generalized Simplex variants of I, then II (Presolve, Postsolve)
			1. generalized version of I.1, using Maros' CF-2 (p13-17), but without translation, calling it CF-2b. (Handwritten p29).
			2. 1. with variable transformation for less lohi calculations, the proper CF-2 formulation
			3. 2. with infeasibility tolerance, type-0 handling and Harris ratio. (CSTM p178)

	LU Factorization
		I. Suhl-suhl (CTSM p136), and original [LULP]
=#

module Jad

module lp_I_1
	using Dcd
	using Lp

	type Dat{T}
		prob::Lp.Cf0_problem
		sol::Lp.Solution
		maxit::Int
		n::Int
		m::Int
		status::Symbol

		# Using CTSM notation.
		# LPFE notation translation.
		# αq - ΔxB, β - xB*, q - j, i - p, t - θ

		iB::Array{Int, 1}
		iR::Array{Int, 1}
		B::Array{T, 2}
		Binv::Array{T, 2}
		β::Array{T, 1}
		π::Array{T, 1}
		cBT::Array{T, 2}
		dJ::Array{T, 1}
		αq::Array{T, 1}
		z::T
		zero::T

		dcd::Dcd.Session
		it::Int
		q::Int
		p::Int
		θ::T

		Dat() = new()
	end

	function construct_dat(T::DataType) return Dat{T}() end

	function fill_dat{T}(prob::Lp.Cf0_problem, dat::Dat{T})
		n = prob.n; m = prob.m;
		dat.maxit = get(prob.params, "maxit", 0)
		dat.prob = prob
		dat.n = n
		dat.m = m
		dat.iB = Array(Int, m)
		dat.iR = Array(Int, n)
		dat.B = Array(T, (m,m))
		#dat.Binv = Array(T, (m,m))
		dat.β = Array(T, m)
		#dat.π = Array(T, m+n)
		#dat.cBT = Array(T, m)
		dat.dJ = Array(T, n)
		#dat.αq = Array(T, n)
		dat.zero = zero(T)
	end

	function warn{T}(dat::Dat{T}, warning::Symbol)
		println("Warning:", warning)
	end

	function succeed{T}(dat::Dat{T})
		dat.status = :Ended
		sol = dat.sol
		sol.iter = dat.it
		sol.solved = true
		sol.status = :Optimal
		sol.z = dat.z
		sol.x = zeros(dat.prob.conv.t, dat.n)
		for i = 1:length(dat.iB)
			xi = dat.iB[i]
			if (xi <= dat.n) sol.x[xi] = dat.β[i]; end
		end
	end

	function fail{T}(dat::Dat{T}, status::Symbol)
		dat.status = :Ended
		sol = dat.sol
		sol.iter = dat.it
		sol.solved = false
		sol.status = status
	end

	function sel_cols{T}(src::Matrix{T}, dst::Matrix{T}, cols::Vector{Int})
		for i = 1:length(cols)
			c = cols[i]
			for r = 1:size(dst)[1]
				dst[r,i] = src[r,c]
			end
		end
	end

	function set_basis_logical{T}(dat::Dat{T})
		dat.iB = [i for i in dat.n+1:dat.n+dat.m]
		dat.iR = [i for i in 1:dat.n]
	end

	function comp_B_R{T}(dat::Dat{T}) sel_cols(dat.prob.A, dat.B, dat.iB) end
	function comp_Binv{T}(dat::Dat{T}) dat.Binv = inv(dat.B) end
	function comp_cBT{T}(dat::Dat{T}) dat.cBT = transpose(dat.prob.c[dat.iB]) end
	function init_β{T}(dat::Dat{T}) dat.β = dat.Binv * (dat.prob.b) end
	function update_β{T}(dat::Dat{T}, p::Int, θ::T) dat.β = dat.β - (θ*dat.αq); dat.β[p] = θ; end
	function init_z{T}(dat::Dat{T}) dat.z = dot(reshape(dat.cBT, length(dat.cBT)), dat.β); end
	function update_z{T}(dat::Dat{T}, q::Int, θ::T) dat.z = dat.z + θ * dat.dJ[q]; end
	function comp_π{T}(dat::Dat{T}) dat.π = reshape(dat.cBT * dat.Binv, length(dat.cBT)); end
	function calc_dj{T}(dat::Dat{T}, j::Int) i = dat.iR[j]; dj = dat.prob.c[i] - dot(dat.π, dat.prob.A[:,i]); return dj; end
	function comp_dJ{T}(dat::Dat{T}) for j = 1:dat.n dat.dJ[j] = calc_dj(dat, j) end; end
	function comp_αq{T}(dat::Dat{T}, q::Int) aq = dat.iR[q]; dat.αq = dat.Binv * (dat.prob.A[:,aq]); end
	function check_optimal_dJ{T}(dat::Dat{T}) return all( dj->(dj >= dat.zero), dat.dJ ); end
	function check_feasible_β{T}(dat::Dat{T}) return all( β->(β >= dat.zero), dat.β ); end

	function pivot_iB_iR{T}(dat::Dat{T}, q::Int, p::Int) dat.iR[q], dat.iB[p] = dat.iB[p], dat.iR[q] end

	function price_full_dantzig{T}(dat::Dat{T})
		# [CTSM].p187
		min_i = 0; min_dj = dat.zero;
		for i = 1:length(dat.dJ)
			dj = dat.dJ[i]
			if (dj < dat.zero && dj <= min_dj)
				if (dj == min_dj)
					warn(dat, :Degen)
				end
				min_i = i; min_dj = dj;
			end
		end
		return min_i
	end

	function chuzro{T}(dat::Dat{T})
		min_i = 0; min_ratio = dat.zero;
		for i = 1:length(dat.αq)
			if (dat.αq[i] > dat.zero)
				ratio = dat.β[i] / dat.αq[i]
				if (min_i == 0 || ratio <= min_ratio)
					if (ratio == min_ratio)
						warn(dat, :Degen)
					end
					min_i = i; min_ratio = ratio;
				end
			end
		end
		return (min_i, min_ratio)
	end

	function solve_init{T}(dat::Dat{T}, sol::Lp.Solution{T})
		dat.sol = sol
		dat.dcd = sol.dcd
		Dcd.@it(dcd)
		dat.it = 0
		dat.status = :Iterating
	end

	function solve_step0{T}(dat::Dat{T})
		set_basis_logical(dat)
		Dcd.@set(dcd, "iB", dat.iB)
		#Initializations
		comp_cBT(dat); comp_B_R(dat); comp_Binv(dat);
		init_β(dat); init_z(dat);
		Dcd.@set(dcd, "β0", dat.β)
	end

	function solve_step_phaseI{T}(dat::Dat{T})
		if (check_feasible_β(dat) == false)
			warn(degen, :PhaseI)
			fail(dat, sol, :Infeasible)
		end
	end

	function solve_iter{T}(dat::Dat{T})
		return (dat.status == :Iterating && (dat.maxit == 0 || dat.it < dat.maxit))
	end

	function solve_next{T}(dat::Dat{T})
		dat.it = dat.it + 1
	end

	function solve_step1{T}(dat::Dat{T})
		Dcd.@it(dcd)
		#Step 1
		Dcd.@set(dcd, "B", dat.B); Dcd.@set(dcd, "Binv", dat.Binv);
		Dcd.@set(dcd, "cBT", dat.cBT);
		comp_π(dat); Dcd.@set(dcd, "π", dat.π);
	end

	function solve_step_price{T}(dat::Dat{T})
		comp_dJ(dat); Dcd.@set(dcd, "dJ", dat.dJ);
		if check_optimal_dJ(dat) succeed(dat); return; end
		dat.q = price_full_dantzig(dat); Dcd.@set(dcd, "q", dat.q);
	end

	function solve_step3{T}(dat::Dat{T})
		comp_αq(dat, dat.q); Dcd.@set(dcd, "αq", dat.αq);
	end

	function solve_step_chuzro{T}(dat::Dat{T})
		dat.p, dat.θ = chuzro(dat); Dcd.@set(dcd, "p", dat.p); Dcd.@set(dcd, "θ", dat.θ);
		if (dat.p == 0) fail(dat, :Unbounded); return; end
		Dcd.@set(dcd, "pivot", (dat.q, dat.p))
	end

	function solve_step5{T}(dat::Dat{T})
		pivot_iB_iR(dat, dat.q, dat.p)
	end

	function solve_update{T}(dat::Dat{T})
		update_β(dat, dat.p, dat.θ); Dcd.@set(dcd, "β", dat.β);
		update_z(dat, dat.q, dat.θ); Dcd.@set(dcd, "z", dat.z);
		comp_cBT(dat); comp_B_R(dat); comp_Binv(dat);
	end

	function solve_dat{T}(dat::Dat{T}, sol::Lp.Solution{T})
		# [CTSM].p33,p30, but with naive basis reinversion instead of update.
		solve_init(dat, sol)
		solve_step0(dat)
		solve_step_phaseI(dat)

		while solve_iter(dat)
			solve_step1(dat)
			solve_step_price(dat);
			if (dat.status != :Iterating) break; end
			solve_step3(dat)
			solve_step_chuzro(dat); if (dat.status != :Iterating) break; end
			solve_step5(dat)
			solve_update(dat)
			solve_next(dat)
		end
		if (dat.status == :Iterating) fail(dat, :Maxit); end
	end

end

module lp_III_1
	using Dcd
	using Lp
	using Jad.lp_I_1

	type Dat{T}
		Base::Module
		prob::Lp.Cf2b_problem{T}
		bdat::lp_I_1.Dat{T}
		is_hi::Vector{Bool}

		p_hi::Bool

		Dat() = new()
	end

	function construct_dat(T::DataType) return Dat{T}() end

	function fill_dat{T}(prob::Lp.Cf2b_problem{T}, dat::Dat{T})
		dat.Base = lp_I_1
		dat.prob = prob
		dat.bdat = dat.Base.construct_dat(T)
		Base.fill_dat(prob.base, dat.bdat)
		dat.is_hi = zeros(Bool, dat.prob.n+dat.prob.m) # all variables at low.
	end

	function init_β{T}(dat::Dat{T}, xlo::Vector{T})
		R = dat.base.prob.A[:,dat.iR]
		dat.β = dat.Binv * (dat.prob.b - R*xlo)
	end

	function init_z{T}(dat::Dat{T}, xlo::Vector{T})
		cR = dat.prob.c[dat.iR]
		zR = dot(cR, xlo)
		zB = dot(reshape(dat.cBT, length(dat.cBT)), dat.β)
		dat.z = zB + zR
	end

	function solve_step0{T}(sdat::Dat{T})
		Base = dat.Base; dat = sdat.bdat;
		Base.set_basis_logical(dat)
		Dcd.@set(dcd, "iB", dat.iB)
		#Initializations
		Base.comp_cBT(dat); Base.comp_B_R(dat); Base.comp_Binv(dat);
		xlo = dat.prob.l[dat.iR]
		init_β(sdat, xlo); init_z(sdat, xlo);
		Dcd.@set(dcd, "β0", dat.β)
	end

	function price_full_dantzig{T}(sdat::Dat{T})
		# [CTSM].p187
		Base = sdat.Base; dat = sdat.bdat;
		min_i = 0; min_dj = dat.zero;
		for i = 1:length(dat.dJ)
			dj = dat.dJ[i]
			if (dj < dat.zero && dj <= min_dj && (sdat.is_hi[dat.iB[dj]] == false) )
				if (dj == min_dj)
					Base.warn(dat, :Degen)
				end
				min_i = i; min_dj = dj;
			end
		end
		return min_i
	end

	function solve_step_price{T}(sdat::Dat{T})
		Base = sdat.Base; dat = sdat.bdat;
		Base.comp_dJ(dat); Dcd.@set(dcd, "dJ", dat.dJ);
		if Base.check_optimal_dJ(dat) Base.succeed(dat); return; end
		dat.q = price_full_dantzig(sdat); Dcd.@set(dcd, "q", dat.q);
	end

	function chuzro{T}(dat::Dat{T})
		# Handwritten notes (p29)
		glob_i = 0; glob_max_θ = dat.zero; glob_hi = false;
		for i = 1:length(dat.αq)
			set_hi = false
			if (dat.αq[i] > dat.zero)
				lo = dat.prob.l[dat.iB[i]]
				max_θ = (dat.β[i] - lo) / dat.αq[i]
				set_hi = false
			elseif (dat.αq[i] < dat.zero)
				hi = dat.prob.h[dat.iB[i]]
				max_θ = (dat.β[i] - hi) / dat.αq[i]
				set_hi = true
			else
				continue
			end
			if (glob_i == 0 || max_θ >= glob_max_θ)
				if (max_θ == globa_max_θ) warn(dat, :Degen); end
				glob_i = i; glob_max_θ = max_θ; glob_hi = set_hi;
			end
		end

		return (glob_i, glob_max_θ, glob_hi)
	end

	function solve_step_chuzro{T}(dat::Dat{T})
		dat.p, dat.θ, dat.p_hi = chuzro(dat); Dcd.@set(dcd, "p", dat.p); Dcd.@set(dcd, "θ", dat.θ);
		if (dat.p == 0) fail(dat, :Unbounded); return; end
		Dcd.@set(dcd, "pivot", (dat.q, dat.p))
	end

	function solve_update{T}(dat::Dat{T})
		dat.Base.solve_update(dat.bdat)
		dat.is_hi[dat.p] = dat.p_hi
	end

	function solve_dat{T}(sdat::Dat{T}, sol::Lp.Solution{T})
		# lp_I_1, with generalized forumulation (w/o translation)
		Base = sdat.Base;	bdat = sdat.bdat;
		Base.solve_init(bdat, sol)
		solve_step0(dat)
		Base.solve_step_phaseI(bdat)

		while solve_iter(bdat)
			Base.solve_step1(bdat)
			solve_step_price(dat); if (bdat.status != :Iterating) break; end
			Base.solve_step3(bdat)
			solve_step_chuzro(dat); if (bdat.status != :Iterating) break; end
			Base.solve_step5(bdat)
			solve_update(dat)
			Base.solve_next(bdat)
		end
		if (dat.status == Iterating) Base.fail(bdat, :Maxit); end
	end

end

module lu_I
	#function factorize{}(B::Array{T, 2})
	#end
end

#=
module lp_I_2
	using Dcd
	using Lp

	type Dat{T}
		prob::Lp.Cf0_problem
		maxit::Int
		n::Int
		m::Int
		# Using CTSM notation.
		iB::Array{Int, 1}
		iR::Array{Int, 1}
		B::Array{T, 2}
		BL::Array{T, 2}
		BU::Array{T, 2}
		BP::Array{T, 1}
		β::Array{T, 1}
		π::Array{T, 1}
		cBT::Array{T, 2}
		dJ::Array{T, 1}
		αq::Array{T, 1}
		z::T
		zero::T

		# LPFE notation translation.
		# αq - ΔxB, β - xB*, q - j, i - p, t - θ

		Dat() = new()
	end

	function fill_dat{T}(prob::Lp.Cf0_problem, dat::Dat{T})
		n = prob.n; m = prob.m;
		dat.maxit = get(prob.params, "maxit", 0)
		dat.prob = prob
		dat.n = n
		dat.m = m
		dat.iB = Array(Int, m)
		dat.iR = Array(Int, n)
		dat.B = Array(T, (m,m))
		#dat.Binv = Array(T, (m,m))
		dat.β = Array(T, m)
		#dat.π = Array(T, m+n)
		#dat.cBT = Array(T, m)
		dat.dJ = Array(T, n)
		#dat.αq = Array(T, n)
		dat.zero = zero(T)
	end

	function sel_cols{T}(src::Matrix{T}, dst::Matrix{T}, cols::Vector{Int})
		for i = 1:length(cols)
			c = cols[i]
			for r = 1:size(dst)[1]
				dst[r,i] = src[r,c]
			end
		end
	end

	function set_basis_logical{T}(dat::Dat{T})
		dat.iB = [i for i in dat.n+1:dat.n+dat.m]
		dat.iR = [i for i in 1:dat.n]
	end

	function comp_B_R{T}(dat::Dat{T}) sel_cols(dat.prob.A, dat.B, dat.iB) end
	function comp_BLU{T}(dat::Dat{T}) dat.BL, dat.BU, dat.BP = lu(dat.B) end
	function comp_cBT{T}(dat::Dat{T}) dat.cBT = transpose(dat.prob.c[dat.iB]) end
	function init_β{T}(dat::Dat{T}) dat.β = dat.Binv * (dat.prob.b) end
	function update_β{T}(dat::Dat{T}, p::Int, θ::T) dat.β = dat.β - (θ*dat.αq); dat.β[p] = θ; end
	function init_z{T}(dat::Dat{T}) dat.z = dot(reshape(dat.cBT, length(dat.cBT)), dat.β); end
	function update_z{T}(dat::Dat{T}, q::Int, θ::T) dat.z = dat.z + θ * dat.dJ[q]; end
	function comp_π{T}(dat::Dat{T}) dat.π = reshape(dat.cBT * dat.Binv, length(dat.cBT)); end
	function calc_dj{T}(dat::Dat{T}, j::Int) i = dat.iR[j]; dj = dat.prob.c[i] - dot(dat.π, dat.prob.A[:,i]); return dj; end
	function comp_dJ{T}(dat::Dat{T}) for j = 1:dat.n dat.dJ[j] = calc_dj(dat, j) end; end
	function comp_αq{T}(dat::Dat{T}, q::Int) aq = dat.iR[q]; dat.αq = dat.Binv * (dat.prob.A[:,aq]); end
	function check_optimal_dJ{T}(dat::Dat{T}) return all( dj->(dj >= dat.zero), dat.dJ ); end
	function check_feasible_β{T}(dat::Dat{T}) return all( β->(β >= dat.zero), dat.β ); end

	function pivot_iB_iR{T}(dat::Dat{T}, q::Int, p::Int) dat.iR[q], dat.iB[p] = dat.iB[p], dat.iR[q] end

	function price_full_dantzig{T}(dat::Dat{T})
		# [CTSM].p187
		min_i = 0; min_dj = dat.zero;
		for i = 1:length(dat.dJ)
			dj = dat.dJ[i]
			if (dj < dat.zero && dj <= min_dj)
				if (dj == min_dj)
					println("Warning: degeneracy.")
				end
				min_i = i; min_dj = dj;
			end
		end
		return min_i
	end

	function chuzro{T}(dat::Dat{T})
		min_i = 0; min_ratio = dat.zero;
		for i = 1:length(dat.αq)
			if (dat.αq[i] > dat.zero)
				ratio = dat.β[i] / dat.αq[i]
				if (min_i == 0 || ratio <= min_ratio)
					if (ratio == min_ratio)
						println("Warning: degeneracy.")
					end
					min_i = i; min_ratio = ratio;
				end
			end
		end
		return (min_i, min_ratio)
	end

	function succeed{T}(dat::Dat{T})
		sol.solved = true
		sol.status = :Optimal
		sol.z = dat.z
		sol.x = eval(parse( "zeros($(dat.prob.type_s), $(dat.n))" ))
		for i = 1:length(dat.iB)
			xi = dat.iB[i]
			if (xi <= dat.n) sol.x[xi] = dat.β[i]; end
		end
	end

	function fail(dat, sol::Lp.Solution, status::Symbol) sol.solved = false; sol.status = status; end

	function solve_dat{T}(dat::Dat{T}, sol::Lp.Solution{T})
		# [CTSM].p33,p30, but with naive basis reinversion instead of update.
		dcd = sol.dcd
		Dcd.@it(dcd)
		it = 0
		#Step 0
		set_basis_logical(dat)
		Dcd.@set(dcd, "iB", dat.iB)
		#Initializations
		comp_cBT(dat); comp_B_R(dat); comp_Binv(dat);
		init_β(dat); init_z(dat);
		Dcd.@set(dcd, "β0", dat.β)
		#todo phaseI
		if (check_feasible_β(dat) == false) println("Warning: phaseI."); fail(dat, sol, :Infeasible); return; end

		while(dat.maxit == 0 || it < dat.maxit)
			Dcd.@it(dcd)
			#Step 1
			Dcd.@set(dcd, "B", dat.B); Dcd.@set(dcd, "Binv", dat.Binv);
			Dcd.@set(dcd, "cBT", dat.cBT);
			comp_π(dat); Dcd.@set(dcd, "π", dat.π);
			#Step 2
			comp_dJ(dat); Dcd.@set(dcd, "dJ", dat.dJ);
			if check_optimal_dJ(dat) succeed(dat, sol); return; end
			q = price_full_dantzig(dat); Dcd.@set(dcd, "q", q);
			#Step 3
			comp_αq(dat, q); Dcd.@set(dcd, "αq", dat.αq);
			#Step 4
			p,θ = chuzro(dat); Dcd.@set(dcd, "p", p); Dcd.@set(dcd, "θ", θ);
			if (p == 0) fail(dat, sol, :Unbounded); return; end
			Dcd.@set(dcd, "pivot", (q, p))
			#Step 5
			pivot_iB_iR(dat, q, p)
			#Updates
			update_β(dat, p, θ); Dcd.@set(dcd, "β", dat.β);
			update_z(dat, q, θ); Dcd.@set(dcd, "z", dat.z);
			comp_cBT(dat); comp_B_R(dat); comp_Binv(dat);
			it = it + 1
		end
		fail(dat, sol, :Maxit)
	end

end
=#

end #module Jad

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
using JuMP
using LightGraphs, SimpleWeightedGraphs

function single_commodity_flow!(model, graph :: SimpleWeightedGraph, k :: Int)
    n::Int = nv(graph)
    es::Vector{SimpleWeightedEdge} = collect(edges(graph))
    sort!(es, by=dst)
    sort!(es, by=src, alg=MergeSort)

    @variables(model, begin
        f[i=1:n, j=1:n; has_edge(graph, i, j)], (lower_bound=0)
        z[i=2:n], (binary=true, lower_bound = 0, upper_bound = 1)
    end)

    alpha::Int = 1

    for i in 2:n
        if has_edge(graph, alpha, i)
            fix(variable_by_name(model, "y[$i,$alpha]"), 0, force=true) # Constraint (14)
            fix(f[i, alpha], 0, force=true) # Constraint (15)
            set_upper_bound(f[alpha,i], k)
        end
        for j in 2:n
            if has_edge(graph, i, j)
                set_upper_bound(f[i,j], k-1) # Constraint (16)
                set_upper_bound(f[j,i], k-1) # Constraint (16)
            end
        end
        @constraint(model, sum(variable_by_name(model, "y[$j,$i]") for j in 1:n if has_edge(graph, j, i)) == z[i]) # Constraint (17)
        #@constraint(model, sum(variable_by_name(model, "y[$i,$j]") for j in 2:n if has_edge(graph, j, i)) <= length(neighbors(graph, i)) - 1)
    end

    for i in 2:n
        @constraint(model,
        (sum(f[dst(e),i] for e in es if src(e) == i) + sum(f[src(e),i] for e in es if dst(e) == i)
        - (sum(f[i,dst(e)] for e in es if src(e) == i && dst(e) != alpha) + sum(f[i,src(e)] for e in es if dst(e) == i && src(e) != alpha)) == z[i])) # Constraint (18)
    end
    @constraint(model, source_flow, sum(f[alpha,j] for j = 2:n if has_edge(graph, alpha, j)) == k) # Constraint (19)
    @constraint(model, one_real_root, sum(variable_by_name(model, "y[$alpha,$j]") for j in 2:n if has_edge(graph, alpha, j)) == 1) # Constraint (20)

    for i in 1:n
        for j in 1:n
            if has_edge(graph, i, j)
                y_ij = variable_by_name(model, "y[$i,$j]")
                if i != 1 && j != 1
                    @constraint(model,
                        f[i,j] <= (k-1) * y_ij # Constraint (21)
                    )
                elseif i == 1
                    @constraint(model,
                        f[i,j] <= k * y_ij # Constraint (22)
                    )
                end
            end
        end
    end

    @constraint(model, sum(z[i] for i = 2:n) == k) # Constraint (23)

end

max_stars_repo_name: Entze/jkMST

max_issues_repo_name: Entze/jkMST

max_forks_repo_name: Entze/jkMST


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct IEEEST <: PSS
        input_code::Int
        remote_bus_control::Int
        A1::Float64
        A2::Float64
        A3::Float64
        A4::Float64
        A5::Float64
        A6::Float64
        T1::Float64
        T2::Float64
        T3::Float64
        T4::Float64
        T5::Float64
        T6::Float64
        Ks::Float64
        Ls_lim::Tuple{Float64, Float64}
        Vcu::Float64
        Vcl::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        states_types::Vector{StateTypes}
        internal::InfrastructureSystemsInternal
    end

IEEE Stabilizing Model PSS.

# Arguments
- `input_code::Int`: Code input for stabilizar, validation range: `(1, 6)`, action if invalid: `error`
- `remote_bus_control::Int`: Remote  Bus number for control.
- `A1::Float64`: Filter coefficient, validation range: `(0, nothing)`, action if invalid: `warn`
- `A2::Float64`: Filter coefficient, validation range: `(0, nothing)`, action if invalid: `warn`
- `A3::Float64`: Filter coefficient, validation range: `(0, nothing)`
- `A4::Float64`: Filter coefficient, validation range: `(0, nothing)`, action if invalid: `warn`
- `A5::Float64`: Filter coefficient, validation range: `(0, nothing)`, action if invalid: `warn`
- `A6::Float64`: Filter coefficient, validation range: `(0, nothing)`, action if invalid: `warn`
- `T1::Float64`: Time constant, validation range: `(0, 10)`, action if invalid: `warn`
- `T2::Float64`: Time constant, validation range: `(0, 10)`, action if invalid: `warn`
- `T3::Float64`: Time constant, validation range: `(0, 10)`, action if invalid: `warn`
- `T4::Float64`: Time constant, validation range: `(0, 10)`, action if invalid: `warn`
- `T5::Float64`: Time constant, validation range: `(0, 10)`, action if invalid: `warn`
- `T6::Float64`: Time constant, validation range: `(eps(), 2.0)`, action if invalid: `error`
- `Ks::Float64`: Proportional gain, validation range: `(0, nothing)`, action if invalid: `warn`
- `Ls_lim::Tuple{Float64, Float64}`: PSS output limits for regulator output `(Ls_min, Ls_max)`
- `Vcu::Float64`: Cutoff limiter upper bound, validation range: `(0, 1.25)`, action if invalid: `warn`
- `Vcl::Float64`: Cutoff limiter lower bound, validation range: `(0, 1.0)`, action if invalid: `warn`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states are:
	x_p1: 1st filter integration,
	x_p2: 2nd filter integration,
	x_p3: 3rd filter integration,
	x_p4: 4rd filter integration,
	x_p5: T1/T2 lead-lag integrator,
	x_p6: T3/T4 lead-lag integrator,
	:x_p7 last integer,
- `n_states::Int`: IEEEST has 7 states
- `states_types::Vector{StateTypes}`: IEEEST has 7 differential states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct IEEEST <: PSS
    "Code input for stabilizar"
    input_code::Int
    "Remote  Bus number for control."
    remote_bus_control::Int
    "Filter coefficient"
    A1::Float64
    "Filter coefficient"
    A2::Float64
    "Filter coefficient"
    A3::Float64
    "Filter coefficient"
    A4::Float64
    "Filter coefficient"
    A5::Float64
    "Filter coefficient"
    A6::Float64
    "Time constant"
    T1::Float64
    "Time constant"
    T2::Float64
    "Time constant"
    T3::Float64
    "Time constant"
    T4::Float64
    "Time constant"
    T5::Float64
    "Time constant"
    T6::Float64
    "Proportional gain"
    Ks::Float64
    "PSS output limits for regulator output `(Ls_min, Ls_max)`"
    Ls_lim::Tuple{Float64, Float64}
    "Cutoff limiter upper bound"
    Vcu::Float64
    "Cutoff limiter lower bound"
    Vcl::Float64
    ext::Dict{String, Any}
    "The states are:
	x_p1: 1st filter integration,
	x_p2: 2nd filter integration,
	x_p3: 3rd filter integration,
	x_p4: 4rd filter integration,
	x_p5: T1/T2 lead-lag integrator,
	x_p6: T3/T4 lead-lag integrator,
	:x_p7 last integer,"
    states::Vector{Symbol}
    "IEEEST has 7 states"
    n_states::Int
    "IEEEST has 7 differential states"
    states_types::Vector{StateTypes}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function IEEEST(input_code, remote_bus_control, A1, A2, A3, A4, A5, A6, T1, T2, T3, T4, T5, T6, Ks, Ls_lim, Vcu, Vcl, ext=Dict{String, Any}(), )
    IEEEST(input_code, remote_bus_control, A1, A2, A3, A4, A5, A6, T1, T2, T3, T4, T5, T6, Ks, Ls_lim, Vcu, Vcl, ext, [:x_p1, :x_p2, :x_p3, :x_p4, :x_p5, :x_p6, :x_p7], 7, [StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Differential], InfrastructureSystemsInternal(), )
end

function IEEEST(; input_code, remote_bus_control, A1, A2, A3, A4, A5, A6, T1, T2, T3, T4, T5, T6, Ks, Ls_lim, Vcu, Vcl, ext=Dict{String, Any}(), states=[:x_p1, :x_p2, :x_p3, :x_p4, :x_p5, :x_p6, :x_p7], n_states=7, states_types=[StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Hybrid, StateTypes.Differential], internal=InfrastructureSystemsInternal(), )
    IEEEST(input_code, remote_bus_control, A1, A2, A3, A4, A5, A6, T1, T2, T3, T4, T5, T6, Ks, Ls_lim, Vcu, Vcl, ext, states, n_states, states_types, internal, )
end

# Constructor for demo purposes; non-functional.
function IEEEST(::Nothing)
    IEEEST(;
        input_code=1,
        remote_bus_control=0,
        A1=0,
        A2=0,
        A3=0,
        A4=0,
        A5=0,
        A6=0,
        T1=0,
        T2=0,
        T3=0,
        T4=0,
        T5=0,
        T6=0,
        Ks=0,
        Ls_lim=(0.0, 0.0),
        Vcu=0,
        Vcl=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`IEEEST`](@ref) `input_code`."""
get_input_code(value::IEEEST) = value.input_code
"""Get [`IEEEST`](@ref) `remote_bus_control`."""
get_remote_bus_control(value::IEEEST) = value.remote_bus_control
"""Get [`IEEEST`](@ref) `A1`."""
get_A1(value::IEEEST) = value.A1
"""Get [`IEEEST`](@ref) `A2`."""
get_A2(value::IEEEST) = value.A2
"""Get [`IEEEST`](@ref) `A3`."""
get_A3(value::IEEEST) = value.A3
"""Get [`IEEEST`](@ref) `A4`."""
get_A4(value::IEEEST) = value.A4
"""Get [`IEEEST`](@ref) `A5`."""
get_A5(value::IEEEST) = value.A5
"""Get [`IEEEST`](@ref) `A6`."""
get_A6(value::IEEEST) = value.A6
"""Get [`IEEEST`](@ref) `T1`."""
get_T1(value::IEEEST) = value.T1
"""Get [`IEEEST`](@ref) `T2`."""
get_T2(value::IEEEST) = value.T2
"""Get [`IEEEST`](@ref) `T3`."""
get_T3(value::IEEEST) = value.T3
"""Get [`IEEEST`](@ref) `T4`."""
get_T4(value::IEEEST) = value.T4
"""Get [`IEEEST`](@ref) `T5`."""
get_T5(value::IEEEST) = value.T5
"""Get [`IEEEST`](@ref) `T6`."""
get_T6(value::IEEEST) = value.T6
"""Get [`IEEEST`](@ref) `Ks`."""
get_Ks(value::IEEEST) = value.Ks
"""Get [`IEEEST`](@ref) `Ls_lim`."""
get_Ls_lim(value::IEEEST) = value.Ls_lim
"""Get [`IEEEST`](@ref) `Vcu`."""
get_Vcu(value::IEEEST) = value.Vcu
"""Get [`IEEEST`](@ref) `Vcl`."""
get_Vcl(value::IEEEST) = value.Vcl
"""Get [`IEEEST`](@ref) `ext`."""
get_ext(value::IEEEST) = value.ext
"""Get [`IEEEST`](@ref) `states`."""
get_states(value::IEEEST) = value.states
"""Get [`IEEEST`](@ref) `n_states`."""
get_n_states(value::IEEEST) = value.n_states
"""Get [`IEEEST`](@ref) `states_types`."""
get_states_types(value::IEEEST) = value.states_types
"""Get [`IEEEST`](@ref) `internal`."""
get_internal(value::IEEEST) = value.internal

"""Set [`IEEEST`](@ref) `input_code`."""
set_input_code!(value::IEEEST, val) = value.input_code = val
"""Set [`IEEEST`](@ref) `remote_bus_control`."""
set_remote_bus_control!(value::IEEEST, val) = value.remote_bus_control = val
"""Set [`IEEEST`](@ref) `A1`."""
set_A1!(value::IEEEST, val) = value.A1 = val
"""Set [`IEEEST`](@ref) `A2`."""
set_A2!(value::IEEEST, val) = value.A2 = val
"""Set [`IEEEST`](@ref) `A3`."""
set_A3!(value::IEEEST, val) = value.A3 = val
"""Set [`IEEEST`](@ref) `A4`."""
set_A4!(value::IEEEST, val) = value.A4 = val
"""Set [`IEEEST`](@ref) `A5`."""
set_A5!(value::IEEEST, val) = value.A5 = val
"""Set [`IEEEST`](@ref) `A6`."""
set_A6!(value::IEEEST, val) = value.A6 = val
"""Set [`IEEEST`](@ref) `T1`."""
set_T1!(value::IEEEST, val) = value.T1 = val
"""Set [`IEEEST`](@ref) `T2`."""
set_T2!(value::IEEEST, val) = value.T2 = val
"""Set [`IEEEST`](@ref) `T3`."""
set_T3!(value::IEEEST, val) = value.T3 = val
"""Set [`IEEEST`](@ref) `T4`."""
set_T4!(value::IEEEST, val) = value.T4 = val
"""Set [`IEEEST`](@ref) `T5`."""
set_T5!(value::IEEEST, val) = value.T5 = val
"""Set [`IEEEST`](@ref) `T6`."""
set_T6!(value::IEEEST, val) = value.T6 = val
"""Set [`IEEEST`](@ref) `Ks`."""
set_Ks!(value::IEEEST, val) = value.Ks = val
"""Set [`IEEEST`](@ref) `Ls_lim`."""
set_Ls_lim!(value::IEEEST, val) = value.Ls_lim = val
"""Set [`IEEEST`](@ref) `Vcu`."""
set_Vcu!(value::IEEEST, val) = value.Vcu = val
"""Set [`IEEEST`](@ref) `Vcl`."""
set_Vcl!(value::IEEEST, val) = value.Vcl = val
"""Set [`IEEEST`](@ref) `ext`."""
set_ext!(value::IEEEST, val) = value.ext = val
"""Set [`IEEEST`](@ref) `states_types`."""
set_states_types!(value::IEEEST, val) = value.states_types = val


text:
using JSON3, StructTypes, Test, Bonsai
using HTTP.Messages
using Sockets
using Sockets: InetAddr
import StructTypes: StructType

include("http.jl")
include("handlers.jl")
include("middleware.jl")
include("router.jl")
include("io.jl")
include("app.jl")
include("openapi.jl")

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
using OpenAIGymAPI

remote_base = "http://127.0.0.1:5000"
client = GymClient(remote_base)
print(client)

# Create environment
env_id = "CartPole-v0"
instance_id = env_create(client, env_id)
print(instance_id)

# List all environments
all_envs = env_list_all(client)
print(all_envs)

# Set up agent
action_space_info = env_action_space_info(client, instance_id)
print(action_space_info)
agent = action_space_info["n"] # perform same action every time

# Run experiment, with monitor
outdir = "/tmp/random-agent-results"
env_monitor_start(client, instance_id, outdir, force = true, resume = false)

episode_count = 100
max_steps = 200
for i in 1:episode_count
  ob = env_reset(client, instance_id)
  done = false
  j = 1
  while j <= 200 && !done
    action = env_action_space_sample(client, instance_id)
    results = env_step(client, instance_id, action, render = true)
    done = results["done"]
    j = j + 1
  end
end

# Dump result info to disk
env_monitor_close(client, instance_id)


text:
module MeshTools

using Meshes

end


text:
#=
    Contents: General utilities.
=#

"Returns true if the application is using secure cookies."
function using_secure_cookies()
    haskey(config, :securecookie)
end

function using_clientside_sessions()
    config[:session][:datastore] == :cookie
end

"Returns: nonce and newform, where newform is form with the nonce inserted as a hidden field."
function add_nonce_to_form(form::AbstractString)
    nonce   = base64encode(csrng(config[:session][:id_length]))
    replc   = "<input type='hidden' id='' name='nonce' value=$nonce></form>"
    newform = replace(form, "</form>", replc)
    nonce,  newform
end

"""
Inserts a hidden field into the form that contains a nonce.
Records the nonce in the server-side session object, form_id => nonce.
"""
function add_nonce_to_form(form_id::AbstractString, form::AbstractString, session_id::AbstractString)
    nonce, newform = add_nonce_to_form(form)
    session_set!(session_id, "forms", form_id, nonce)    # Store nonce in session
    newform
end

"""
Inserts a hidden field into the form that contains a nonce.
Records the nonce in the client-side session object, form_id => nonce.
"""
function add_nonce_to_form(form_id::AbstractString, form::AbstractString, session::Dict)
    nonce, newform = add_nonce_to_form(form)
    session["forms"][form_id] = nonce
    newform
end

"Cryptographically secure RNG"
function csrng(numbytes::Integer)
    entropy = MbedTLS.Entropy()
    rng     = MbedTLS.CtrDrbg()
    MbedTLS.seed!(rng, entropy)
    rand(rng, numbytes)
end

"""
Convert number to byte array.
The implementation is a modified version of that found in Stack Overflow question 3076680.
"""
function num_to_bytearray(x)
    io = IOBuffer()
    write(io, x)
    seekstart(io)
    readbytes(io)
end

#=
"""
Returns true is password adheres to the following formatting rules:
  R1) Length at least 8 characters
  R2) At least 5 unique characters
  R3) At least 1 upper case letter
  R4) At least 1 lower case letter
  R5) At least 1 number
  R6) At least 1 special character
"""
function password_is_permissible(password)
    # Determine the truth values of rules 1 to 6.
    uniq = unique(password)
    n    = length(uniq)
    r1   = length(password) >= 8
    r2   = n >= 5
    r3   = false
    r4   = false
    r5   = false
    r6   = false
    for i = 1:n
	c     = Int(uniq[i])    # Integer representation of Char
	if c >= 65 && c <= 90
	    r3 = true
	elseif c >= 97 && c <= 122
	    r4 = true
	elseif c >= 48 && c <= 57
	    r5 = true
	else
	    r6 = true
	end
    end

    # Calculate result
    result = false
    if r1 && r2 && r3 && r4 && r5 && r6
	result = true
    end
    result
end
=#

# EOF


text:
@testset "Test using MSA files" begin

    msa_types = (
        Matrix{Residue},
        NamedResidueMatrix{Array{Residue,2}},
        MultipleSequenceAlignment,
        AnnotatedMultipleSequenceAlignment
        )

    pf09645_sto = joinpath(pwd(), "data", "PF09645_full.stockholm")
    gaoetal2011 = joinpath(pwd(), "data", "Gaoetal2011.fasta")

    gaoetal_msas = [ read(gaoetal2011, FASTA, T) for T in msa_types ]
    pfam_msas    = [ read(pf09645_sto, Stockholm, T) for T in msa_types ]

    @testset "getindex" begin

        residues = permutedims(hcat(res"DAWAEE",
                                    res"DAWAEF",
                                    res"DAWAED",
                                    res"DAYCMD",
                                    res"DAYCMT",
                                    res"DAYCMT"), [2,1])

        for index in [ 2, (2,:), (:,2), (2,2), (2:4, 2:4), (2,[3,4,5]),
                      (:, [1,2,3,4,5,6] .< 4 ), ([1,2,3,4,5,6] .< 4, :),
                      (:,:), ([1,2,3,4,5,6] .< 4, [1,2,3,4,5,6] .< 4) ]

            for msa in gaoetal_msas
                @test msa[index...] == residues[index...]
            end
        end
    end

    @testset "setindex! and copy" begin

        for msa in gaoetal_msas
            copy_msa = copy(msa)
            deepcopy_msa = deepcopy(msa)

            for (index, value) in [ ((1), Residue('H')),
                                    ((:,1), res"HHHHHH") ]

                deepcopy_msa[index...] = value
                @test deepcopy_msa[index...] == value
                @test msa[index...] != value

                copy_msa[index...] = value
                @test copy_msa[index...] == value
                @test msa[index...] != value
            end

            for (index, value) in [ ( 1 , Residue('H')),
                                    ( 4 , Residue('X')),
                                    ( : , res"HHHHHH" ) ]

                seq = copy(getsequence(msa, 4))
                seq[1, index] = value
                @test seq[1, index] == value
                @test msa[4, index] != value # since seq is a copy
            end
        end
    end

    @testset "Size" begin

        for aln in gaoetal_msas
            @test size(aln) == (6,6)
            @test length(aln) == 36
            @test ncolumns(aln) == 6
            @test nsequences(aln) == 6
            @test ncolumns(getsequence(aln, 4)) == 6
        end

        for aln in pfam_msas
            @test size(aln) == (4,110)
            @test length(aln) == 440
            @test ncolumns(aln) == 110
            @test nsequences(aln) == 4
        end
    end

    @testset "AnnotatedAlignedSequence and AlignedSequence" begin

        @testset "Creation" begin
            seq_types = (   Matrix{Residue},
                            NamedResidueMatrix{Array{Residue,2}},
                            AlignedSequence,
                            AnnotatedAlignedSequence    )

            for i in 1:length(msa_types)
                M = msa_types[i]
                S = seq_types[i]
                msa = pfam_msas[i]

                if M != Matrix{Residue}
                    for id in [ "C3N734_SULIY/1-95","H2C869_9CREN/7-104",
                                "Y070_ATV/2-70",    "F112_SSV1/3-112"   ]
                        annseq = getsequence(msa, id)
                        @test msa[id,:] == vec(annseq) # Sequences are matrices
                        @test isa(annseq, S)
                    end
                end

                for seq in 1:4
                    annseq = getsequence(msa, seq)
                    @test msa[seq,:] == vec(annseq) # Sequences are matrices
                    @test isa(annseq, S)
                end
            end
        end

        @testset "Annotations" begin
            msa = read(pf09645_sto, Stockholm)

            @test   getannotcolumn(msa, "SS_cons") ==
                    getannotcolumn(getsequence(msa, 4), "SS_cons")
            @test   getannotresidue(msa, "F112_SSV1/3-112", "SS") ==
                    getannotresidue(getsequence(msa, 4), "F112_SSV1/3-112", "SS")
            @test   getannotfile(msa) == getannotfile(getsequence(msa, 4))
            @test   getannotsequence(msa,"F112_SSV1/3-112","DR") ==
                    getannotsequence(getsequence(msa, 4), "F112_SSV1/3-112", "DR")
        end
    end

    @testset "Print" begin
        pfam = read(pf09645_sto, Stockholm)
        gao  = read(gaoetal2011, FASTA)

        @test stringsequence(gao, 4) == "DAYCMD"
        @test stringsequence(pfam, 1) == stringsequence(pfam, "C3N734_SULIY/1-95")

        for T in (Stockholm, FASTA, Raw)
            buffer = IOBuffer()
            print(buffer, pfam, T)
            @test parse(String(take!(buffer)), T) == pfam

            print(buffer, gao, T)
            @test parse(String(take!(buffer)), T) == gao
        end
    end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to help screen.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    _help!(pargerd::Pager)

Open a new pager with the help.

"""
function _help!(pagerd::Pager)
    # Unpack values.
    term = pagerd.term
    buf  = pagerd.buf

    if get(term.out_stream, :color, true)::Bool
        _b  = string(crayon"bold")
        _d  = string(Crayon(reset = true))
        _cb = string(crayon"cyan bold")
        _c  = string(crayon"cyan")
        _g  = string(crayon"dark_gray")
        _r  = string(crayon"red bold")
        _y  = string(crayon"yellow bold")
    else
        _b = ""
        _d = ""
        _c = ""
        _g = ""
        _r = ""
        _y = ""
    end

    # Get the current key bindings.

    # General
    kb_help         = _getkb(:help)
    kb_quit         = _getkb(:quit)
    kb_quit_eot     = _getkb(:quit_eot)
    kb_toggle_ruler = _getkb(:toggle_ruler)

    # Movement
    kb_up        = _getkb(:up)
    kb_down      = _getkb(:down)
    kb_left      = _getkb(:left)
    kb_right     = _getkb(:right)
    kb_fastup    = _getkb(:fastup)
    kb_fastdown  = _getkb(:fastdown)
    kb_fastleft  = _getkb(:fastleft)
    kb_fastright = _getkb(:fastright)
    kb_pageup    = _getkb(:pageup)
    kb_pagedown  = _getkb(:pagedown)
    kb_hpageup   = _getkb(:halfpageup)
    kb_hpagedown = _getkb(:halfpagedown)
    kb_bol       = _getkb(:bol)
    kb_eol       = _getkb(:eol)
    kb_home      = _getkb(:home)
    kb_end       = _getkb(:end)

    # Searching
    kb_search         = _getkb(:search)
    kb_next_match     = _getkb(:next_match)
    kb_previous_match = _getkb(:previous_match)
    kb_quit_search    = _getkb(:quit_search)

    # Freezing data
    kb_change_freeze     = _getkb(:change_freeze)
    kb_change_title_rows = _getkb(:change_title_rows)

    help_str =
"""
  $(_cb)TerminalPager.jl $(PKG_VERSION)$(_d)

  The pager can execute some type of actions as shown in the following. The key
  bindings of each actions can be changed using the function $(_c)set_keybinding$(_d).

  Some actions are only available if a feature is enabled. The set of enabled
  features can be selected using the keyword $(_c)features$(_d) when calling the
  pager.

$(_b)                                     General$(_d)
$(_y)  :help$(_d)
    Show this screen.
$(_c)    Keybindings: $(kb_help)$(_d)
$(_g)    This action requires the feature :help.
$(_y)  :quit$(_d)
    Quit the pager.
$(_c)    Keybindings: $(kb_quit)$(_d)
$(_y)  :quit_eot$(_d)
    This is an special quit action design for the $(_c)END OF TRANSMISSION (^D)$(_d)
    keycode. If we are in a search operation, then it quits the search. If not,
    then it quits the pager.
$(_c)    Keybindings: $(kb_quit_eot)$(_d)
$(_y)  :toggle_ruler$(_d)
    Toggle the vertical ruler.
$(_c)    Keybindings: $(kb_toggle_ruler)$(_d)

$(_b)                                    Movement$(_d)
$(_y)  :up$(_d)
    Move the display one line up.
$(_c)    Keybindings: $(kb_up)$(_d)
$(_y)  :down$(_d)
    Move the display one line down.
$(_c)    Keybindings: $(kb_down)$(_d)
$(_y)  :left$(_d)
    Move the display one column to the left.
$(_c)    Keybindings: $(kb_left)$(_d)
$(_y)  :right$(_d)
    Move the display one column to the right.
$(_c)    Keybindings: $(kb_right)$(_d)
$(_y)  :fastup$(_d)
    Move the display five lines up.
$(_c)    Keybindings: $(kb_fastup)$(_d)
$(_y)  :fastdown$(_d)
    Move the display five lines down.
$(_c)    Keybindings: $(kb_fastdown)$(_d)
$(_y)  :fastleft$(_d)
    Move the display ten columns to the left.
$(_c)    Keybindings: $(kb_fastleft)$(_d)
$(_y)  :fastright$(_d)
    Move the display ten columns to the right.
$(_c)    Keybindings: $(kb_fastright)$(_d)
$(_y)  :pageup$(_d)
    Move the display one page up (a page has the same size of the view).
$(_c)    Keybindings: $(kb_pageup)$(_d)
$(_y)  :pagedown$(_d)
    Move the display one page down (a page has the same size of the view).
$(_c)    Keybindings: $(kb_pagedown)$(_d)
$(_y)  :halfpageup$(_d)
    Move the display half page up (a page has the same size of the view).
$(_c)    Keybindings: $(kb_hpageup)$(_d)
$(_y)  :halfpagedown$(_d)
    Move the display half page down (a page has the same size of the view).
$(_c)    Keybindings: $(kb_hpagedown)$(_d)
$(_y)  :bol$(_d)
    Move the display to the first column.
$(_c)    Keybindings: $(kb_bol)$(_d)
$(_y)  :eol$(_d)
    Move the display to show the last column.
$(_c)    Keybindings: $(kb_eol)$(_d)
$(_y)  :home$(_d)
    Move the display to the first line.
$(_c)    Keybindings: $(kb_home)$(_d)
$(_y)  :end$(_d)
    Move the display to show the last line.
$(_c)    Keybindings: $(kb_end)$(_d)

$(_b)                                   Searching$(_d)
$(_y)  :search$(_d)
    Request a regex in the command line and highlight all the matches.
$(_c)    Keybindings: $(kb_search)$(_d)
$(_y)  :next_match$(_d)
    Go to the next match of the search.
$(_c)    Keybindings: $(kb_next_match)$(_d)
$(_y)  :previous_match$(_d)
    Go to the previous match of the search.
$(_c)    Keybindings: $(kb_previous_match)$(_d)
$(_y)  :quit_search$(_d)
    Quit searching, removing all the highlights (only during search mode).
$(_c)    Keybindings: $(kb_quit_search)$(_d)

$(_b)                                 Freezing data$(_d)
$(_g)  These action requires the feature :change_freeze.
$(_y)  :change_freeze$(_d)
    Two values will be requested in the command line. The first is the number of
    columns and the second is the number of rows that will be frozen. If the
    value is equal or lower than 0, then no row or column will be frozen.
$(_c)    Keybindings: $(kb_change_freeze)$(_d)
$(_y)  :change_title_rows$(_d)
    Define the number of rows within the frozen rows that will be considered as
    titles. In this case, these rows will not scroll horizontally.
$(_c)    Keybindings: $(kb_change_title_rows)$(_d)
"""

    _pager!(pagerd.term, help_str; hashelp = false)

    return nothing
end

################################################################################
#                              Private functions
################################################################################

function _getkb(action::Symbol)
    kb = [_kbtostr(k) for (k, v) in _keybindings if v == action]
    num_kb = length(kb)

    str = ""

    @inbounds for i = 1:num_kb
        str *= kb[i]
        i != num_kb && (str *= ", ")
    end

    return str
end

function _kbtostr(kb::Tuple{Union{Symbol, String}, Bool, Bool, Bool})
    if kb[1] == " "
        str = "space"
    else
        str = string(kb[1])
    end

    kb[2] && (str = "ALT " * str)
    kb[3] && (str = "CTRL " * str)
    kb[4] && (str = "SHIFT " * str)

    return str
end


text:
module HigherOrderNetworks
using MatrixNetworks
using SparseArrays
using StatsBase
using LinearAlgebra

include("Utilities.jl")
include("Hypergraphs.jl")
export Hypergraphs, add_node!, remove_node!, add_edge!, remove_edge!

include("Hyperkron.jl")
export hyperkron_graph, kron_params

include("HypergraphConversions.jl")
export to_adjacency_matrix, hypergraph_to_matrixnetwork, hyper_to_bipartite, bipartite_to_hyper, is_bipartite

include("RandomCliqueCovers.jl")
export rcc

include("RandomHypergraphs.jl")
export stub_matching, pairwise_reshuffle!, MCMC, MCMC_s, MCMC_v

include("Cliques.jl")
export kcliques

include("ManageData.jl")
export load_data

end # module


text:
# this runs in your environment and assumes that all your scripts are up to date etc.
# do not run this if you have not previously run `serve()` and checked everything was fine.

Pkg.activate()

using JuDoc, Pkg, Logging, NodeJS, Literate

##################################################################
## PART 1: generate all scripts and notebooks and push everything
# on the gh-pages branch (yes we're abusing the system)
##################################################################

Logging.disable_logging(Logging.LogLevel(1500))

ifiles = joinpath.("scripts", readdir("scripts"))

nbpath = joinpath("generated", "notebooks")
scpath = joinpath("generated", "scripts")

isdir(nbpath) || mkpath(nbpath)
isdir(scpath) || mkpath(scpath)

ACTIVATE = """
    # Before running this, please make sure to activate and instantiate the environment
    # corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Manifest.toml)
    # so that you get an environment which matches the one used to generate the tutorials:
    #
    # ```julia
    # cd("MLJTutorials") # cd to folder with the *.toml
    # using Pkg; Pkg.activate("."); Pkg.instantiate()
    # ```

    """

preproc(s) = ACTIVATE * s

# Remove lines that end with `# hide`
postproc(s) = replace(s, r"(^|\n).*?#(\s)*?(?i)hide"=>s"\1")

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

get_fn(fp) = splitext(splitdir(fp)[2])[1]

for file in ifiles
   # Generate annotated notebooks
   Literate.notebook(file, nbpath,
                      preprocess=preproc, postprocess=postproc,
                      execute=false, documenter=false)

   # Generate annotated scripts
   Literate.script(file, scpath,
                    postprocess=preproc ∘ postproc,
                    keep_comments=true, documenter=false)

   # Generate stripped scripts (no comments)
   Literate.script(file, scpath, name=get_fn(file)*"-raw",
                    postprocess=preproc ∘ postproc,
                    keep_comments=false, documenter=false)
end

JS_GHP = """
    var ghpages = require('gh-pages');
    ghpages.publish('generated/', function(err) {});
    """

run(`$(nodejs_cmd()) -e $JS_GHP`)

##################################################################
## PART 2
# optimize (compress and pre-render) the website and publish it
# to the master branch this time.
##################################################################

publish()

max_stars_repo_name: vdayanand/MLJTutorials

max_issues_repo_name: vdayanand/MLJTutorials

max_forks_repo_name: vdayanand/MLJTutorials


text:
@testset "debuginfo" begin

DEBUG_METADATA_VERSION()

Context() do ctx
      mod = parse(LLVM.Module,  """
          define void @foo() !dbg !5 {
          top:
            ret void, !dbg !7
          }

          define void @bar() {
          top:
            ret void
          }

          !llvm.module.flags = !{!0, !1}
          !llvm.dbg.cu = !{!2}

          !0 = !{i32 2, !"Dwarf Version", i32 4}
          !1 = !{i32 1, !"Debug Info Version", i32 3}
          !2 = distinct !DICompileUnit(language: DW_LANG_C89, file: !3, producer: "julia", isOptimized: true, runtimeVersion: 0, emissionKind: FullDebug, enums: !4)
          !3 = !DIFile(filename: "REPL[1]", directory: ".")
          !4 = !{}
          !5 = distinct !DISubprogram(name: "foo", linkageName: "foo", scope: null, file: !3, line: 1, type: !6, isLocal: false, isDefinition: true, scopeLine: 1, isOptimized: true, unit: !2)
          !6 = !DISubroutineType(types: !4)
          !7 = !DILocation(line: 1, scope: !5)"""; ctx)

    foo = functions(mod)["foo"]

    if LLVM.version() >= v"8.0"
        sp = LLVM.get_subprogram(foo)
        @test sp !== nothing

        bar = functions(mod)["bar"]
        @test LLVM.get_subprogram(bar) === nothing
        LLVM.set_subprogram!(bar, sp)
        @test LLVM.get_subprogram(bar) == sp
      end

    bb = entry(foo)
    inst = first(instructions(bb))

    @test !isempty(metadata(inst))
    strip_debuginfo!(mod)
    @test isempty(metadata(inst))
end

end


text:
module prompt

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
module GreekAnalytical

__precompile__(true)

using Distributions

export blackScholesValue, blackScholesDelta, blackScholesGamma, blackScholesVega, blackScholesTheta
export blackScholesVanna, blackScholesVolga, blackScholesSpeed, blackScholesCharm, blackScholesColor, blackScholesZomma
export blackScholesRho, blackScholesRhoFutures, blackScholesTouchProb, blackScholesDriftlessTheta

std_normal = Normal(0, 1)

function blackScholesValue(isCall::Bool, spot::T, strike::T, τ::T, r::T, b::T, σ::T) where {T}
    sign = isCall ? 1 : -1

    d1 = (log(spot / strike) + (b + 0.5 * σ^2) * τ) / (σ * sqrt(τ))
    d2 = d1 - σ * sqrt(τ)

    sign * (spot * exp((b-r) * τ) * cdf(std_normal, sign * d1) - strike * exp(-r * τ) * cdf(std_normal, sign * d2))
end

function blackScholesDelta(isCall::Bool, spot::T, strike::T, τ::T, r::T, b::T, σ::T) where {T}
    d1 = (log(spot / strike) + (b + 0.5 * σ^2) * τ) / (σ * sqrt(τ))
    if isCall
        exp((b-r) * τ) * cdf(std_normal, d1)
    else
        exp((b-r) * τ) * (cdf(std_normal, d1) - 1.0)
    end
end

function blackScholesGamma(isCall::Bool, spot::T, strike::T, τ::T, r::T, b::T, σ::T) where {T}
    d1 = (log(spot / strike) + (b + 0.5 * σ^2) * τ) / (σ * sqrt(τ))
    pdf(std_normal, d1) * exp((b-r) * τ) / spot / σ / sqrt(τ)
end

function blackScholesVega(isCall::Bool, spot::T, strike::T, τ::T, r::T, b::T, σ::T) where {T}
    d1 = (log(spot / strike) + (b + 0.5 * σ^2) * τ) / (σ * sqrt(τ))
    spot * exp((b-r) * τ) * pdf(std_normal, d1) * sqrt(τ)
end

function blackScholesTheta(isCall::Bool, spot::T, strike::T, τ::T, r::T, b::T, σ::T) where {T}
    sign = isCall ? 1 : -1

    d1 = (log(spot / strike) + (b + 0.5 * σ^2) * τ) / (σ * sqrt(τ))
    d2 = d1 - σ * sqrt(τ)

    term1 = -spot * exp((b-r) * τ) * pdf(std_normal, d1) * σ / 2 / sqrt(τ)
    term2 = -sign * (b-r) * spot * exp((b-r) * τ) * cdf(std_normal, sign * d1)
    term3 = -sign * r * strike * exp(-r * τ) * cdf(std_normal, sign * d2)

    term1 + term2 + term3
end

function blackScholesVanna(isCall::Bool, spot::T, strike::T, τ::T, r::T, b::T, σ::T) where {T}
    d1 = (log(spot / strike) + (b + 0.5 * σ^2) * τ) / (σ * sqrt(τ))
    d2 = d1 - σ * sqrt(τ)

    -d2 * pdf(std_normal, d1) / σ * exp((b-r) * τ)
end

function blackScholesCharm(isCall::Bool, spot::T, strike::T, τ::T, r::T, b::T, σ::T) where {T}
    sign = isCall ? 1 : -1

    d1 = (log(spot / strike) + (b + 0.5 * σ^2) * τ) / (σ * sqrt(τ))
    d2 = d1 - σ * sqrt(τ)

    term1 = pdf(std_normal, d1) * (b / σ / sqrt(τ) - d2 / 2 / τ)
    term2 = sign * (b - r) * cdf(std_normal, sign * d1)

    -exp((b-r) * τ) * (term1 + term2)
end

function blackScholesVolga(isCall::Bool, spot::T, strike::T, τ::T, r::T, b::T, σ::T) where {T}
    d1 = (log(spot / strike) + (b + 0.5 * σ^2) * τ) / (σ * sqrt(τ))
    d2 = d1 - σ * sqrt(τ)

    blackScholesVega(isCall, spot, strike, τ, r, b, σ) * d1 * d2 / σ
end

function blackScholesSpeed(isCall::Bool, spot::T, strike::T, τ::T, r::T, b::T, σ::T) where {T}
    d1 = (log(spot / strike) + (b + 0.5 * σ^2) * τ) / (σ * sqrt(τ))

    -1.0 * blackScholesGamma(isCall, spot, strike, τ, r, b, σ) / spot * (1 + d1 / σ / sqrt(τ))
end

function blackScholesZomma(isCall::Bool, spot::T, strike::T, τ::T, r::T, b::T, σ::T) where {T}
    d1 = (log(spot / strike) + (b + 0.5 * σ^2) * τ) / (σ * sqrt(τ))
    d2 = d1 - σ * sqrt(τ)

    blackScholesGamma(isCall, spot, strike, τ, r, b, σ) * (d1 * d2 - 1.0) / σ
end

function blackScholesColor(isCall::Bool, spot::T, strike::T, τ::T, r::T, b::T, σ::T) where {T}
    d1 = (log(spot / strike) + (b + 0.5 * σ^2) * τ) / (σ * sqrt(τ))
    d2 = d1 - σ * sqrt(τ)

    blackScholesGamma(isCall, spot, strike, τ, r, b, σ) * (r - b + b * d1 / σ / sqrt(τ) + (1 - d1 * d2) / 2 / τ)
end

function blackScholesDriftlessTheta(isCall::Bool, spot::T, strike::T, τ::T, r::T, b::T, σ::T) where {T}
    d1 = (log(spot / strike) + (b + 0.5 * σ^2) * τ) / (σ * sqrt(τ))
    println(pdf(std_normal, d1))
    -spot * pdf(std_normal, d1) * σ / 2 / sqrt(τ)
end

function blackScholesRho(isCall::Bool, spot::T, strike::T, τ::T, r::T, b::T, σ::T) where {T}
    sign = isCall ? 1 : -1

    d1 = (log(spot / strike) + (b + 0.5 * σ^2) * τ) / (σ * sqrt(τ))
    d2 = d1 - σ * sqrt(τ)

    sign * τ * strike * exp(-r * τ) * cdf(std_normal, sign * d2)
end

function blackScholesRhoFutures(isCall::Bool, spot::T, strike::T, τ::T, r::T, b::T, σ::T) where{T}
    -1.0 * τ * blackScholesValue(isCall, spot, strike, τ, r, b, σ)
end

function blackScholesTouchProb(isCall::Bool, spot::T, strike::T, τ::T, r::T, b::T, σ::T) where {T}
    if isCall
        if strike <= spot
            return 1.0
        end
    else
        if strike >= spot
            return 1.0
        end
    end

    sign = isCall ? -1 : 1

    μ = (b - σ^2 / 2) / (σ * σ)
    λ = sqrt(μ^2 + 2 * r / (σ * σ))
    z = log(strike / spot) / σ / sqrt(τ) + λ * σ * sqrt(τ)

    (strike / spot) ^ (μ + λ) * cdf(std_normal, sign * z) + (strike / spot) ^ (μ - λ) * cdf(std_normal, sign * (z - 2 * λ * σ * sqrt(τ)))
end

end

max_stars_repo_name: Stochastic-Adventure/BlackScholesGreeks

max_issues_repo_name: Stochastic-Adventure/BlackScholesGreeks

max_forks_repo_name: Stochastic-Adventure/BlackScholesGreeks


text:
function starting_points!(
  pt0::Point{T},
  fd::Abstract_QM_FloatData{T},
  id::QM_IntData,
  itd::IterData{T},
  spd::StartingPointData{T},
) where {T <: Real}
  itd.Qx = mul!(itd.Qx, Symmetric(fd.Q, fd.uplo), pt0.x)
  fd.uplo == :U ? mul!(itd.ATy, fd.A, pt0.y) : mul!(itd.ATy, fd.A', pt0.y)
  spd.dual_val .= itd.Qx .- itd.ATy .+ fd.c
  pt0.s_l = spd.dual_val[id.ilow]
  pt0.s_u = -spd.dual_val[id.iupp]

  # check distance to bounds δ for x, s_l and s_u
  itd.x_m_lvar .= @views pt0.x[id.ilow] .- fd.lvar[id.ilow]
  itd.uvar_m_x .= @views fd.uvar[id.iupp] .- pt0.x[id.iupp]
  if id.nlow == 0
    δx_l1, δs_l1 = zero(T), zero(T)
  else
    δx_l1 = max(-T(1.5) * minimum(itd.x_m_lvar), T(1.e-2))
    δs_l1 = max(-T(1.5) * minimum(pt0.s_l), T(1.e-4))
  end
  if id.nupp == 0
    δx_u1, δs_u1 = zero(T), zero(T)
  else
    δx_u1 = max(-T(1.5) * minimum(itd.uvar_m_x), T(1.e-2))
    δs_u1 = max(-T(1.5) * minimum(pt0.s_u), T(1.e-4))
  end

  # correct components that to not respect the bounds
  itd.x_m_lvar .+= δx_l1
  itd.uvar_m_x .+= δx_u1
  spd.s0_l1 .= pt0.s_l .+ δs_l1
  spd.s0_u1 .= pt0.s_u .+ δs_u1
  xs_l1, xs_u1 = dot(spd.s0_l1, itd.x_m_lvar), dot(spd.s0_u1, itd.uvar_m_x)
  if id.nlow == 0
    δx_l2, δs_l2 = zero(T), zero(T)
  else
    δx_l2 = δx_l1 + xs_l1 / sum(spd.s0_l1) / 2
    δs_l2 = δs_l1 + xs_l1 / sum(itd.x_m_lvar) / 2
  end
  if id.nupp == 0
    δx_u2, δs_u2 = zero(T), zero(T)
  else
    δx_u2 = δx_u1 + xs_u1 / sum(spd.s0_u1) / 2
    δs_u2 = δs_u1 + xs_u1 / sum(itd.uvar_m_x) / 2
  end
  δx = max(δx_l2, δx_u2)
  δs = max(δs_l2, δs_u2)
  pt0.x[id.ilow] .+= δx
  pt0.x[id.iupp] .-= δx
  pt0.s_l .= pt0.s_l .+ δs
  pt0.s_u .= pt0.s_u .+ δs

  # deal with the compensation phaenomenon in x if irng != []
  update_rngbounds!(pt0.x, id.irng, fd.lvar, fd.uvar, T(1e-4))

  # verify bounds
  @assert all(pt0.x .> fd.lvar) && all(pt0.x .< fd.uvar)
  @assert all(pt0.s_l .> zero(T)) && all(pt0.s_u .> zero(T))

  # update itd
  update_IterData!(itd, pt0, fd, id, false)
end

function update_rngbounds!(x::Vector{T}, irng, lvar, uvar, ϵ) where {T <: Real}
  @inbounds @simd for i in irng
    if lvar[i] >= x[i]
      x[i] = lvar[i] + ϵ
    end
    if x[i] >= uvar[i]
      x[i] = uvar[i] - ϵ
    end
    if (lvar[i] < x[i] < uvar[i]) == false
      x[i] = (lvar[i] + uvar[i]) / 2
    end
  end
end


text:
# "in this version, set data and smoothness terms using arrays.
#  grid neighborhood structure is assumed.#
#
using GCoptimization

width = 10
height = 5
num_pixels = width * height
num_labels = 7
result = zeros(num_pixels)    # "stores result of optimization"

immutable ForDataFn
    data::Vector{Cint}
    numLab::Cint
end

function dataFn(p::GCOSiteID, l::GCOLabelID, data::Ptr{Void})::GCOEnergyTermType
    myData = unsafe_pointer_to_objref(data)
    numLab = myData.numLab
    return myData.data[p*numLab+l+1]
end

smoothFn(p1::GCOSiteID, p2::GCOSiteID, l1::GCOLabelID, l2::GCOLabelID)::GCOEnergyTermType = (l1-l2)*(l1-l2) <= 4 ? (l1-l2)*(l1-l2) : 4

# "first set up the array for data costs"
data = Vector{Cint}(num_pixels * num_labels)
for i = 0:num_pixels-1, l = 0:num_labels-1
    if i < 25
        if l == 0
            data[i*num_labels+l+1] = 0
        else
            data[i*num_labels+l+1] = 10
        end
    else
        if l == 5
            data[i*num_labels+l+1] = 0
        else
            data[i*num_labels+l+1] = 10
        end
    end
end

gco = GCoptimizationGridGraph(width, height, num_labels)

# "set up the needed data to pass to function for the data costs"
toFn = ForDataFn(data, num_labels)
setDataCost(gco, dataFn, pointer_from_objref(toFn))

# "smoothness comes from function pointer"
setSmoothCost(gco, smoothFn)

println("Before optimization energy is ", compute_energy(gco))
expansion(gco, 2)    # "run expansion for 2 iterations."
println("After optimization energy is ", compute_energy(gco))

for i = 0:num_pixels-1
    result[i+1] = whatLabel(gco, i)
end

showall(result)


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 05-2019
# Chebyshev Polynomials: Gauss Lobatto grid
# See Boyd F.8
#--------------------------------------------------------------------

export collocation, derivative, integral

function naturalcollocation(S::ChebyshevGL{Tag, N, T}, i::Int)::T where {Tag, N, T}
    @assert i <= N
    return  cospi(T(i-1)/T(N-1))
end

function naturalderivative(S::ChebyshevGL{Tag, N, T}, i::Int, j::Int)::T where {Tag, N, T}
    @assert i <= N && j <= N
	if i==j==1
        return  (2T(N-1)^2 + 1)/T(6)
	elseif i==j==N
        return -(2T(N-1)^2 + 1)/T(6)
	elseif i==j
        return - naturalcollocation(S,j)/(2(1-naturalcollocation(S,j)^2))
	else
		ci = (i == 1 || i == N) ?  2 : 1
		cj = (j == 1 || j == N) ?  2 : 1
		s  = (i + j) % 2 != 0   ? -1 : 1
        return (T(ci)/T(cj))*(T(s)/(naturalcollocation(S, i) - naturalcollocation(S,j)))
	end
end

function naturalintegral(S::ChebyshevGL{Tag, N, T}, i::Int)::T where {Tag, N, T}
    @assert i <= N
    W = T(0)
	for j in 1:N
        w = (j == 1 ? 1 : (j-1)%2 == 0 ? T(2)/T((1-(j-1)^2)) : 0)
        l = (i == 1 || i == N ? (T(1)/T(N-1))*cospi(T(i-1)*T(j-1)/T(N-1)) : (T(2)/T(N-1))*cospi(T(i-1)*T(j-1)/T(N-1)))
        W = W + w*l
	end
	return W
end

function collocation(space::S, i::Int)::T where {S <: Cardinal{Tag, N, T}} where {Tag, N, T}
    return naturalcollocation(space,i)*T(space.max - space.min)/2 + T(space.max + space.min)/2
end

function derivative(space::S, i::Int, j::Int)::T where {S <: Cardinal{Tag, N, T}} where {Tag, N, T}
    return naturalderivative(space,i,j)*(2/T(space.max - space.min))
end

function integral(space::S, i::Int)::T where {S <: Cardinal{Tag, N, T}} where {Tag, N, T}
    return naturalintegral(space,i)*T(space.max - space.min)/T(2)
end


text:
mutable struct GFRegenerateState
    prev_trace::DynamicDSLTrace
    trace::DynamicDSLTrace
    selection::Selection
    weight::Float64
    visitor::AddressVisitor
    params::Dict{Symbol,Any}
end

function GFRegenerateState(gen_fn, args, prev_trace,
                           selection, params)
    visitor = AddressVisitor()
    GFRegenerateState(prev_trace, DynamicDSLTrace(gen_fn, args), selection,
        0., visitor, params)
end

function traceat(state::GFRegenerateState, dist::Distribution{T},
                 args, key) where {T}
    local prev_retval::T
    local retval::T

    # check that key was not already visited, and mark it as visited
    visit!(state.visitor, key)

    # check for previous choice at this key
    has_previous = has_choice(state.prev_trace, key)
    if has_previous
        prev_choice = get_choice(state.prev_trace, key)
        prev_retval = prev_choice.retval
        prev_score = prev_choice.score
    end

    # check whether the key was selected
    in_selection = key in state.selection

    # get return value
    if has_previous && in_selection
        retval = random(dist, args...)
    elseif has_previous
        retval = prev_retval
    else
        retval = random(dist, args...)
    end

    # compute logpdf
    score = logpdf(dist, retval, args...)

    # update weight
    if has_previous && !in_selection
        state.weight += score - prev_score
    end

    # add to the trace
    add_choice!(state.trace, key, retval, score)

    retval
end

function traceat(state::GFRegenerateState, gen_fn::GenerativeFunction{T,U},
                 args, key) where {T,U}
    local prev_retval::T
    local trace::U
    local retval::T

    # check that key was not already visited, and mark it as visited
    visit!(state.visitor, key)

    # check whether the key was selected
    subselection = state.selection[key]

    # get subtrace
    has_previous = has_call(state.prev_trace, key)
    if has_previous
        prev_call = get_call(state.prev_trace, key)
        prev_subtrace = prev_call.subtrace
        get_gen_fn(prev_subtrace) === gen_fn || gen_fn_changed_error(key)
        (subtrace, weight, _) = regenerate(
            prev_subtrace, args, map((_) -> UnknownChange(), args), subselection)
    else
        (subtrace, weight) = generate(gen_fn, args, EmptyChoiceMap())
    end

    # update weight
    state.weight += weight

    # add to the trace
    add_call!(state.trace, key, subtrace)

    # get return value
    retval = get_retval(subtrace)

    retval
end

function splice(state::GFRegenerateState, gen_fn::DynamicDSLFunction,
                args::Tuple)
    prev_params = state.params
    state.params = gen_fn.params
    retval = exec(gen_fn, state, args)
    state.params = prev_params
    retval
end

function regenerate_delete_recurse(prev_trie::Trie{Any,ChoiceOrCallRecord},
                             visited::EmptySelection)
    noise = 0.
    for (key, choice_or_call) in get_leaf_nodes(prev_trie)
        if !choice_or_call.is_choice
            noise += choice_or_call.noise
        end
    end
    for (key, subtrie) in get_internal_nodes(prev_trie)
        noise += regenerate_delete_recurse(subtrie, EmptySelection())
    end
    noise
end

function regenerate_delete_recurse(prev_trie::Trie{Any,ChoiceOrCallRecord},
                             visited::DynamicSelection)
    noise = 0.
    for (key, choice_or_call) in get_leaf_nodes(prev_trie)
        if !(key in visited) && !choice_or_call.is_choice
            noise += choice_or_call.noise
        end
    end
    for (key, subtrie) in get_internal_nodes(prev_trie)
        subvisited = visited[key]
        noise += regenerate_delete_recurse(subtrie, subvisited)
    end
    noise
end

function regenerate(trace::DynamicDSLTrace, args::Tuple, argdiffs::Tuple,
                    selection::Selection)
    gen_fn = trace.gen_fn
    state = GFRegenerateState(gen_fn, args, trace, selection, gen_fn.params)
    retval = exec(gen_fn, state, args)
    set_retval!(state.trace, retval)
    visited = state.visitor.visited
    state.weight -= regenerate_delete_recurse(trace.trie, visited)
    (state.trace, state.weight, UnknownChange())
end


text:
using SpecialFunctions

function abramm1(z::Complex{Float64}, ifexpon=1)
    local w:: Complex{Float64} = complex(0.0, 0.0)
    ccall((:abramm1_,:abram), Nothing,
          (Ptr{Complex{Float64}},Ptr{Complex{Float64}},Ptr{Int64}),
          Ref(Complex{Float64}(z)), Ref{Complex{Float64}}(w), Ref{Int64}(ifexpon))
    w
end

function abram0(z::Complex{Float64}, ifexpon=1)
    local w:: Complex{Float64} = complex(0.0, 0.0)
    ccall((:abram0_,:abram), Nothing,
          (Ptr{Complex{Float64}},Ptr{Complex{Float64}},Ptr{Int64}),
          Ref(Complex{Float64}(z)), Ref{Complex{Float64}}(w), Ref{Int64}(ifexpon))
    w
end

function abram1(z::Complex{Float64}, ifexpon=1)
    local w:: Complex{Float64} = complex(0.0, 0.0)
    ccall((:abram1_,:abram), Nothing,
          (Ptr{Complex{Float64}},Ptr{Complex{Float64}},Ptr{Int64}),
          Ref(Complex{Float64}(z)), Ref{Complex{Float64}}(w), Ref{Int64}(ifexpon))
    w
end

function abram2(z::Complex{Float64}, ifexpon=1)
    local w:: Complex{Float64} = complex(0.0, 0.0)
    ccall((:abram2_,:abram), Nothing,
          (Ptr{Complex{Float64}},Ptr{Complex{Float64}},Ptr{Int64}),
          Ref(Complex{Float64}(z)), Ref{Complex{Float64}}(w), Ref{Int64}(ifexpon))
    w
end

function test()

ntest = 1000000
z = rand(ntest) + 1im*rand(ntest)
z = z * 10

println("n = ", ntest)

println("abramm1:")
@time w = abramm1.(z)
@time w = abramm1.(z)
@time w = abramm1.(z)
println("isfinite(w) = ", all(isfinite.(w)))

println("abram0:")
@time w = abram0.(z)
@time w = abram0.(z)
@time w = abram0.(z)
println("isfinite(w) = ", all(isfinite.(w)))

println("abram1:")
@time w = abram1.(z)
@time w = abram1.(z)
@time w = abram1.(z)
println("isfinite(w) = ", all(isfinite.(w)))

println("abram2:")
@time w = abram2.(z)
@time w = abram2.(z)
@time w = abram2.(z)
println("isfinite(w) = ", all(isfinite.(w)))

println("exp:")
@time w = exp.(z)
@time w = exp.(z)
@time w = exp.(z)
println("isfinite(w) = ", all(isfinite.(w)))

println("done.")

end

test()
test()
test()

max_stars_repo_name: zgimbutas/abramowitz

max_issues_repo_name: zgimbutas/abramowitz

max_forks_repo_name: zgimbutas/abramowitz


text:
### Abstract parameter NStates

abstract ParameterNState{S<:ValueSupport, F<:VariateForm} <: VariableNState{F}

typealias MarkovChain ParameterNState

typealias DiscreteParameterNState{F<:VariateForm} ParameterNState{Discrete, F}
typealias ContinuousParameterNState{F<:VariateForm} ParameterNState{Continuous, F}

typealias UnivariateParameterNState{S<:ValueSupport} ParameterNState{S, Univariate}
typealias MultivariateParameterNState{S<:ValueSupport} ParameterNState{S, Multivariate}
typealias MatrixvariateParameterNState{S<:ValueSupport} ParameterNState{S, Matrixvariate}

diagnostics(nstate::ParameterNState) =
  Dict(zip(nstate.diagnostickeys, Any[nstate.diagnosticvalues[i, :][:] for i = 1:size(nstate.diagnosticvalues, 1)]))

Base.copy!{S<:ValueSupport, F<:VariateForm}(nstate::ParameterNState{S, F}, state::ParameterState{S, F}, i::Integer) =
  nstate.copy(nstate, state, i)