text:
using Test, Random, InteractiveUtils, SparseArrays, LinearAlgebra, StatsBase, Distributions

using BayesianNonparametricStatistics

Base.Libc.srand(123)
e = Base.MathConstants.e
sqrttwo = sqrt(2.0)

@testset "BayesianNonparametricStatistics package" begin
    @testset "model.jl" begin
        @test supertype(SDEModel) == AbstractModel
        @test supertype(AbstractModel) == Any
        for σ in -2.0:0.1:2.0
            for beginvalue in -2.0:0.1:2.0
                for endtime in 0.1:0.1:2.0
                    for Δ in 0.01:0.01:0.99
                        M = SDEModel(σ, beginvalue, endtime, Δ)
                        @test M.σ == σ
                        @test M.beginvalue == beginvalue
                        @test M.endtime == endtime
                        @test M.Δ==Δ
                    end
                end
            end
        end
        σ = identity
        M = SDEModel(σ, 0.0,10.0, 0.01)
        @test M.σ == identity
        σ(x) = 1.0
        M = SDEModel(σ, 0.0,10.0,0.01)
        @test M.σ == σ

        σ2(x) = 2+sin(x)
        M = SDEModel(σ2, 1.3, 101.0, 0.5)
        @test M.σ == σ2
        @test M.beginvalue == 1.3
        @test M.endtime == 101.0
        @test M.Δ == 0.5

        @test_throws ArgumentError SDEModel(sin, 0.0, -1.0, 0.1)
        @test_throws ArgumentError SDEModel(sin, 0.0, 1.0, -0.1)
        @test_throws ArgumentError SDEModel(sin, 0.0, -1.0, -0.1)
        @test_throws ArgumentError SDEModel(1.0, 0.0, 0.0, 0.1)
        @test_throws ArgumentError SDEModel(1.0, 0.0, 10.0, 0.0)
        @test_throws ArgumentError SDEModel(1.0, 0.0, 0.0, 0.0)
    end # perfect

    #Test functions are correct, complete and the same in both versions.
    @testset "samplepath.jl" begin
        @test supertype(SamplePath) == AbstractSamplePath

        @test supertype(AbstractSamplePath) == Any

        @test BayesianNonparametricStatistics.isstrictlyincreasing(0:0.1:1)
        @test !BayesianNonparametricStatistics.isstrictlyincreasing(0:-0.1:-1.0)
        @test !BayesianNonparametricStatistics.isstrictlyincreasing([0.0, 0.1, 0.2, 0.1, 0.3, 0.4])
        @test !BayesianNonparametricStatistics.isstrictlyincreasing([0.0, -0.1, 0.0, 0.2, 0.3])
        @test !BayesianNonparametricStatistics.isstrictlyincreasing([1, 2, 3, 4, -10, 11, 12])
        @test BayesianNonparametricStatistics.isstrictlyincreasing([1, 2, 3, 4, 500, 600, 2000, 10000])
        @test BayesianNonparametricStatistics.isstrictlyincreasing([1.0, 2.0, 2.00000001, 2.00000002])
        @test !BayesianNonparametricStatistics.isstrictlyincreasing([2.0,2.0,2.0,2.0,2.0])
        @test !BayesianNonparametricStatistics.isstrictlyincreasing([2.0,2.0,2.1,2.2,2.3])
        @test !BayesianNonparametricStatistics.isstrictlyincreasing([2.0,2.1,2.2,2.2,2.3])
        @test !BayesianNonparametricStatistics.isstrictlyincreasing([2.1,2.2,2.3,2.4,2.4])
        @test !BayesianNonparametricStatistics.isstrictlyincreasing([2.0,2.0, 3.0, 3.0001])

        X = SamplePath(0.0:0.1:1.0, sin)
        @test X.timeinterval == 0.0:0.1:1.0
        @test X.samplevalues == sin.(0.0:0.1:1.0)

        X = SamplePath(0.0:0.1:1.0, 1.0:11.0)

        @test X.timeinterval == 0.0:0.1:1.0
        @test X.samplevalues == 1.0:11.0

        X = SamplePath(0.0:0.1:1.0, x->2*x)
        @test X.timeinterval == 0.0:0.1:1.0
        @test X.samplevalues == 0.0:0.2:2.0

        @test step(X) == 0.1
        @test length(X) == 11

        X = SamplePath([1., 2.], [1.1, 2.1])
        @test X.timeinterval == [1., 2.]
        @test X.samplevalues == [1.1, 2.1]

        @test length(X) == 2

        X = SamplePath([1.0, 3.0, 5.0, 9.0], x->x^2)
        @test X.timeinterval == [1.0, 3.0, 5.0, 9.0]
        @test X.samplevalues == [1.0, 9.0, 25.0, 81.0]

        X = SamplePath(collect(-2.0:0.01:2.0), sin)
        X.timeinterval == -2.0:0.01:2.0
        X.samplevalues == sin.(-2.0:0.01:2.0)

        X = SamplePath(-2.0:0.01:2.0, sin)
        X.timeinterval == -2.0:0.01:2.0
        X.samplevalues == sin.(-2.0:0.01:2.0)

        X = SamplePath(0.0:0.1:2π, sin)
        @test length(X) == length(0.0:0.1:2π)

        @test_throws DimensionMismatch SamplePath(0.0:2.0, [1.0, 2.0])
        @test_throws DimensionMismatch SamplePath(0.0:2.0, [1.0, 2.0, 4.0, 5.0])

        @test_throws ArgumentError SamplePath(0.0:-0.1:-0.1, [0.0, 2.0])

        @test_throws DimensionMismatch SamplePath([0.1, 2.0],[3.0])
        @test_throws DimensionMismatch SamplePath([0.1, 2.0],Float64[])
        @test_throws DimensionMismatch SamplePath(Float64[],[3.0])

        @test_throws ArgumentError SamplePath([0.1, 2.0, 1.0], [1.0,2.0,3.0])
        @test_throws ArgumentError SamplePath(0.0:-1.0:-10.0, rand(11))

        @test step(SamplePath(0.0:0.2:1.0, sin)) == 0.2
        @test step(SamplePath(0.0:10.0:1000.0, 0.0:100.0)) == 10.0

        for i in 1:10
            @test length(SamplePath(1.0:i, rand(i))) == i
        end

        @test length(SamplePath([1.0, 2.0, 4.0, 5.0, 16.0, 30.0, 34.0, 100.0], x->x^2)) == 8
        @test length(SamplePath([1.0, 2.0, 4.0, 5.0, 16.0, 30.0, 34.0, 100.0], [90.0, 56.0, 49.0, 25.0, 6.0, 10.0, 19.0, 1000000.0])) == 8

        t = 0.0:0.1:2π
        X = SamplePath(t, sin)
        @test [value for value in X] == X.samplevalues
        @test try [value for value in X]
            true
        catch
            false
        end
        @test sum(X) ≈ sum(X.samplevalues)
        @test [value for value in Iterators.Reverse(X)] == reverse(X.samplevalues)
        @test sum(Iterators.Reverse(X)) ≈ sum(X)
        @test maximum(X) == maximum(X.samplevalues)
        @test minimum(X) == minimum(X.samplevalues)
        @test extrema(X) == extrema(X.samplevalues)

        @test firstindex(X) == 1
        @test lastindex(X) == length(X.samplevalues)

        Y = SamplePath([0.1, 0.2], [1.0, 2.0])
        @test eltype(typeof(X)) == Float64
        @test eltype(typeof(Y)) == Float64

        @test X[begin] == sin(0.0)
        @test X[end] == sin(t[end])

        X = SamplePath(0:100.0, 0:100.0)
        @test sum(X) == 5050
        @test X[begin] == 0.0
        @test X[end] == 100.0
        @test [value for value in Iterators.Reverse(X)][1] == 100.0
    end # Perfect

    #Test functions are correct, complete and the same in both versions.
    @testset "basisfunctions.jl" begin
        numberofbasisfunctionstested = 200
        numberoffaberschauderlevelstested = convert(Int64, round(log2(numberofbasisfunctionstested), RoundUp)) - 1
        numberoffourierfunctionstested = numberofbasisfunctionstested

        # Exceptions

        @test_throws AssertionError fourier(-1)
        @test_throws AssertionError fourier(-10)

        #The Fourier functions are orthogonal.
        t = 0.0:1/(4*ceil(numberoffourierfunctionstested/2)):1.0
        tprecise = 0.0:10.0^-1/(4*ceil(numberoffourierfunctionstested/2)):1.0
        for k in 0:numberoffourierfunctionstested-1
            for l in k+1:numberoffourierfunctionstested
                phikvalues = fourier(k).(tprecise[1:end-1])
                philvalues = fourier(l).(tprecise[1:end-1])
                @test abs(sum(step(tprecise)*phikvalues.*philvalues)) < 10.0^-10
            end
        end
        # And the square integrates to one.
        for k in 0:numberoffourierfunctionstested
            values = fourier(k).(tprecise)
            @test abs(step(tprecise)*sum(values.*values)-1.0) < 10.0^-2
        end

        # They are one-periodic.
        for k in 0:numberoffourierfunctionstested
            for y in t
                value = fourier(k)(y)
                for m in -10.0:10.0
                    valuetwo = fourier(k)(y+m)
                    @test abs(value - valuetwo) < 10.0^-10
                end
            end
        end

        # fourier(0) is constant 1.0.
        for x in -10.0:0.01:10.0
            @test fourier(0)(x) == 1.0
        end
        # # For k≥1 they have minimum -sqrt(2) and maximum sqrt(2).
        for k in 1:numberoffourierfunctionstested
            if k % 2 == 1
                t = 0.0:0.1*1/(2k+2):1.0
            else
                t = 0.0:0.1*1/k:1.0
            end
            values = fourier(k).(t)
            @test abs(minimum(values)+sqrttwo) < 10.0^-5
            @test abs(maximum(values)-sqrttwo) < 10.0^-5
        end
        # # For k≥1:
        # # Functions \phi_k with k odd are zero when x=m/(k+1), m=0,...,k.
        # # Functions \phi_k with k even are zero when x=(1+2m)/(2k), m=0,...,k-1.
        for k in 1:2:numberoffourierfunctionstested
            for m in 0:k
                @test abs(fourier(k)(m/(k+1))) < 10.0^-5
            end
        end
        #
        for k in 2:2:numberoffourierfunctionstested
            for m in 0:k-1
                @test abs(fourier(k)((1+2m)/(2k))) < 10.0^-5
            end
        end
        #
        # # When k is odd, then \phi_k(x)=\sqrt(2) when x = (4m+1)/(2k+2), m=0,..,(k-1)/2.
        for k in 1:2:numberoffourierfunctionstested
            for m in 0:convert(Int64,(k-1)/2)
                @test abs(fourier(k)((4m+1)/(2k+2))-sqrt(2)) < 10.0^-5
            end
        end
        #
        # When k is odd, then \phi_k(x)= -sqrt(2) when x = (3+4m)/(2k+2), m=0,..,(k-1)/2.
        for k in 1:2:numberoffourierfunctionstested
            for m in 0:convert(Int64, (k-1)/2)
                @test abs(fourier(k)((3+4m)/(2k+2))+sqrt(2)) < 10.0^-5
            end
        end
        #
        # When k is even, then \phi_k(x)=sqrt(2) when x = 2m/k, m=0,...,k/2-1.
        for k in 2:2:numberoffourierfunctionstested
            for m in 0:convert(Int64, k/2)-1
                @test abs(fourier(k)(2m/k)-sqrt(2)) < 10.0^-5
            end
        end
        #
        # when k is even, then \phi_k(x)=-sqrt(2) when x=(2m+1)/k, m = 0,...,k/2-1.
        for k in 2:2:numberoffourierfunctionstested
            for m in 0:convert(Int64, k/2)-1
                @test abs(fourier(k)((2m+1)/k)+sqrt(2)) < 10.0^-5
            end
        end

        # tests for Faber-Schauder functions
        #
        @test_throws AssertionError faberschauder(-1, 1)
        @test_throws AssertionError faberschauder(-3, 1)
        @test_throws AssertionError faberschauder(1, 0)
        @test_throws AssertionError faberschauder(1, -2)
        @test_throws AssertionError faberschauder(1, 3)

        for j in 0:numberoffaberschauderlevelstested
            if j>0
                @test_throws AssertionError faberschauder(-j, 1)
            end
            @test_throws AssertionError faberschauder(j, 0)
            @test_throws AssertionError faberschauder(j, -1)
            @test_throws AssertionError faberschauder(j, 2^j + 1)
            @test_throws AssertionError faberschauder(j, 2^j + 2)
        end

        # # Below we test Faber-Schauder functions on their mathematical properities, up
        # # to level numberoffaberschauderlevelstested, defined above.
        t = 0.:2.0^(-numberoffaberschauderlevelstested-2):1.0

        # Orthogonality tests

        # Faber-Schauder functions are level-wise orthogonal, for levels j≥1.
        for j in 1:numberoffaberschauderlevelstested
            for k in 1:2^j-1
                for l in k+1:2^j
                    psijkvalues = faberschauder(j,k).(t[1:end-1])
                    psijlvalues = faberschauder(j,l).(t[1:end-1])
                    @test abs(step(t)*sum(psijkvalues.*psijlvalues))<10.0^-3
                end
            end
        end

        # ψ_{j,k} is independent with ψ_{j+r, 1},..., ψ_{j+r, (k-1)*2^r} and with
        # ψ_{j+r, k*2^r+1},...,ψ_{j+r, 2^(j+r)}
        for j in 1:numberoffaberschauderlevelstested
            for k in 1:2^j
                for r in 0:(numberoffaberschauderlevelstested-j)
                    for ℓ in 1:(k-1)*2^r
                        psijkvalues = faberschauder(j,k).(t[1:end-1])
                        psijplusrlvalues = faberschauder(j+r,ℓ).(t[1:end-1])
                        @test abs(step(t)*sum(psijkvalues.*psijplusrlvalues))<10.0^-3
                    end
                    for ℓ in k*2^r + 1:2^(j+r)
                        psijkvalues = faberschauder(j,k).(t[1:end-1])
                        psijplusrlvalues = faberschauder(j+r,ℓ).(t[1:end-1])
                        @test abs(step(t)*sum(psijkvalues.*psijplusrlvalues))<10.0^-3
                    end
                end
            end
        end

        # Integral

        t = 0.:2.0^(-numberoffaberschauderlevelstested-2):1.0
        # They integrate to 2^(-j-1)
        psionevalues = faberschauderone.(t[1:end-1])
        @test abs(step(t)*sum(psionevalues)-2.0^-1) < 10.0^-3

        for j in 1:numberoffaberschauderlevelstested
            for k in 1:2^j
                psijkvalues = faberschauder(j,k).(t[1:end-1])
                @test abs(step(t)*sum(psijkvalues)-2.0^(-j-1)) < 10.0^-3
            end
        end

        # Maxima and minima.

        t = 0.:2.0^(-numberoffaberschauderlevelstested-3):1.0
        # They have global minimum 0.0 and global maximum 1.0
        @test maximum(faberschauderone.(t)) == 1.0
        @test minimum(faberschauderone.(t)) == 0.0

        # psi_1 has maximum 1 at 0 and 1.
        @test faberschauderone(0.0) == 1.0
        @test faberschauderone(1.0) == 1.0

        # psi_1 has minimum 0 at 0.5.
        @test faberschauderone(0.5) == 0.0

        # Also ψ_{j,k} has global minimum 0 and global maximum 1.
        for j in 0:numberoffaberschauderlevelstested
            for k in 1:2^j
                values = faberschauder(j,k).(t)
                @test minimum(values) == 0.0
                @test maximum(values) == 1.0
            end
        end

        # psi_{j,k} has maximum 1.0 at 2.0^(-j-1)+(k-1)*2.0^-j
        for j in 0:numberoffaberschauderlevelstested
            for k in 1:2^j
                @test faberschauder(j,k)(2.0^(-j-1)+(k-1)*2.0^-j) == 1.0
            end
        end

        #psi_{j,k} has minimum 0.0 at (k-1)*2.0^-j and k*2.0^-j
        for j in 0:numberoffaberschauderlevelstested
            for k in 1:2^j
                @test faberschauder(j,k)((k-1)*2.0^-j) == 0.0
                @test faberschauder(j,k)(k*2.0^-j) == 0.0
            end
        end

        # Values at certain places.

        #psi_{j,k} are zero for 0.0≤x≤(k-1)*2.0^-j and k*2.0^-j ≤x≤1.0
        x = 0.:2.0^(-numberoffaberschauderlevelstested-3):1.0
        for j in 1:numberoffaberschauderlevelstested
            for k in 1:2^j
                t_left = 0.:2.0^(-numberoffaberschauderlevelstested-3):(k-1)*2.0^-j
                t_right = k*2.0^-j:2.0^(-numberoffaberschauderlevelstested-3):1.0
                ψjk = faberschauder(j,k)
                @test all(ψjk.(t_left) .== 0.0)
                @test all(ψjk.(t_right) .== 0.0)
            end
        end

        # ψ_1 is 1/2 at 0.25 and 0.75
        # ψ_{j,k} is 1/2 at (k-1+1/4)*2.0^-j and at (k-1+3/4)*2^-j.
        @test faberschauderone(0.25) == 1/2
        @test faberschauderone(0.75) == 1/2
        for j in 0:numberoffaberschauderlevelstested
            for k in 1:2^j
                ψjk = faberschauder(j,k)
                @test ψjk((k-1+1/4)*2.0^-j) == 1/2
                @test ψjk((k-1+3/4)*2.0^-j) == 1/2
            end
        end

        # ψ_1 is 1/4 at 3/8 and 5/8
        # ψ_{j,k} is 1/4 at (k-1+1/8)*2.0^-j and at (k-1+7/8)*2^-j.
        @test faberschauderone(3/8) == 1/4
        @test faberschauderone(5/8) == 1/4
        for j in 0:numberoffaberschauderlevelstested
            for k in 1:2^j
                ψjk = faberschauder(j,k)
                @test ψjk((k-1+1/8)*2.0^-j) == 1/4
                @test ψjk((k-1+7/8)*2.0^-j) == 1/4
            end
        end

        # ψ_1 is 3/4 at 1/8 and 7/8
        # ψ_{j,k} is 1/4 at (k-1+3/8)*2.0^-j and at (k-1+5/8)*2^-j.
        @test faberschauderone(1/8) == 3/4
        @test faberschauderone(7/8) == 3/4
        for j in 0:numberoffaberschauderlevelstested
            for k in 1:2^j
                ψjk = faberschauder(j,k)
                @test ψjk((k-1+3/8)*2.0^-j) == 3/4
                @test ψjk((k-1+5/8)*2.0^-j) == 3/4
            end
        end

        # Periodicity

        # They are one-periodic.
        x = 0.:2.0^(-numberoffaberschauderlevelstested-3):1.0
        for y in x
            value = faberschauderone(y)
            for k in -10.0:10.0
                valuetwo = faberschauderone(y+k)
                @test value == valuetwo
            end
        end

        for j in 0:numberoffaberschauderlevelstested
            for k in 1:2^j
                for y in x
                    value = faberschauder(j,k)(y)
                    for m in -10.:10.
                        valuetwo = faberschauder(j,k)(y+m)
                        @test value == valuetwo
                    end
                end
            end
        end

        # End of tests for mathematical properities of Faber-Schauder functions.
    end # Perfect.

    #Test functions are correct, complete and the same in both versions.
    @testset "SDE.jl" begin
        @test supertype(SDE) == AbstractSDE

        @test supertype(AbstractSDE) == Any

        @test SDE in subtypes(AbstractSDE)

        model = SDEModel(1.0, 0.0, 10.0, 0.001)
        sde = SDE(sin, model)
        @test sde.b == sin
        @test sde.model.beginvalue == 0.0
        @test sde.model.σ == 1.0
        @test sde.model.endtime == 10.0
        @test sde.model.Δ == 0.001

        model = SDEModel(abs, -2.0, 1098.0, 0.02)
        sde = SDE(cos, model)
        @test sde.b == cos
        @test sde.model.beginvalue == -2.0
        @test sde.model.σ == abs
        @test sde.model.endtime == 1098.0
        @test sde.model.Δ == 0.02

        # The following should represent a Brownian motion.
        model = SDEModel(1.0, 0.0, 1.0, 0.001)
        sde = SDE(x->0.0, model)

        lengthvector = 1000
        x = Vector{Float64}(undef, lengthvector)

        for k in 1:lengthvector
            X = rand(sde)
            x[k] = X[end]
        end

        @test abs(StatsBase.mean(x)) < 0.1
        @test 0.9 < StatsBase.var(x) < 1.1

        X = rand(sde)

        @test length(X) == length(0.0:sde.model.Δ:sde.model.endtime)
        @test X.timeinterval == 0.0:sde.model.Δ:sde.model.endtime

        model = SDEModel(identity, 1.0, 1.0, 0.001)
        geometricbrownianmotion = SDE(x->0.5*x, model)

        lengthvector = 1000_000
        x = Vector{Float64}(undef, lengthvector)
        for k in 1:lengthvector
            X = rand(geometricbrownianmotion)
            x[k] = X.samplevalues[end]
        end

        @test abs(e^(0.5)-StatsBase.mean(x)) < 0.1

        X = rand(geometricbrownianmotion)

        @test length(X) == length(0.0:geometricbrownianmotion.model.Δ:geometricbrownianmotion.model.endtime)
        @test X.timeinterval ==  0.0:geometricbrownianmotion.model.Δ:geometricbrownianmotion.model.endtime

        sde = SDE(x->0.0, SDEModel(1.0, 0.0, 1.0, 0.001))

        for k in 1:lengthvector
            X = rand(sde)
            x[k] = X.samplevalues[end]
        end

        @test abs(StatsBase.mean(x)) < 0.1
        @test 0.9 < StatsBase.var(x) < 1.1

        X = rand(sde)
        @test X.timeinterval == 0.0:sde.model.Δ:sde.model.endtime
        @test length(X) == length(0.0:sde.model.Δ:sde.model.endtime)

        ornsteinuhlenbeckprocess = SDE(x-> -x, SDEModel(1.0, 2.72, 1.0, 0.001))

        for k in 1:lengthvector
            X = rand(ornsteinuhlenbeckprocess)
            x[k] = X.samplevalues[end]
        end

        @test abs(StatsBase.mean(x) - 1) < 0.1

        X = rand(ornsteinuhlenbeckprocess)

        @test X.timeinterval == 0.0:ornsteinuhlenbeckprocess.model.Δ:ornsteinuhlenbeckprocess.model.endtime
        @test length(X) == length(0.0:ornsteinuhlenbeckprocess.model.Δ:ornsteinuhlenbeckprocess.model.endtime)

        f = BayesianNonparametricStatistics.calculatenextsamplevalue
        σ = 1.0
        Δ = 1.0
        prevXval = 1.0
        b = identity
        BMincrement = 1.0
        @test f(prevXval, σ, b, Δ, BMincrement) == 3.0

        σ = 1.0
        Δ = 1.0
        prevXval = 1.0
        b = x->1.0
        BMincrement = 1.0
        @test f(prevXval, σ, b, Δ, BMincrement) == 3.0

        σ = x -> 4.0
        Δ = 3.0
        prevXval = 0.0
        b = x->3.0
        BMincrement = 4.0
        @test f(prevXval, σ, b, Δ, BMincrement) == 5.0^2

        # End tests "sde.jl".
    end # Perfect

    #Test functions are correct, complete and the same in both versions.
    @testset "gaussianvectors.jl" begin
        m = [10.0, -3.5]
        A = [1.0 3.0; 5.0 7.0]
        B = [1.0 1.0 1.0; 2.0 2.0 2.0 ; 3.0 4.0 5.0]
        C = m
        GV = GaussianVector(A)
        @test mean(GV) == GV.mean == [0.0, 0.0]
        @test var(GV) == GV.var == A * A'
        @test GaussianVector(C).mean ≈ [0.0, 0.0]
        GV = GaussianVector(m, A)
        @test mean(GV) ==  GV.mean == m
        @test var(GV) == GV.var == A * A'
        @test length(GV) == GV.length == 2
        @test GV.numberofGaussians == 2
        @test_throws DimensionMismatch GaussianVector([0.0], A)
        @test_throws DimensionMismatch GaussianVector([3.4, 2.2, 9.5], A)
        @test_throws DimensionMismatch GaussianVector([0.0], [1.0, 2.0])
        @test_throws DimensionMismatch GaussianVector([0.0, 0.0], [2.0])
        GV = GaussianVector([0.0, 0.0], [1.0, 2.0])
        @test mean(GV) == GV.mean == [0.0, 0.0]
        @test var(GV) == GV.var == [1.0, 2.0] * [1.0, 2.0]'
        @test length(GV) == GV.length == 2
        @test GV.numberofGaussians == 1
        randomvector = rand(GV)
        @test 2 * randomvector[1] ≈ randomvector[2]
        m = [1.0, 2.0]
        Σ = [1.0 0.0; 0.0 2.0]
        n = 100000
        GV = GaussianVector(m, Σ)
        nrandomdraws = Array{Float64, 2}(undef, n, 2)
        for k in 1:n
            nrandomdraws[k, :] = rand(GV)
        end
        mempirical = StatsBase.mean(nrandomdraws, dims = 1)
        vempirical = StatsBase.var(nrandomdraws, dims = 1)
        @test norm([1.0 2.0]-mempirical) < 0.1
        @test norm(vempirical-[1.0 4.0]) < 0.1
    end # Skip

    #Test functions are correct, complete and the same in both versions.
    @testset "gaussianprocess.jl" begin
        @test Set([GaussianProcess, FaberSchauderExpansionWithGaussianCoefficients]) ⊆ Set(subtypes(AbstractGaussianProcess))

        @test supertype(AbstractGaussianProcess) == Any

        @test supertype(GaussianProcess) == AbstractGaussianProcess
        @test supertype(FaberSchauderExpansionWithGaussianCoefficients) == AbstractGaussianProcess

        # Tests for GaussianProcess

        d = MvNormal([1.0,2.0, 3.0, 4.0])

        @test_throws DimensionMismatch GaussianProcess([fourier(k) for k in 1:2], d)

        @test_throws DimensionMismatch GaussianProcess([fourier(k) for k in 1:5], d)

        X = GaussianProcess([sinpi, cospi], MvNormal([1.0, 2.0]))

        @test X.basis == [sinpi, cospi]
        @test typeof(X.distribution) <: AbstractMvNormal
        @test cov(X.distribution) == [1.0 0.0; 0.0 4.0]
        @test mean(X.distribution) == [0.0,0.0]
        @test length(X) == 2
        @test length(X) == 2
        @test typeof(rand(X)) <: Function

        X = GaussianProcess([sinpi, cospi], MvNormal([1.0,1.0]))
        n = 100000
        a = sinpi(1/4)
        x = Vector{Float64}(undef, n)
        for k in 1:n
            f = rand(X)
            x[k] = f(1/4)
        end

        @test abs(StatsBase.mean(x))<0.01
        @test abs(StatsBase.var(x)-2*a^2)<0.1

        # Tests for calculateboundssupport

        f = BayesianNonparametricStatistics.calculateboundssupport

        for higestlevel in 1:10
            leftbounds, rightbounds = f(higestlevel)
            leftboundsnext, rightboundsnext = f(higestlevel+1)
            @test length(leftbounds) == 2^(higestlevel+1)
            @test length(rightbounds) == 2^(higestlevel+1)
            @test leftbounds[1] == 0.0
            @test rightbounds[1] == 1.0
            @test leftboundsnext[1:length(leftbounds)] == leftbounds
            @test rightboundsnext[1:length(rightbounds)] == rightbounds
            @test leftbounds[end-2^higestlevel+1] == 0.0
            @test rightbounds[end] == 1.0
            for k in 1:2^higestlevel
                @test leftbounds[end - 2^higestlevel + k] == (k-1)*2.0^(-higestlevel)
                @test rightbounds[end - 2^higestlevel + k] == k*2.0^(-higestlevel)
            end
        end

        # Tests for createFaberSchauderBasisUpToLevelHigestLevel

        f = BayesianNonparametricStatistics.createFaberSchauderBasisUpToLevelHigestLevel

        numberoflevelstested=10
        for i in 1:numberoflevelstested
            v = f(i)
            vnext = f(i+1)
            @test length(v) == 2^(i+1)
            @test v[1] == faberschauderone
            @test vnext[1:length(v)] == v
            @test length(unique(v)) == length(v)
            # the elements are callable.
            @test try [item(0.5) for item in v]
                true
            catch
                false
            end
            grid = 0.0:2.0^(-numberoflevelstested-3):1.0
            @test v[1].(grid) == faberschauderone.(grid)
            index = 2
            for j in 0:i
                for k in 1:2^j
                    @test v[index].(grid) == faberschauder(j,k).(grid)
                    index += 1
                end
            end
        end

        # Tests for FaberSchauderExpansionWithGaussianCoefficients type

        # HIER VERDER!!!

        x = 0.0:0.1:1.0
        y = BayesianNonparametricStatistics.sumoffunctions([sin, cos], [1., 1.]).(x)
        z = map(x-> sin(x)+cos(x), x)
        @test y ≈ z

        A = sparse(Diagonal([1.0, 1.0]))

        distribution = GaussianVector(A)

        @test_throws AssertionError FaberSchauderExpansionWithGaussianCoefficients(1,distribution)

        @test_throws AssertionError FaberSchauderExpansionWithGaussianCoefficients(Vector{Float64}(undef, 0))

        Π = FaberSchauderExpansionWithGaussianCoefficients(0,distribution)
        @test length(Π) == 2

        f = BayesianNonparametricStatistics.createvectorofstandarddeviationsfromstandarddeviationsperlevel

        standarddeviationsperlevel = 1.0:4.0

        @test f(standarddeviationsperlevel) == [1.0, 1.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, 4.0]

        @test_throws AssertionError f(Float64[])

        for k in 1:10
            vectorofstandarddeviations = f(1.0:k)
            @test length(vectorofstandarddeviations) == 2^k
            @test vectorofstandarddeviations[1] == 1.0
            @test vectorofstandarddeviations[2] == 1.0
            @test vectorofstandarddeviations[end] == k
            for ℓ in 1:k
                @test all(vectorofstandarddeviations[1+2^(ℓ-1):2^ℓ] .== ℓ)
            end
        end

        @test_throws AssertionError BayesianNonparametricStatistics.sumoffunctions([sin], [1.0,2.0])

        vectoroffunctions = [sin, x -> cos(x+π/2)]
        vectorofscalars = [1.0, 1.0]
        f = BayesianNonparametricStatistics.sumoffunctions(vectoroffunctions, vectorofscalars)
        x = 0.0:0.01:1.0
        y = f.(x)
        @test maximum(abs.(y)) < 0.001

        n = 10000
        x = Vector{Float64}(undef, n)

        for k in 1:n
            f = rand(Π)
            x[k] = f(0.5)
        end

        @test abs(StatsBase.mean(x))<0.1
        @test abs(StatsBase.var(x)-1.0) <0.1

        for k in 1:n
            f = rand(Π)
            x[k] = f(0.0)
        end

        @test abs(StatsBase.mean(x))<0.1
        @test abs(StatsBase.var(x)-1.0) <0.1

        for k in 1:n
            f = rand(Π)
            x[k] = f(1.0)
        end

        @test abs(StatsBase.mean(x))<0.1
        @test abs(StatsBase.var(x)-1.0) <0.1

        distribution = GaussianVector([1.0, -2.0], A)
        Π = FaberSchauderExpansionWithGaussianCoefficients(0,distribution)

        @test length(Π) == 2

        n = 10000
        x = Vector{Float64}(undef, n)

        for k in 1:n
            f = rand(Π)
            x[k] = f(0.5)
        end

        @test abs(StatsBase.mean(x)+2.0)<0.1
        @test abs(StatsBase.var(x)-1.0) <0.1

        for k in 1:n
            f = rand(Π)
            x[k] = f(0.0)
        end

        @test abs(StatsBase.mean(x)-1.0)<0.1
        @test abs(StatsBase.var(x)-1.0) <0.1

        for k in 1:n
            f = rand(Π)
            x[k] = f(1.0)
        end

        @test abs(StatsBase.mean(x)-1.0)<0.1
        @test abs(StatsBase.var(x)-1.0) <0.1

        Π = FaberSchauderExpansionWithGaussianCoefficients([1.0, sqrt(2)])

        n = 10000
        x = Vector{Float64}(undef, n)
        for k in 1:n
            f = rand(Π)
            x[k] = f(0.25)
        end

        @test abs(StatsBase.mean(x)) < 0.1
        @test abs(StatsBase.var(x) - 2.5) < 0.1

        A = SparseArrays.sparse(1.0*I, 4, 4)

        distribution = GaussianVector([1.0,1.0,1.0,1.0], A)

        Π = FaberSchauderExpansionWithGaussianCoefficients(1, distribution)

        n = 10000
        x = Vector{Float64}(undef, n)
        for k in 1:n
            f = rand(Π)
            x[k] = f(0.25)
        end

        @test abs(StatsBase.mean(x)-2.0) < 0.1
        @test abs(StatsBase.var(x)-1.5) < 0.1
    end #testset "gaussianprocess.jl"

    #Test functions are correct, complete and the same in both versions.
    @testset "calculateposterior.jl" begin
        # Test for two cases whether calculatedependentfsfunctions
        # calculates the right dependency structure of the Faber-Schauder basis.
        A = [
        true true true true true true true true;
        true true true true true true true true;
        true true true false true true false false;
        true true false true false false true true;
        true true true false true false false false;
        true true true false false true false false;
        true true false true false false true false;
        true true false true false false false true
        ]

        A = Symmetric(A)
        f = BayesianNonparametricStatistics.calculatedependentfaberschauderfunctions
        lengthvectors, rowindices, columnindices = f(2)

        @test lengthvectors == length(rowindices) == length(columnindices)

        B = falses(8, 8)
        for i in 1:lengthvectors
            B[rowindices[i], columnindices[i]] = true
            B[columnindices[i], rowindices[i]] = true
        end
        for i in 1:8
            B[i,i] = true
        end

        @test A == B

        A = [
        true true true true true true true true true true true true true true true true; #1
        true true true true true true true true true true true true true true true true; #2
        true true true false true true false false true true true true false false false false; #3
        true true false true false false true true false false false false true true true true; #4
        true true true false true false false false true true false false false false false false; #5
        true true true false false true false false false false true true false false false false; #6
        true true false true false false true false false false false false true true false false; #7
        true true false true false false false true false false false false false false true true; #8
        true true true false true false false false true false false false false false false false; #9
        true true true false true false false false false true false false false false false false; #10
        true true true false false true false false false false true false false false false false; #11
        true true true false false true false false false false false true false false false false; #12
        true true false true false false true false false false false false true false false false; #13
        true true false true false false true false false false false false false true false false; #14
        true true false true false false false true false false false false false false true false; #15
        true true false true false false false true false false false false false false false true #16
        ]

        A = Symmetric(A)

        B = falses(16, 16)
        lengthvectors, rowindices, columnindices = f(3)

        @test lengthvectors == length(rowindices) == length(columnindices)

        for i in 1:lengthvectors
            B[rowindices[i], columnindices[i]] = true
            B[columnindices[i], rowindices[i]] = true
        end
        for i in 1:16
            B[i,i] = true
        end

        @test A == B

        f = BayesianNonparametricStatistics.calculategirsanovmatrixelement

        @test_throws DimensionMismatch f(collect(1.0:2.0), collect(1.0:3.0),
            collect(1.0:2.0), collect(1.0:2.0))

        @test_throws DimensionMismatch f(collect(1.0:2.0), collect(1.0:3.0),
            collect(1.0:2.0), 3.0)

        @test_throws DimensionMismatch f(collect(1.0:2.0), collect(1.0:3.0),
            2.5, collect(1.0:2.0))

        @test_throws DimensionMismatch f(collect(1.0:2.0), collect(1.0:3.0),
            0.9, 1.2)

        @test_throws DimensionMismatch f(collect(1.0:3.0), collect(1.0:2.0),
            collect(1.0:2.0), collect(1.0:2.0))

        @test_throws DimensionMismatch f(collect(1.0:3.0), collect(1.0:2.0),
            collect(1.0:2.0), 3.0)

        @test_throws DimensionMismatch f(collect(1.0:3.0), collect(1.0:2.0),
            2.5, collect(1.0:2.0))

        @test_throws DimensionMismatch f(collect(1.0:3.0), collect(1.0:2.0),
            0.9, 1.2)

        @test_throws DimensionMismatch f(collect(1.0:2.0), collect(1.0:2.0),
            collect(1.0:3.0), collect(1.0:2.0))

        @test_throws DimensionMismatch f(collect(1.0:2.0), collect(1.0:2.0),
            collect(1.0:3.0), collect(1.0:3.0))

        @test_throws DimensionMismatch f(collect(1.0:2.0), collect(1.0:2.0),
            collect(1.0:3.0), 1.2)

        @test f(collect(1.0:100.0),ones(Float64,100),1.0,1.0) == 5050.0

        @test f(collect(1.0:100.0),2.0*ones(Float64,100),1.0,1.0) == 10100.0

        @test f(1.0:100.0, 1.0:100.0, 1.0:100.0, ones(Float64, 100)) == 100.0

        @test f(1.0:100.0, 1.0:100.0, 1.0:100.0, 1.0) == 100.0

        @test f(1.0:100.0, 1.0:100.0, 1.0, 1.0 ./ (1.0:100.0)) == 5050.0

        @test f(Float64[], Float64[], Float64[], Float64[]) == 0.0

        @test f(Float64[], Float64[], Float64[], 10.0) == 0.0

        @test f(Float64[], Float64[], 2.5, Float64[]) == 0.0

        @test f(Float64[], Float64[], 2.0, 3.0) == 0.0

        samplevalueindices = vcat(trues(100), falses(100))

        @test_throws BoundsError f(samplevalueindices, samplevalueindices, 1.0:100.0, 1.0:100.0,
            1.0, 1.0)

        @test_throws BoundsError f(samplevalueindices, samplevalueindices, 1.0:200.0, 1.0:100.0,
            1.0, 1.0)

        @test_throws BoundsError f(samplevalueindices, samplevalueindices, 1.0:200.0, 1.0:200.0,
            1.0:100.0, 1.0)

        @test_throws BoundsError f(samplevalueindices, samplevalueindices, 1.0:200.0, 1.0:200.0,
            1.0:200.0, 1.0:100.0)

        @test_throws BoundsError f(samplevalueindices, samplevalueindices, 1.0:300.0, 1.0:300.0,
            1.0, 1.0)

        @test_throws BoundsError f(samplevalueindices, samplevalueindices, 1.0:200.0, 1.0:300.0,
            1.0, 1.0)

        @test_throws BoundsError f(samplevalueindices, samplevalueindices, 1.0:200.0, 1.0:200.0,
            1.0:300.0, 1.0)

        @test_throws BoundsError f(samplevalueindices, samplevalueindices, 1.0:200.0, 1.0:200.0,
            1.0:200.0, 1.0:300.0)

        @test f(samplevalueindices, samplevalueindices, 1.0:200.0, ones(Float64, 200), 1.0, 1.0) ==
            5050.0

        @test f(samplevalueindices, samplevalueindices, 1.0:200.0, ones(Float64, 200),
            ones(Float64,200), 1.0) == 5050.0

        @test f(samplevalueindices, samplevalueindices, 1.0:200.0, ones(Float64, 200), 1.0,
            ones(Float64,200)) == 5050.0

        @test f(samplevalueindices, samplevalueindices, 1.0:200.0, ones(Float64, 200),
            ones(Float64, 200), ones(Float64,200)) == 5050.0

        @test f(samplevalueindices, samplevalueindices, ones(Float64, 200), 1.0:200.0,
            ones(Float64, 200), ones(Float64,200)) == 5050.0

        @test f(samplevalueindices, samplevalueindices, ones(Float64, 200), 1.0:200.0, 1.0, 1.0) ==
            5050.0

        @test f(samplevalueindices, samplevalueindices, ones(Float64, 200), 1.0:200.0,
            ones(Float64,200), 1.0) == 5050.0

        @test f(samplevalueindices, samplevalueindices, ones(Float64, 200), 1.0:200.0, 1.0,
            ones(Float64,200)) == 5050.0

        f = BayesianNonparametricStatistics.calculateΔt

        @test f(0.0:0.1:100.0) == 0.1

        @test f(0.0:0.1:-0.1) == 0.1

        @test f([0.0, 1.0, 100.0]) == [1.0, 99.0]

        f = BayesianNonparametricStatistics.calculategirsanovmatrix

        # bla.

        f = BayesianNonparametricStatistics.calculategirsanovvectorelement

        @test f(1.0:100.0, ones(Float64, 100), ones(Float64,100)) == 5050.0

        @test f(1.0:100.0, ones(Float64, 100), 1.0) == 5050.0

        @test_throws DimensionMismatch f(1.0:100.0, 1.0:50.0, 2.0)

        @test_throws DimensionMismatch f(1.0:66.0, 1.0:100.0, 2.0)

        @test_throws DimensionMismatch f(1.0:100.0, 1.0:50.0, 1.0:50.0)

        @test_throws DimensionMismatch f(1.0:50.0, 1.0:100.0, 1.0:50.0)

        @test_throws DimensionMismatch f(1.0:100.0, 1.0:100.0, 1.0:50.0)

        @test f(1.0:10.0, 1.0:10.0, 0.0) == Inf

        @test f(1.0:2.0, 1.0:2.0, [0.0, 1.0]) == Inf

        samplevalueindices = vcat(trues(100), falses(100))

        @test f(samplevalueindices, collect(1.0:200.0), collect(1.0:200.0).^2, 1.0) == sum(k^3
            for k in 1:100)

        @test f(samplevalueindices, collect(1.0:200.0), collect(1.0:200.0).^2, collect(1.0/k
            for k in 1.0:200.0)) == sum(k^5 for k in 1:100)

        @test_throws BoundsError f(samplevalueindices, 1.0:100.0, 1.0:100.0,
            1.0)

        @test_throws BoundsError f(samplevalueindices, 1.0:100.0, 1.0:100.0,
            1.0:100.0)

        f = BayesianNonparametricStatistics.calculategirsanovvector

        # lengthvector = 2
        samplevalues = [1.0, 2.0]
        ψXt = [[1.0], [1.0]]
        σXt = [1.0]
        @test f(2, samplevalues, ψXt, σXt) == [1.0, 1.0]

        # lengthvector = 2
        samplevalues = [1.0, 2.0]
        ψXt = [[1.0], [1.0]]
        σXt = 1.0
        @test f(2, samplevalues, ψXt, σXt) == [1.0, 1.0]

        samplevalueindices = [BitArray([true, false]), BitArray([false, true])]

        # lengthvector = 2
        samplevalues = [1.0, 2.0, 3.0]
        ψXt = [[1.0, 1.0], [1.0, 1.0]]
        σXt = [1.0, 1.0]
        @test f(2, samplevalueindices, samplevalues, ψXt, σXt) == [1.0, 1.0]

        # lengthvector = 2
        samplevalues = [1.0, 2.0, 3.0]
        ψXt = [[1.0, 1.0], [1.0, 1.0]]
        σXt = 1.0
        @test f(2, samplevalueindices, samplevalues, ψXt, σXt) == [1.0, 1.0]

        f = BayesianNonparametricStatistics.calculateσXt

        @test f(1.0, 1.0:10.0) == 1.0

        @test f(identity, 1.0:10.0) == collect(1.0:10.0)

        α = 0.5
        Π = FaberSchauderExpansionWithGaussianCoefficients([2^(α*j) for j in 1:5])

        model = SDEModel(1.0, 0.0, 10000.0, 0.01)
        sde = SDE(x->0.0, model)
        X = rand(sde)

        postΠ = calculateposterior(Π, X, model)

        n = 1000
        y = 0.0:0.01:1.0
        x = Array{Float64}(undef, length(y),n)
        for k in 1:n
            f = rand(postΠ)
            x[:,k] = f.(y)
        end
        @test maximum(abs.(StatsBase.mean(x, dims = 2))) < 0.1
        @test maximum(abs.(mean(postΠ).(y))) < 0.1

        Π = GaussianProcess([fourier(k) for k in 1:40], GaussianVector(sparse(Diagonal([k^(-1.0) for k in 1:40]))))

        model = SDEModel(1.0, 0.0, 10000.0, 0.01)
        sde = SDE(x->0.0, model)
        X = rand(sde)

        postΠ = calculateposterior(Π, X, model)

        n = 1000
        y = 0.0:0.01:1.0
        x = Array{Float64}(undef, length(y),n)
        for k in 1:n
            f = rand(postΠ)
            x[:,k] = f.(y)
        end

        @test maximum(abs.(StatsBase.mean(x, dims = 2))) < 0.1
        @test maximum(abs.(mean(postΠ).(y))) < 0.1
    end # testset calculateposterior.jl
