text:
import StaticArrays: SVector

@testset "Impulse operations" begin

    # test adding analytic impulses
        ω_vec = 0.0:0.1:1.01
        t_vec = 0.0:0.1:1.01

        gauss = GaussianImpulse(maximum(ω_vec))
        dirac = FreqDiracImpulse(ω_vec[1])
        impulse = gauss +  dirac * 2.0

        @test_throws(DomainError,[1]*dirac)

        @test all(impulse.in_freq.(ω_vec) .== gauss.in_freq.(ω_vec) + 2.0 .* dirac.in_freq.(ω_vec))
        @test all(impulse.in_time.(t_vec) .== gauss.in_time.(t_vec) + 2.0 .* dirac.in_time.(t_vec))

    # test a shifted discrete impulse
        t_vec = 0.0:0.1:10.01
        ω_vec = t_to_ω(t_vec)

        gauss = GaussianImpulse(maximum(ω_vec) / 2)

        discrete_impulse = continuous_to_discrete_impulse(gauss, t_vec, ω_vec; t_shift = 4.2, ω_shift = maximum(ω_vec) / 2);

        time_response = frequency_to_time(discrete_impulse.in_freq, ω_vec, t_vec);
        freq_response = time_to_frequency(real.(discrete_impulse.in_time), t_vec, ω_vec);

        @test maximum(abs.(time_response - discrete_impulse.in_time)) < 1e-2
        @test maximum(abs.(freq_response - discrete_impulse.in_freq)) < 1e-2

    # compare with the analytic shifts
        gauss_shift = GaussianImpulse(maximum(ω_vec) / 2; t_shift = 4.2, ω_shift = maximum(ω_vec) / 2)
        discrete_impulse2 = continuous_to_discrete_impulse(gauss_shift, t_vec, ω_vec);

        @test maximum(abs.(discrete_impulse2.in_time - discrete_impulse.in_time)) < 1e-10
        @test maximum(abs.(discrete_impulse2.in_freq - discrete_impulse.in_freq)) < 1e-10
end

@testset "Time Result" begin
    sound_p = Acoustic(.1, 0.1 + 0.0im,2)
    particles = [Particle(sound_p,Sphere([10.5,0.0], .5))]

    sound_sim = Acoustic(1., 1. + 0.0im,2)
    source = plane_source(sound_sim, [0.0,0.0], [1.0,0.0], 1.)
    sim = FrequencySimulation(particles, source)

    ω_vec = 0.0:0.01:5.01
    @test LinRange(ω_vec) ≈ t_to_ω(ω_to_t(ω_vec)) # only exact for length(ω_vec) = even number

    # invertability of dft
        x_vec = [ [0.0,0.0], [3.0,0.0]]
        ω_vec = 0.0:0.1:1.01
        t_vec = ω_to_t(ω_vec)
        simres = run(sim, x_vec, ω_vec)
        # choose an impulse which does nothing and so can be inverted
        discrete_impulse = DiscreteTimeDiracImpulse(0.0, t_vec, ω_vec)
        timres = frequency_to_time(simres; method=:dft, discrete_impulse = discrete_impulse)
        simres2 = time_to_frequency(timres; method=:dft, discrete_impulse = discrete_impulse)

        @test norm(field(simres) - field(simres2)) / norm(field(simres)) < 1e-14
        timres2 = frequency_to_time(simres2; method=:dft, impulse = TimeDiracImpulse(0.0));
        @test norm(field(timres2) - field(timres))/norm(field(timres)) < 1e-14

    # test impulse consistency by seeing if we get the same result when going from freq to time, and when going from time to freq.
        impulse = GaussianImpulse(maximum(ω_vec))
        timres2 = frequency_to_time(simres; method=:dft, impulse = impulse)
        simres2 = time_to_frequency(timres; method=:dft, impulse = impulse)
        # invert without an impulse
        timres3 = frequency_to_time(simres2; method=:dft, discrete_impulse = discrete_impulse)
        @test norm(field(timres2) - field(timres3))/norm(field(timres2)) < 1e-14

    # Compare dft and trapezoidal integration
        ω_vec = 0.0:0.1:100.01 # need a high frequency to match a delta impluse function!
        simres = run(sim, x_vec, ω_vec)
        timres1 = frequency_to_time(simres; method=:trapezoidal, impulse = TimeDiracImpulse(0.0))
        timres2 = frequency_to_time(simres; method=:dft, discrete_impulse = DiscreteTimeDiracImpulse(0.0, ω_to_t(simres.ω)))
        @test norm(field(timres1) - field(timres2))/norm(field(timres1)) < 0.02

        ω_vec = 0.0:0.0001:2.01 # need a high sampling to match a delta impluse function!
        t_vec = 0.:0.5:20.

        impulse = GaussianImpulse(maximum(ω_vec))

        simres = run(sim, x_vec, ω_vec)
        timres1 = frequency_to_time(simres; t_vec = t_vec, method=:trapezoidal, impulse = impulse)
        timres2 = frequency_to_time(simres; t_vec = t_vec, method=:dft, impulse = impulse)
        @test norm(field(timres1) - field(timres2))/norm(field(timres1)) < 2e-5
        # plot(timres1.t, [field(timres1)[1,:]-field(timres2)[1,:]])
end

@testset "Time shift signal" begin

    sound_sim = Acoustic(2., 5.2 + 0.0im,2)
    # suppose we measure in time the signal
    t = 0.1:0.01:30.
    t0 = 5.
    impulse_func(t) = exp( -10000.0 .* (t - t0)^4)
    measured_response = impulse_func.(t)
    plot(t, measured_response)

    ω = t_to_ω(t)
    measured_f = time_to_frequency(measured_response, t, ω)

    # at the point
    source_x = [-3.0,2.0]
    source_direction = [1.0,1.0]
    x_measure =  source_x + 4.0*source_direction

    # then can we reverse-engineer what plane wave source (of the form below) would create this measured_response
    travel_time = norm(x_measure - source_x)/real(sound_sim.c)
    impulse_time = t .- travel_time;
    impulse_response = impulse_func.(t);

    amp0 = 1.0
    source = plane_source(sound_sim, source_x, source_direction, amp0)
    sim = FrequencySimulation(source)
    simres = run(sim, [x_measure, x_measure + source_direction], ω)
    short_time = 1.:0.00005:1.5
    timres = frequency_to_time(simres; t_vec = short_time)
    i1 = findmax(abs.(field(timres)[1,:]))[2]
    i2 = findmax(abs.(field(timres)[2,:]))[2]
    # test wave peak propagation speed
    @test (short_time[i2] - short_time[i1] - norm(source_direction/sound_sim.c)) < 1e-4*norm(source_direction/sound_sim.c)

    amps = measured_f./field(simres)[1,:]

    # predict the source impulse
    predict_impulse = frequency_to_time(amps, ω, impulse_time[300:600])
    # correct mistake from ω=0
    predict_impulse = predict_impulse .- mean(predict_impulse[1:2])

    @test norm(predict_impulse - impulse_response[300:600]) < 1e-12*norm(impulse_response[300:600])
end


text:
using ApproxFun, LinearAlgebra, SpecialFunctions, FastTransforms, Test
    import ApproxFun: testbandedoperator, testraggedbelowoperator, InterlaceOperator, testspace,
                        testbandedbelowoperator, testbandedblockbandedoperator, testfunctional, factor
## This includes extra tests that are too time consuming for Travis

include("runtests.jl")

@time include("FractionalTest.jl")

@testset "Full Operator" begin
    @time for M in (Multiplication(Fun(CosSpace(),[1.]),CosSpace()),
                    Multiplication(Fun(CosSpace(),[1.]),SinSpace()),
                    Multiplication(Fun(SinSpace(),[1.]),SinSpace()),
                    Multiplication(Fun(SinSpace(),[1.]),CosSpace()),
                    Derivative(SinSpace()),Derivative(CosSpace()))
          testbandedoperator(M)
    end

    S=Chebyshev()
    @time for io in (
            [Dirichlet(S);Derivative(Chebyshev());lneumann(S)],
            [Dirichlet(S);Derivative(Chebyshev())+Fun(cos);lneumann(S)],
            [Dirichlet(S);Derivative(Chebyshev())],
            [Dirichlet(S);Derivative(Chebyshev())+Fun(cos)],
            [Derivative(Chebyshev());Dirichlet(S)],
            [Derivative(Chebyshev())+Fun(cos);Dirichlet(S)])
        testraggedbelowoperator(io)
    end

    @testset "nullspace" begin
        d=ChebyshevInterval()
        D=Derivative(d)
        A=D^2-I
        @time κ=nullspace(A)
        @test length(κ) == 2
    end

    @testset "Newton iteration bug" begin
        S=Chebyshev(0..7)
        ω=2π
        N = u->Any[Fun(u(0.)-0.1);Fun(u(ω)-u(0.));Fun(u'(ω)-u'(0.));u''+u+u^3]

        u=0.1Fun(cos,S)
        @test u(0.0) ≈ 0.1

        D=Derivative(S)

        Z=ApproxFun.ZeroOperator(ApproxFun.ConstantSpace())

        testfunctional(Evaluation(S,0))
        testfunctional(Evaluation(S,ω)-Evaluation(S,0))
        testfunctional(Evaluation(S,ω,1)-Evaluation(S,0,1))

        A=[Z                      Evaluation(S,0);
                             u'(ω)    Evaluation(S,ω)-Evaluation(S,0);
                             u''(ω)   Evaluation(S,ω,1)-Evaluation(S,0,1);
                              0         D^2+I+3u^2]

        ApproxFun.backend_testinfoperator(A)

        DS=WeightedJacobi(0.1+1,0.2+1)
        D=Derivative(DS)[2:end,:]
        @time ApproxFun.testbandedoperator(D)
    end
end

@testset "Full ODE" begin
    @testset "Null space" begin
        d=ChebyshevInterval()
        D=Derivative(d)
        A=D^2-I
        @time κ=nullspace(A)
        @test length(κ) == 2

        c=[κ(0.);κ'(0.)]\[exp(0.),exp(0.)]
        u=(κ*c)[1]

        @test u(1.0) ≈ ℯ
    end

    d=(-50..5.)
    x=Fun(identity,d)
    D=Derivative(d)
    @time u=nullspace(D^2-x)
    c=[u(leftendpoint(d)); u(rightendpoint(d))]\[airyai.(endpoints(d))...]
    @test norm((u*c)[1]-Fun(airyai,d))<10000eps()

    ## constant forcing
    d = Interval(0.,50.)
    D = Derivative(d)
    t = Fun(identity,d)

    F = D^2 +.5D + I

    A= [ 0    ldirichlet(d);
         0    lneumann(d);
         0    rdirichlet(d);
        -1    F; ]

    @time u,x=A\[1.,0.,2.,0.]

    @test norm(F*x-u)<1000eps()

    @testset "Bessel" begin
        @time for ν in (1.,0.5,2.,3.5)
            println("        ν = $ν")
            S=JacobiWeight(-ν,0.,Chebyshev(0..1))
            D=Derivative(S)
            x=Fun(identity,domain(S))
            L=(x^2)*D^2+x*D+(x^2-ν^2);
            u=\([rdirichlet(S);rneumann(S);L],[bessely(ν,1.),.5*(bessely(ν-1.,1.)-bessely(ν+1.,1.)),0];
                        tolerance=1E-10)
            @test ≈(u(.1),bessely(ν,.1);atol=eps(1000000.)*max(abs(u(.1)),1))
            u=Fun(x->bessely(ν,x),S)
            @test ≈(u(.1),bessely(ν,.1);atol=eps(10000.)*max(abs(u(.1)),1))
            u=Fun(x->besselj(ν,x),S)
            @test ≈(u(.1),besselj(ν,.1);atol=eps(10000.)*max(abs(u(.1)),1))
        end

        @time for ν in (1.,0.5,0.123,3.5)
            println("        ν = $ν")
            S=JacobiWeight(ν,0.,Chebyshev(0..1))
            D=Derivative(S)
            x=Fun(identity,domain(S))
            L=(x^2)*D^2+x*D+(x^2-ν^2);

            u=\([rdirichlet(S);rneumann(S);L],[besselj(ν,1.),.5*(besselj(ν-1.,1.)-besselj(ν+1.,1.)),0];
                        tolerance=1E-10)
            @test ≈(u(.1),besselj(ν,.1);atol=eps(1000000.)*max(abs(u(.1)),1))
            u=Fun(x->besselj(ν,x),S)
            @test ≈(u(.1),besselj(ν,.1);atol=eps(10000.)*max(abs(u(.1)),1))
        end
    end
end

@testset "Full Jacobi" begin
    sp = Jacobi(.5,2.124)
    f = Fun(exp,sp)
    sp2 = Jacobi(1.5,2.124)
    M = Multiplication(f,sp2)
    @time testbandedoperator(M)

    ## Legendre conversions
    testspace(Ultraspherical(1); haslineintegral=false)
    testspace(Ultraspherical(2); haslineintegral=false)
    # minpoints is a tempory fix a bug
    @time testspace(Ultraspherical(1//2); haslineintegral=false, minpoints=2)
    @test norm(Fun(exp,Ultraspherical(1//2))-Fun(exp,Jacobi(0,0))) < 100eps()

    C=Conversion(Jacobi(0,0),Chebyshev())
    @time testbandedbelowoperator(C)
    @test norm(C*Fun(exp,Jacobi(0,0))  - Fun(exp)) < 100eps()

    C=Conversion(Ultraspherical(1//2),Chebyshev())
    @time testbandedbelowoperator(C)
    @test norm(C*Fun(exp,Ultraspherical(1//2))  - Fun(exp)) < 100eps()

    C=Conversion(Chebyshev(),Ultraspherical(1//2))
    @time testbandedbelowoperator(C)
    @test norm(C*Fun(exp)-Fun(exp,Legendre())) < 100eps()

    C=Conversion(Chebyshev(),Jacobi(0,0))
    @time testbandedbelowoperator(C)
    @test norm(C*Fun(exp)  - Fun(exp,Jacobi(0,0))) < 100eps()

    C=Conversion(Chebyshev(),Jacobi(1,1))
    @time testbandedbelowoperator(C)
    @test norm(C*Fun(exp) - Fun(exp,Jacobi(1,1))) < 100eps()

    C=Conversion(Ultraspherical(1//2),Ultraspherical(1))
    @time testbandedbelowoperator(C)

    λ1 = ApproxFun.order(domainspace(C))
    λ2 = ApproxFun.order(rangespace(C))

    # test against version that doesn't use lgamma
    Cex = Float64[(if j ≥ k && iseven(k-j)
            gamma(λ2)*(k-1+λ2)/(gamma(λ1)*gamma(λ1-λ2))*
                (gamma((j-k)/2+λ1-λ2)/gamma((j-k)/2+1))*
                (gamma((k+j-2)/2+λ1)/gamma((k+j-2)/2+λ2+1))
        else
            0.0
        end) for k=1:20,j=1:20]

    @test norm(Cex - C[1:20,1:20]) < 100eps()

    @test norm(C*Fun(exp,Ultraspherical(1//2))-Fun(exp,Ultraspherical(1))) < 100eps()

    C=Conversion(Jacobi(0,0),Ultraspherical(1))
    testbandedbelowoperator(C)
    @test norm(C*Fun(exp,Jacobi(0,0))-Fun(exp,Ultraspherical(1))) < 100eps()

    C=Conversion(Ultraspherical(1),Jacobi(0,0))
    testbandedbelowoperator(C)
    @test norm(C*Fun(exp,Ultraspherical(1))-Fun(exp,Jacobi(0,0))) < 100eps()
end

@testset "Full multivariate" begin
    ## ProductFun
    u0   = ProductFun((x,y)->cos(x)+sin(y) +exp(-50x.^2-40(y-0.1)^2)+.5exp(-30(x+0.5)^2-40(y+0.2)^2))

    @test values(u0)-values(u0|>LowRankFun)|>norm < 1000eps()
    @test chebyshevtransform(values(u0))-coefficients(u0)|>norm < 100eps()

    ##TODO: need to do adaptive to get better accuracy
    @test sin(u0)(.1,.2)-sin(u0(.1,.2))|>abs < 10e-4

    F = LowRankFun((x,y)->hankelh1(0,10abs(y-x)),Chebyshev(1.0..2.0),Chebyshev(Segment(1.0im,2.0im)))

    @test F(1.5,1.5im) ≈ hankelh1(0,10abs(1.5im-1.5))

    ## Periodic
    f=LowRankFun((x,y)->cos(x)*sin(y),PeriodicSegment(),PeriodicSegment())
    @test f(.1,.2) ≈ cos(.1)*sin(.2)

    f=LowRankFun((x,y)->cos(cos(x)+sin(y)),PeriodicSegment(),PeriodicSegment())
    @test f(.1,.2) ≈ cos(cos(.1)+sin(.2))
    @test norm(Float64[cos(cos(x)+sin(y)) for x=ApproxFun.vecpoints(f,1),y=ApproxFun.vecpoints(f,2)]-values(f))<10000eps()

    f=ProductFun((x,y)->cos(cos(x)+sin(y)),PeriodicSegment()^2)
    @test f(.1,.2) ≈ cos(cos(.1)+sin(.2))
    x,y=points(f)
    @test norm(Float64[cos(cos(x[k,j])+sin(y[k,j])) for k=1:size(f,1),j=1:size(f,2)]-values(f))<10000eps()

    d=PeriodicSegment()^2
    f=ProductFun((x,y)->exp(-10(sin(x/2)^2+sin(y/2)^2)),d)
    @test (transpose(f)-f|>coefficients|>norm)< 1000eps()

    @testset "Functional*Fun" begin
        d=ChebyshevInterval()
        B=ldirichlet(d)
        f=ProductFun((x,y)->cos(cos(x)*sin(y)),d^2)

        @test norm(B*f-Fun(y->cos(cos(-1)*sin(y)),d))<20000eps()
        @test norm(f*B-Fun(x->cos(cos(x)*sin(-1)),d))<20000eps()
    end

    @testset "matrix" begin
        f=Fun((x,y)->[exp(x*cos(y));cos(x*sin(y));2],ChebyshevInterval()^2)
        @test f(0.1,0.2) ≈ [exp(0.1*cos(0.2));cos(0.1*sin(0.2));2]

        f=Fun((x,y)->[exp(x*cos(y)) cos(x*sin(y)); 2 1],ChebyshevInterval()^2)
        @test f(0.1,0.2) ≈ [exp(0.1*cos(0.2)) cos(0.1*sin(0.2));2 1]
    end

    @testset "Cauchy fun" begin
        f = Fun((x,y)->1/(2π*(x^2+y^2+1)^(3/2)),Line()^2)
        @test f(0.1,0.2) ≈ 1/(2π*(0.1^2+0.2^2+1)^(3/2))

        #TODO: improve tolerance
        f = LowRankFun((x,y)->1/(2π*(x^2+y^2+1)^(3/2)),JacobiWeight(2.,2.,Line())^2)
        @test ≈(f(0.1,0.2),1/(2π*(0.1^2+0.2^2+1)^(3/2));atol=1E-4)
    end

    @testset "2d derivative (issue #346)" begin
        d = Chebyshev()^2
        f = Fun((x,y) -> sin(x) * cos(y), d)
        C=Conversion(Chebyshev()⊗Chebyshev(),Ultraspherical(1)⊗Ultraspherical(1))
        @test (C*f)(0.1,0.2) ≈ f(0.1,0.2)
        Dx = Derivative(d, [1,0])
        f = Fun((x,y) -> sin(x) * cos(y), d)
        fx = Fun((x,y) -> cos(x) * cos(y), d)
        @test (Dx*f)(0.2,0.3) ≈ fx(0.2,0.3)
        Dy = Derivative(d, [0,1])
        fy = Fun((x,y) -> -sin(x) * sin(y), d)
        @test (Dy*f)(0.2,0.3) ≈ fy(0.2,0.3)
        L=Dx+Dy
        testbandedblockbandedoperator(L)

        @test (L*f)(0.2,0.3) ≈ (fx(0.2,0.3)+fy(0.2,0.3))

        B=ldirichlet(factor(d,1))⊗ldirichlet(factor(d,2))
        @test Number(B*f) ≈ f(-1.,-1.)

        B=Evaluation(factor(d,1),0.1)⊗ldirichlet(factor(d,2))
        @test Number(B*f) ≈ f(0.1,-1.)

        B=Evaluation(factor(d,1),0.1)⊗Evaluation(factor(d,2),0.3)
        @test Number(B*f) ≈ f(0.1,0.3)

        B=Evaluation(d,(0.1,0.3))
        @test Number(B*f) ≈ f(0.1,0.3)
    end
end

include("FullPDETest.jl")
println("Speed tests")
include("SpeedTest.jl")
include("SpeedODETest.jl")
include("SpeedPDETest.jl")

include("ReadmeTest.jl")

ext: jl


text:
using Documenter, CoupledFields

makedocs(
    modules = [CoupledFields]
)

deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo   = "github.com/Mattriks/CoupledFields.jl.git",
    julia  = "release"
)


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

struct NetworkSecurityGroupsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Creates or updates a network security group in the specified resource group.
Param: resourceGroupName::String (required)
Param: networkSecurityGroupName::String (required)
Param: parameters::NetworkSecurityGroup (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkSecurityGroup
"""
function _swaggerinternal_networkSecurityGroupsCreateOrUpdate(_api::NetworkSecurityGroupsApi, resourceGroupName::String, networkSecurityGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PUT", NetworkSecurityGroup, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkSecurityGroups/{networkSecurityGroupName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkSecurityGroupName", networkSecurityGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkSecurityGroupsCreateOrUpdate(_api::NetworkSecurityGroupsApi, resourceGroupName::String, networkSecurityGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkSecurityGroupsCreateOrUpdate(_api, resourceGroupName, networkSecurityGroupName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkSecurityGroupsCreateOrUpdate(_api::NetworkSecurityGroupsApi, response_stream::Channel, resourceGroupName::String, networkSecurityGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkSecurityGroupsCreateOrUpdate(_api, resourceGroupName, networkSecurityGroupName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Deletes the specified network security group.
Param: resourceGroupName::String (required)
Param: networkSecurityGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: Nothing
"""
function _swaggerinternal_networkSecurityGroupsDelete(_api::NetworkSecurityGroupsApi, resourceGroupName::String, networkSecurityGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkSecurityGroups/{networkSecurityGroupName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkSecurityGroupName", networkSecurityGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkSecurityGroupsDelete(_api::NetworkSecurityGroupsApi, resourceGroupName::String, networkSecurityGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkSecurityGroupsDelete(_api, resourceGroupName, networkSecurityGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkSecurityGroupsDelete(_api::NetworkSecurityGroupsApi, response_stream::Channel, resourceGroupName::String, networkSecurityGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkSecurityGroupsDelete(_api, resourceGroupName, networkSecurityGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets the specified network security group.
Param: resourceGroupName::String (required)
Param: networkSecurityGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: expand::String
Return: NetworkSecurityGroup
"""
function _swaggerinternal_networkSecurityGroupsGet(_api::NetworkSecurityGroupsApi, resourceGroupName::String, networkSecurityGroupName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkSecurityGroup, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkSecurityGroups/{networkSecurityGroupName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkSecurityGroupName", networkSecurityGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_param(_ctx.query, "\$expand", expand)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkSecurityGroupsGet(_api::NetworkSecurityGroupsApi, resourceGroupName::String, networkSecurityGroupName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_networkSecurityGroupsGet(_api, resourceGroupName, networkSecurityGroupName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkSecurityGroupsGet(_api::NetworkSecurityGroupsApi, response_stream::Channel, resourceGroupName::String, networkSecurityGroupName::String, api_version::String, subscriptionId::String; expand=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_networkSecurityGroupsGet(_api, resourceGroupName, networkSecurityGroupName, api_version, subscriptionId; expand=expand, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all network security groups in a resource group.
Param: resourceGroupName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkSecurityGroupListResult
"""
function _swaggerinternal_networkSecurityGroupsList(_api::NetworkSecurityGroupsApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkSecurityGroupListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkSecurityGroups", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkSecurityGroupsList(_api::NetworkSecurityGroupsApi, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkSecurityGroupsList(_api, resourceGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkSecurityGroupsList(_api::NetworkSecurityGroupsApi, response_stream::Channel, resourceGroupName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkSecurityGroupsList(_api, resourceGroupName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets all network security groups in a subscription.
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkSecurityGroupListResult
"""
function _swaggerinternal_networkSecurityGroupsListAll(_api::NetworkSecurityGroupsApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", NetworkSecurityGroupListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Network/networkSecurityGroups", ["azure_auth"])
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkSecurityGroupsListAll(_api::NetworkSecurityGroupsApi, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkSecurityGroupsListAll(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkSecurityGroupsListAll(_api::NetworkSecurityGroupsApi, response_stream::Channel, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkSecurityGroupsListAll(_api, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Updates a network security group tags.
Param: resourceGroupName::String (required)
Param: networkSecurityGroupName::String (required)
Param: parameters::TagsObject (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: NetworkSecurityGroup
"""
function _swaggerinternal_networkSecurityGroupsUpdateTags(_api::NetworkSecurityGroupsApi, resourceGroupName::String, networkSecurityGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "PATCH", NetworkSecurityGroup, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/networkSecurityGroups/{networkSecurityGroupName}", ["azure_auth"], parameters)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "networkSecurityGroupName", networkSecurityGroupName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function networkSecurityGroupsUpdateTags(_api::NetworkSecurityGroupsApi, resourceGroupName::String, networkSecurityGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkSecurityGroupsUpdateTags(_api, resourceGroupName, networkSecurityGroupName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function networkSecurityGroupsUpdateTags(_api::NetworkSecurityGroupsApi, response_stream::Channel, resourceGroupName::String, networkSecurityGroupName::String, parameters, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_networkSecurityGroupsUpdateTags(_api, resourceGroupName, networkSecurityGroupName, parameters, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export networkSecurityGroupsCreateOrUpdate, networkSecurityGroupsDelete, networkSecurityGroupsGet, networkSecurityGroupsList, networkSecurityGroupsListAll, networkSecurityGroupsUpdateTags


text:
using eFEM, Test

# mesh input from gmsh
mesh1 = Mesh("test-data/step.msh")
mesh2 = FluidMesh("test-data/step.msh")


text:
function optimize(f, ranges, ps::ParticleSwarm, iter)
    fields = getproperty.(ranges, :field)
    state = PSO.initialize(ranges, ps)
    for _ in 1:iter
        PSO.retrieve!(state, ps)
        measurements = [
            f(zip(fields, params)) for params in zip(state.parameters...)
        ]
        PSO.pbest!(state, measurements, ps)
        PSO.gbest!(state)
        PSO.move!(ps.rng, state, ps.w, ps.c1, ps.c2)
    end
    min, particle = findmin(state.pbest)
    return min, state.pbest_X[particle, :]
end

function optimize(f, ranges, ps::AdaptiveParticleSwarm, iter)
    fields = getproperty.(ranges, :field)
    state = PSO.initialize(ranges, ps)
    phase = nothing
    c1, c2 = ps.c1, ps.c2
    for _ in 1:iter
        PSO.retrieve!(state, ps)
        measurements = [
            f(zip(fields, params)) for params in zip(state.parameters...)
        ]
        PSO.pbest!(state, measurements, ps)
        PSO.gbest!(state)
        factor = PSO._evolutionary_factor(state.X, argmin(state.pbest))
        phase = PSO._evolutionary_phase(factor, phase)
        w, c1, c2 = PSO._adapt_parameters(ps.rng, c1, c2, factor, phase)
        PSO.move!(ps.rng, state, w, c1, c2)
    end
    min, particle = findmin(state.pbest)
    return min, state.pbest_X[particle, :]
end

for PS in (ParticleSwarm, AdaptiveParticleSwarm)
    @testset "Optimize Ackley with $PS" begin
        function ackley(x; a=20, b=0.2, c=2π)
            d = length(x)
            return -a * exp(-b * sqrt(sum(x.^2) / d)) -
                exp(sum(cos.(c .* x)) / d) + a + ℯ
        end

        @testset "Numeric Ackley" begin
            r1 = range(Float64, :r1; lower=-2.0, upper=2.0)
            r2 = range(Float64, :r2; lower=-2.0, upper=2.0)
            ps = PS(n_particles=10, rng=StableRNG(1234))
            min, params = optimize([r1, r2], ps, 1_500_000) do pairs
                x = [param for (field, param) in pairs]
                ackley(x)
            end
            @test all(isapprox.(params, 0; atol=3)) # analytical solution
        end

        @testset "Integer Ackley" begin
            r1 = range(Int, :r1; lower=-20, upper=20)
            r2 = range(Int, :r2; lower=-20, upper=20)
            ps = PS(n_particles=10, rng=StableRNG(1234))
            min, params = optimize([r1, r2], ps, 250) do pairs
                x = [param for (field, param) in pairs]
                ackley(x)
            end
            @test all(round.(params) .== 0) # analytical solution
        end

        @testset "Nominal Ackley" begin
            vals = shuffle(StableRNG(1234), -2:2)
            r1 = range(Int, :r1; values=vals)
            r2 = range(Int, :r2; values=vals)
            ps = PS(n_particles=10, rng=StableRNG(1234))
            min, params = optimize([r1, r2], ps, 100) do pairs
                x = [param for (field, param) in pairs]
                ackley(x)
            end
            true_params = [1., 0., 0., 0., 0., 1., 0., 0., 0., 0.] # analytical solution
            @test all(isapprox.(params, true_params; atol=1e-3))
        end

        @testset "Mixed Ackley" begin
            vals = shuffle(StableRNG(1234), -2:2)
            r1 = range(Float64, :r1; lower=-2.0, upper=2.0)
            r2 = range(Int, :r2; lower=-20, upper=20)
            r3 = range(Int, :r3; values=vals)
            ps = PS(n_particles=20, prob_shift=0.15, rng=StableRNG(1234))
            min, params = optimize([r1, r2, r3], ps, 20000) do pairs
                x = [param for (field, param) in pairs]
                ackley(x)
            end
            # Compare with analytical solution
            @test isapprox(params[1], 0; atol=1e-3)
            @test round(params[2]) == 0
            @test argmax(params[3:7]) == 1
        end
    end
end


text:
export ConvectionDiffusion

"""
    struct ConvectionDiffusion{N} <: Application{N}
Defines linear convection-diffusion equation u_t + c⋅∇u - dΔu = 0
Default Dirichlet BC (element neighbor == -1): u = 0.0
Default Neumann BC (element neighbor == -2): ∇u⋅n = 0.0
Must initialize struct with constants c,d
"""
struct ConvectionDiffusion{N} <: Application{N}
    c    # Convection constant
    d    # Diffusion constant
    bc_d # Dirichlet BC
    bc_n # Neumann BC
    ConvectionDiffusion{N}(c,d,bc_d=0.0,bc_n=0.0) where {N} = new{N}(c,d,bc_d,bc_n)
    function ConvectionDiffusion{N}(param_dict::Dict{String,Any}) where {N}
        c = param_dict["c"]
        d = param_dict["d"]
        bc_d = get!(param_dict, "bc_d", 0.0)
        bc_n = get!(param_dict, "bc_n", 0.0)
        new{N}(c,d, bc_d,bc_n)
    end
end
is_second_order(app::ConvectionDiffusion) = true
"""
    flux_c!(flux, app::ConvectionDiffusion, u)
Computes the convective flux for convection-diffusion equation: f_c(u) = cu
"""
flux_c!(flux, app::ConvectionDiffusion, u) = (flux .= u*app.c')
"""
    flux_d!(flux, app::ConvectionDiffusion, u, Du)
Computes the diffusion flux for convection-diffusion equation: f_d(u,∇u)=d∇u
"""
flux_d!(flux, app::ConvectionDiffusion, u, Du) where {N} = (flux .= app.d.*Du)

function boundary_flux_c!(flux, app::ConvectionDiffusion, u, bc, normal_k)
    flux_k = similar(u)
    if bc == -1 # Dirichlet BC: u = bc_d
        boundary_u = fill(app.bc_d, size(u))
        flux_c!(flux_k, app, boundary_u)
    elseif bc == -2 # Neumann BC: ∇u⋅n = bc_n
        flux_c!(flux_k, app, u)
    end
    flux .= flux_k.*normal_k
end

function boundary_flux_d!(flux, app::ConvectionDiffusion, u,Du, bc, normal_k)
    fluxes = similar(Du)
    if bc == -1 # Dirichlet BC: u = bc_d
        boundary_u = fill(app.bc_d,size(u))
        flux_d!(fluxes, app, boundary_u, Du)
    elseif bc == -2 # Neumann BC: ∇u⋅n = bc_n
        boundary_Du = fill(app.bc_n,size(Du))
        flux_d!(fluxes, app, u, boundary_Du)
    end
    flux .= fluxes.*normal_k
end

function boundary_flux_l!(flux, app::ConvectionDiffusion, u, bc, normal_k)
    flux_k = similar(u)
    if bc == -1 # Dirichlet BC: u = bc_d
        boundary_u = fill(app.bc_d, size(u))
        flux_l!(flux_k, app, boundary_u)
    elseif bc == -2 # Neumann BC: ∇u⋅n = bc_n
        flux_l!(flux_k, app, u)
    end
    flux .= flux_k.*normal_k
end

max_stars_repo_name: NoseKnowsAll/DGToolkit

max_issues_repo_name: NoseKnowsAll/DGToolkit

max_forks_repo_name: NoseKnowsAll/DGToolkit


text:
export num_grover_step, inference_oracle, GroverIter, groverblock, groveriter, prob_match_oracle

"""
    inference_oracle([nbit::Int,] locs::Vector{Int}) -> ControlBlock

A simple inference oracle, e.g. inference([-1, -8, 5]) is a control block that flip the bit if values of bits on position [1, 8, 5] match [0, 0, 1].
"""
inference_oracle(locs::Vector{Int}) = control(locs[1:end-1], abs(locs[end]) => (locs[end]>0 ? Z : chain(phase(π), Z)))
inference_oracle(nbit::Int, locs::Vector{Int}) = inference_oracle(locs)(nbit)

"""
    target_space(oracle) -> Vector{Bool}

Return a mask, that disired subspace of an oracle are masked true.
"""
function target_space(nbit::Int, oracle)
    r = register(ones(ComplexF64, 1<<nbit))
    apply!(r, oracle)
    real(statevec(r)) .< 0
end

prob_inspace(psi::DefaultRegister, ts) = norm(statevec(psi)[ts])^2

"""
    prob_match_oracle(psi, oracle) -> Float64

Return the probability that `psi` matches oracle.
"""
prob_match_oracle(psi::DefaultRegister, oracle) = prob_inspace(psi, target_space(nqubits(psi), oracle))

"""
    num_grover_step(psi::DefaultRegister, oracle) -> Int

Return number of grover steps needed to match the oracle.
"""
num_grover_step(psi::DefaultRegister, oracle) = _num_grover_step(prob_match_oracle(psi, oracle))

_num_grover_step(prob::Real) = Int(round(pi/4/sqrt(prob)))-1

"""
    GroverIter{N, T}

    GroverIter(oracle, ref::ReflectBlock{N, T}, psi::DefaultRegister, niter::Int)

an iterator that perform Grover operations step by step.
An Grover operation consists of applying oracle and Reflection.
"""
struct GroverIter{N, T}
    psi::DefaultRegister
    oracle
    ref::ReflectBlock{N, T}
    niter::Int
end

groveriter(psi::DefaultRegister, oracle, ref::ReflectBlock{N, T}, niter::Int) where {N, T} = GroverIter{N, T}(psi, oracle, ref, niter)
groveriter(psi::DefaultRegister, oracle, niter::Int) = groveriter(psi, oracle, ReflectBlock(psi |> copy), niter)
groveriter(psi::DefaultRegister, oracle) = groveriter(psi, oracle, ReflectBlock(psi |> copy), num_grover_step(psi, oracle))

function Base.iterate(it::GroverIter, st=1)
    if it.niter + 1 == st
        nothing
    else
        apply!(it.psi, it.oracle)
        apply!(it.psi, it.ref), st+1
    end
end

Base.length(it::GroverIter) = it.niter

"""
    groverblock(oracle, ref::ReflectBlock{N, T}, niter::Int=-1)
    groverblock(oracle, psi::DefaultRegister, niter::Int=-1)

Return a ChainBlock/Sequential as Grover Iteration, the default `niter` will stop at the first optimal step.
"""
function groverblock(oracle::MatrixBlock{N, T}, ref::ReflectBlock{N, T}, niter::Int=-1) where {N, T}
    if niter == -1 niter = num_grover_step(ref.psi, oracle) end
    chain(N, chain(oracle, ref) for i = 1:niter)
end

function groverblock(oracle, ref::ReflectBlock{N, T}, niter::Int=-1) where {N, T}
    if niter == -1 niter = num_grover_step(ref.psi, oracle) end
    sequence(sequence(oracle, ref) for i = 1:niter)
end

groverblock(oracle, psi::DefaultRegister, niter::Int=-1) = groverblock(oracle, ReflectBlock(psi |> copy), niter)


text:
module TestLearningNetworkMachines

using Test
using ..Models
using ..TestUtilities
using MLJBase
using Tables
using StableRNGs
rng = StableRNG(616161)

# A dummy clustering model:
mutable struct DummyClusterer <: Unsupervised
    n::Int
end
DummyClusterer(; n=3) = DummyClusterer(n)
function MLJBase.fit(model::DummyClusterer, verbosity::Int, X)
    Xmatrix = Tables.matrix(X)
    n = min(size(Xmatrix, 2), model.n)
    centres = Xmatrix[1:n, :]
    levels = categorical(1:n)
    report = (centres=centres,)
    fitresult = levels
    return fitresult, nothing, report
end
MLJBase.transform(model::DummyClusterer, fitresult, Xnew) =
    selectcols(Xnew, 1:length(fitresult))
MLJBase.predict(model::DummyClusterer, fitresult, Xnew) =
    [fill(fitresult[1], nrows(Xnew))...]

@testset "wrapping a learning network in a machine" begin

    N = 20
    X = (a = rand(N), b = categorical(rand("FM", N)))

    # unsupervised:
    Xs = source(X)
    W = transform(machine(OneHotEncoder(), Xs), Xs)
    clust = DummyClusterer(n=2)
    m = machine(clust, W)
    yhat = predict(m, W)
    Wout = transform(m, W)

    mach = machine(Unsupervised(), Xs; predict=yhat, transform=Wout)
    @test mach.args == (Xs, )
    @test mach.args[1] == Xs
    fit!(mach, force=true)

    report(mach)
    fitted_params(mach)

    # supervised
    y = rand("ab", N) |> categorical;
    ys = source(y)
    mm = machine(ConstantClassifier(), W, ys)
    yhat = predict(mm, W)
    @test_throws Exception machine(predict=yhat)
    mach = machine(Probabilistic(), Xs, ys; predict=yhat)
    @test mach.model isa Probabilistic
    @test_throws ArgumentError machine(Probabilistic(), Xs, ys)

    # supervised - predict_mode
    fit!(mach)
    @test predict_mode(mach, X) == mode.(predict(mach, X))
    @test predict(mach, rows=1:2) == predict(mach, rows=:)[1:2]

    # evaluate a learning machine
    evaluate!(mach, measure=LogLoss())

    # supervised - predict_median, predict_mean
    X, y = make_regression(20)
    Xs = source(X); ys = source(y)
    mm = machine(ConstantRegressor(), Xs, ys)
    yhat = predict(mm, Xs)
    mach = machine(Probabilistic(), Xs, ys; predict=yhat) |> fit!
    @test predict_mean(mach, X) ≈ mean.(predict(mach, X))
    @test predict_median(mach, X) ≈ median.(predict(mach, X))

end

# build a learning network:
x1 = map(n -> mod(n,3), rand(rng, UInt8, 100)) |> categorical;
x2 = randn(rng, 100);
X = (x1=x1, x2=x2);
y = x2.^2;

Xs = source(X)
ys = source(y)
z = log(ys)
stand = UnivariateStandardizer()
standM = machine(stand, z)
u = transform(standM, z)
hot = OneHotEncoder()
hotM = machine(hot, Xs)
W = transform(hotM, Xs)
knn = KNNRegressor()
knnM = machine(knn, W, u)
oak = DecisionTreeRegressor()
oakM = machine(oak, W, u)
uhat = 0.5*(predict(knnM, W) + predict(oakM, W))
zhat = inverse_transform(standM, uhat)
yhat = exp(zhat)

@testset "replace method for learning network machines" begin

    fit!(yhat)

    # test nested reporting:
    r = MLJBase.report(yhat)
    d = r.report_given_machine
    ms = machines(yhat)
    @test ms == r.machines |> reverse
    @test all(mach -> report(mach) == d[mach], ms)

    hot2 = deepcopy(hot)
    knn2 = deepcopy(knn)

    # duplicate a learning network machine:
    mach  = machine(Deterministic(), Xs, ys; predict=yhat)
    mach2 = replace(mach, hot=>hot2, knn=>knn2,
                    ys=>source(ys.data);
                    empty_unspecified_sources=true)
    ss = sources(glb(mach2.fitresult...))
    @test isempty(ss[1])
    mach2 = @test_logs((:warn, r"No replacement"),
                       replace(mach, hot=>hot2, knn=>knn2,
                               ys=>source(ys.data)))
    yhat2 = mach2.fitresult.predict
    fit!(mach, verbosity=0)
    fit!(mach2, verbosity=0)
    @test predict(mach, X) ≈ predict(mach2, X)

    @test mach2.args[1]() == Xs()
    @test mach2.args[2]() == ys()

    ## EXTRA TESTS FOR TRAINING SEQUENCE

    # pickout the newly created machines:
    standM2 = machines(yhat2, stand) |> first
    oakM2 = machines(yhat2, oak) |> first
    knnM2 = machines(yhat2, knn) |> first
    hotM2 = machines(yhat2, hot) |> first

    @test_mach_sequence(fit!(yhat2, force=true),
                        [(:train, standM2), (:train, hotM2),
                         (:train, knnM2), (:train, oakM2)],
                        [(:train, hotM2), (:train, standM2),
                         (:train, knnM2), (:train, oakM2)],
                        [(:train, standM2), (:train, hotM2),
                         (:train, oakM2), (:train, knnM2)],
                        [(:train, hotM2), (:train, standM2),
                         (:train, oakM2), (:train, knnM2)])

    @test length(MLJBase.machines(yhat)) == length(MLJBase.machines(yhat2))
    @test MLJBase.models(yhat) == MLJBase.models(yhat2)
    @test sources(yhat) == sources(yhat2)
    @test MLJBase.tree(yhat) == MLJBase.tree(yhat2)
    @test yhat() ≈ yhat2()

    # this change should trigger retraining of all machines except the
    # univariate standardizer:
    hot2.drop_last = true
    @test_mach_sequence(fit!(yhat2),
                        [(:skip, standM2), (:update, hotM2),
                         (:train, knnM2), (:train, oakM2)],
                        [(:update, hotM2), (:skip, standM2),
                         (:train, knnM2), (:train, oakM2)],
                        [(:skip, standM2), (:update, hotM2),
                         (:train, oakM2), (:train, knnM2)],
                        [(:update, hotM2), (:skip, standM2),
                         (:train, oakM2), (:train, knnM2)])
end

mutable struct DummyComposite <: DeterministicComposite
    stand1
    stand2
end

@testset "issue 377" begin
    stand = Standardizer()
    model = DummyComposite(stand, stand)

    Xs = source()
    mach1 = machine(model.stand1, Xs)
    X1 = transform(mach1, Xs)
    mach2 = machine(model.stand2, X1)
    X2 = transform(mach2, X1)

    mach = machine(Unsupervised(), Xs; transform=X2)
    @test_logs((:error, r"The hyper"),
               @test_throws(ArgumentError,
                            MLJBase.network_model_names(model, mach)))
end

end

true


text:
# Test the downstream general `solve` in DifferenceEquations.jl
# Currently not included in the tests. Just in casae it is needed in the future.

using DifferentiableStateSpaceModels, DifferenceEquations, LinearAlgebra, Test, Zygote
using Distributions, DistributionsAD
using DifferentiableStateSpaceModels.Examples

@testset "Sequence Simulation, 1st order" begin
    m = @include_example_module(Examples.rbc_observables)
    p_f = (ρ = 0.2, δ = 0.02, σ = 0.01, Ω_1 = 0.1)
    p_d = (α = 0.5, β = 0.95)

    c = SolverCache(m, Val(1), p_d)
    sol = generate_perturbation(m, p_d, p_f; cache = c)

    T = 9
    eps_value = [[0.22], [0.01], [0.14], [0.03], [0.15], [0.21], [0.22], [0.05], [0.18]]
    x0 = zeros(m.n_x)

    # General simulation
    problem = StateSpaceProblem(
        DifferentiableStateSpaceModels.dssm_evolution,
        DifferentiableStateSpaceModels.dssm_volatility,
        DifferentiableStateSpaceModels.dssm_observation,
        x0,
        (0, T),
        sol,
        noise = eps_value
    )
    simul = solve(problem, NoiseConditionalFilter())
    @test simul.z[2:end] ≈
          [[-0.0014843113235688628, 0.0],
           [-0.001672969226342977, -0.013660616212663025],
           [-0.0025907902434120613, -0.016424018091334355],
           [-0.002808352075911724, -0.02507880927301923],
           [-0.003749820417827369, -0.027731843802628737],
           [-0.00514124775142971, -0.036595970922849726],
           [-0.006607765317710925, -0.05006822400565087],
           [-0.006885970365182856, -0.06457803532896965],
           [-0.00789049635407183, -0.06822941930528102]]

    # inference
    @inferred solve(problem, NoiseConditionalFilter())
end

function likelihood_test_joint_first_general(p_d_input, p_f, ϵ, x0, m, tspan, z)
    p_d = (α = p_d_input[1], β = p_d_input[2])
    sol = generate_perturbation(m, p_d, p_f, Val(1))
    problem = StateSpaceProblem(
        DifferentiableStateSpaceModels.dssm_evolution,
        DifferentiableStateSpaceModels.dssm_volatility,
        DifferentiableStateSpaceModels.dssm_observation,
        x0,
        tspan,
        sol,
        noise = ϵ,
        obs_noise = sol.D,
        observables = z
    )
    return solve(problem, NoiseConditionalFilter()).loglikelihood
end

@testset "Gradients, generate_perturbation + likelihood, 1st order" begin
    m = @include_example_module(Examples.rbc_observables)
    p_f = (ρ = 0.2, δ = 0.02, σ = 0.01, Ω_1 = 0.1)
    p_d = (α = 0.5, β = 0.95)
    p_d_input = [0.5, 0.95]
    ϵ = [[0.22], [0.01], [0.14], [0.03], [0.15], [0.21], [0.22], [0.05], [0.18]]
    z = [[-0.6949847708598687, -0.8456988740809867],
         [-0.7804117657996692, 0.07781473603479207],
         [-1.1363021614363802, -2.41253450179418],
         [-0.2140813001516194, -0.10914617826240575],
         [-1.0365874981404577, 0.9869373465251516],
         [-0.7321498641416826, 0.012293325072265942],
         [-0.054809260599132194, -1.8233591236618099],
         [0.5407452466493482, -0.9773559802938866],
         [1.3968232347532277, -2.139194998843768]]
    x0 = zeros(m.n_x)
    tspan = (0, length(z))

    res = gradient((p_d_input, ϵ) -> likelihood_test_joint_first_general(p_d_input, p_f, ϵ, x0, m, tspan, z), p_d_input, ϵ)
    @test res[1] ≈ [303.7133186356109, 553.6149537473261]
    @test res[2] ≈ [[40.62454806083384], [39.38899479341156], [25.095297618483304],
                    [26.06697625612332], [33.10959536324157], [31.484308705831474],
                    [19.172319198105615], [11.464791870737214], [-0.9477420442978448]]

end

@testset "Sequence Simulation, 2nd order" begin
    m = @include_example_module(Examples.rbc_observables)
    p_f = (ρ = 0.2, δ = 0.02, σ = 0.01, Ω_1 = 0.1)
    p_d = (α = 0.5, β = 0.95)

    c = SolverCache(m, Val(2), p_d)
    sol = generate_perturbation(m, p_d, p_f, Val(2); cache = c)

    T = 9
    eps_value = [[0.22], [0.01], [0.14], [0.03], [0.15], [0.21], [0.22], [0.05], [0.18]]
    x0 = zeros(m.n_x)

    # General simulation
    problem = StateSpaceProblem(
        DifferentiableStateSpaceModels.dssm_evolution,
        DifferentiableStateSpaceModels.dssm_volatility,
        DifferentiableStateSpaceModels.dssm_observation,
        [x0; x0],
        (0, T),
        sol,
        noise = eps_value
    )
    simul = solve(problem, NoiseConditionalFilter())
    @test simul.z[2:end] ≈ [[-0.0014120420256672264, -7.824904812715083e-5],
                            [-0.001607843339241866, -0.013798593509356017],
                            [-0.0025317633821568975, -0.016632915260522855],
                            [-0.002755836083102567, -0.02534820495624617],
                            [-0.0037026560860619877, -0.028065833613511802],
                            [-0.005097998299327853, -0.036982697654476426],
                            [-0.006567668012302603, -0.05049239849879983],
                            [-0.006851208817373075, -0.06503101829364497],
                            [-0.007859486666066144, -0.06873403795558215]]

    @inferred solve(problem, NoiseConditionalFilter())
end

function likelihood_test_joint_second(p_d_input, p_f, ϵ, x0, m, tspan, z)
    p_d = (α = p_d_input[1], β = p_d_input[2])
    sol = generate_perturbation(m, p_d, p_f, Val(2))
    problem = StateSpaceProblem(
        DifferentiableStateSpaceModels.dssm_evolution,
        DifferentiableStateSpaceModels.dssm_volatility,
        DifferentiableStateSpaceModels.dssm_observation,
        x0,
        tspan,
        sol,
        noise = ϵ,
        obs_noise = sol.D,
        observables = z
    )
    return solve(problem, NoiseConditionalFilter()).loglikelihood
end

@testset "Gradients, generate_perturbation + likelihood, 2nd order" begin
    m = @include_example_module(Examples.rbc_observables)
    p_f = (ρ = 0.2, δ = 0.02, σ = 0.01, Ω_1 = 0.1)
    p_d = (α = 0.5, β = 0.95)
    p_d_input = [0.5, 0.95]
    ϵ = [[0.22], [0.01], [0.14], [0.03], [0.15], [0.21], [0.22], [0.05], [0.18]]
    z = [[-0.6949847708598687, -0.8456988740809867],
         [-0.7804117657996692, 0.07781473603479207],
         [-1.1363021614363802, -2.41253450179418],
         [-0.2140813001516194, -0.10914617826240575],
         [-1.0365874981404577, 0.9869373465251516],
         [-0.7321498641416826, 0.012293325072265942],
         [-0.054809260599132194, -1.8233591236618099],
         [0.5407452466493482, -0.9773559802938866],
         [1.3968232347532277, -2.139194998843768]]
    x0 = zeros(2 * m.n_x)
    tspan = (0, length(z))

    res = gradient((p_d_input, ϵ) -> likelihood_test_joint_second(p_d_input, p_f, ϵ, x0, m, tspan, z), p_d_input, ϵ)
    @test res[1] ≈ [305.5874661276336, 559.166700806099]
    @test res[2] ≈ [[40.5141940179588], [39.32706019833505], [25.02785099195666],
                    [26.010688843169483], [33.01985483763039], [31.381238099783715],
                    [19.106378855992403], [11.441562042277948], [-0.9454627257067805]]

    # inferred
    # @inferred likelihood_test_joint_second(p_d_input, p_f, ϵ, x0, m, tspan, z)
    # @inferred likelihood_test_joint_second_sol(p_d_input, p_f, ϵ, x0, m, tspan, z)
end


text:
function comparefiles(reference, candidates)
    refc3d = @test_nowarn readc3d(reference; missingpoints=false)
    for candfn in candidates
        file = @test_nowarn readc3d(candfn; missingpoints=false)

        @testset "Parameters equivalency between files" begin
            @testset "Compare groups with $(file.name)" begin
                @test keys(file.groups) ⊆ keys(refc3d.groups)
            end
            for grp in keys(refc3d.groups)
                @testset "Compare :$(refc3d.groups[grp].symname) parameters with $(file.name)" begin
                    @test keys(file.groups[grp].params) ⊆ keys(refc3d.groups[grp].params)
                end
                @testset "Compare the :$(refc3d.groups[grp].symname) parameters" begin
                    for param in keys(refc3d.groups[grp].params)
                        if eltype(refc3d.groups[grp].params[param].data) <: Number
                            if grp == :POINT && param == :SCALE
                                @test abs.(refc3d.groups[grp].params[param].data) ≈ abs.(file.groups[grp].params[param].data)
                            elseif grp == :POINT && param == :DATA_START && any(basename(candfn) .== ("TESTBPI.c3d", "TESTCPI.c3d", "TESTDPI.c3d"))
                                @test file.groups[grp].params[param].data == 20
                            else
                                @test refc3d.groups[grp].params[param].data ≈ file.groups[grp].params[param].data
                            end
                        else
                            @test reduce(*,refc3d.groups[grp].params[param].data .== file.groups[grp].params[param].data)
                        end
                    end
                end
            end
        end

        @testset "Data equivalency between file types" begin
            @testset "Ensure data equivalency between $(refc3d.name) and $(file.name)" begin
                for sig in keys(refc3d.point)
                    @testset "$sig" begin
                        @test haskey(file.point,sig)
                        @test all(isapprox.(refc3d.point[sig], file.point[sig]; atol=0.3))
                    end
                end
                for sig in keys(refc3d.analog)
                    @testset "$sig" begin
                        @test haskey(file.analog,sig)
                        @test all(isapprox.(refc3d.analog[sig], file.analog[sig]; atol=0.3))
                    end
                end
            end
        end
    end
end

@testset "Comparing different files with identical data" begin
    sample01 = ( ("sample01", "Eb015pr.c3d"),
                 ("sample01", "Eb015pi.c3d"),
                 ("sample01", "Eb015vr.c3d"),
                 ("sample01", "Eb015vi.c3d"),
                 ("sample01", "Eb015sr.c3d"),
                 ("sample01", "Eb015si.c3d"))

    sample02 = ( ("sample02", "pc_real.c3d"),
                 ("sample02", "pc_int.c3d"),
                 ("sample02", "dec_real.c3d"),
                 ("sample02", "dec_int.c3d"),
                 ("sample02", "sgi_real.c3d"),
                 ("sample02", "sgi_int.c3d"))

    sample08 = ( ("sample08", "EB015PI.c3d"),
                 ("sample08", "TESTAPI.c3d"),
                 ("sample08", "TESTBPI.c3d"),
                 ("sample08", "TESTCPI.c3d"),
                 ("sample08", "TESTDPI.c3d"))

    @testset "sample01 files" begin
        comparefiles(joinpath(datadir, sample01[1]...), broadcast(x -> joinpath(datadir, x...), sample01[2:end]))
    end
    @testset "sample02 files" begin
        comparefiles(joinpath(datadir, sample02[1]...), broadcast(x -> joinpath(datadir, x...), sample02[2:end]))
    end
    @testset "sample08 files" begin
        comparefiles(joinpath(datadir, sample08[1]...), broadcast(x -> joinpath(datadir, x...), sample08[2:end]))
    end
end


text:
module TwoBoxModel
#=
This module serves to load the 2-box model matrix and grid
from the Sarmiento and Gruber book (2006) book.

┌────────────────────────────────┐
│                        Surface │
├────────────────────────────────┤
│                                │
│                                │
│                                │
│                                │
│                           Deep │
└────────────────────────────────┘

=#

using LinearAlgebra, SparseArrays
using Unitful
using Reexport
@reexport using OceanGrids            # To store the grid
using ..CirculationGeneration
CG = CirculationGeneration

build_wet3D() = trues(1, 1, 2)

function build_grid()
    wet3D = build_wet3D()
    vdeep = 1.26e18u"m^3"     # Sarmiento and Gruber, 2006
    vtot = 1.292e18u"m^3"     # Archer et al., 2000
    area_ocean = 349e12u"m^2" # Archer et al., 2000
    # Infer other values
    depth_surface = (vtot-vdeep) / area_ocean
    depth_ocean = vtot / area_ocean
    # Build DY3d
    elon = [0,360] * u"°"
    elat = [-90, 90] * u"°"
    edepth = [0u"m", depth_surface, depth_ocean]
    return OceanGrid(elon, elat, edepth, wet3D)
end

function build_T(grid)
    # From Archer et al. [2000]
    v3D = array_of_volumes(grid)
    nb = length(v3D)
    # ν = 38 Sv (Sarmiento and Gruber, 2008) (1 Sv = 10⁶ m³/s)
    T  = CG.T_diffusion(38e6u"m^3/s", 1, 2, v3D, nb)
    return T
end

"""
    load

Returns `grd` and `T` (in that order).
"""
function load()
    print("Creating 2-box model")
    grid = build_grid()
    T = build_T(grid)
    println(" ✔")
    @info """
            You are about to use the 2-box model of Sarmiento and Gruber [2006].
            If you use it for research, please cite:

            - Sarmiento, Jorge L., and Nicolas Gruber. Ocean biogeochemical dynamics. Princeton University Press, 2006.

            You can find the corresponding BibTeX entry in the CITATION.bib file
            at the root of the AIBECS.jl package repository.
            (Look for the "Sarmiento_Gruber_2006" key.)
            """
    return grid, T
end

end

export TwoBoxModel


text:
using Test
using CrypticCrosswords.SemanticSimilarity
using Random
using WordNet
using ProgressMeter

# include("caches.jl")
include("metrics.jl")


text:
# Julia wrapper for header: rocfft-version.h
# Automatically generated using Clang.jl

# Julia wrapper for header: rocfft.h
# Automatically generated using Clang.jl

function rocfft_setup()
    @check ccall((:rocfft_setup, librocfft), rocfft_status, ())
end

function rocfft_cleanup()
    @check ccall((:rocfft_cleanup, librocfft), rocfft_status, ())
end

function rocfft_plan_create(plan, placement, transform_type, precision, dimensions,
                            lengths, number_of_transforms, description)
    @check ccall((:rocfft_plan_create, librocfft), rocfft_status,
                 (Ptr{rocfft_plan}, rocfft_result_placement, rocfft_transform_type,
                  rocfft_precision, Csize_t, Ptr{Csize_t}, Csize_t, rocfft_plan_description),
                 plan, placement, transform_type, precision, dimensions, lengths,
                 number_of_transforms, description)
end

function rocfft_execute(plan, in_buffer, out_buffer, info)
    @check ccall((:rocfft_execute, librocfft), rocfft_status,
                 (rocfft_plan, Ptr{Ptr{Cvoid}}, Ptr{Ptr{Cvoid}}, rocfft_execution_info),
                 plan, in_buffer, out_buffer, info)
end

function rocfft_plan_destroy(plan)
    @check ccall((:rocfft_plan_destroy, librocfft), rocfft_status,
                 (rocfft_plan,),
                 plan)
end

function rocfft_plan_description_set_data_layout(description, in_array_type,
                                                 out_array_type, in_offsets, out_offsets,
                                                 in_strides_size, in_strides, in_distance,
                                                 out_strides_size, out_strides, out_distance)
    @check ccall((:rocfft_plan_description_set_data_layout, librocfft), rocfft_status,
                 (rocfft_plan_description, rocfft_array_type, rocfft_array_type,
                  Ptr{Csize_t}, Ptr{Csize_t}, Csize_t, Ptr{Csize_t}, Csize_t, Csize_t,
                  Ptr{Csize_t}, Csize_t),
                 description, in_array_type, out_array_type, in_offsets, out_offsets,
                 in_strides_size, in_strides, in_distance, out_strides_size, out_strides,
                 out_distance)
end

function rocfft_get_version_string(buf, len)
    @check ccall((:rocfft_get_version_string, librocfft), rocfft_status,
                 (Cstring, Csize_t),
                 buf, len)
end

function rocfft_plan_get_work_buffer_size(plan, size_in_bytes)
    @check ccall((:rocfft_plan_get_work_buffer_size, librocfft), rocfft_status,
                 (rocfft_plan, Ptr{Csize_t}),
                 plan, size_in_bytes)
end

function rocfft_plan_get_print(plan)
    @check ccall((:rocfft_plan_get_print, librocfft), rocfft_status,
                 (rocfft_plan,),
                 plan)
end

function rocfft_plan_description_create(description)
    @check ccall((:rocfft_plan_description_create, librocfft), rocfft_status,
                 (Ptr{rocfft_plan_description},),
                 description)
end

function rocfft_plan_description_destroy(description)
    @check ccall((:rocfft_plan_description_destroy, librocfft), rocfft_status,
                 (rocfft_plan_description,),
                 description)
end

function rocfft_execution_info_create(info)
    @check ccall((:rocfft_execution_info_create, librocfft), rocfft_status,
                 (Ptr{rocfft_execution_info},),
                 info)
end

function rocfft_execution_info_destroy(info)
    @check ccall((:rocfft_execution_info_destroy, librocfft), rocfft_status,
                 (rocfft_execution_info,),
                 info)
end

function rocfft_execution_info_set_work_buffer(info, work_buffer, size_in_bytes)
    @check ccall((:rocfft_execution_info_set_work_buffer, librocfft), rocfft_status,
                 (rocfft_execution_info, Ptr{Cvoid}, Csize_t),
                 info, work_buffer, size_in_bytes)
end

function rocfft_execution_info_set_stream(info, stream)
    @check ccall((:rocfft_execution_info_set_stream, librocfft), rocfft_status,
                 (rocfft_execution_info, Ptr{Cvoid}),
                 info, stream)
end
# Julia wrapper for header: rocfft-export.h
# Automatically generated using Clang.jl


text:
include("settings.jl")

makedocs(sitename = sitename; settings...)

deploydocs(repo = "github.com/ITensor/ITensors.jl.git")


text:
"""
    getValidDatabases()

Return a `DataFrame` object with two columns, `DatabaseDisplayName` and `DatabaseId`.
The dataframe containing all supported databases. Values from the column
`DatabaseId` can be used as input parameter `targetDb` and `sourceDb`
in several functions contained in the GeneNameGenie module, such as the function
`convertFromTo`.
"""
function getValidDatabases()
   x = Neo4j.cypherQuery(NEO4J_CONNECTION, "CALL rcsi.params.getValidDatabases() YIELD value RETURN value.DatabaseDisplayName AS DatabaseDisplayName, value.DatabaseId AS DatabaseId");
   return x;
end

"""
    getValidAttributes()

Return a `DataFrame` object with a single column, `Parameters`. Values from
the `Parameters` column can be used as input parameter `attributes` in the
`convertFromToExtended` function, for example `ensg:region`, which will return
the genomic region of the corresponding Ensembl gene.
"""
function getValidAttributes()
   x = Neo4j.cypherQuery(NEO4J_CONNECTION, "CALL rcsi.params.getValidAttributes() YIELD value RETURN value AS Parameters");
   return x;
end

"""
    getValidMirnaMetadataValues()

Return a `DataFrame` object with a single column, `Parameters`. Values from
the `Parameters` column can be used as input parameter `metadata` in the
`convertToCurrentMirbaseVersion` function, for example the parameter `region`
returns the genomic region of the corresponding precursor miRNA gene.
"""
function getValidMirnaMetadataValues()
   x = Neo4j.cypherQuery(
         NEO4J_CONNECTION,
         "CALL rcsi.params.getValidMirnaMetadataValues() YIELD value RETURN value AS Parameter");
   return x;
end

"""
    getCurrentMirbaseVersion()

Return the current miRBase release version information (the highest supported version).
"""
function getCurrentMirbaseVersion()
   x = Neo4j.cypherQuery(
         NEO4J_CONNECTION,
         "MATCH (db :MirnaDB {id: 'miRBase_mature_name'}) RETURN db.release AS release;")[1,:release];
   # mirbaseVersion = parse(match(r"\d+(\.?\d+)", split(x, ",")[1]).match);
   # mirbaseVersion = split(x, ",");
   return x;
end

"""
    getEnsemblVersion()

Return the Ensembl DB version of GeneNameGenie.
"""
function getEnsemblVersion()
   x = Neo4j.cypherQuery(
         NEO4J_CONNECTION,
         "MATCH (db :EnsemblDB {id: 'ArrayExpress'}) RETURN db.release AS release;")[1,:release];
   # mirbaseVersion = parse(match(r"\d+(\.?\d+)", split(x, ",")[1]).match);
   # mirbaseVersion = split(x, ",");
   return x;
end

max_stars_repo_name: StefanHaunsberger/GeneNameGenieJ

max_issues_repo_name: StefanHaunsberger/GeneNameGenieJ

max_forks_repo_name: StefanHaunsberger/GeneNameGenieJ


text:
###############################################################################
#
# Diff function to minimize by ConstrainedRootSolvers
# Useful for nocturnal stomatal conductance model
#
###############################################################################
"""
    dRdE(clayer::CanopyLayer{FT},
         envir::AirLayer{FT},
         LAI::FT
    ) where {FT<:AbstractFloat}

Calculate the marginal decrease of respiration rate, given
- `clayer` [`CanopyLayer`](@ref) type of struct
- `envir` `AirLayer` type struct
- `LAI` Total leaf area index of whole canopy
"""
function dRdE(
            clayer::CanopyLayer{FT},
            envir::AirLayer{FT}
) where {FT<:AbstractFloat}
    @unpack T = clayer;

    return dTdE(clayer, envir) * clayer.ps.PSM.r_d * clayer.ps.PSM.TD_R.ΔHA / GAS_R(FT) / T^2
end

"""
    dTdE(clayer::CanopyLayer{FT},
         envir::AirLayer{FT}
    ) where {FT<:AbstractFloat}

Calculate the margian decrease in leaf temperature, given
- `clayer` [`CanopyLayer`](@ref) type of struct
- `envir` `AirLayer` type struct
- `LAI` Total leaf area index of whole canopy
"""
function dTdE(
            clayer::CanopyLayer{FT},
            envir::AirLayer{FT}
) where {FT<:AbstractFloat}
    @unpack T, tLAI, width = clayer;
    @unpack wind = envir;

    lambda = latent_heat_vapor(T) * M_H₂O(FT);
    gbe    = FT(0.189 * sqrt(wind/(0.72*width)));
    emis   = FT(0.97);
    denom  = 2 * CP_D_MOL(FT) * gbe +
             4 / tLAI * K_STEFAN(FT) * emis * T^3;
    ∂T∂E   = lambda / denom;

    return ∂T∂E
end

"""
    dΘdE(clayer::CanopyLayer{FT},
         sm::OSMWang{FT}(),
         g_sw::FT
    ) where {FT<:AbstractFloat}

Calculate the margian carbon cost related to nighttime transpiration, given
- `clayer` [`CanopyLayer`](@ref) type of struct
- `sm` [`OSMWang`](@ref) type stomatal model
- `g_sw` Given leaf level stomatal conductance
"""
function dΘdE(
            clayer::CanopyLayer{FT},
            sm::OSMWang{FT},
            g_sw::FT
) where {FT<:AbstractFloat}
    @unpack APAR_m, ec, envir_m, ff = clayer;
    @unpack P_AIR, p_H₂O = envir_m
    clayer.ps_m.apar = APAR_m;

    # calculate g_lc and a_net using memory clayer and envir
    g_sc = g_sw / FT(1.6);
    g_lc = 1 / (1/g_sc + 1/clayer.ps_m.g_CO₂_b);
    g_bw = clayer.ps_m.g_CO₂_b * FT(1.35);
    g_lw = 1 / (1/g_sw + 1/g_bw);
    flow = g_lw * max(1, clayer.ps_m.p_H₂O_sat - p_H₂O) / P_AIR;
    leaf_photosynthesis!(clayer.ps_m, envir_m, GCO₂Mode(), g_lc);

    return clayer.ps_m.PSM.a_net / (ec - flow) * ff
end

"""
    nocturnal_diff!(
                x::FT,
                clayer::CanopyLayer{FT},
                envir::AirLayer{FT},
                sm::OSMWang{FT}()
    ) where {FT<:AbstractFloat}

Calculate the difference between marginal gain and risk of nighttime
    transpiration, given
- `x` Given leaf level stomatal conductance
- `clayer` [`CanopyLayer`](@ref) type of struct
- `envir` `AirLayer` type struct
- `sm` [`OSMWang`](@ref) type stomatal model
"""
function nocturnal_diff!(
            x::FT,
            clayer::CanopyLayer{FT},
            envir::AirLayer{FT},
            sm::OSMWang{FT}
) where {FT<:AbstractFloat}
    ∂R∂E = dRdE(clayer, envir);
    ∂Θ∂E = dΘdE(clayer, sm, x);

    return ∂R∂E - ∂Θ∂E
end


text:
using Documenter, PDSampler

try
    # refreshes all the examples
    include("readexamples.jl")
catch e
    # probably on travis CI
    println(e)
end

makedocs(
    modules = [PDSampler],
    format = :html,
    sitename = "PDSampler.jl",
    authors = "Thibaut Lienart",
    pages = Any[
        "Introduction" => "index.md",
        "About PDSampler"   => "aboutpdsampler.md",
        "Examples" => Any[
            "Global BPS" => "examples/ex_gbps1.md",
            "Local BPS"  => "examples/ex_lbps1.md"
        ],
        "Technical Documentation" => Any[
            "Code structure" => "techdoc/structure.md",
            "Core tools"     => "techdoc/coretools.md",
            "Models"         => "techdoc/models.md",
            "Global sampler" => "techdoc/global.md",
            "Local sampler"  => "techdoc/local.md",
#            "Types"          => "techdoc/types.md",
        ],
        "Contributing" => Any[
            "New example" => "contributing/addingexample.md"
            "New feature" => "contributing/addingfeature.md"
        ]
    ],
    assets = String[
                "assets/partial.css",
                "assets/BPS.svg",
                "asserts/truncatedgaussian.png"],
)

deploydocs(
    repo   = "github.com/alan-turing-institute/PDSampler.jl.git",
    deps   = Deps.pip("mkdocs", "pygments"),
    target = "build",
    julia  = "0.5",
    make = nothing
)


text:
import MPI, HDF5, RIOPA
import Test: @testset, @test, @test_throws

@testset "main" begin
    @test RIOPA.main(["hello"]) == 0
    if HDF5.has_parallel()
        @test_exists_and_rm("hello.h5")
    else
        @test_exists_and_rm("hello_0.h5")
    end
    filename = "testcase-temp-config.yaml"
    @test RIOPA.main(["-c", filename, "generate-config"]) == 0
    @test_exists_and_rm(filename)
end


text:
using JLD2
using Oceananigans.Fields: location, XFaceField, YFaceField, ZFaceField, CenterField
using Oceananigans.Grids: Face, Center

location_guide = Dict(:u => (Face, Center, Center),
    :v => (Center, Face, Center),
    :w => (Center, Center, Face))

# If the location isn't correct in `observations`, `interior` won't be computed correctly
# in InverseProblems.transpose_model_output
function infer_location(field_name)
    if field_name in keys(location_guide)
        return location_guide[field_name]
    else
        return (Center, Center, Center)
    end
end

function legacy_data_field_time_serieses(path, field_names, times)

    # Build a grid, assuming it's a 1D RectilinearGrid
    file = jldopen(path)

    old_Nz = file["grid/Nz"]
    Hz = file["grid/Hz"]
    Lz = file["grid/Lz"]

    αg = file["buoyancy/gravitational_acceleration"] * file["buoyancy/equation_of_state/α"]
    Qᵇ = file["parameters/boundary_condition_θ_top"] * αg
    Qᵘ = file["parameters/boundary_condition_u_top"]
    u_bottom = file["parameters/boundary_condition_u_bottom"]
    b_bottom = file["parameters/boundary_condition_θ_bottom"] * αg

    close(file)

    b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ), bottom = GradientBoundaryCondition(b_bottom))
    u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ), bottom = GradientBoundaryCondition(u_bottom))

    grid = RectilinearGrid(size = old_Nz, halo = Hz, z = (-Lz, 0), topology = (Flat, Flat, Bounded))

    boundary_conditions = (; u = u_bcs, b = b_bcs)

    field_time_serieses = NamedTuple(name => FieldTimeSeries(path, string(name);
                                                             times = times,
                                                             location = infer_location(name),
                                                             grid = grid,
                                                             boundary_conditions = boundary_conditions)
                                    for name in field_names)

    return field_time_serieses
end


text:
module Classifiers
export LogisticClassifier, centerscale!, add_designmatrix!, fit!, designmatrix, predict, score, traintestsplit

using ..Classification: Optimizer, GDContext, SGDContext, σ, accuracy
import ..Classification: fit!, addvalidationset!
using Random: shuffle
using Statistics

mutable struct LogisticClassifier
    optimizer::Optimizer
    designmatrix::Matrix{<:Real}
    β::Vector{Float32}
    mean::Matrix{<:Real}
    scale::Matrix{<:Real}
    function LogisticClassifier(optimizer::V) where V<:Optimizer
        new(optimizer, Matrix{Float32}(undef, 0, 0),
            Float64[], Matrix{Float32}(undef, 0, 0),
            Matrix{Float32}(undef, 0, 0))
    end
end
function LogisticClassifier(context::GDContext)
    LogisticClassifier(GradientDescent(context))
end

function LogisticClassifier(context::SGDContext)
    LogisticClassifier(StochasticGradientDescent(context))
end

function centerscale!(classifier::LogisticClassifier)
    classifier.mean = mean(classifier.designmatrix, dims=1)
    classifier.scale = std(classifier.designmatrix, dims=1)
end

function designmatrix(classifier::LogisticClassifier, X::Matrix{Float64}; center=true,
                       scale=true, addintercept=true)::Matrix{Float64}
    mat = copy(X)
    center && mat .-= classifier.mean
    scale  && mat ./= classifier.scale
    if addintercept
        [ones(typeof(mat[1]), size(mat, 1)) mat]
    else
        mat
    end
end

function add_designmatrix!(classifier::LogisticClassifier, X::Matrix{T}; center=true,
                           scale=true, addintercept=true) where T<:Real
    classifier.designmatrix = copy(X)
    centerscale!(classifier)
    center && classifier.designmatrix .-= classifier.mean
    scale  && classifier.designmatrix ./= classifier.scale
    if addintercept
        classifier.designmatrix = [ones(typeof(classifier.designmatrix[1]), size(classifier.designmatrix, 1)) classifier.designmatrix]::Matrix{T}
    end
    classifier.β = zeros(Float64, size(classifier.designmatrix, 2))
end

function fit!(classifier::LogisticClassifier, X::Matrix{<:Real}, y::Vector{<:Real})
    add_designmatrix!(classifier, X)
    fit!(classifier, y)
end

function addvalidationset!(classifier::LogisticClassifier, X, y, earlystopping=false)
    X = designmatrix(classifier, X)
    addvalidationset!(classifier.optimizer, X, y, earlystopping)
end

function fit!(classifier::LogisticClassifier, y::Vector{<:Real})
    fit!(classifier.optimizer,
         classifier.β,
         classifier.designmatrix,
         y)
end

function predict(classifier::LogisticClassifier, X)
    X = designmatrix(classifier, X)
    y = (σ.(X*classifier.β) .> 0.5) .|> Int
end

function score(classifier::LogisticClassifier, X, y)
    ŷ = predict(classifier, X)
    accuracy(ŷ, y)
end

function traintestsplit(X, y, splitratio, portion=nothing)
    if portion === nothing
        portion = size(X, 1)
    end
    if portion > size(X, 1)
        @warn "Data set too small; using all"
        portion = size(X, 1)
    end
    M = shuffle(1:length(y))
    traintest = splitratio*portion |> x -> round(Int, x)
    Xtrain = X[M, :][1:traintest, :]
    Ytrain = y[M][1:traintest]
    Xtest  = X[M, :][traintest+1:portion, :]
    Ytest  = y[M][traintest+1:portion];
    (Xtrain, Ytrain), (Xtest, Ytest)
end

end


text:
function f(x, y)
    2x+4y
end


text:
using IndexableBitVectors
import IndexableBitVectors
using Humanize

srand(1024)

function bench{T}(::Type{T}, len::Int, ord::Vector{Int}; r=0.5)
    b = convert(T, rand(len) .> (1.0 - r))
    b[1]
    rank1(b, 1)
    select1(b, 1)

    n = 0

    gc()
    t0 = time_ns()
    for i in ord; n += b[i]; end
    t1 = time_ns()
    @printf "%s (access) : %10.3f ns/op\n" T (t1 - t0) / length(ord)

    gc()
    t0 = time_ns()
    for i in ord; n += rank1(b, i); end
    t1 = time_ns()
    @printf "%s (rank1)  : %10.3f ns/op\n" T (t1 - t0) / length(ord)

    gc()
    t0 = time_ns()
    for i in ord; n += select1(b, i); end
    t1 = time_ns()
    @printf "%s (select1): %10.3f ns/op\n" T (t1 - t0) / length(ord)

    return sizeof(b)
end

let
    random = false
    r = 0.5
    while !isempty(ARGS)
        arg = shift!(ARGS)
        if arg == "--random" || arg == "-r"
            random = true
        elseif arg == "--sparsity" || arg == "-s"
            r = parse(Float64, shift!(ARGS))
        else
            error("$arg")
        end
    end

    for p in [6, 12, 18, 24, 30]
        len = 2^p
        ord = collect(1:len)
        if random
            ord = ord[rand(1:length(ord), 10_000)]
        end
        println("length: $len bits = $(datasize(div(len, 8), style=:bin))")
        for t in [CompactBitVector, SucVector, CSucVector, RRR, RRRNP]
            size = bench(t, len, ord, r=r)
            println("sizeof(bitvector) = $(datasize(size, style=:bin))")
        end
        println()
    end
end


text:
#==============================================================================#
# OpenSCADPolyhedron.jl
# December 2014
#==============================================================================#

# This module helps generate complex polyhedron commands for OpenSCAD.
# The key function is generatePolyhedron, which takes an initial set of points
# and faces, a set of transformations, and a dictionary of symmetries. All the
# transforms are applied to the initial base polyhedron, resulting in a new set
# of points and faces. The symmetry dictionary is used to eradicate duplicate
# points, so that the resulting polyhedron can be described as simply as
# possible.

# Background
# The polyhedron command in OpenSCAD looks like this:
# polyhedron(points = [[-1.0, 0.0, 0.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 1.0]],
#         faces = [[0, 1, 2], [0, 3, 1], [0, 2, 3], [1, 3, 2]])
# Note that the faces must be specified clocwise, when viewing the final shape from
# the outside.

#-------------------------------------------------------------------------------
module OpenSCADPolyhedron
import Base.append!

export Point3, Polyhedron, stringify, generatePolyhedron, expandSymmetries

#-------------------------------------------------------------------------------
type Point3{T<:Real}
  x::Array{T,1}
end

# This object can easily be turned into a polyhedron command in OpenSCAD.
type Polyhedron{T<:Real}
  points::Array{Point3{T},1}
  faces::Array{Array{Int,1},1}

  Polyhedron{T<:Real} (points::Array{Point3{T},1}) = new(points, Array{Int,1}[])
  Polyhedron{T<:Real} (points::Array{Point3{T},1}, faces::Array{Array{Int,1},1}) = new(points, faces)
end
append!{T} (poly::Polyhedron, point::Point3{T}) = append!(poly.points, Point3{T}[point])
append!{T} (poly::Polyhedron, points::Array{Point3{T},1}) = append!(poly.points, points)
append! (p::Polyhedron, face::Array{Int,1}) = append!(p.faces, Array{Int,1}[face])

# These methods are used to generate the OpenSCAD command strings.
stringify (p::Point3) = "[$(p.x[1]), $(p.x[2]), $(p.x[3])]"
stringify (x::Int) = string(x-1)
function stringify{T} (collection::Array{T})
  strings = [ stringify(x) for x in collection ]
  return "[$(join(strings, ", "))]"
end
stringify (poly::Polyhedron) = "polyhedron(points = $(stringify(poly.points)), faces = $(stringify(poly.faces)));"

#-------------------------------------------------------------------------------
# This object is a helper object used by generatePolyhedron.
# For each transformation i and point j, it stores the index of the resulting
# point in the points array of the polyhedron object being built.
# This allows duplicate points to not be recorded in the polyhedron object,
# since the same index can be recorded multiple times in this object.
type PointIndex
  npoints::Int
  pointids::Array{Int}

  PointIndex (npoints::Int, ntransforms::Int) = new(npoints, Array(Int, npoints*ntransforms))
end
packindex (p::PointIndex, i::Int, j::Int) = p.npoints*(i-1) + j
getindex (p::PointIndex, i::Int) = p.pointids[i]
getindex (p::PointIndex, i::Int, j::Int) = p[packindex(p, i, j)]
setindex! (p::PointIndex, x::Int, i::Int) = p.pointids[i] = x
setindex! (p::PointIndex, x::Int, i::Int, j::Int) = p[packindex(p, i, j)] = x

#-------------------------------------------------------------------------------
# This is the critical function. It generates a new polyhedron that is composed
# of one copy of the base polyhedron (p) for each transform. The symmetries
# dictionary is used to link these polyhedra together, so that their shared
# vertices (and hence edges) need not be recorded multiple times in the final
# command.
function generatePolyhedron{T} (p::Polyhedron{T}, transforms::Array{Function}, symmetries::Dict{(Int, Int), (Int, Int)})
  nt = length(transforms)
  np = length(p.points)
  nf = length(p.faces)
  pointids = PointIndex(np, nt)
  result = Polyhedron{T}(Point3{T}[])
  for i in 1:nt
    for j in 1:np
      if haskey(symmetries, (i, j))
        pointids[i, j] = pointids[symmetries[(i, j)]...]
      else
        append!(result, transforms[i](p.points[j]))
        pointids[i, j] = length(result.points)
      end
    end
    for k in 1:nf
      face::Array{Int,1} = p.faces[k]
      for j in 1:length(face)
      end
      append!(result, [ pointids[i, face[j]] for j in 1:length(face) ])
    end
  end
  return result
end

#-------------------------------------------------------------------------------
# This function allows more complicated shapes to be built up from simpler ones
function expandSymmetries (fundamentalSymmetries::Dict{(Int, Int), (Int, Int)}, fundamentalFace::Array{Int,1}, layers::Set{Array{Int,1}})
  symmetries = Dict{(Int, Int), (Int, Int)}()
  for layer in layers
    # build point lookup table
    lookup = Dict{Int, Int}()
    for i::Int in 1:length(fundamentalFace)
      lookup[fundamentalFace[i]] = layer[i]
    end

    # build new symmetries
    for ((t1, p1), (t2, p2)) in fundamentalSymmetries
      symmetries[(t1, lookup[p1])] = (t2, lookup[p2])
    end
  end

  return symmetries
end

end

max_stars_repo_name: erikstrand/OpenSCADPolyhedron

max_issues_repo_name: erikstrand/OpenSCADPolyhedron

max_forks_repo_name: erikstrand/OpenSCADPolyhedron


text:
function setCollisionSmoothingRadius(shape, collisionSmoothingRadius) end

# function setCollisionSmoothingRadius(shape::Sphere) end # not defined

# function setCollisionSmoothingRadius(shape::Ellipsoid) end # not defined

function setCollisionSmoothingRadius(shape::Box, collisionSmoothingRadius)
    @assert(collisionSmoothingRadius >= 0.0)
    collisionSmoothingRadius2 = min(collisionSmoothingRadius, 0.1*min(shape.lengthX, shape.lengthY, shape.lengthZ)) # at most 10% of the smallest edge length
    if collisionSmoothingRadius2 < collisionSmoothingRadius
      println("collisionSmoothingRadius in Box changed from $collisionSmoothingRadius to $collisionSmoothingRadius2")
    end
    @assert(shape.lengthX >= 2*collisionSmoothingRadius2)
    @assert(shape.lengthY >= 2*collisionSmoothingRadius2)
    @assert(shape.lengthZ >= 2*collisionSmoothingRadius2)
    return collisionSmoothingRadius2
end

function setCollisionSmoothingRadius(shape::Cylinder, collisionSmoothingRadius)
    @assert(collisionSmoothingRadius >= 0.0)
    collisionSmoothingRadius = min(collisionSmoothingRadius, 0.1*min(shape.diameter, shape.length))  # at most 10% of the smallest edge length
    return collisionSmoothingRadius
end

function setCollisionSmoothingRadius(shape::Cone, collisionSmoothingRadius)
    @assert(collisionSmoothingRadius >= 0.0)
    collisionSmoothingRadius2 = min(collisionSmoothingRadius, 0.1*min(shape.diameter, shape.length)) # at most 10% of the smallest edge length
    return collisionSmoothingRadius2
end

function setCollisionSmoothingRadius(shape::Capsule, collisionSmoothingRadius)
    @assert(collisionSmoothingRadius >= 0.0)
    collisionSmoothingRadius2 = min(collisionSmoothingRadius, 0.1*min(shape.diameter, shape.length)) # at most 10% of the smallest edge length
    return collisionSmoothingRadius2
end

function setCollisionSmoothingRadius(shape::Beam, collisionSmoothingRadius)
    @assert(collisionSmoothingRadius >= 0.0)
    collisionSmoothingRadius2 = min(collisionSmoothingRadius, 0.1*min(shape.length, shape.width, shape.thickness)) # at most 10% of the smallest edge length
    return collisionSmoothingRadius2
end


text:
module CUDA

using Compat
import Compat.String
importall ..JuCUDA

include(joinpath(Pkg.dir("JuCUDA"),"lib/7.5/libcuda.jl"))
include(joinpath(Pkg.dir("JuCUDA"),"lib/7.5/libcuda_types.jl"))

@compat if is_windows()
    const libcuda = Libdl.find_library(["nvcuda"])
else
    const libcuda = Libdl.find_library(["libcuda"])
end
isempty(libcuda) && error("CUDA library cannot be found.")

function check_curesult(status)
    status == CUDA_SUCCESS && return nothing
    warn("CUDA error triggered from:")
    Base.show_backtrace(STDOUT, backtrace())
    p = Ptr{UInt8}[0]
    cuGetErrorString(status, p)
    throw(bytestring(p[1]))
end

function version()
    p = Cint[0]
    cuDriverGetVersion(p)
    Int(p[1])
end

include("context.jl")
include("device.jl")
include("stream.jl")
include("pointer.jl")
include("module.jl")
include("function.jl")
include("execution.jl")

const contexts = CUcontext[]
const streams = CuStream[]

info("Initializing CUDA...")
cuInit(0)
for dev = 0:ndevices()-1
    p = CUcontext[0]
    cuCtxCreate(p, 0, dev)
    push!(contexts, p[1])
end
setdevice(0)
info("CUDA driver version: $(version())")

end


text:
push!(LOAD_PATH, "/users/creanj/julialib_fork/PUMI.jl")
using PumiInterface

dmg_name = ".null"
smb_name = "tri8l.smb"

tmp, num_Entities, m_ptr, mshape_ptr = apf.init2(dmg_name, smb_name, 1)

apf.writeVtkFiles("output_pre", m_ptr)

#=
println("about to create field")
f_ptr = apf.createPackedField(m_ptr, "testfield", 4)
println("created field")
comps = [1.0, 2, 3, 4]

resetVertIt()
vert_ptr = getVert()  # get a vertex
apf.setComponents(f_ptr, vert_ptr, 0, comps)
println("set components")
comps_retrieved = [0.0, 0, 0, 0]
apf.getComponents(f_ptr, vert_ptr, 0, comps_retrieved)
println("comps_retrieved = ", comps_retrieved)
=#

function func2(entity_ptr, m_ptr, u_::Ptr{T}) where T
# determine the mesh size at a vertex
# this function can have arbitrary arguments, but can only access data that is passed in as an argument (cannot violate trival closure condition)

  u = pointer_to_array(u_, 9)  # pointer_to_array(p, dims), where p is the pointer, dims can either be an integer that is the number of elements for a 1D array, or an array (tuple?) that tells the number of elements in each dimension of the array.  This enables the C code that is behind pointer_to_array to make sense of the data laid out in memory
  coords = zeros(3,1)
  apf.getVertCoords(entity_ptr, coords, 3, 1)
  x_coords = coords[1]
  frac = (x_coords + 1.5)/2.0
  h_value = 1.25*frac
  println("u = \n", u)

  return h_value  # a Cdouble
end

function func3(entity_ptr, r_::Ptr{T}, h_::Ptr{T}, m_ptr, u_::Ptr{T}) where T
# an anisotropic function
# populates h with the desired mesh size in all three dimension

  println("entered func3")
  # load arrays from C pointers
  r = pointer_to_array(r_, (3,3))  # REMEMBER: r is transposed when it is passed back to C
  println("r = \n", r)
  h = pointer_to_array(h_, 3)
  println("h = \n", h)
  u = pointer_to_array(u_, 9)
  println("u = \n", u)

  r[1,1] = 1.0
  r[2,2] = 1.0
  r[3,3] = 1.0

  println("in julia, r = ", r)

  h[1] = 0.5
  h[2] = 1.0
  h[3] = 1.0

return nothing
end

u = ones(Float64, 9)

#=
cfunc2 = cfunction(func2, Cdouble, (Ptr{Void}, Ptr{Void}, Ptr{Float64}))
println("typeof(cfunc2 = ", typeof(cfunc2))
apf.createIsoFunc(m_ptr, cfunc2, u)
runIsoAdapt(m_ptr)
=#
cfunc3 = cfunction(func3, Void, (Ptr{Void}, Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Ptr{Float64}))

createAnisoFunc(m_ptr, cfunc3, u)
runAnisoAdapt(m_ptr)

apf.writeVtkFiles("output_post", m_ptr)
tmp, num_Entities, m_ptr, mshape_ptr = apf.init2(dmg_name, smb_name, 1, 0)  # re-apf.initilize the same mesh


text:
module Turkie

using Makie: Makie, Figure, Scene, Point2f0
using Makie: barplot!, lines!, scatter! # Plotting tools
using Makie: Observable, Node, lift, on # Observable tools
using Makie: recordframe! # Recording tools
using Makie.MakieLayout # Layouting tool
using Colors, ColorSchemes # Colors tools
using KernelDensity # To be able to give a KDE
using OnlineStats # Estimators
using Turing: DynamicPPL.VarInfo, DynamicPPL.Model, Inference._params_to_array

using MCMCChains

export TurkieCallback
export addIO!, record

# Uses the colorblind scheme of seaborn by default
const std_colors = ColorSchemes.seaborn_colorblind

name(s::Symbol) = name(Val(s))
name(::Val{T}) where {T} = string(T)
name(s::OnlineStat) = string(nameof(typeof(s)))

include(joinpath("live_sampling", "turkie_callback.jl"))
include(joinpath("chain_plotting", "makie_plots.jl"))

end


text:
# Starting example file for Julia number types

# Julia defines a set of specific sub-types for numbers
# Signed integers: Int, Int8, Int16, Int32, Int64, BigInt
a = Int16(1)
println(typeof(a))

b = Int(2000)
println(typeof(b))

# Unsigned integers: UInt, UInt8, UInt16, UInt32, UInt64
c = UInt16(1)
println(typeof(c))

d = UInt(2000)
println(typeof(d))

# typemax() and typemin() will provide max and min values
println(typemax(Int8))
println(typemax(Int32))
println(typemax(Int64))
println(typemin(Int64))

# Use the WORDSIZE property to see what type of system this is
println(Sys.WORD_SIZE)

# Trying to assign a number too large for the type will fail and give an error
# CANNOT do: a = UInt(300)

# Special values represent Infinity and not-a-number
println(1.0 / 0.0)
println(1.0 / Inf)
println(0.0 / 0.0)

# zero() and one() functions produce values for a given type
println(zero(Float64))
println(one(UInt32))

max_stars_repo_name: scaotravis/Julia-Course

max_issues_repo_name: scaotravis/Julia-Course

max_forks_repo_name: scaotravis/Julia-Course


text:
__precompile__()

module MultivectorAlgebras

import Base: +, -, *, /, \, real, isreal, conj, inv, abs, abs2, show, zero, one, iszero

"Error when finding the inverse of zero."
const ZeroInverse = "inverse of zero"

"Error when finding the inverse of a zero divisor."
const ZeroDivisorInverse = "inverse of zero divisor"

"Error when the denominator in a quotient is zero."
const ZeroDenominator = "denominator is zero"

"Error when the denominator in a quotient is a zero divisor."
const ZeroDivisorDenominator = "denominator is zero divisor"

include("AbstractMultivector.jl")
include("Multivector1.jl")
include("Multivector2.jl")
include("Multivector3.jl")
include("Multivector4.jl")

export AbstractMultivector, asarray, iszerodivisor, commutator, associator, jacobiator, crossratioL, crossratioR, möbiusL, möbiusR
export Multivector1, unreal, cloak, dagger, star, selfstar, antiselfstar, random
export Multivector2
export Multivector3
export Multivector4

end


text:
# TODO: Type-stability of HDF5 functions should be checked.
# TODO: Implement proper tests.
# TODO: the passthrough function can be more efficient.
# TODO: Skip iter in collectdata it is `nothing` in the generated function.
# TODO: Writers close their remote channel to communicate that they are finished. Should
#       implement automatic reopening of the channel when the Writer is re-used.
# TODO: Check that HDF5Writer has the correct supported datatypes.

module AgentBased

using HDF5
using DataFrames

include("utils.jl")
include("base.jl")
include("writer.jl")
include("writer_hdf5.jl")
include("writer_dataframe.jl")
include("runbatch.jl")

export Reporter, Collector, Writer, HDF5Writer, DataFrameWriter, runbatch, collectdata,
       flushdata

end # module AgentBased


text:
# Packages required for testing
using Test
using Random
using LinearAlgebra
using Logging
using UUIDs
using SatelliteDynamics

# Package Under Test
using SatelliteTasking
using SatelliteTasking.SatellitePlanning
using SatelliteTasking.Analysis

# Set logging level
global_logger(SimpleLogger(stderr, Logging.Debug))

# Fix randomness during tests
Random.seed!(0)

# Check equality of two arrays
@inline function array_isapprox(x::AbstractArray{F},
                  y::AbstractArray{F};
                  rtol::F=sqrt(eps(F)),
                  atol::F=zero(F)) where {F<:AbstractFloat}

    # Easy check on matching size
    if length(x) != length(y)
        return false
    end

    for (a,b) in zip(x,y)
        @test isapprox(a,b, rtol=rtol, atol=atol)
    end
end

# Check if array equals a single value
@inline function array_isapprox(x::AbstractArray{F},
                  y::F;
                  rtol::F=sqrt(eps(F)),
                  atol::F=zero(F)) where {F<:AbstractFloat}

    for a in x
        @test isapprox(a, y, rtol=rtol, atol=atol)
    end
end

@inline function array_isapprox(x::AbstractArray{F}, y::AbstractArray{F}) where {F<:Integer}
    # Easy check on matching size
    if length(x) != length(y)
        return false
    end

    for (a,b) in zip(x,y)
        @test a == b
    end
end

# Check if array equals a single value
@inline function array_isapprox(x::AbstractArray{F}, y::F) where {F<:Integer}
    for a in x
        @test a == y
    end
end

@time @testset "SatelliteTasking Package Tests" begin
    testdir = joinpath(dirname(@__DIR__), "test")

end


text:
function solve()
    x, y = [parse(Int, x) for x in split(readline())]
    y % 2 == 0 && 2x ≤ y ≤ 4x ? "Yes" : "No"
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using Test
using RegionTrees
using StaticArrays
using NBInclude
using LinearAlgebra

include("hyperrectangle.jl")
include("cell.jl")
include("twosarray.jl")
include("asdfs.jl")

# const notebooks_dir = joinpath(dirname(@__FILE__), "..", "examples")
# @nbinclude(joinpath(notebooks_dir, "demo.ipynb"))
# @nbinclude(joinpath(notebooks_dir, "adaptive_distances.ipynb"))


text:
module BadukGoWeiqiTools

include("scraping.jl")
include("namesdb.jl")
include("headtohead.jl")

end


text:
using CSV
using YAML
using Colors
using DataFrames

"""
    palette_csv2yaml(palette_csv_path::AbstractString, palette_yaml_path::AbstractString)

This function makes a yaml of generator categorization based on a csv input.

# Arguments
- `palette_csv_path::AbstractString`: the csv filepath
- `palette_yaml_path::AbstractString`: the yaml save path

# Example
palette_csv2yaml("deps/NREL-palette.csv", "report_templates/color-palette.yaml")

"""
function palette_csv2yaml(
    palette_csv_path::AbstractString,
    palette_yaml_path::AbstractString,
)
    palette = DataFrame(CSV.read(palette_csv_path))

    if !("Over generation" in palette.Technology)
        push!(
            palette,
            Dict(
                Symbol("Technology") => "Over generation",
                Symbol("R:") => 12,
                Symbol("G:") => 244,
                Symbol("B:") => 252,
                Symbol("HEX") => "0CF4FC",
            ),
        )
    end

    palette_yaml = Dict()
    for row in eachrow(palette)
        RGB = "rgba($(row[Symbol("R:")]), $(row[Symbol("G:")]), $(row[Symbol("B:")]), 1)"
        #color = Colors.RGBA(row[Symbol("R:")]/255, row[Symbol("B:")]/255, row[Symbol("G:")]/255, 1)
        palette_yaml[row.Technology] = Dict("RGB" => RGB, "order" => getfield(row, :row))
    end

    YAML.write_file(palette_yaml_path, palette_yaml)
end


text:
module GuidedProposals

    using DiffEqCallbacks, OrdinaryDiffEq, LinearAlgebra
    using StaticArrays, Random, RecursiveArrayTools, Trajectories
    using DiffusionDefinition
    using ObservationSchemes
    const DD = DiffusionDefinition
    const OBS = ObservationSchemes
    import DiffusionDefinition: set_parameters!, same_entries

    include("utility_functions.jl")
    include("containers.jl")
    include("best_containers.jl")
    include("ode_solver_general.jl")
    include("guided_proposals.jl")
    include("log_likelihood.jl")
    include("sampling.jl")
    include("reparameterizations.jl")
    include("ode_solver_hfc.jl")
    include("ode_solver_mlmu.jl")
    include("ode_solver_pnu.jl")

    export GuidProp, H, F, c, recompute_guiding_term!, build_guid_prop
    export loglikhd, loglikhd_obs
    export forward_guide!, backward_filter!
    export guid_prop_for_blocking, set_obs!
    export standard_guid_prop_time_transf
end


text:
function help(v)
    i = 1
    while true
        if v[i] == 1
            v[v[i + 3] + 1] = v[v[i + 1] + 1] + v[v[i + 2] + 1]
        elseif v[i] == 2
            v[v[i + 3] + 1] = v[v[i + 1] + 1] * v[v[i + 2] + 1]
        else return v[1]
        end
        i += 4
    end
    v[1]
end

function parsev()
    v = [1,0,0,3,1,1,2,3,1,3,4,3,1,5,0,3,2,1,9,19,1,10,19,23,2,9,23,27,1,6,27,31,2,31,9,35,1,5,35,39,1,10,39,43,1,10,43,47,2,13,47,51,1,10,51,55,2,55,10,59,1,9,59,63,2,6,63,67,1,5,67,71,1,71,5,75,1,5,75,79,2,79,13,83,1,83,5,87,2,6,87,91,1,5,91,95,1,95,9,99,1,99,6,103,1,103,13,107,1,107,5,111,2,111,13,115,1,115,6,119,1,6,119,123,2,123,13,127,1,10,127,131,1,131,2,135,1,135,5,0,99,2,14,0,0]

    for i in 0:99, j in 0:99
        z = deepcopy(v)
        z[2] = i
        z[3] = j
        if help(z) == 19690720
            return 100 * i + j
        end
    end
end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
function allergic_to(score, allergen)

end

function allergy_list(score)

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    NeighborSearchMethod

A method for searching neighbors in a spatial object given a reference point.
"""
abstract type NeighborSearchMethod end

"""
    object(method)

Return the spatial object containing all possible neighbors.
"""
object(method::NeighborSearchMethod) = method.object

"""
    search(xₒ, method, mask=nothing)

Return neighbors of coordinates `xₒ` using `method` and a `mask` over
the spatial object.
"""
function search end

"""
    search(ind, method, mask=nothing)

Return neighbors of index `ind` in spatial object using `method` and a `mask`.
"""
search(ind::Int, method::NeighborSearchMethod; mask=nothing) =
  search(coordinates(object(method), ind), method; mask=mask)

"""
    BoundedNeighborSearchMethod

A method for searching neighbors with the property that the number of neighbors
is bounded above by a known constant (e.g. k-nearest neighbors).
"""
abstract type BoundedNeighborSearchMethod <: NeighborSearchMethod end

"""
    maxneighbors(method)

Return the maximum number of neighbors obtained with `method`.
"""
function maxneighbors end

"""
    search!(neighbors, xₒ, method, mask)

Update `neighbors` of coordinates `xₒ` using `method` and `mask`,
and return number of neighbors found.
"""
function search! end

"""
    search!(neighbors, ind, method, mask)

Update `neighbors` of index `ind` in spatial object using `method` and `mask`,
and return number of neighbors found.
"""
search!(neighbors, ind::Int, method::BoundedNeighborSearchMethod; mask=nothing) =
  search!(neighbors, coordinates(object(method), ind), method; mask=mask)

function search(xₒ::AbstractVector, method::BoundedNeighborSearchMethod; mask=nothing)
  neighbors = Vector{Int}(undef, maxneighbors(method))
  nneigh = search!(neighbors, xₒ, method; mask=mask)
  view(neighbors, 1:nneigh)
end

#------------------
# IMPLEMENTATIONS
#------------------
include("neighborsearch/neighborhood.jl")
include("neighborsearch/knearest.jl")
include("neighborsearch/kball.jl")
include("neighborsearch/bounded.jl")