end # testset.


text:
# coupling.jl

"""
Implementation of Frohlich's Hamiltonian and coupling parameter α applied to the Polaron model.

See Frohlich 1952:
https://www.jstor.org/stable/99166

"""

"""
frohlich_α(ϵ_optic::Float64, ϵ_static::Float64, freq::Float64, m_eff::Float64)

    Calculates the Frohlich alpha parameter, for a given dielectric constant, frequency (f) of phonon in Hertz, and effective mass (in units of the bare electron mass).

    See Equation 2.37.

"""
function frohlich_α(ϵ_optic, ϵ_static, freq, m_eff)

    ω = 2 * π * freq

    α = 0.5 / (4 * π * ϵ_0) *
       (1 / ϵ_optic - 1 / ϵ_static) *
       (eV^2 / (ħ * ω)) *
       sqrt(2 * m_e * m_eff * ω / ħ)

    return α
end

ext: jl


text:
function surround(string ; border = :default, padding = 0)

 ve, ho, ul, ur, dl, dr =
   border == :round ? ("\u2502","\u2500","\u256d","\u256e","\u2570","\u256f") :
   border == :bold  ? ("\u2503","\u2501","\u250F","\u2513","\u2517","\u251b") :
   border == :double? ("\u2551","\u2550","\u2554","\u2557","\u255a","\u255d") :
   border == :dotted? ("\u254e","\u254c","\u250c","\u2510","\u2514","\u2518") :
   border == :cross ? ("\u2502","\u2500","\u253c","\u253c","\u253c","\u253c") :
                      ("\u2502","\u2500","\u250c","\u2510","\u2514","\u2518")

 println(ul, ho^(length(string) + 2padding),  ur, "\n",
         ve, " "^padding, string," "^padding, ve, "\n",
         dl, ho^(length(string) + 2padding),  dr)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
Constitutive driver.
"""
function driver(stress,
                ats,
                matstat,
                matpar,
                time_history::Vector,
                strain_history::Matrix,
                stress_history::Matrix,
                strain_control::Vector{Bool};
                xtol::Real = 0.0,
                ftol::Real = 1e-5,
                iterations::Integer = 100,
                method::Symbol = :trust_region,
                show_trace::Bool = false)

    ɛ_dim, nɛ = size(strain_history)
    σ_dim, nσ = size(stress_history)
    nt = length(time_history)

    nɛ == nt || throw(ArgumentError("number of strains must be equal to number of time steps"))
    nσ == nt || throw(ArgumentError("number of stesses must be equal to number of time steps"))
    ɛ_dim == σ_dim || throw(ArgumentError("stresses and strains must have equal number of components"))

    mp = matpar
    ms = matstat
    # Stress control indices
    sc = !strain_control

    matstats = [ms]
    εs = strain_history[:,1]
    σs = stress_history[:,1]
    ɛ_sol = zeros(sum(sc))
    ms_grad = deepcopy(ms)

    i = 1
    t_prev = time_history[1]
    ɛ_sol = zeros(sum(sc))
    ∆εₙ₋₁ = zeros(sum(sc))
    ∆ε₀ = zeros(sum(sc))
    ∆ε = zeros(sum(sc))
    for t in time_history[2:end]
        i += 1
        dt = t - t_prev
        # Guess total strain
        ɛ = strain_history[:, i]
        σ = stress_history[:, i]
        ms_new = deepcopy(ms)
        if !(sum(sc) == 0)
            copy!(∆ε₀, ∆εₙ₋₁)

            function f(∆ε_red)
                dε = zeros(ɛ_dim)
                dε[sc] = ∆ε_red
                σ_res, ms_new = stress(ɛ + dε, dt, mp, ms)
                σ_res[sc] - σ[sc]
            end

            function g(∆ε_red)
                dε = zeros(ɛ_dim)
                dε[sc] = ∆ε_red
                grad = ats(ɛ + dε, dt, mp, ms_new)
                grad[sc, sc]
            end
            try
                res = nlsolve(not_in_place(f, g), ∆ε₀; xtol=xtol, ftol=ftol,
                              iterations=iterations, method=method, show_trace=show_trace)
                if !converged(res)
                    i -= 1
                    warn(NON_CONV_MESSAGE)
                    @goto ret
                end
                copy!(∆ε, res.zero)
            catch e
                isa(e, MaterialNonConvergenceError) || rethrow(e)
                i -= 1
                warn(MAT_NON_CONV_MESSAGE)
                @goto ret
            end
        end

        copy!(∆εₙ₋₁, ∆ε)

        ɛ[sc] += ∆ε

        # Update σ and material status at converged point
        try
            σ, ms = stress(ɛ, dt, mp, ms)
            push!(matstats, ms)
            append!(σs, σ)
            append!(εs, ɛ)
        catch e
            println(e)
            isa(e, MaterialNonConvergenceError) || rethrow(e)
            i -= 1
            warn(MAT_NON_CONV_MESSAGE)
            @goto ret
        end

        t_prev = t
    end # time

    @label ret
    ɛs_mat = reshape(εs, (ɛ_dim, i))
    σs_mat = reshape(σs, (σ_dim, i))
    return ɛs_mat, σs_mat, matstats
end


text:
# Send shifts to API
function set_shifts(shifts, task, organization)
    path = string(
        "/api/v1/organizations/",
        task["organization_id"],
        "/locations/",
        task["location_id"],
        "/roles/",
        task["role_id"],
        "/schedules/",
        task["schedule_id"],
        "/shifts/",
    )
    for e in keys(shifts)
        for shift in shifts[e]
            if StaffJoy.is_unassigned_shift(e)
                data = {
                    "user_id" => 0, # This line could be omitted
                    "day" => shift["day"],
                    "start" => shift["start"] + TIME_OFFSET,
                    "length" => shift["length"],
                }
            else
                data = {
                    "user_id" => e, # Flask can convert strings to int
                    "day" => shift["day"],
                    "start" => shift["start"] + TIME_OFFSET,
                    "length" => shift["length"],
                }
            end
            results = Requests.post(build_uri(config, path), json=data)
            if results.status != 201
                code = results.status
                message = results.data
                Logging.err("Unable to create shift $data - status %code - message $message")
                error()
            end
            Logging.info("Created shift - $data")
        end
    end
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
export AmplModel, AmplException,
       reset!,
       write_sol, amplmodel_finalize, varscale, lagscale, conscale,
       obj, grad, grad!,
       cons, cons!, jth_con, jth_congrad, jth_congrad!, jth_sparse_congrad,
       jac_coord, jac, jprod, jprod!, jtprod, jtprod!,
       jth_hprod, jth_hprod!, ghjvprod, ghjvprod!,
       hess_coord, hess, hprod, hprod!

# Convenience macro.
macro asl_call(func, args...)
  args = map(esc,args)
  quote
    ccall(($(esc(func)), libasl), $(args...))
  end
end

type AmplException
  msg :: String
end

macro check_ampl_model()
  esc(:(nlp.__asl == C_NULL && throw(AmplException("Uninitialized AMPL model"))))
end

type AmplModel <: AbstractNLPModel
  meta  :: NLPModelMeta;     # Problem metadata.
  __asl :: Ptr{Void};        # Pointer to internal ASL structure. Do not touch.

  counters :: Counters       # Evaluation counters
  safe :: Bool               # Always evaluate the objective before the Hessian.
  cvar :: Array{Int32}       # x[cv(i)] is complementary to the constraints cv(i)

  function AmplModel(stub :: AbstractString; safe :: Bool=false)
    asl = @compat @asl_call(:asl_init, Ptr{Void}, (Ptr{UInt8},), stub);
    asl == C_NULL && error("Error allocating ASL structure")

    minimize = @asl_call(:asl_objtype, Int32, (Ptr{Void},), asl) == 0;
    islp = @asl_call(:asl_islp, Int32, (Ptr{Void},), asl) != 0;

    nlo = @compat Int(@asl_call(:asl_nlo, Int32, (Ptr{Void},), asl));

    nvar = @compat Int(@asl_call(:asl_nvar, Int32, (Ptr{Void},), asl));
    ncon = @compat Int(@asl_call(:asl_ncon, Int32, (Ptr{Void},), asl));

    x0   = @compat unsafe_wrap(Array, @asl_call(:asl_x0,   Ptr{Float64}, (Ptr{Void},), asl),
                            (nvar,), false)
    y0   = @compat unsafe_wrap(Array, @asl_call(:asl_y0,   Ptr{Float64}, (Ptr{Void},), asl),
                            (ncon,), false)

    lvar = @compat unsafe_wrap(Array, @asl_call(:asl_lvar, Ptr{Float64}, (Ptr{Void},), asl),
                            (nvar,), false)
    uvar = @compat unsafe_wrap(Array, @asl_call(:asl_uvar, Ptr{Float64}, (Ptr{Void},), asl),
                            (nvar,), false)

    nzo = @compat Int(@asl_call(:asl_nzo, Int32, (Ptr{Void},), asl))
    nbv = @compat Int(@asl_call(:asl_nbv, Int32, (Ptr{Void},), asl))
    niv = @compat Int(@asl_call(:asl_niv, Int32, (Ptr{Void},), asl))
    nlvb = @compat Int(@asl_call(:asl_nlvb, Int32, (Ptr{Void},), asl))
    nlvo = @compat Int(@asl_call(:asl_nlvo, Int32, (Ptr{Void},), asl))
    nlvc = @compat Int(@asl_call(:asl_nlvc, Int32, (Ptr{Void},), asl))
    nlvbi = @compat Int(@asl_call(:asl_nlvbi, Int32, (Ptr{Void},), asl))
    nlvci = @compat Int(@asl_call(:asl_nlvci, Int32, (Ptr{Void},), asl))
    nlvoi = @compat Int(@asl_call(:asl_nlvoi, Int32, (Ptr{Void},), asl))
    nwv = @compat Int(@asl_call(:asl_nwv, Int32, (Ptr{Void},), asl))
#Tangi add's
    n_cc = @compat Int(@asl_call(:asl_n_cc, Int32, (Ptr{Void},), asl))
    nlcc = @compat Int(@asl_call(:asl_nlcc, Int32, (Ptr{Void},), asl))
    ndcc = @compat Int(@asl_call(:asl_ndcc, Int32, (Ptr{Void},), asl))
    nzlb = @compat Int(@asl_call(:asl_nzlb, Int32, (Ptr{Void},), asl))
#Fin Tangi add's
    lcon = @compat unsafe_wrap(Array, @asl_call(:asl_lcon, Ptr{Float64}, (Ptr{Void},), asl),
                            (ncon,), false)
    ucon = @compat unsafe_wrap(Array, @asl_call(:asl_ucon, Ptr{Float64}, (Ptr{Void},), asl),
                            (ncon,), false)
#Tangi add's
    if n_cc>0
    #not sure what to do if nzlb!=0
    cvar = @compat unsafe_wrap(Array, @asl_call(:asl_cvar, Ptr{Int32}, (Ptr{Void},), asl),
                            (ncon,), false)
    @show n_cc nlcc ndcc nzlb cvar nvar ncon
    else
    cvar=[]
    end
#Fin Tangi add's
    nlnet = @compat Int(@asl_call(:asl_lnc, Int32, (Ptr{Void},), asl))
    nnnet = @compat Int(@asl_call(:asl_nlnc, Int32, (Ptr{Void},), asl))
    nnln = @compat(Int(@asl_call(:asl_nlc,  Int32, (Ptr{Void},), asl))) - nnnet
    nlin = ncon - nnln - nnnet

    nln  = 1 : nnln
    nnet = nnln+1 : nnln+nnnet
    lnet = nnln+nnnet+1 : nnln+nnnet+nlnet
    lin  = nnln+nnnet+nlnet+1 : ncon

    nnzj = @compat Int(@asl_call(:asl_nnzj, Int32, (Ptr{Void},), asl))
    nnzh = @compat Int(@asl_call(:asl_nnzh, Int32, (Ptr{Void},), asl))

    meta = NLPModelMeta(nvar, x0=x0, lvar=lvar, uvar=uvar,
                        nlo=nlo, nnzo=nzo,
                        ncon=ncon, y0=y0, lcon=lcon, ucon=ucon,
                        nnzj=nnzj, nnzh=nnzh,
                        nbv=nbv, niv=niv,
                        nlvb=nlvb, nlvo=nlvo, nlvc=nlvc,
                        nlvbi=nlvbi, nlvci=nlvci, nlvoi=nlvoi, nwv=nwv,
                        lin=lin, nln=nln, nnet=nnet, lnet=lnet,
                        nlin=nlin, nnln=nnln, nnet=nnet, nlnet=nlnet,
                        minimize=minimize, islp=islp, name=stub)

    nlp = new(meta, asl, Counters(), safe, cvar)

    finalizer(nlp, amplmodel_finalize)
    return nlp
  end

end

# Import methods we override.
import Base.show, Base.print
import NLPModels.reset!
import NLPModels.varscale, NLPModels.lagscale, NLPModels.conscale
import NLPModels.obj, NLPModels.grad, NLPModels.grad!
import NLPModels.cons, NLPModels.cons!, NLPModels.jth_con
import NLPModels.jth_congrad, NLPModels.jth_congrad!, NLPModels.jth_sparse_congrad
import NLPModels.jac_coord, NLPModels.jac
import NLPModels.jprod, NLPModels.jtprod, NLPModels.jprod!, NLPModels.jtprod!
import NLPModels.jth_hprod, NLPModels.jth_hprod!
import NLPModels.ghjvprod, NLPModels.ghjvprod!
import NLPModels.hess_coord, NLPModels.hess, NLPModels.hprod, NLPModels.hprod!
import NLPModels.NLPtoMPB

# Methods associated to AmplModel instances.

"Reset evaluation counters in `nlp`."
function reset!(nlp :: AmplModel)
  reset!(nlp.counters)
  return nlp
end

"Write message `msg` along with primal and dual variables `x` and `y` to file."
function write_sol(nlp :: AmplModel, msg :: String, x :: Vector{Float64}, y :: Vector{Float64})
  @check_ampl_model
  length(x) == nlp.meta.nvar || error("x must have length $(nlp.meta.nvar)")
  length(y) == nlp.meta.ncon || error("y must have length $(nlp.meta.ncon)")

  @compat @asl_call(:asl_write_sol, Void,
                    (Ptr{Void}, Ptr{UInt8}, Ptr{Float64}, Ptr{Float64}),
                     nlp.__asl, msg,        x,            y)
end

function amplmodel_finalize(nlp :: AmplModel)
  if nlp.__asl == C_NULL
    return
  end
  @asl_call(:asl_finalize, Void, (Ptr{Void},), nlp.__asl)
  nlp.__asl = C_NULL
end

# Displaying AmplModel instances.

function show(io :: IO, nlp :: AmplModel)
  @check_ampl_model
  show(io, nlp.meta)
end

function print(io :: IO, nlp :: AmplModel)
  @check_ampl_model
  print(io, nlp.meta)
end

# Scaling AmplModel instances.

"Scale the vector of variables by the vector `s`."
function varscale(nlp :: AmplModel, s :: Vector{Float64})
  @check_ampl_model
  length(s) >= nlp.meta.nvar || error("s must have length at least $(nlp.meta.nvar)")

  err = Cint[0]
  @asl_call(:asl_varscale, Void, (Ptr{Void}, Ptr{Float64}, Ptr{Cint}), nlp.__asl, s, err)
  err[1] == 0 || throw(AmplException("Error while scaling variables"))
end

"""Set the scaling factor σ in the Lagrangian:
    L(x,y) = f(x) + σ ∑ yi ci(x).
"""
function lagscale(nlp :: AmplModel, σ :: Float64)
  @check_ampl_model
  err = Cint[0]
  @asl_call(:asl_lagscale, Void, (Ptr{Void}, Float64, Ptr{Cint}), nlp.__asl, σ, err)
  err[1] == 0 || throw(AmplException("Error while scaling Lagrangian"))
end

"Scale the vector of constraints by the vector `s`."
function conscale(nlp :: AmplModel, s :: Vector{Float64})
  @check_ampl_model
  length(s) >= nlp.meta.ncon || error("s must have length at least $(nlp.meta.ncon)")

  err = Cint[0]
  @asl_call(:asl_conscale, Void, (Ptr{Void}, Ptr{Float64}, Ptr{Cint}), nlp.__asl, s, err)
  err[1] == 0 || throw(AmplException("Error while scaling constraints"))
end

# Evaluating objective, constraints and derivatives.

"Evaluate the objective function of `nlp` at `x`."
function obj(nlp :: AmplModel, x :: Vector{Float64})
  @check_ampl_model
  length(x) >= nlp.meta.nvar || error("x must have length at least $(nlp.meta.nvar)")

  err = Cint[0]
  f = @asl_call(:asl_obj, Float64, (Ptr{Void}, Ptr{Float64}, Ptr{Cint}), nlp.__asl, x, err)
  nlp.counters.neval_obj += 1
  err[1] == 0 || throw(AmplException("Error while evaluating objective"))
  return f
end

"Evaluate the gradient of the objective function at `x`."
function grad(nlp :: AmplModel, x :: Vector{Float64})
  g = Array{Float64}(nlp.meta.nvar)
  return grad!(nlp, x, g)
end

"Evaluate the gradient of the objective function at `x` in place."
function grad!(nlp :: AmplModel, x :: Vector{Float64}, g :: Vector{Float64})
  @check_ampl_model
  length(x) >= nlp.meta.nvar || error("x must have length at least $(nlp.meta.nvar)")

  err = Cint[0]
  @asl_call(:asl_grad, Ptr{Float64},
            (Ptr{Void}, Ptr{Float64}, Ptr{Float64}, Ptr{Cint}),
             nlp.__asl, x,            g,            err)
  nlp.counters.neval_grad += 1
  err[1] == 0 || throw(AmplException("Error while evaluating objective gradient"))
  return g
end

"Evaluate the constraints at `x`."
function cons(nlp :: AmplModel, x :: Vector{Float64})
  c = Array{Float64}(nlp.meta.ncon)
  return cons!(nlp, x, c)
end

"Evaluate the constraints at `x` in place."
function cons!(nlp :: AmplModel, x :: Vector{Float64}, c :: Vector{Float64})
  @check_ampl_model
  length(x) >= nlp.meta.nvar || error("x must have length at least $(nlp.meta.nvar)")

  err = Cint[0]
  @asl_call(:asl_cons, Void,
            (Ptr{Void}, Ptr{Float64}, Ptr{Float64}, Ptr{Cint}),
             nlp.__asl, x,            c,            err)
  nlp.counters.neval_cons += 1
  err[1] == 0 || throw(AmplException("Error while evaluating constraints"))
  return c
end

"Evaluate the `j`-th constraint at `x`."
function jth_con(nlp :: AmplModel, x :: Vector{Float64}, j :: Int)
  @check_ampl_model
  (1 <= j <= nlp.meta.ncon)  || error("expected 0 ≤ j ≤ $(nlp.meta.ncon)")
  length(x) >= nlp.meta.nvar || error("x must have length at least $(nlp.meta.nvar)")

  err = Cint[0]
  cj = @asl_call(:asl_jcon, Float64,
                 (Ptr{Void}, Ptr{Float64}, Int32, Ptr{Cint}),
                  nlp.__asl, x,            j-1,   err)
  nlp.counters.neval_jcon += 1
  err[1] == 0 || throw(AmplException("Error while evaluating $j-th constraint"))
  return cj
end

"Evaluate the `j`-th constraint gradient at `x`."
function jth_congrad(nlp :: AmplModel, x :: Vector{Float64}, j :: Int)
  g = Array{Float64}(nlp.meta.nvar)
  return jth_congrad!(nlp, x, j, g)
end

"Evaluate the `j`-th constraint gradient at `x` in place."
function jth_congrad!(nlp :: AmplModel, x :: Vector{Float64}, j :: Int, g :: Vector{Float64})
  @check_ampl_model
  (1 <= j <= nlp.meta.ncon)  || error("expected 0 ≤ j ≤ $(nlp.meta.ncon)")
  length(x) >= nlp.meta.nvar || error("x must have length at least $(nlp.meta.nvar)")

  err = Cint[0]
  @asl_call(:asl_jcongrad, Ptr{Float64},
            (Ptr{Void}, Ptr{Float64}, Ptr{Float64}, Int32, Ptr{Cint}),
             nlp.__asl, x,            g,            j-1,   err)
  nlp.counters.neval_jgrad += 1
  err[1] == 0 || throw(AmplException("Error while evaluating $j-th constraint gradient"))
  return g
end

"Evaluate the `j`-th constraint sparse gradient at `x`."
function jth_sparse_congrad(nlp :: AmplModel, x :: Vector{Float64}, j :: Int)
  @check_ampl_model
  (1 <= j <= nlp.meta.ncon)  || error("expected 0 ≤ j ≤ $(nlp.meta.ncon)")
  length(x) >= nlp.meta.nvar || error("x must have length at least $(nlp.meta.nvar)")

  nnz = @asl_call(:asl_sparse_congrad_nnz, Csize_t,
                  (Ptr{Void}, Cint), nlp.__asl, j-1)

  err = Cint[0]
  inds = Array{Int64}(nnz)
  vals = Array{Float64}(nnz)
  @asl_call(:asl_sparse_congrad, Void,
            (Ptr{Void}, Ptr{Float64}, Int32, Ptr{Int64}, Ptr{Float64}, Ptr{Cint}),
             nlp.__asl, x,            j-1,   inds,       vals,         err)
  nlp.counters.neval_jgrad += 1
  err[1] == 0 || throw(AmplException("Error while evaluating $j-th sparse constraint gradient"))
  # Use 1-based indexing.
  return sparsevec(inds+1, vals, nlp.meta.nvar)
end

"Evaluate the constraints Jacobian at `x` in sparse coordinate format."
function jac_coord(nlp :: AmplModel, x :: Vector{Float64})
  @check_ampl_model
  length(x) >= nlp.meta.nvar || error("x must have length at least $(nlp.meta.nvar)")

  err = Cint[0]
  rows = Array{Int64}(nlp.meta.nnzj)
  cols = Array{Int64}(nlp.meta.nnzj)
  vals = Array{Float64}(nlp.meta.nnzj)
  @asl_call(:asl_jac, Void,
            (Ptr{Void}, Ptr{Float64}, Ptr{Int64}, Ptr{Int64}, Ptr{Float64}, Ptr{Cint}),
             nlp.__asl, x,            rows,       cols,       vals,         err)
  nlp.counters.neval_jac += 1
  err[1] == 0 || throw(AmplException("Error while evaluating constraints Jacobian"))
  # Use 1-based indexing.
  return (rows+1, cols+1, vals)
end

"Evaluate the constraints Jacobian at `x` as a sparse matrix."
function jac(nlp :: AmplModel, x :: Vector{Float64})
  @check_ampl_model
  (rows, cols, vals) = jac_coord(nlp, x)
  return sparse(rows, cols, vals, nlp.meta.ncon, nlp.meta.nvar)
end

"""
Evaluate the Jacobian-vector product at `x`.
Warning: Currently building the Jacobian for this.
"""
function jprod(nlp :: AmplModel, x :: Vector{Float64}, v :: Vector{Float64})
  Jv = zeros(nlp.meta.ncon)
  return jprod!(nlp, x, v, Jv)
end

"""
Evaluate the Jacobian-vector product at `x` in place.
Warning: Currently building the Jacobian for this.
"""
function jprod!(nlp :: AmplModel,
                x :: Vector{Float64},
                v :: Vector{Float64},
                Jv :: Vector{Float64})
  nlp.counters.neval_jac -= 1
  nlp.counters.neval_jprod += 1
  Jv[1:nlp.meta.ncon] = jac(nlp, x) * v
  return Jv
end

"""
Evaluate the transposed-Jacobian-vector product at `x`.
Warning: Currently building the Jacobian for this.
"""
function jtprod(nlp :: AmplModel, x :: Vector{Float64}, v :: Vector{Float64})
  Jtv = zeros(nlp.meta.nvar)
  return jtprod!(nlp, x, v, Jtv)
end

"""
Evaluate the transposed-Jacobian-vector product at `x` in place.
Warning: Currently building the Jacobian for this.
"""
function jtprod!(nlp :: AmplModel,
                 x :: Vector{Float64},
                 v :: Vector{Float64},
                 Jtv :: Vector{Float64})
  nlp.counters.neval_jac -= 1
  nlp.counters.neval_jtprod += 1
  Jtv[1:nlp.meta.nvar] = jac(nlp, x)' * v
  return Jtv
end

"Evaluate the product of the Lagrangian Hessian at `(x,y)` with the vector `v`."
function hprod(nlp :: AmplModel,
               x :: Vector{Float64},
               v :: Vector{Float64};
               y :: Vector{Float64} = nlp.meta.y0,
               obj_weight :: Float64 = 1.0)
  hv = Array{Float64}(nlp.meta.nvar);
  return hprod!(nlp, x, v, hv, y=y, obj_weight=obj_weight)
end

"Evaluate the product of the Lagrangian Hessian at `(x,y)` with the vector `v` in place."
function hprod!(nlp :: AmplModel,
                x :: Vector{Float64},
                v :: Vector{Float64},
                hv :: Vector{Float64};
                y :: Vector{Float64} = nlp.meta.y0,
                obj_weight :: Float64 = 1.0)
  # Note: x is in fact not used in hprod.
  @check_ampl_model
  length(x) >= nlp.meta.nvar || error("x must have length at least $(nlp.meta.nvar)")
  length(y) >= nlp.meta.ncon || error("y must have length at least $(nlp.meta.ncon)")
  length(v) >= nlp.meta.nvar || error("v must have length at least $(nlp.meta.nvar)")

  if nlp.safe
    _ = obj(nlp, x) ; nlp.counters.neval_obj -= 1
    _ = cons(nlp, x) ; nlp.counters.neval_cons -= 1
  end
  @asl_call(:asl_hprod, Ptr{Float64},
            (Ptr{Void}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Float64),
             nlp.__asl, y,            v,            hv,           obj_weight);
  nlp.counters.neval_hprod += 1
  return hv
end

"""Evaluate the product of the `j`-th constraint Hessian at `x` with the vector `v`.
The objective Hessian is used if `j=0`.
"""
function jth_hprod(nlp :: AmplModel,
                   x :: Vector{Float64}, v :: Vector{Float64}, j :: Int)
  hv = Array{Float64}(nlp.meta.nvar)
  return jth_hprod!(nlp, x, v, j, hv)
end

"""Evaluate the product of the `j`-th constraint Hessian at `x` with the vector `v` in place.
The objective Hessian is used if `j=0`.
"""
function jth_hprod!(nlp :: AmplModel,
                    x :: Vector{Float64}, v :: Vector{Float64},
                    j :: Int, hv :: Vector{Float64})
# Note: x is in fact not used in hprod.
  @check_ampl_model
  length(x) >= nlp.meta.nvar || error("x must have length at least $(nlp.meta.nvar)")
  length(v) >= nlp.meta.nvar || error("v must have length at least $(nlp.meta.nvar)")
  (1 <= j <= nlp.meta.ncon)  || error("expected 0 ≤ j ≤ $(nlp.meta.ncon)")

  if nlp.safe
    if j == 0
      _ = obj(nlp, x) ; nlp.counters.neval_obj -= 1
    else
      _ = cons(nlp, x) ; nlp.counters.neval_cons -= 1
    end
  end
  @asl_call(:asl_hvcompd, Ptr{Float64},
            (Ptr{Void}, Ptr{Float64}, Ptr{Float64}, Int),
             nlp.__asl, v,            hv,           j-1);
  nlp.counters.neval_jhprod += 1
  return hv
end

"""Compute the vector of dot products `(g, Hj*v)`
where `Hj` is the Hessian of the `j`-th constraint at `x`.
"""
function ghjvprod(nlp :: AmplModel,
                  x :: Vector{Float64}, g :: Vector{Float64}, v :: Vector{Float64})
  gHv = Array{Float64}(nlp.meta.ncon);
  return ghjvprod!(nlp, x, g, v, gHv)
end

"""Compute the vector of dot products `(g, Hj*v)` in place
where `Hj` is the Hessian of the `j`-th constraint at `x`.
"""
function ghjvprod!(nlp :: AmplModel,
                   x :: Vector{Float64}, g :: Vector{Float64},
                   v :: Vector{Float64}, gHv :: Vector{Float64})
  # Note: x is in fact not used.
  @check_ampl_model
  length(x) >= nlp.meta.nvar || error("x must have length at least $(nlp.meta.nvar)")
  length(g) >= nlp.meta.nvar || error("g must have length at least $(nlp.meta.nvar)")
  length(v) >= nlp.meta.nvar || error("v must have length at least $(nlp.meta.nvar)")

  if nlp.safe
    _ = cons(nlp, x) ; nlp.counters.neval_cons -= 1
  end
  @asl_call(:asl_ghjvprod, Ptr{Float64},
            (Ptr{Void}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}),
             nlp.__asl, g,            v,            gHv);
  nlp.counters.neval_hprod += nlp.meta.ncon
end

"""Evaluate the Lagrangian Hessian at `(x,y)` in sparse coordinate format.
Only the lower triangle is returned.
"""
function hess_coord(nlp :: AmplModel,
                    x :: Vector{Float64};
                    y :: Vector{Float64} = nlp.meta.y0,
                    obj_weight :: Float64 = 1.0)
  # Note: x is in fact not used.
  @check_ampl_model
  length(x) >= nlp.meta.nvar || error("x must have length at least $(nlp.meta.nvar)")
  length(y) >= nlp.meta.ncon || error("y must have length at least $(nlp.meta.ncon)")

  if nlp.safe
    _ = obj(nlp, x) ; nlp.counters.neval_obj -= 1
    _ = cons(nlp, x) ; nlp.counters.neval_cons -= 1
  end

  rows = Array{Int64}(nlp.meta.nnzh)
  cols = Array{Int64}(nlp.meta.nnzh)
  vals = Array{Float64}(nlp.meta.nnzh)
  @asl_call(:asl_hess, Void,
            (Ptr{Void}, Ptr{Float64}, Float64,    Ptr{Int64}, Ptr{Int64}, Ptr{Float64}),
             nlp.__asl, y,            obj_weight, rows,       cols,       vals)
  nlp.counters.neval_hess += 1
  # Use 1-based indexing.
  # Swap rows and cols to obtain the lower triangle.
  return (cols+1, rows+1, vals)
end

"""Evaluate the Lagrangian Hessian at `(x,y)` as a sparse matrix.
Only the lower triangle is returned.
"""
function hess(nlp :: AmplModel,
              x :: Vector{Float64};
              y :: Vector{Float64} = nlp.meta.y0,
              obj_weight :: Float64 = 1.0)
  @check_ampl_model
  (rows, cols, vals) = hess_coord(nlp, x, y=y, obj_weight=obj_weight);
  return sparse(rows, cols, vals, nlp.meta.nvar, nlp.meta.nvar)
end


text:
# ---
# title: 1386. Cinema Seat Allocation
# id: problem1386
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Greedy
# link: <https://leetcode.com/problems/cinema-seat-allocation/description/>
# hidden: true
# ---
#
# ![](https://assets.leetcode.com/uploads/2020/02/14/cinema_seats_1.png)
#
# A cinema has `n` rows of seats, numbered from 1 to `n` and there are ten seats
# in each row, labelled from 1 to 10 as shown in the figure above.
#
# Given the array `reservedSeats` containing the numbers of seats already
# reserved, for example, `reservedSeats[i] = [3,8]` means the seat located in
# row **3** and labelled with **8**  is already reserved.
#
# _Return the maximum number of four-person groups  you can assign on the cinema
# seats._ A four-person group occupies four adjacent seats **in one single
# row**. Seats across an aisle (such as [3,3] and [3,4]) are not considered to
# be adjacent, but there is an exceptional case on which an aisle split a four-
# person group, in that case, the aisle split a four-person group in the middle,
# which means to have two people on each side.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/02/14/cinema_seats_3.png)
#
#
#
#     Input: n = 3, reservedSeats = [[1,2],[1,3],[1,8],[2,6],[3,1],[3,10]]
#     Output: 4
#     Explanation: The figure above shows the optimal allocation for four groups, where seats mark with blue are already reserved and contiguous seats mark with orange are for one group.
#
#
# **Example 2:**
#
#
#
#     Input: n = 2, reservedSeats = [[2,1],[1,8],[2,6]]
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: n = 4, reservedSeats = [[4,3],[1,4],[4,6],[1,7]]
#     Output: 4
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 10^9`
#   * `1 <= reservedSeats.length <= min(10*n, 10^4)`
#   * `reservedSeats[i].length == 2`
#   * `1 <= reservedSeats[i][0] <= n`
#   * `1 <= reservedSeats[i][1] <= 10`
#   * All `reservedSeats[i]` are distinct.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
@doc raw"""
    Niblack <: AbstractImageBinarizationAlgorithm
    Niblack(; window_size = 7, bias = 0.2)

    binarize([T,] img, f::Niblack)
    binarize!([out,] img, f::Niblack)

Applies Niblack adaptive thresholding [1] under the assumption that the input
image is textual.

# Output

Return the binarized image as an `Array{Gray{T}}` of size `size(img)`. If
`T` is not specified, it is inferred from `out` and `img`.

# Details

The input image is binarized by varying the threshold across the image, using
a modified version of Niblack's algorithm [2]. A threshold ``T`` is defined for
each pixel based on the mean ``m`` and standard deviation ``s`` of the
intensities of neighboring pixels in a window around it. This threshold is given
by

```math
T(x,y) = m(x,y) + k \cdot s(x,y),
```

where ``k`` is a user-defined parameter weighting the influence of the standard
deviation on the value of ``T``.

Note that Niblack's algorithm is highly sensitive to variations in the gray
values of background pixels, which often exceed local thresholds and appear as
artifacts in the binarized image. The [`Sauvola`](@ref) algorithm included in
this package implements an attempt to address this issue [2].

# Arguments

## `img`

An image which is binarized according to a per-pixel adaptive
threshold into background (0) and foreground (1) pixel values.

## `bias::Real`  (denoted by ``k`` in the publication)

A user-defined biasing parameter on threshold. This can take negative values.
Larger `bias` encourages more black pixels in the output.

## `window_size::Integer` (denoted by ``w`` in the publication)

The threshold for each pixel is a function of the distribution of the intensities
of all neighboring pixels in a square window around it. The side length of this
window is ``2w + 1``, with the target pixel in the center position.

If not specified, `window_size` is `7`.

# Example

Binarize the "cameraman" image in the `TestImages` package.

```julia
using TestImages, ImageBinarization

img = testimage("cameraman")
img₀₁ = binarize(img, Niblack())
```

# References

[1] Wayne Niblack (1986). *An Introduction to Image Processing*. Prentice-Hall, Englewood Cliffs, NJ: 115-16.
[2] J. Sauvola and M. Pietikäinen (2000). "Adaptive document image binarization". *Pattern Recognition* 33 (2): 225-236. [doi:10.1016/S0031-3203(99)00055-2](https://doi.org/10.1016/S0031-3203(99)00055-2)
"""
struct Niblack <: AbstractImageBinarizationAlgorithm
    window_size::Int
    bias::Float32

    function Niblack(window_size::Integer, bias::Real)
        window_size < 0 && throw(ArgumentError("window_size should be non-negative."))
        new(window_size, bias)
    end
end

Niblack(; window_size::Integer=7, bias::Real=0.2) = Niblack(window_size, bias)

function (f::Niblack)(out::GenericGrayImage,
                      img::GenericGrayImage)

    size(out) == size(img) || throw(ArgumentError("out and img should have the same shape, instead they are $(size(out)) and $(size(img))"))

    k = f.bias
    img_raw = channelview(img)
    I = integral_image(img_raw)
    I² = integral_image(img_raw.^2)

    function threshold(pixel::CartesianIndex{2})
        row₀, col₀, row₁, col₁ = get_window_bounds(img, pixel, f.window_size)
        m = μ_in_window(I, row₀, col₀, row₁, col₁)
        s = σ_in_window(I², m, row₀, col₀, row₁, col₁)
        return m + (k * s)
    end

    @simd for pixel in CartesianIndices(img)
        out[pixel] = img[pixel] <= threshold(pixel) ? 0 : 1
    end

    return out
end

(f::Niblack)(out::GenericGrayImage, img::AbstractArray{<:Color3}) =
    f(out, of_eltype(Gray, img))


text:
export TraceTestPencil

"""
    TraceTestPencil(F, p₀, v) <: AbstractSystem

The system `H(x,t) := [F(x, p₀+x[end]v); l ⋅[x[1:end-1];1] * λ + τ*γ)`.
"""
struct TraceTestPencil{S<:AbstractSystem,T} <: AbstractSystem
    F::TraceTestSystem{S,T}
    l::Vector{T}
    γ::T
end

function TraceTestPencil(F::TraceTestSystem, l)
    TraceTestPencil(F, l, cis(2π * rand()))
end

struct TraceTestPencilCache{SC<:TraceTestSystemCache} <: AbstractSystemCache
    cache::SC
end

function cache(F::TraceTestPencil, x, τ)
    c = cache(F.F, x, F.l)
    TraceTestPencilCache(c)
end

Base.size(F::TraceTestPencil) = size(F.F)

function evaluate!(u, F::TraceTestPencil, x, τ, c::TraceTestPencilCache)
    λ = x[end]
    evaluate!(u, F.F, x, F.l, c.cache)
    u[end] = u[end] * λ + τ[1] * F.γ
    u
end
function evaluate(F::TraceTestPencil, x, τ, c::TraceTestPencilCache)
    evaluate!(evaluate(F.F, x, F.l, c.cache), F, x, τ, c)
end

function jacobian!(U, F::TraceTestPencil, x, τ, c::TraceTestPencilCache)
    λ = x[end]
    jacobian!(U, F.F, x, F.l, c.cache)
    for j = 1:length(x)-1
        U[end, j] *= λ
    end
    U[end, end] = F.l[end]
    for i = 1:length(x)-1
        U[end, end] += F.l[i] * x[i]
    end
    U
end
function jacobian(F::TraceTestPencil, x, τ, c::TraceTestPencilCache)
    jacobian!(jacobian(F.F, x, F.l, c.cache), F, x, τ, c)
end

function evaluate_and_jacobian!(u, U, F::TraceTestPencil, x, τ, c::TraceTestPencilCache)
    λ = x[end]
    evaluate_and_jacobian!(u, U, F.F, x, F.l, c.cache)
    u[end] = u[end] * λ + τ[1] * F.γ
    for j = 1:length(x)-1
        U[end, j] *= λ
    end
    U[end, end] = F.l[end]
    for i = 1:length(x)-1
        U[end, end] += F.l[i] * x[i]
    end
    nothing
end

function differentiate_parameters!(u, F::TraceTestPencil, x, τ, c::TraceTestPencilCache)
    u .= zero(eltype(u))
    u[end, 1] = F.γ
    u
end

function differentiate_parameters(F::TraceTestPencil, x, τ, c::TraceTestPencilCache)
    u = similar(c.cache.u, size(F)[1], 1)
    differentiate_parameters!(u, F, x, τ, c)
end


text:
using FlagColors
using Test

@testset "FlagColors.jl" begin
    # Write your tests here.
end