text:
using AcousticFeatures
using AcousticFeatures: name

using AlphaStableDistributions
using BenchmarkTools
using Distributions
using LazyWAVFiles
using LinearAlgebra
using SignalAnalysis
using Test
using WAV

tmpdir = mktempdir()
fs = 100_000
N = 100_000
A = 1.0
frequency = 10_000
t = (0:N-1)./fs

@testset "AcousticFeatures" begin

    @testset "Energy" begin
        @info "Testing Energy"

        x = A.*sin.(2π*frequency*t)
        @test Score(Energy(), x)[1] ≈ (A^2)/2
        winlens = [1_000, 10_000, 1_001, 10_001]
        noverlaps = [0, 100, 500]
        for winlen in winlens, noverlap in noverlaps
            subseq = Subsequence(x, winlen, noverlap)
            sc = Score(Energy(), x; winlen=winlen, noverlap=noverlap)
            spart = sc[(sc.axes[1] .> subseq.winlen÷2) .& (sc.axes[1] .< length(x)-subseq.winlen÷2)]
            @test all(isapprox.(spart, repeat([(A^2)/2], length(spart)), atol=0.001))
        end
        WAV.wavwrite(x[1:N÷2], joinpath(tmpdir, "1.wav"), Fs=fs)
        WAV.wavwrite(x[N÷2:end], joinpath(tmpdir, "2.wav"), Fs=fs)
        dfile = DistributedWAVFile(tmpdir)
        for winlen in winlens, noverlap in noverlaps
            subseq = Subsequence(dfile, winlen, noverlap)
            sc = Score(Energy(), dfile; winlen=winlen, noverlap=noverlap)
            spart = sc[(sc.axes[1] .> subseq.winlen÷2) .& (sc.axes[1] .< length(x)-subseq.winlen÷2)]
            @test all(isapprox.(spart, repeat([(A^2)/2], length(spart)), atol=0.001))
        end

        @test name(Energy()) == ["Energy"]
    end

    @testset "Myriad" begin
        @info "Testing Myriad"

        α = 1.9999
        scale = 1.0
        x = rand(AlphaStable(α=α, scale=scale), N)
        d = fit(AlphaStable, x)
        sqKscale = myriadconstant(d.α, d.scale)
        @test Score(Myriad(), x)[1]/N ≈ (log((d.α/(2-d.α+eps()))*(d.scale^2))) atol=0.1
        winlens = [1_000, 10_000, 1_001, 10_001]
        noverlaps = [0, 100, 500]
        for winlen in winlens, noverlap in noverlaps
            subseq = Subsequence(x, winlen, noverlap)
            sc = Score(Myriad(sqKscale), x; winlen=winlen, noverlap=noverlap)
            spart = sc[(sc.axes[1] .> subseq.winlen÷2) .& (sc.axes[1] .< length(x)-subseq.winlen÷2)]
            @test all(isapprox.(spart./subseq.winlen, repeat([(log((d.α/(2-d.α+eps()))*(d.scale^2)))], length(spart)), atol=0.1))
        end
        WAV.wavwrite(x[1:N÷2], joinpath(tmpdir, "1.wav"), Fs=fs)
        WAV.wavwrite(x[N÷2:end], joinpath(tmpdir, "2.wav"), Fs=fs)
        dfile = DistributedWAVFile(tmpdir)
        for winlen in winlens, noverlap in noverlaps
            subseq = Subsequence(dfile, winlen, noverlap)
            sc = Score(Myriad(sqKscale), dfile; winlen=winlen, noverlap=noverlap)
            spart = sc[(sc.axes[1] .> subseq.winlen÷2) .& (sc.axes[1] .< length(x)-subseq.winlen÷2)]
            @test all(isapprox.(spart./subseq.winlen, repeat([(log((d.α/(2-d.α+eps()))*(d.scale^2)))], length(spart)), atol=0.1))
        end

        @test name(Myriad()) == ["Myriad"]
    end

    # @testset "VMyriad" begin
    #     @info "Testing VMyriad"

    #     # α = 1.3
    #     # identitymatrix = zeros(5, 5)
    #     # identitymatrix[diagind(identitymatrix)] .= 1.0
    #     # d = AlphaSubGaussian(α=α, n=N)
    #     # x = rand(d)
    #     # d̂ = fit(AlphaSubGaussian, x, 4)

    #     α = 1.3
    #     d = AlphaSubGaussian(;α=α, n=N)
    #     x = rand(d)
    #     d̂ = fit(AlphaSubGaussian, x, 4)

    #     f0 = VMyriad(vmyriadconstant(d̂.α, d̂.R)...)
    #     f1 = VMyriad(vmyriadconstant(1.8, d̂.R)...)
    #     identitymatrix = zeros(5, 5)
    #     identitymatrix[diagind(identitymatrix)] .= 1.0
    #     f2 = VMyriad(vmyriadconstant(1.8, identitymatrix)...)

    #     @test Score(f0, x).s[1] ≈ Score(VMyriad(vmyriadconstant(d̂.α, d̂.R)...), x).s[1]
    #     @test Score(f0, x).s[1] < Score(f1, x).s[1]
    #     @test Score(f0, x).s[1] < Score(f2, x).s[1]
    # end

    @testset "FrequencyContours" begin
        @info "Testing FrequencyContours"

        duration = N/fs
        f11 = 10_000; f21 = 50_000
        f12 = 1_000; f22 = 20_000
        x1 = real(samples(chirp(f11, f21, duration, fs)))+real(samples(chirp(f12, f22, duration, fs)))
        x2 = real(samples(chirp(f11, f21, duration, fs)))
        n = 512
        nv = 256
        tnorm = 1.0
        fd = 1000.0
        minhprc = 99.0
        minfdist = 1000.0
        mintlen = 0.05
        sc1 = Score(FrequencyContours(fs, n, nv, tnorm, fd, minhprc, minfdist, mintlen), x1)
        sc2 = Score(FrequencyContours(fs, n, nv, tnorm, fd, minhprc, minfdist, mintlen), x2)
        @test sc1[1] > sc2[1]
        winlens = [10_000, 10_001]
        noverlaps = [0, 100, 500]
        for winlen in winlens, noverlap in noverlaps
            subseq = Subsequence(x1, winlen, noverlap)
            sc1 = Score(FrequencyContours(fs, n, nv, tnorm, fd, minhprc, minfdist, mintlen), x1, winlen=winlen, noverlap=noverlap)
            sc2 = Score(FrequencyContours(fs, n, nv, tnorm, fd, minhprc, minfdist, mintlen), x2, winlen=winlen, noverlap=noverlap)
            spart1 = sc1[(sc1.axes[1] .> subseq.winlen÷2) .& (sc1.axes[1] .< length(x1)-subseq.winlen÷2)]
            spart2 = sc2[(sc2.axes[1] .> subseq.winlen÷2) .& (sc2.axes[1] .< length(x1)-subseq.winlen÷2)]
            @test all(isless.(spart2, spart1))
        end
        tmpdir1 = mktempdir()
        tmpdir2 = mktempdir()
        WAV.wavwrite(x1[1:N÷2], joinpath(tmpdir1, "1.wav"), Fs=fs)
        WAV.wavwrite(x1[N÷2:end], joinpath(tmpdir1, "2.wav"), Fs=fs)
        WAV.wavwrite(x2[1:N÷2], joinpath(tmpdir2, "1.wav"), Fs=fs)
        WAV.wavwrite(x2[N÷2:end], joinpath(tmpdir2, "2.wav"), Fs=fs)
        dfile1 = DistributedWAVFile(tmpdir1)
        dfile2 = DistributedWAVFile(tmpdir2)
        for winlen in winlens, noverlap in noverlaps
            subseq = Subsequence(dfile1, winlen, noverlap)
            sc1 = Score(FrequencyContours(fs, n, nv, tnorm, fd, minhprc, minfdist, mintlen), dfile1, winlen=winlen, noverlap=noverlap)
            sc2 = Score(FrequencyContours(fs, n, nv, tnorm, fd, minhprc, minfdist, mintlen), dfile2, winlen=winlen, noverlap=noverlap)
            spart1 = sc1[(sc1.axes[1] .> subseq.winlen÷2) .& (sc1.axes[1] .< length(dfile1)-subseq.winlen÷2)]
            spart2 = sc2[(sc2.axes[1] .> subseq.winlen÷2) .& (sc2.axes[1] .< length(dfile1)-subseq.winlen÷2)]
            @test all(isless.(spart2, spart1))
        end

        @test name(FrequencyContours(fs, n, nv, tnorm, fd, minhprc, minfdist, mintlen)) == ["Frequency Contours"]
    end

    @testset "SoundPressureLevel" begin
        @info "Testing SoundPressureLevel"

        x = A.*sin.(2π*frequency*t)
        x = pressure(x, 0.0, 0.0)
        sc = Score(SoundPressureLevel(), x)
        @test sc[1] ≈ 20*log10(1/sqrt(2))
        winlens = [1_000, 10_000, 1_001, 10_001]
        noverlaps = [0, 100, 500]
        for winlen in winlens, noverlap in noverlaps
            subseq = Subsequence(x, winlen, noverlap)
            sc = Score(SoundPressureLevel(), x; winlen=winlen, noverlap=noverlap)
            spart = sc[(sc.axes[1] .> subseq.winlen÷2) .& (sc.axes[1] .< length(x)-subseq.winlen÷2)]
            @test all(isapprox.(spart, repeat([20*log10(1/sqrt(2))], length(spart)), atol=0.01))
        end
        WAV.wavwrite(x[1:length(x)÷2], joinpath(tmpdir, "1.wav"), Fs=fs)
        WAV.wavwrite(x[length(x)÷2:end], joinpath(tmpdir, "2.wav"), Fs=fs)
        dfile = DistributedWAVFile(tmpdir)
        for winlen in winlens, noverlap in noverlaps
            subseq = Subsequence(dfile, winlen, noverlap)
            sc = Score(SoundPressureLevel(), dfile; winlen=winlen, noverlap=noverlap)
            spart = sc[(sc.axes[1] .> subseq.winlen÷2) .& (sc.axes[1] .< length(x)-subseq.winlen÷2)]
            @test all(isapprox.(spart, repeat([20*log10(1/sqrt(2))], length(spart)), atol=0.01))
        end

        @test name(SoundPressureLevel()) == ["SPL"]
    end

    @testset "ImpulseStats" begin
        @info "Testing ImpulseStats"

        trueindices = [101, 2254, 5322, 8888]
        x = zeros(N)
        x[trueindices] .= 10.0
        x += 0.1 .* randn(N)

        sc1 = Score(ImpulseStats(fs), x)
        sc2 = Score(ImpulseStats(fs, 10, 1e-3), x)
        @test sc1[1, 1] == sc2[1, 1] == length(trueindices)
        truetimeintervals = diff(trueindices)
        @test sc1[1, 2] == sc2[1, 2] == mean(truetimeintervals)/fs
        @test sc1[1, 3] == sc2[1, 3] == var(truetimeintervals)/fs

        @test name(ImpulseStats(fs)) == ["Nᵢ", "μᵢᵢ", "varᵢᵢ"]
    end

    @testset "SymmetricAlphaStableStats" begin
        @info "Testing SymmetricAlphaStableStats"

        α = 1.6
        scale = 2.0
        d = AlphaStable(α=α, scale=scale)
        x = rand(d, N)
        sc = Score(SymmetricAlphaStableStats(), x)
        @test sc[1,1] ≈ α atol=0.1
        @test sc[1,2] ≈ scale atol=0.1

        @test name(SymmetricAlphaStableStats()) == ["α", "scale"]
    end

    @testset "Entropy" begin
        @info "Testing Entropy"

        x = A.*sin.(2π*6250*t)
        sc = Score(Entropy(fs, 256, 128), x)
        @test sc[1] ≈ 1.0 atol=1e-2
        @test sc[2] ≈ 0.0 atol=1e-2
        @test sc[3] ≈ 0.0 atol=1e-2

        @test name(Entropy(fs, 256, 128)) == ["Temporal Entropy","Spectral Entropy","Entropy Index"]
    end

    @testset "ZeroCrossingRate" begin
        @info "Testing ZeroCrossingRate"
        x = [100.0, 1.0, -2.0, 2.0, -100, 0.0, 10.0]
        sc = Score(ZeroCrossingRate(), x)
        @test sc[1] == 4 / length(x)

        @test name(ZeroCrossingRate()) == ["ZCR"]
    end

    @testset "SpectralCentroid" begin
        @info "Testing SpectralCentroid"

        x = A.*sin.(2π*6250*t)
        sc = Score(SpectralCentroid(fs), x)
        @test sc[1] ≈ 6250 atol=0.0001

        @test name(SpectralCentroid(fs)) == ["Spectral Centroid"]
    end

    @testset "SpectralFlatness" begin
        @info "Testing SpectralFlatness"

        x = A.*sin.(2π*6250*t)
        sc = Score(SpectralFlatness(), x)
        @test sc[1] ≈ 0.0 atol=0.0001

        x = randn(N)
        scnormal = Score(SpectralFlatness(), x)
        @test scnormal[1] > sc[1]

        @test name(SpectralFlatness()) == ["Spectral Flatness"]
    end

    @testset "PermutationEntropy" begin
        @info "Testing PermutationEntropy"

        x = [4,7,9,10,6,11,3]
        m = 3
        τ = 1
        norm1 = false
        norm2 = true
        sc1 = Score(PermutationEntropy(m, τ, norm1), x)
        sc2 = Score(PermutationEntropy(m, τ, norm2), x)
        sc3 = Score(PermutationEntropy(m), x)
        @test sc1[1] ≈ 1.5219 atol=0.0001
        @test sc2[1] ≈ 0.5887 atol=0.0001
        @test sc2[1] == sc3[1]

        @test name(PermutationEntropy(m)) == ["Permutation Entropy"]
    end

    @testset "PSD" begin
        @info "Testing PSD"

        freq = 3000
        fs = 96000
        x = cw(freq, 0.1, fs) |> real |> samples
        sc = Score(PSD(fs, 64, 32), x)
        @test sc.axes[2][argmax(sc)[2]] == "PSD-$(round(freq; digits=1))Hz"
    end

    @testset "AcousticComplexityIndex" begin
        @info "Testing AcousticComplexityIndex"

        freq1, freq2 = 1000, 48000
        fs = 96000
        x1 = chirp(freq1, freq2, 5.0, fs) |> real |> samples
        x2 = cw(freq1, 5.0, fs) |> real |> samples
        x1 += randn(length(x1))
        x2 += randn(length(x2))
        sc1 = Score(AcousticComplexityIndex(fs, 1048, 0, 30), x1)
        sc2 = Score(AcousticComplexityIndex(fs, 1048, 0, 30), x2)
        @test sc1[1] > sc2[1]
    end

    @testset "Score" begin
        @info "Testing Score"

        @test_throws ArgumentError Score(Energy(), randn(1000); winlen=1001)
        f = Energy()
        @test Score(f, randn(100000))[1] ≈ f(randn(100000))[1] atol=0.1
    end

    @testset "Subsequences" begin
        @info "Testing Subsequences"

        x = [1, 2, 3, 4, 5, 6, 7]
        tmpdir = mktempdir()
        a, b = [1, 2, 3, 4], [5, 6, 7]
        WAV.wavwrite(a, joinpath(tmpdir, "1.wav"), Fs=100)
        WAV.wavwrite(b, joinpath(tmpdir, "2.wav"), Fs=100)
        dfile = DistributedWAVFile(tmpdir)

        winlen = 3
        noverlap = 1
        subseq1 = [[0, 1, 2],
                    [2, 3, 4],
                    [4, 5, 6],
                    [6, 7, 0]]
        subseq2 = Subsequence(x, winlen, noverlap)
        subseq3 = Subsequence(dfile, winlen, noverlap)
        @test length(subseq1) == length(subseq2) == length(subseq3)
        for (s1, s2, s3) in zip(subseq1, subseq2, subseq3)
            @test s1 == s2 == s3
        end

        winlen = 3
        noverlap = 0
        subseq1 = [[0, 1, 2],
        [3, 4, 5],
        [6, 7, 0]]
        subseq2 = Subsequence(x, winlen, noverlap)
        subseq3 = Subsequence(dfile, winlen, noverlap)
        @test length(subseq1) == length(subseq2) == length(subseq3)
        for (i, (subseq1, subseq2, subseq3)) in enumerate(zip(subseq1, subseq2, subseq3))
            @test subseq1 == subseq2 == subseq3
            @test subseq1[i] == subseq2[i] == subseq3[i]
        end

        winlen = 4
        noverlap = 1
        subseq1 = [[0, 1, 2, 3],
                    [3, 4, 5, 6],
                    [6, 7, 0, 0]]
        subseq2 = Subsequence(x, winlen, noverlap)
        subseq3 = Subsequence(dfile, winlen, noverlap)
        @test length(subseq1) == length(subseq2) == length(subseq3)
        for (i, (subseq1, subseq2, subseq3)) in enumerate(zip(subseq1, subseq2, subseq3))
            @test subseq1 == subseq2 == subseq3
            @test subseq1[i] == subseq2[i] == subseq3[i]
        end

        winlen = 4
        noverlap = 1
        subseq1 = [[1, 1, 2, 3],
                    [3, 4, 5, 6],
                    [6, 7, 7, 7]]
        subseq2 = Subsequence(x, winlen, noverlap; padtype=:replicate)
        subseq3 = Subsequence(dfile, winlen, noverlap; padtype=:replicate)
        @test length(subseq1) == length(subseq2) == length(subseq3)
        for (i, (subseq1, subseq2, subseq3)) in enumerate(zip(subseq1, subseq2, subseq3))
            @test subseq1 == subseq2 == subseq3
            @test subseq1[i] == subseq2[i] == subseq3[i]
        end

        winlen = 4
        noverlap = 1
        subseq1 = [[7, 1, 2, 3],
                    [3, 4, 5, 6],
                    [6, 7, 1, 2]]
        subseq2 = Subsequence(x, winlen, noverlap; padtype=:circular)
        subseq3 = Subsequence(dfile, winlen, noverlap; padtype=:circular)
        @test length(subseq1) == length(subseq2) == length(subseq3)
        for (i, (subseq1, subseq2, subseq3)) in enumerate(zip(subseq1, subseq2, subseq3))
            @test subseq1 == subseq2 == subseq3
            @test subseq1[i] == subseq2[i] == subseq3[i]
        end

        winlen = 4
        noverlap = 1
        subseq1 = [[1, 1, 2, 3],
                    [3, 4, 5, 6],
                    [6, 7, 7, 6]]
        subseq2 = Subsequence(x, winlen, noverlap; padtype=:symmetric)
        subseq3 = Subsequence(dfile, winlen, noverlap; padtype=:symmetric)
        @test length(subseq1) == length(subseq2) == length(subseq3)
        for (i, (subseq1, subseq2, subseq3)) in enumerate(zip(subseq1, subseq2, subseq3))
            @test subseq1 == subseq2 == subseq3
            @test subseq1[i] == subseq2[i] == subseq3[i]
        end

        winlen = 4
        noverlap = 1
        subseq1 = [[2, 1, 2, 3],
                    [3, 4, 5, 6],
                    [6, 7, 6, 5]]
        subseq2 = Subsequence(x, winlen, noverlap; padtype=:reflect)
        subseq3 = Subsequence(dfile, winlen, noverlap; padtype=:reflect)
        @test length(subseq1) == length(subseq2) == length(subseq3)
        for (i, (subseq1, subseq2, subseq3)) in enumerate(zip(subseq1, subseq2, subseq3))
            @test subseq1 == subseq2 == subseq3
            @test subseq1[i] == subseq2[i] == subseq3[i]
        end

    end

    @testset "Utils" begin
        @info "Testing Utils"

        x = [1, 2, 3, 4, 5, 6, 7]
        Nnorm = 3
        xfilt = x - [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 6.0]
        @test spectrumflatten(x, Nnorm) == xfilt

        x = [[1 2 3 4 5 6 7];
             [8 9 10 11 12 13 14];
             [15 16 17 18 19 20 21]]
        xfiltrow = x - [[1.0 2.0 3.0 4.0 5.0 6.0 6.0];
                        [8.0 9.0 10.0 11.0 12.0 13.0 13.0];
                        [15.0 16.0 17.0 18.0 19.0 20.0 20.0]]
        xfiltcol = x - [[1.0 2.0 3.0 4.0 5.0 6.0 7.0];
                        [8.0 9.0 10.0 11.0 12.0 13.0 14.0];
                        [8.0 9.0 10.0 11.0 12.0 13.0 14.0]]
        @test spectrumflatten(x, Nnorm) == xfiltrow
        @test spectrumflatten(x, Nnorm; dims=1) == xfiltcol

        nbits = 16
        vref = 1.0
        xvolt = vref.*real(samples(cw(64, 1, 512)))
        xbit = xvolt*(2^(nbits-1))
        sensitivity = 0.0
        gain = 0.0
        p1 = pressure(xvolt, sensitivity, gain)
        p2 = pressure(xbit, sensitivity, gain, voltparams=(nbits, vref))
        @test p1 == p2

        p = [1,2,3,4,5,6,7]
        @test AcousticFeatures.ordinalpatterns(p,3,1) == [1.0]
        @test AcousticFeatures.ordinalpatterns(p,3,2) == [1.0]
        p = [1,2,1,2,1,2,1]
        @test AcousticFeatures.ordinalpatterns(p,3,1) == [0.6,0.4]
        @test AcousticFeatures.ordinalpatterns(p,3,2) == [1.0]
    end

    @testset "Benchmarks" begin
        @info "Benchmarks"

        path = mktempdir()
        y = sin.((0:99999999)/48000*2pi*440);
        wavwrite(y, joinpath(path, "test1.wav"), Fs=48000)

        dfile = DistributedWAVFile(path)
        subseqdf = Subsequence(dfile, 96000, 0)

        t = @belapsed $subseqdf[1]
        @test t < 0.1
        t = @belapsed $subseqdf[10]
        @test t < 0.1
        t = @belapsed $subseqdf[100]
        @test t < 0.1

        filepath = joinpath(path, "test1.wav")
        lfile = LazyWAVFile(filepath)
        subseqlf = Subsequence(lfile, 96000, 0)

        t = @belapsed $subseqlf[1]
        @test t < 0.01
        t = @belapsed $subseqlf[10]
        @test t < 0.01
        t = @belapsed $subseqlf[100]
        @test t < 0.01
    end

end


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

mutable struct IoK8sApiCoreV1Node <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    spec::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeSpec } # spec name: spec
    status::Any # spec type: Union{ Nothing, IoK8sApiCoreV1NodeStatus } # spec name: status

    function IoK8sApiCoreV1Node(;apiVersion=nothing, kind=nothing, metadata=nothing, spec=nothing, status=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1Node, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiCoreV1Node, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiCoreV1Node, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiCoreV1Node, Symbol("spec"), spec)
        setfield!(o, Symbol("spec"), spec)
        validate_property(IoK8sApiCoreV1Node, Symbol("status"), status)
        setfield!(o, Symbol("status"), status)
        o
    end
end # type IoK8sApiCoreV1Node

const _property_map_IoK8sApiCoreV1Node = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("spec")=>Symbol("spec"), Symbol("status")=>Symbol("status"))
const _property_types_IoK8sApiCoreV1Node = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("spec")=>"IoK8sApiCoreV1NodeSpec", Symbol("status")=>"IoK8sApiCoreV1NodeStatus")
Base.propertynames(::Type{ IoK8sApiCoreV1Node }) = collect(keys(_property_map_IoK8sApiCoreV1Node))
Swagger.property_type(::Type{ IoK8sApiCoreV1Node }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiCoreV1Node[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1Node }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1Node[property_name]

function check_required(o::IoK8sApiCoreV1Node)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1Node }, name::Symbol, val)
end

ext: jl


text:
# function to set y annotations with equal displacement
function yd!(y,d)
  y[1] = y[1] - d
  return y[1]
end

# function to draw a rectangle
rectangle(xmin, xmax, ymin, ymax) = Shape(xmin .+ [0,xmax-xmin,xmax-xmin,0], ymin .+ [0,0,ymax-ymin,ymax-ymin])


text:
"""
Module providing functionality to numerically discretize the (thermoacoustic) Helmholtz equation by first, second, and hermitian-order finite elements.
"""
module Helmholtz
import SparseArrays, LinearAlgebra, ProgressMeter
import FFTW: fft
using ..Meshutils, ..NLEVP #TODO:check where find_smplx is introduced to scope
import ..Meshutils: get_line_idx
import ..NLEVP: generate_1_gz
include("Meshutils_exports.jl")
include("NLEVP_exports.jl")
include("./FEM/FEM.jl")
include("./FEM/helmholtz_getters.jl")
include("shape_sensitivity.jl")
include("Bloch.jl")
export discretize

##
function outer(ii,aa,jj,bb)
    #TODO: preallocation
    mm=ComplexF64[]
    II=UInt32[]
    JJ=UInt32[]
    for (a,i) in zip(aa,ii)
        for (b,j) in zip(bb,jj)
            push!(mm,a*b)
            push!(II,i)
            push!(JJ,j)
        end
    end
    #TODO: conversion to array
    return II,JJ,mm
end

"""
    L=discretize(mesh, dscrp, C; order=:lin, b=:__none__, mass_weighting=true,source=false, output=true)

Discretize the Helmholtz equation using the mesh `mesh`.

# Arguments
- `mesh::Mesh`: tetrahedral mesh
- `dscrp::Dict `: dictionary containing information on where to apply physical constraints. Such as standard wave propagation, boundary conditions, flame responses, etc.
- `C:Array`: array defining the speed of sound. If `length(C)==length(mesh.tetrahedra)` the speed of sound is constant along one tetrahedron. If `length(C)==size(mesh.points,2)` the speed of sound is linearly interpolated between the vertices of the mesh.
- `order::Symbol = :lin`: optional paramater to select between first (`order==:lin` the default), second (`order==:quad`),or hermitian-order (`order==:herm`) finite elements.
- `b::Symbol=:__none__`: optional parameter defining the Bloch wave number. If `b=:__none__` (the default) no Blochwave formalism is applied.
- `mass_weighting=true`: optional parameter if true mass matrix is used as weighting matrix for householder, otherwise this matrix is not set.
- `source::Bool=false`: optional parameter to toggle the return of a source vector (experimental)
- `output::Bool=false': optional parameter to toggle live progress report.

# Returns
- `L::LinearOperatorFamily`: parametereized discretization of the specified Helmholtz equation.
- `rhs::LinearOperatorFamily`: parameterized discretization of the source vector. Only returned if `source==true`.  (experimental)
"""
function discretize(mesh::Mesh, dscrp, C; order=:lin, b=:__none__, mass_weighting=true, source=false, output=true)
    triangles,tetrahedra,dim=aggregate_elements(mesh,order)
    N_points=size(mesh.points,2)

    if length(C)==length(mesh.tetrahedra)
        C_tet=C
        if length(mesh.tri2tet)!=0 && mesh.tri2tet[1]==0xffffffff #TODO: Initialize as empty instead with sentinel value
            link_triangles_to_tetrahedra!(mesh)
        end
        C_tri=C[mesh.tri2tet]
    elseif length(C)==size(mesh.points,2)
        C_tet=Array{Float64,1}[] #TODO: Preallocation
        C_tri=Array{Float64,1}[]
        for tet in tetrahedra
            push!(C_tet,C[tet[1:4]])
        end
        for tri in triangles
            push!(C_tri,C[tri[1:3]])
        end
    end

    #initialize linear operator family...
    L=LinearOperatorFamily(["ω","λ"],complex([0.,Inf]))
    #...and source vector
    rhs=LinearOperatorFamily(["ω"],complex([0.,]))

    if b!=:__none__
        bloch=true
        naxis=mesh.dos.naxis
        nxbloch=mesh.dos.nxbloch
        nsector=naxis+mesh.dos.nxsector
        naxis_ln=mesh.dos.naxis_ln+N_points
        nsector_ln=mesh.dos.naxis_ln+mesh.dos.nxsector_ln+N_points
        Δϕ=2*pi/mesh.dos.DOS
        exp_plus(z,k)=exp_az(z,Δϕ*1.0im,k) #check whether this is a closure
        exp_minus(z,k)=exp_az(z,-Δϕ*1.0im,k)
        bloch_filt=zeros(ComplexF64,mesh.dos.DOS)
        bloch_filt[1]=(1.0+0.0im)/mesh.dos.DOS
        bloch_filt=fft(bloch_filt)
        bloch_filt=generate_Σy_exp_ikx(bloch_filt)
        anti_bloch_filt=generate_1_gz(bloch_filt)
        bloch_exp_plus=generate_gz_hz(bloch_filt,exp_plus)
        bloch_exp_minus=generate_gz_hz(bloch_filt,exp_minus)
        txt_plus="*exp(i$(b)2π/$(mesh.dos.DOS))"
        txt_minus="*exp(-i$(b)2π/$(mesh.dos.DOS))"
        txt_filt="*δ($b)"
        txt_filt_plus=txt_filt*txt_plus
        txt_filt_minus=txt_filt*txt_minus
        #txt_filt="*1($b)"
        #bloch_filt=pow0
        L.params[b]=0.0+0.0im
        if order==:lin
            dim-=mesh.dos.nxbloch
        elseif order==:quad
            dim-=mesh.dos.nxbloch+mesh.dos.nxbloch_ln
        elseif order==:herm
            dim-=4*mesh.dos.nxbloch #TODO: check blochify for hermitian elements
        end
        #N_points_dof=N_points-mesh.dos.nxbloch
    else
        bloch=false
        naxis=nsector=0
    end
    #wrapper for FEM constructors TODO: use multipledispatch and move to FEM.jl
    function stiff(J,c)
        if order==:lin
            if length(c)==1
                return -c^2*s43nv1nu1(J)
            elseif length(c)==4
                return -s43nv1nu1cc1(J,c)
            end
        elseif order==:quad
            if length(c)==1
                return -c^2*s43nv2nu2(J)
            elseif length(c)==4
                return -s43nv2nu2cc1(J,c)
            end
        elseif order==:herm
            if length(c)==1
                return -c^2*s43nvhnuh(J)
            elseif length(c)==4
                return -s43nvhnuhcc1(J,c)
            end
        end
    end
    function mass(J)
        if order==:lin
            return s43v1u1(J)
        elseif order==:quad
            return s43v2u2(J)
        elseif order==:herm
            return s43vhuh(J)#massh(J)
        end
    end

    function bound(J,c)
        if order==:lin
            if length(c)==1
                return c*s33v1u1(J)
            elseif length(c)==3
                return s33v1u1c1(J,c)
            end
        elseif order==:quad
            if length(c)==1
                return c*s33v2u2(J)
            elseif length(c)==3
                return s33v2u2c1(J,c)
            end
        elseif order==:herm
            if length(c)==1
                return c*s33vhuh(J)
            elseif length(c)==3
                return s33vhuhc1(J,c)
            end
        end
    end

    function volsrc(J)
        if order==:lin
            return s43v1(J)
        elseif order==:quad
            return s43v2(J)
        elseif order==:herm
            return s43vh(J)
        end
    end

    function gradsrc(J,n,x)
        if order==:lin
            return s43nv1rx(J,n,x)
        elseif order==:quad
            return s43nv2rx(J,n,x)
        elseif order==:herm
            return s43nvhrx(J,n,x)
        end
    end

    function wallsrc(J,c)
        if length(c)==1
            if order==:lin
                return c*s33v1(J)
            elseif order==:quad
                return c*s33v2(J)
            elseif order==:herm
                return c*s33vh(J)
            end
        else
            if order==:lin
                return s33v1c1(J,c)
            elseif order==:quad
                return s33v2c1(J,c)
            elseif order==:herm
                return s33vhc1(J,c)
            end
        end
    end

    #prepare progressbar
    if output
        n_task=0
        for (domain,(type,data)) in dscrp
            if type==:interior
                task_factor=2
            else
                task_factor=1
            end
                n_task+=task_factor*length(mesh.domains[domain]["simplices"])
        end

        p = ProgressMeter.Progress(n_task,desc="Discretize... ", dt=.5,
             barglyphs=ProgressMeter.BarGlyphs('|','█', ['▁' ,'▂' ,'▃' ,'▄' ,'▅' ,'▆', '▇'],' ','|',),
             barlen=20)

    end

    ## build discretization matrices and vectors from domain definitions
    for (domain,(type,data)) in dscrp
        if type==:interior
            make=[:M,:K]
            stiff_func=()
            stiff_arg=()
            stiff_txt=""

        elseif type==:mass
            make=[:M]

        elseif type==:stiff
            make=[:K]
            stiff_func,stiff_arg,stiff_txt=data
            for args in stiff_arg
                for arg in args
                    L.params[arg]=0.0
                end
            end

        elseif type in (:admittance,:speaker)
            make=[]
            if type==:speaker
                append!(make,[:m])
                speak_sym,speak_val = data[1:2]
                rhs.params[speak_sym]=speak_val
                data = data[3:end]
            end

            if length(data)>0
                append!(make,[:C])
                if length(data)==2
                    adm_sym,adm_val=data

                    adm_txt="ω*"*string(adm_sym)
                    if adm_sym ∉ keys(L.params)
                        L.params[adm_sym]=adm_val
                        if type==:speaker
                            rhs.params[adm_sym]=adm_val
                        end
                    end
                    boundary_func=(pow1,pow1,)
                    boundary_arg=((:ω,),(adm_sym,),)
                    boundary_txt=adm_txt
                elseif length(data)==1
                    boundary_func=(generate_z_g_z(data[1]),)
                    boundary_arg=((:ω,),)
                    boundary_txt="ω*Y(ω)"
                elseif length(data)==4
                    Ass,Bss,Css,Dss = data
                    adm_txt="ω*C_s(iωI-A)^{-1}B"
                    func_z_stsp = generate_z_g_z(generate_stsp_z(Ass,Bss,Css,Dss))
                    boundary_func=(func_z_stsp,)
                    boundary_arg=((:ω,),)
                    boundary_txt=adm_txt
                end
            end

        elseif type==:flame #TODO: unified interface with flameresponse and normalize n_ref
            make=[:Q]
            isntau=false
            if length(data)==9 ## ref_idx is not specified by the user
                gamma,rho,nglobal,x_ref,n_ref,n_sym,tau_sym,n_val,tau_val=data
                ref_idx=0
                isntau=true
            elseif length(data)==10## ref_idx is specified by the user
                gamma,rho,nglobal,ref_idx,x_ref,n_ref,n_sym,tau_sym,n_val,tau_val=data
                isntau=true
            elseif length(data)==6 #custom FTF
                gamma,rho,nglobal,x_ref,n_ref,FTF=data
                ref_idx=0
                flame_func = (FTF,)
                flame_arg = ((:ω,),)
                if applicable(FTF,:ω)
                    flame_txt=FTF(:ω)
                else
                    flame_txt = "FTF(ω)"
                end
            elseif length(data)==5 #plain FTF
                gamma,rho,nglobal,x_ref,n_ref=data
                ref_idx=0
                L.params[:FTF]=0.0
                flame_func = (pow1,)
                flame_arg = ((:FTF,),)
                flame_txt = "FTF"
                gamma,rho,nglobal,x_ref,n_ref=data
            else
             println("Error: Data length does not match :flame option!")
            end

            nlocal=(gamma-1)/rho*nglobal/compute_size!(mesh,domain)
            if isntau
                if n_sym ∉ keys(L.params)
                    L.params[n_sym]=n_val
                end
                if tau_sym ∉ keys(L.params)
                    L.params[tau_sym]=tau_val
                end
                flame_func=(pow1,exp_delay,)
                flame_arg=((n_sym,),(:ω,tau_sym))
                flame_txt="$(string(n_sym))*exp(-iω$(string(tau_sym)))"
            end

            if ref_idx==0
                ref_idx=find_tetrahedron_containing_point(mesh,x_ref)
                #TODO: check whethere ref_idx is still 0
            end
            if ref_idx ∈ mesh.domains[domain]["simplices"]
                println("Warning: your reference point is inside the domain of heat release. (short-circuited FTF!)")
            end

        elseif type==:flameresponse
            make=[:Q]
            gamma,rho,nglobal,x_ref,n_ref,eps_sym,eps_val=data
            ref_idx=0
            nlocal=(gamma-1)/rho*nglobal/compute_size!(mesh,domain)
            if eps_sym ∉ keys(L.params)
                L.params[eps_sym]=eps_val
            end
            flame_func=(pow1,)
            flame_arg=((eps_sym,),)
            flame_txt="$(string(eps_sym))"
            if ref_idx==0
                ref_idx=find_tetrahedron_containing_point(mesh,x_ref)
            end

        elseif type==:fancyflame
            make=[:Q]
            gamma,rho,nglobal,x_ref,n_ref,n_sym,tau_sym, a_sym, n_val, tau_val, a_val=data
            ref_idx=0
            nlocal=(gamma-1)/rho*nglobal/compute_size!(mesh,domain)

            if typeof(n_val)<:Number #TODO: sanity check that other lists are same length
                if n_sym ∉ keys(L.params)
                    L.params[n_sym]=n_val
                end
                if tau_sym ∉ keys(L.params)
                    L.params[tau_sym]=tau_val
                end
                if a_sym ∉ keys(L.params)
                    L.params[a_sym]=a_val
                end

                flame_func=(pow1,exp_az2mzit,)
                flame_arg=((n_sym,),(:ω,tau_sym,a_sym,))
                flame_txt="$(string(n_sym))* exp($(string(a_sym))ω^2-iω$(string(tau_sym)))"
            else
                flame_arg=(:ω,)
                flame_txt=""
                for (ns, ts, as, nv, tv, av) in zip(n_sym, tau_sym, a_sym, n_val, tau_val,a_val)
                    L.params[ns]=nv
                    L.params[ts]=tv
                    L.params[as]=av
                    flame_arg=(flame_arg...,ns,ts,as,)
                    flame_txt*="[$(string(ns))* exp($(string(as))ω^2-iω$(string(ts)))+"
                end
                flame_txt=flame_txt[1:end-1]*"]"
                flame_arg=(flame_arg,)
                flame_func=(Σnexp_az2mzit,)
            end
            if ref_idx==0
                ref_idx=find_tetrahedron_containing_point(mesh,x_ref)
            end

        else
            make=[]
        end

        for opr in make
            V=ComplexF64[]
            I=UInt32[] #TODO: consider preallocation
            J=UInt32[]
            if opr==:M
                matrix=true #sentinel value to toggle assembley of matrix (true) or vector (false)
                for smplx in tetrahedra[mesh.domains[domain]["simplices"]]
                    CT=CooTrafo(mesh.points[:,smplx[1:4]])
                    ii,jj=create_indices(smplx)
                    vv=mass(CT)
                    append!(V,vv[:])
                    append!(I,ii[:])
                    append!(J,jj[:])
                    if output
                        ProgressMeter.next!(p)
                    end
                end
                func=(pow2,)
                arg=((:ω,),)
                txt="ω^2"
                mat="M"
            elseif opr==:K
                matrix=true
                for (smplx,c) in zip(tetrahedra[mesh.domains[domain]["simplices"]],C_tet[mesh.domains[domain]["simplices"]])
                    CT=CooTrafo(mesh.points[:,smplx[1:4]])
                    ii,jj=create_indices(smplx)
                    #println("#############")
                    #println("$CT")
                    vv=stiff(CT,c)
                    append!(V,vv[:])
                    append!(I,ii[:])
                    append!(J,jj[:])
                    if output
                        ProgressMeter.next!(p)
                    end
                end
                func=stiff_func
                arg=stiff_arg
                txt=stiff_txt
                #V=-V
                mat="K"
            elseif opr==:C
                matrix=true
                for (smplx,c) in zip(triangles[mesh.domains[domain]["simplices"]],C_tri[mesh.domains[domain]["simplices"]])
                    CT=CooTrafo(mesh.points[:,smplx[1:3]])
                    ii,jj=create_indices(smplx)
                    vv=bound(CT,c)
                    append!(V,vv[:])
                    append!(I,ii[:])
                    append!(J,jj[:])
                    if output
                        ProgressMeter.next!(p)
                    end
                end
                V.*=-1im
                func=boundary_func # func=(pow1,pow1,)
                arg=boundary_arg  # arg=((:ω,),(adm_sym,),)
                txt=boundary_txt  # txt=adm_txt
                mat="C"
            elseif opr==:Q
                matrix=true
                S=ComplexF64[]
                G=ComplexF64[]
                for smplx in tetrahedra[mesh.domains[domain]["simplices"]]
                    CT=CooTrafo(mesh.points[:,smplx[1:4]])
                    mm=volsrc(CT)
                    append!(S,mm[:])
                    append!(I,smplx[:])
                    if output
                        ProgressMeter.next!(p)
                    end
                end
                smplx=tetrahedra[ref_idx]
                CT=CooTrafo(mesh.points[:,smplx[1:4]])
                mm=gradsrc(CT,n_ref,x_ref)
                append!(G,mm[:])
                append!(J,smplx[:])
                G=-nlocal.*G
                I,J,V=outer(I,S,J,G)
                func =flame_func
                arg=flame_arg
                txt=flame_txt
                mat="Q"
            elseif opr==:m
                matrix=false
                for (smplx,c) in zip(triangles[mesh.domains[domain]["simplices"]],C_tri[mesh.domains[domain]["simplices"]])
                    CT=CooTrafo(mesh.points[:,smplx[1:3]])
                    ii=smplx[:]
                    vv=wallsrc(CT,c)
                    append!(V,vv[:])
                    append!(I,ii[:])
                    if output
                        ProgressMeter.next!(p)
                    end
                end
                V./=1im
                func=(boundary_func...,pow1,)
                arg=(boundary_arg...,(speak_sym,),)
                txt="speaker"
                mat="m"
            end

            if matrix
                #assemble discretization matrices in sparse format
                if bloch
                    for (i,j,v,f,a,t) in zip(blochify(I,J,V,naxis,nxbloch,nsector,naxis_ln,nsector_ln,N_points)...,[(),(exp_plus,),(exp_minus,),(bloch_filt,),(bloch_exp_plus,),(bloch_exp_minus,),],[(), ((b,),), ((b,),),((b,),),((b,),),((b,),),], ["", txt_plus, txt_minus, txt_filt, txt_filt_plus, txt_filt_minus,])
                        M =SparseArrays.sparse(i,j,v,dim,dim)
                        push!(L,Term(M,(func...,f...),(arg...,a...),txt*t,mat))
                    end
                else
                    M =SparseArrays.sparse(I,J,V,dim,dim)
                    push!(L,Term(M,func,arg,txt,mat))
                end
            else
                #assemble discretization vectors in sparse format
                M=SparseArrays.sparsevec(I,V,dim)
                push!(rhs,Term(M,func,arg,txt,mat))
            end

        end
    end

    if mass_weighting||bloch
        V=ComplexF64[]
        I=UInt32[] #IDEA: consider preallocation
        J=UInt32[]
        for smplx in tetrahedra
            CT=CooTrafo(mesh.points[:,smplx[1:4]])
            ii,jj=create_indices(smplx)
            vv=mass(CT)
            append!(V,vv[:])
            append!(I,ii[:])
            append!(J,jj[:])
        end
    end
    if bloch

        #IDEA: implement switch to select weighting matrix

        IJV=blochify(I,J,V,naxis,nxbloch,nsector,naxis_ln,nsector_ln,N_points,axis=false)#TODO check whether this matrix needs more blochfication... especially because tetrahedra is not periodic
        I=[IJV[1][1]..., IJV[1][2]..., IJV[1][3]...]
        J=[IJV[2][1]..., IJV[2][2]..., IJV[2][3]...]
        V=[IJV[3][1]..., IJV[3][2]..., IJV[3][3]...]
        M=SparseArrays.sparse(I,J,-V,dim,dim)

        if 0<naxis #modify axis dof for essential boundary condition when blochwave number !=0
            DI=1:naxis
            DV=ones(ComplexF64,naxis) #NOTE: make this more compact
            if order==:quad
                DI=vcat(DI,(N_points+1:naxis_ln).-nxbloch)
                DV=vcat(DV,ones(ComplexF64,naxis_ln-N_points))
            end
            for idx=1:naxis
                DV[idx]=1/M[idx,idx]
            end
            if order==:quad
                for (idx,ii) in enumerate((N_points+1:naxis_ln).-nxbloch)
                    DV[idx+naxis]=1/M[ii,ii]
                end
            end
            DM=SparseArrays.sparse(DI,DI,DV,dim,dim)
            push!(L,Term(DM,(anti_bloch_filt,),((:b,),),"(1-δ(b))","D"))
        end
    end

    if mass_weighting&&!bloch
        M=SparseArrays.sparse(I,J,-V,dim,dim)
    end
    push!(L,Term(M,(pow1,),((:λ,),),"-λ","__aux__"))

    if source #experimental return mode returning the source term
        return L,rhs
    else #classic return mode
        return L
    end
end

end #module Helmholtz
##
# module lin
# export assemble_volume_source, assemble_gradient_source, assemble_mass_operator, assemble_stiffness_operator, assemble_boundary_mass_operator
# include("FEMlin.jl")
# end
# module quad20
# export assemble_volume_source, assemble_gradient_source, assemble_mass_operator, assemble_stiffness_operator, assemble_boundary_mass_operator
# include("FEMquadlin.jl")
# end
#
# module quad
# export assemble_volume_source, assemble_gradient_source, assemble_mass_operator, assemble_stiffness_operator, assemble_boundary_mass_operator
# include("FEMquad.jl")
# end
#
# import .lin
# import .quad
# import .quad20
#
# """
#     L=discretize(mesh, dscrp, C; el_type=1, c_type=0, b=:__none__, mass_weighting=true)
#
# Discretize the Helmholtz equation using the mesh `mesh`.
#
# # Arguments
# - `mesh::Mesh`: tetrahedral mesh
# - `dscrp::Dict `: dictionary containing information on where to apply physical constraints. Such as standard wave propagation, boundary conditions, flame responses, etc.
# - `C:Array`: array defining the speed of sound. If `c_type==0` the speed of sound is constant along one tetrahedron and `length(C)==length(mesh.tetrahedra)`. If `c_type==1` the speed of sound is linearly interpolated between the vertices of the mesh and `length(C)==size(mesh.points,2)`.
# - `el_type = 1`: optional paramater to select between first (`el_type==1` the default) and second order (`el_type==2`) finite elements.
# - `c_type = 1`: optional parameter controlling whether speed of sound is constant on a tetrahedron or linearly interpolated between vertices.
# - `b::Symbol=:__none__`: optional parameter defining the Bloch wave number. If `b=:__none__` (the default) no Blochwave formalism is applied.
# - `mass_weighting=true`: optional parameter if true mass matrix is used as weighting matrix for householder, otherwise this matrix is not set.
#
# # Returns
# - `L::LinearOperatorFamily`: parametereized discretization of the specified Helmholtz equation.
# """
# function discretize(mesh::Mesh, dscrp, C; el_type=1, c_type=0, b=:__none__, mass_weighting=true)
#     if el_type==1 && c_type==0
#         #println("using linear finite elements and locally constant speed of sound")
#         assemble_volume_source, assemble_gradient_source, assemble_mass_operator, assemble_stiffness_operator, assemble_boundary_mass_operator =lin.assemble_volume_source, lin.assemble_gradient_source, lin.assemble_mass_operator, lin.assemble_stiffness_operator, lin.assemble_boundary_mass_operator
#     elseif el_type==2 && c_type==0
#         #println("using linear finite elements and locally constant speed of sound")
#         assemble_volume_source, assemble_gradient_source, assemble_mass_operator, assemble_stiffness_operator, assemble_boundary_mass_operator =quad20.assemble_volume_source, quad20.assemble_gradient_source, quad20.assemble_mass_operator, quad20.assemble_stiffness_operator, quad20.assemble_boundary_mass_operator
#     elseif el_type==2 && c_type==1
#         #println("using quadratic finite elements and locally linear speed of sound")
#         assemble_volume_source, assemble_gradient_source, assemble_mass_operator, assemble_stiffness_operator, assemble_boundary_mass_operator =quad.assemble_volume_source, quad.assemble_gradient_source, quad.assemble_mass_operator, quad.assemble_stiffness_operator, quad.assemble_boundary_mass_operator
#     else
#         println("ERROR: Elements not supported!")
#         return
#     end
#
#     L=LinearOperatorFamily(["ω","λ"],complex([0.,Inf]))
#     N_points=size(mesh.points)[2]
#     dim=deepcopy(N_points)
#     #Prepare Bloch wave analysis
#     if b!=:__none__
#         bloch=true
#         naxis=mesh.dos.naxis
#         nxbloch=mesh.dos.nxbloch
#         nsector=naxis+mesh.dos.nxsector
#         naxis_ln=mesh.dos.naxis_ln+N_points
#         nsector_ln=mesh.dos.naxis_ln+mesh.dos.nxsector_ln+N_points
#         Δϕ=2*pi/mesh.dos.DOS
#         exp_plus(z,k)=exp_az(z,Δϕ*1.0im,k) #check whether this is a closure
#         exp_minus(z,k)=exp_az(z,-Δϕ*1.0im,k)
#         bloch_filt=zeros(ComplexF64,mesh.dos.DOS)
#         bloch_filt[1]=(1.0+0.0im)/mesh.dos.DOS
#         bloch_filt=fft(bloch_filt)
#         bloch_filt=generate_Σy_exp_ikx(bloch_filt)
#         anti_bloch_filt=generate_1_gz(bloch_filt)
#         bloch_exp_plus=generate_gz_hz(bloch_filt,exp_plus)
#         bloch_exp_minus=generate_gz_hz(bloch_filt,exp_minus)
#         txt_plus="*exp(i$(b)2π/$(mesh.dos.DOS))"
#         txt_minus="*exp(-i$(b)2π/$(mesh.dos.DOS))"
#         txt_filt="*δ($b)"
#         txt_filt_plus=txt_filt*txt_plus
#         txt_filt_minus=txt_filt*txt_minus
#         #txt_filt="*1($b)"
#         #bloch_filt=pow0
#         L.params[b]=0.0+0.0im
#         dim-=mesh.dos.nxbloch
#         #N_points_dof=N_points-mesh.dos.nxbloch
#     else
#         bloch=false
#         naxis=nsector=0
#     end
#
#
#     #aggregator to create local elements ( tetrahedra, triangles)
#     if el_type==2
#         triangles=Array{UInt32,1}[] #TODO: preallocation?
#         tetrahedra=Array{UInt32,1}[]
#         tet=Array{UInt32}(undef,10)
#         tri=Array{UInt32}(undef,6)
#         for (idx,smplx) in enumerate(mesh.tetrahedra) #TODO: no enumeration
#             tet[1:4]=smplx[:]
#             tet[5]=get_line_idx(mesh,smplx[[1,2]])+N_points#find_smplx(mesh.lines,smplx[[1,2]])+N_points #TODO: type stability
#             tet[6]=get_line_idx(mesh,smplx[[1,3]])+N_points
#             tet[7]=get_line_idx(mesh,smplx[[1,4]])+N_points
#             tet[8]=get_line_idx(mesh,smplx[[2,3]])+N_points
#             tet[9]=get_line_idx(mesh,smplx[[2,4]])+N_points
#             tet[10]=get_line_idx(mesh,smplx[[3,4]])+N_points
#             push!(tetrahedra,copy(tet))
#         end
#         for (idx,smplx) in enumerate(mesh.triangles)
#             tri[1:3]=smplx[:]
#             tri[4]=get_line_idx(mesh,smplx[[1,2]])+N_points
#             tri[5]=get_line_idx(mesh,smplx[[1,3]])+N_points
#             tri[6]=get_line_idx(mesh,smplx[[2,3]])+N_points
#             push!(triangles,copy(tri))
#         end
#         dim+=size(mesh.lines)[1]
#         if bloch
#             dim-=mesh.dos.nxbloch_ln
#         end
#
#     elseif el_type==1
#             tetrahedra=mesh.tetrahedra
#             triangles=mesh.triangles
#     end
#
#     if c_type==0
#         C_tet=C
#         if mesh.tri2tet[1]==0xffffffff
#             link_triangles_to_tetrahedra!(mesh)
#         end
#         C_tri=C[mesh.tri2tet]
#     elseif c_type==1
#         C_tet=Array{UInt32,1}[] #TODO: Preallocation
#         C_tri=Array{UInt32,1}[]
#         for tet in tetrahedra
#             push!(C_tet,C[tet[1:4]])
#         end
#         for tri in triangles
#             push!(C_tri,C[tri[1:3]])
#         end
#     end
#
#     ## build matrices from domain definitions
#     for (domain,(type,data)) in dscrp
#         if type==:interior
#             make=[:M,:K]
#
#         elseif type==:admittance
#             make=[:C]
#             if length(data)==2
#                 adm_sym,adm_val=data
#
#                 adm_txt="ω*"*string(adm_sym)
#                 if adm_sym ∉ keys(L.params)
#                     L.params[adm_sym]=adm_val
#                 end
#                 boundary_func=(pow1,pow1,)
#                 boundary_arg=((:ω,),(adm_sym,),)
#                 boundary_txt=adm_txt
#             elseif length(data)==1
#                 boundary_func=(generate_z_g_z(data[1]),)
#                 boundary_arg=((:ω,),)
#                 boundary_txt="ω*Y(ω)"
#             elseif length(data)==4
#                 Ass,Bss,Css,Dss = data
#                 adm_txt="ω*C_s(iωI-A)^{-1}B"
#                 func_z_stsp = generate_z_g_z(generate_stsp_z(Ass,Bss,Css,Dss))
#                 boundary_func=(func_z_stsp,)
#                 boundary_arg=((:ω,),)
#                 boundary_txt=adm_txt
#             end
#
#         elseif type==:flame #TODO: unified interface with flameresponse and normalize n_ref
#             make=[:Q]
#             gamma,rho,nglobal,x_ref,n_ref,n_sym,tau_sym,n_val,tau_val=data
#             nlocal=(gamma-1)/rho*nglobal/compute_size!(mesh,domain)
#             if n_sym ∉ keys(L.params)
#                 L.params[n_sym]=n_val
#             end
#             if tau_sym ∉ keys(L.params)
#                 L.params[tau_sym]=tau_val
#             end
#             flame_func=(pow1,exp_delay,)
#             flame_arg=((n_sym,),(:ω,tau_sym))
#             flame_txt="$(string(n_sym))*exp(-iω$(string(tau_sym)))"
#
#         elseif type==:flameresponse
#             make=[:Q]
#             gamma,rho,nglobal,x_ref,n_ref,eps_sym,eps_val=data
#             nlocal=(gamma-1)/rho*nglobal/compute_size!(mesh,domain)
#             if eps_sym ∉ keys(L.params)
#                 L.params[eps_sym]=eps_val
#             end
#             flame_func=(pow1,)
#             flame_arg=((eps_sym,),)
#             flame_txt="$(string(eps_sym))"
#         elseif type==:fancyflame
#             make=[:Q]
#             gamma,rho,nglobal,x_ref,n_ref,n_sym,tau_sym, a_sym, n_val, tau_val, a_val=data
#             nlocal=(gamma-1)/rho*nglobal/compute_size!(mesh,domain)
#             if n_sym ∉ keys(L.params)
#                 L.params[n_sym]=n_val
#             end
#             if tau_sym ∉ keys(L.params)
#                 L.params[tau_sym]=tau_val
#             end
#             if a_sym ∉ keys(L.params)
#                 L.params[a_sym]=a_val
#             end
#
#             flame_func=(pow1,exp_ax2mxit,)
#             flame_arg=((n_sym,),(:ω,tau_sym,a_sym,))
#             flame_txt="$(string(n_sym))* exp($(string(a_sym))ω^2-iω$(string(tau_sym)))"
#         else
#             make=[]
#         end
#
#         for opr in make
#             if opr==:M
#                 I,J,V=assemble_mass_operator(mesh.points,tetrahedra[mesh.domains[domain]["simplices"]])
#                 func=(pow2,)
#                 arg=((:ω,),)
#                 txt="ω^2"
#                 mat="M"
#             elseif opr==:K
#                 I,J,V=assemble_stiffness_operator(mesh.points,tetrahedra[mesh.domains[domain]["simplices"]],C_tet[mesh.domains[domain]["simplices"]])
#                 func=()
#                 arg=()
#                 txt=""
#                 #V=-V
#                 mat="K"
#             elseif opr==:C
#                 I,J,V=assemble_boundary_mass_operator(mesh.points,triangles[mesh.domains[domain]["simplices"]],C_tri[mesh.domains[domain]["simplices"]])
#                 V*=-1im
#                 func=boundary_func # func=(pow1,pow1,)
#                 arg=boundary_arg  # arg=((:ω,),(adm_sym,),)
#                 txt=boundary_txt  # txt=adm_txt
#                 mat="C"
#             elseif opr==:Q
#                 I,S=assemble_volume_source(mesh.points,tetrahedra[mesh.domains[domain]["simplices"]])
#                 ref_idx=find_tetrahedron_containing_point(mesh,x_ref)
#                 J,G=assemble_gradient_source(mesh.points,tetrahedra[ref_idx],x_ref,n_ref)
#                 G=-nlocal*G
#                 #G=-G
#                 I,J,V=outer(I,S,J,G)
#                 func =flame_func
#                 arg=flame_arg
#                 txt=flame_txt
#                 mat="Q"
#             end
#
#             if bloch
#                 for (i,j,v,f,a,t) in zip(blochify(I,J,V,naxis,nxbloch,nsector,naxis_ln,nsector_ln,N_points)...,[(),(exp_plus,),(exp_minus,),(bloch_filt,),(bloch_exp_plus,),(bloch_exp_minus,),],[(), ((b,),), ((b,),),((b,),),((b,),),((b,),),], ["", txt_plus, txt_minus, txt_filt, txt_filt_plus, txt_filt_minus,])
#                     M =SparseArrays.sparse(i,j,v,dim,dim)
#                     push!(L,Term(M,(func...,f...),(arg...,a...),txt*t,mat))
#                 end
#
#             else
#                 M =SparseArrays.sparse(I,J,V,dim,dim)
#                 push!(L,Term(M,func,arg,txt,mat))
#             end
#
#         end
#     end
#
#
#     if mass_weighting||bloch
#         I,J,V=assemble_mass_operator(mesh.points,tetrahedra)
#     end
#     if bloch
#
#         #TODO: implement switch to select weighting matrix
#
#         IJV=blochify(I,J,V,naxis,nxbloch,nsector,naxis_ln,nsector_ln,N_points,axis=false)#TODO check whether this matrix needs more blochfication... especially because tetrahedra is not periodic
#         I=[IJV[1][1]..., IJV[1][2]..., IJV[1][3]...]
#         J=[IJV[2][1]..., IJV[2][2]..., IJV[2][3]...]
#         V=[IJV[3][1]..., IJV[3][2]..., IJV[3][3]...]
#         M=SparseArrays.sparse(I,J,-V,dim,dim)
#
#         if 0<naxis #modify axis dof for essential boundary condition when blochwave number !=0
#             DI=1:naxis
#             DV=ones(ComplexF64,naxis) #TODO: make this more compact
#             if el_type==2
#                 DI=vcat(DI,(N_points+1:naxis_ln).-nxbloch)
#                 DV=vcat(DV,ones(ComplexF64,naxis_ln-N_points))
#             end
#             for idx=1:naxis
#                 DV[idx]=1/M[idx,idx]
#             end
#             if el_type==2
#                 for (idx,ii) in enumerate((N_points+1:naxis_ln).-nxbloch)
#                     DV[idx+naxis]=1/M[ii,ii]
#                 end
#             end
#             DM=SparseArrays.sparse(DI,DI,DV,dim,dim)
#             push!(L,Term(DM,(anti_bloch_filt,),((:b,),),"(1-δ(b))","D"))
#         end
#     end
#
#     if mass_weighting&&!bloch
#         M=SparseArrays.sparse(I,J,-V,dim,dim)
#     end
#     push!(L,Term(M,(pow1,),((:λ,),),"-λ","__aux__"))
#
#     return L
# end
#end#module


text:
#=
# ReLU layer

In this example we look at what the ReLU layer hidden units look like,
for different parameter values.
=#

#=
First load some packages.
=#

import RestrictedBoltzmannMachines as RBMs
using CairoMakie, Statistics
nothing #hide

#=
Now initialize our ReLU layer, with unit parameters spanning an interesting range.
=#

θs = [0; 10]
γs = [5; 10]
layer = RBMs.ReLU([θ for θ in θs, γ in γs], [γ for θ in θs, γ in γs])
nothing #hide

#=
Now we sample our layer to collect some data.
=#

data = RBMs.sample_from_inputs(layer, zeros(size(layer)..., 10^6))
nothing #hide

#=
Let's plot the resulting histogram of the activations of each unit.
We also overlay the analytical PDF.
=#

fig = Figure(resolution=(700,500))
ax = Axis(fig[1,1])
xs = repeat(reshape(range(minimum(data), maximum(data), 100), 1, 1, 100), size(layer)...)
ps = exp.(RBMs.free_energies(layer) .- RBMs.energies(layer, xs))
for (iθ, θ) in enumerate(θs), (iγ, γ) in enumerate(γs)
    hist!(ax, data[iθ, iγ, :], normalization=:pdf, label="θ=$θ, γ=$γ")
    lines!(xs[iθ, iγ, :], ps[iθ, iγ, :], linewidth=2)
end
axislegend(ax)
fig


text:
if Pkg.installed("ImageView") != nothing
    using ImageView

    export triplot

    function triplot(img::Image, shape::Matrix{Float64}, trigs::Matrix{Int64})
        imgc, img2 = view(img)
        for i=1:size(trigs, 1)
            a = (shape[trigs[i, 1], 2], shape[trigs[i, 1], 1])
            b = (shape[trigs[i, 2], 2], shape[trigs[i, 2], 1])
            c = (shape[trigs[i, 3], 2], shape[trigs[i, 3], 1])
            annotate!(imgc, img2, AnnotationLine(a, b))
            annotate!(imgc, img2, AnnotationLine(b, c))
            annotate!(imgc, img2, AnnotationLine(c, a))
        end
        imgc, img2
    end

    triplot{T,N}(mat::Array{T, N}, shape::Matrix{}, trigs::Matrix{Int64}) =
        triplot(convert(Image, mat), shape, trigs)
end


text:
# For defining the ODEs for describing all the molecular dynamics.
# Designed for use with the Sundials ODE solver.
# Args: xdot = vector, [dCa2+/dt, dBCa/dt, dECa/dt, dImCa/dt, dBCa*/dt]
#       x = current state of continuous variables, [Ca2+, BCa, ECa, ImCa, BCa*]
#       t = time
#       params =  Dict, parameters [calcium_baseline, ca_influx, ca_outflux, b_i, f_i, indicator,
#                 b_e, f_e, endogeneous, b_im, f_im, immobile,
#                 excitation_rate, release_rate, previous_BCa*]
# Returns: nothing

function molecularODE(t::Float64, x::Array{Float64}, xdot::Array{Float64}, baseline::Int64,
                      calcium_rate::Float64, indicator::Int64, endogeneous::Int64,
                      immobile::Int64, b_i::Float64, f_i::Float64, b_e::Float64, f_e::Float64,
                      b_im::Float64, f_im::Float64, excitation::Float64, release::Float64)
  B, E, Im = indicator - (x[2] + x[5]), endogeneous - x[3], immobile - x[4]
  # d[Ca2+]/dt = b_i*[BCa] + b_e*[ECa] + b_im*x[ImCa] - f_i*[Ca2+]*[B] - f_e*[Ca2+]*[E] + f_im*x[Ca2+]*[Im] + \beta_{Ca}([Ca2+]_0 - [Ca2+])
  xdot[1] = b_i*x[2] + b_e*x[3] + b_im*x[4] - f_i*x[1]*B - f_e*x[1]*E - f_im*x[1]*Im + calcium_rate * (baseline - x[1])
  # d[BCa]/dt = f_i*[Ca2+]*[B] - b_i*[BCa] + release*[BCa*] - excitation*[BCa]
  xdot[2] = f_i*x[1]*B - b_i*x[2] + release*x[5] - excitation*x[2]
  # d[ECa]/dt = f_e*[Ca2+]*[E] - b_e*[ECa]
  xdot[3] = f_e*x[1]*E - b_e*x[3]
  # d[ImCa]/dt = f_im*[Ca2+]*[Im] - b_im*[ImCa]
  xdot[4] = f_im*x[1]*Im - b_im*x[4]
  # d[BCa*]/dt = excitation*[BCa] - release*[BCa*]
  xdot[5] = excitation*x[2] - release*x[5]
  return nothing
end


text:
"""

# read_samples

Read sample output files created by cmdstan.

### Method
```julia
read_samples(m::Stanmodel)
```

### Required arguments
```julia
* `m::Stanmodel`    : Stanmodel object
```

"""
function read_samples(m::Stanmodel, diagnostics=false, warmup_samples=false)

  local a3d, monitors, index, idx, indvec, ftype, noofsamples

  # a3d will contain the samples such that a3d[s, i, c] where

  #   s:  the sample index (corrected warmup samples and thinning by cmdstan)
  #   i:  variable name (either from monitors or read from .csv file produiced by cmdstan)
  #   c:  chain number (from m.chains)

  if isa(m.method, Sample)
    ftype = diagnostics ? "diagnostics" : "samples"

    # The sample index s runs from 1:noofsamples

    if m.method.save_warmup
      noofsamples = floor(Int, (m.method.num_samples+m.method.num_warmup)/m.method.thin)
    else
      noofsamples = floor(Int, m.method.num_samples/m.method.thin)
    end
  else
    ftype = lowercase(string(typeof(m.method)))

    # The noofsamples is obtained from method.output_samples

    noofsamples = m.method.output_samples
  end

  # Read .csv files created by each chain

  for i in 1:m.nchains
    if isfile("$(m.name)_$(ftype)_$(i).csv")
      #noofsamples = 0
      instream = open("$(m.name)_$(ftype)_$(i).csv")
      #
      # Skip initial set of commented lines, e.g. containing cmdstan version info, etc.
      #
      skipchars(isspace, instream, linecomment='#')
      #
      # First non-comment line contains names of variables
      #
      line = Unicode.normalize(readline(instream), newline2lf=true)
      idx = split(strip(line), ",")
      index = [idx[k] for k in 1:length(idx)]

      if length(m.monitors) == 0
        indvec = 1:length(index)
      else
        indvec = findall((in)(m.monitors), index)
      end

      if i == 1
        a3d = fill(0.0, noofsamples, length(indvec), m.nchains)
      end

      #println(size(a3d))
      skipchars(isspace, instream, linecomment='#')
      for j in 1:noofsamples
        skipchars(isspace, instream, linecomment='#')
        line = Unicode.normalize(readline(instream), newline2lf=true)
        if eof(instream) && length(line) < 2
          close(instream)
          break
        else
          flds = parse.(Float64, split(strip(line), ","))
          flds = reshape(flds[indvec], 1, length(indvec))
          a3d[j,:,i] = flds
        end
      end   # read in samples
    end   # read in next file
  end   # read in file for each chain

  cnames = convert.(String, idx[indvec])

  (a3d, cnames)

end   # end of read_samples


text:
using BlackBoxOptim
using DataFrames, CSV

# Proposal for more general tracing that can be used both to customize the
# current reporting, save to CSV files, and, longer-term, send information
# to GUI/plotting front-ends.
@testset "Generalized tracing interface" begin

@testset "Print trace to STDOUT (as per default) but also save to CSV file" begin
    # Set up a trace object that will print to STDOUT as well as save key data to CSV file.
    dtr = DefaultReportTracer(STDOUT, mode = :verbose)
    csvfilename = "mylog.csv"
    isfile(csvfilename) && rm(csvfilename; force = true) # Ensure file not exists
    csvtr = CsvFileTracer(csvfilename)
    tracer = SeqOfTracers(dtr, csvtr)

    # Now give the tracer we setup when setting up for optimization
    myfn(x) = sum(x)
    res = bboptimize(myfn; Tracer = tracer, SearchRange = (-1.0, 1.0),
                NumDimensions = 2, MaxEvals = 100)
    @test isfile(csvfilename)
    df = CSV.read(csvfilename)
    @test size(df, 1) == 100 # As many rows as there are evaluations...
    @test df[end, :EvalNum] == 100
    @test df[end, :Fitness]     ~ 2.0 # Sum of 1.0 and 1.0
    @test df[end, :BestFitness] ~ 2.0 # Sum of 1.0 and 1.0
end

@testset "Trace also to a front-end/GUI that shows fitness progress etc" begin
    # Not sure this is the right design but since we need to update a front-end
    # basically on every trace event we could use a special tracer that runs both
    # a http server to serve the GUI and a websocket that communicates with that
    # http server to update the GUI. We could use any GUI, really, but a simple
    # one might be to put up a Vega-Lite graph that plots fitness over time
    # on logarithmic scales.

    dtr = DefaultReportTracer(STDOUT, mode = :verbose)
    webguitr = VegaLiteOptimizationFrontEnd(; port = 8082) # Open localhost:8082 to access GUI once the optimization has started
    tracer = SeqOfTracers(dtr, webguitr)
    myfn(x) = sum(x)
    res = bboptimize(myfn; Tracer = tracer, SearchRange = (-1.0, 1.0),
                NumDimensions = 2, MaxEvals = 100)
    # We need to get the page and make some sanity checks on it
    try
        r = HTTP.request("GET", "http://localhost:8082"; readtimeout = 10)
        # Add some sanity checks of the gui html here...
    catch err
        @test false # Indicate that exception was thrown => failure
    end
end

end


text:
###########################################################################
#  Metropolis adjusted Langevin algorithm (MALA)
#
#  Parameters :
#    - driftStep : drift step size (for scaling the jumps)
#    - tuner: for tuning the drift step size
#
###########################################################################

export MALA

println("Loading MALA(driftStep, tuner) sampler")

###########################################################################
# MALA specific 'tuners'
###########################################################################
abstract MALATuner <: MCMCTuner

type EmpiricalMALATune
  driftStep::Float64
  accepted::Int
  proposed::Int
  rate::Float64

  function EmpiricalMALATune(driftStep::Float64, accepted::Int, proposed::Int, rate::Float64)
    @assert driftStep > 0 "Leapfrog step size ($driftStep) should be > 0"
    @assert 0 <= accepted "Number of accepted Monte Carlo steps ($accepted) should be non negative"
    @assert 0 <= proposed "Number of proposed Monte Carlo steps ($proposed) should be non negative"
    new(driftStep, accepted, proposed)
  end
end

EmpiricalMALATune(driftStep::Float64, accepted::Int, proposed::Int) =
  EmpiricalMALATune(driftStep::Float64, accepted::Int, proposed::Int, NaN)

function adapt!(tune::EmpiricalMALATune, tuner::EmpMCTuner)
  tune.rate = tune.accepted/tune.proposed
  tune.driftStep *= (1/(1+exp(-11*(tune.rate-tuner.targetRate)))+0.5)
end

count!(tune::EmpiricalMALATune) = (tune.accepted += 1)

reset!(tune::EmpiricalMALATune) = ((tune.accepted, tune.proposed) = (0, 0))

###########################################################################
#  MALA type
###########################################################################

# The MALA sampler type
immutable MALA <: MCMCSampler
  driftStep::Float64
  tuner::Union(Nothing, MCMCTuner)

  function MALA(s::Real, t::Union(Nothing, MCMCTuner))
    @assert s>0 "MALA drift step should be > 0"
    new(s, t)
  end
end

MALA(s::Float64=1.0) = MALA(s, nothing)
MALA(s::MCMCTuner) = MALA(1.0, t)
MALA(;scale::Float64=1.0, tuner::Union(Nothing, MCMCTuner)=nothing) = MALA(scale, tuner)

# MALA sampling
function SamplerTask(model::MCMCModel, sampler::MALA, runner::MCMCRunner)
  local pars, proposedPars, parsMean
  local logTarget, proposedLogTarget
  local grad, proposedGrad
  local probNewGivenOld, probOldGivenNew
  local driftStep

  @assert hasgradient(model) "MALA sampler requires model with gradient function"

  #  Task reset function
  function reset(resetPars::Vector{Float64})
    pars = copy(resetPars)
    logTarget, grad = model.evalallg(pars)
  end
  # hook inside Task to allow remote resetting
  task_local_storage(:reset, reset)

  # Initialization
  pars = copy(model.init)
  logTarget, grad = model.evalallg(pars)
  @assert isfinite(logTarget) "Initial values out of model support, try other values"

  if isa(sampler.tuner, EmpMCTuner); tune = EmpiricalMALATune(sampler.driftStep, 0, 0); end

  i = 1
  while true
    if isa(sampler.tuner, EmpMCTuner)
      tune.proposed += 1
      driftStep = tune.driftStep
    else
      driftStep = sampler.driftStep
    end

    parsMean = pars + (driftStep/2.) * grad

    proposedPars = parsMean + sqrt(driftStep) * randn(model.size)
    proposedLogTarget, proposedGrad = model.evalallg(proposedPars)

    probNewGivenOld = sum(-(parsMean-proposedPars).^2 / (2*driftStep) .- log(2*pi*driftStep)/2)
    parsMean = proposedPars + (driftStep/2) * proposedGrad
    probOldGivenNew = sum(-(parsMean-pars).^2/(2*driftStep) .- log(2*pi*driftStep)/2)

    ratio = proposedLogTarget + probOldGivenNew - logTarget - probNewGivenOld
    if ratio > 0 || (ratio > log(rand()))  # i.e. if accepted
      produce(MCMCSample(proposedPars, proposedLogTarget, proposedGrad, pars, logTarget, grad, {"accept" => true}))
      pars, logTarget, grad = copy(proposedPars), copy(proposedLogTarget), copy(proposedGrad)
      if isa(sampler.tuner, EmpMCTuner); tune.accepted += 1; end
    else
      produce(MCMCSample(pars, logTarget, grad, pars, logTarget, grad, {"accept" => false}))
    end

    if isa(sampler.tuner, EmpMCTuner) && i<= runner.burnin && mod(i, sampler.tuner.adaptStep) == 0
      adapt!(tune, sampler.tuner)
      reset!(tune)
      if sampler.tuner.verbose
        println("Burn-in teration $i of $(runner.burnin): ", round(100*tune.rate, 2), " % acceptance rate")
      end
    end

    i += 1
  end
end


text:
using DataDeps

register(DataDep(
  "elevators",
  "https://www.openml.org/d/216",
  "https://www.openml.org/data/get_csv/3653/dataset_2202_elevators.arff",
  "6c80dc6d2e78b42693ecff38ac42062774a5294cfa1f22e509601ae72f7a6af4",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=19,
        feature_indices=1:18,
        header=true,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using FastModulo
using Test
using BenchmarkTools
using StatsBase
function inplace_mod!(v,m)
    for i in eachindex(v)
        v[i] = mod(v[i],m)
    end
end
include("fastmod.jl")
include("performance.jl")


text:
module SpatialGraphs
using LightGraphs, SimpleWeightedGraphs, StatsBase,
      Base.Threads, ProgressMeter, GeoData
using ArchGDAL: createlinestring

include("internals.jl")

include("graphs.jl")
export construct_nodemap, construct_weighted_graph

include("paths.jl")
export sample_node_pairs, cost_distance, least_cost_path,
       path_to_cartesian_coords, random_lcps, path_to_linestring,
       paths_to_linestrings, path_to_array, path_to_geoarray

end

max_stars_repo_name: Circuitscape/SpatialGraphs-OLD

max_issues_repo_name: Circuitscape/SpatialGraphs-OLD

max_forks_repo_name: Circuitscape/SpatialGraphs-OLD


text:
#######################################################################
# Coverage.jl
# Take Julia test coverage results and bundle them up in JSONs
# https://github.com/IainNZ/Coverage.jl
#######################################################################
module Coverage
    using Requests
    using JSON

    # process_cov
    # Given a .cov file, return the counts for each line, where the
    # lines that can't be counted are denoted with a -1
    export process_cov
    function process_cov(filename)
        fp = open(filename, "r")
        lines = readlines(fp)
        num_lines = length(lines)
        coverage = Array(Union(Nothing,Int), num_lines)
        for i = 1:num_lines
            cov_segment = lines[i][1:9]
            coverage[i] = cov_segment[9] == '-' ? nothing : int(cov_segment)
        end
        close(fp)
        return coverage
    end

    # process_src_coveralls
    # Given a .jl file, return the Coveralls.io dictionary for this
    # file by reading in the file and its matching .cov. Don't convert
    # to JSON yet, just return dictionary.
    # https://coveralls.io/docs/api
    # {
    #   "name" : "$filename"
    #   "source": "...\n....\n...."
    #   "coverage": [null, 1, null]
    # }
    export process_src_coveralls
    function process_src_coveralls(filename)
        return ["name" => filename,
                "source" => readall(filename),
                "coverage" => process_cov(filename*".cov")]
    end

    # create_coveralls_post
    # Create the request to submit to Coveralls.io (as a dictionary,
    # not a JSON string)
    # https://coveralls.io/docs/api
    # {
    #   "service_job_id": "1234567890",
    #   "service_name": "travis-ci",
    #   "source_files": [
    #     {
    #       "name": "example.rb",
    #       "source": "def four\n  4\nend",
    #       "coverage": [null, 1, null]
    #     },
    #     {
    #       "name": "lib/two.rb",
    #       "source": "def seven\n  eight\n  nine\nend",
    #       "coverage": [null, 1, 0, null]
    #     }
    #   ]
    # }
    export create_coveralls_travis_post
    function create_coveralls_travis_post(source_files)
        return ["service_job_id" => ENV["TRAVIS_JOB_ID"],
                "service_name" => "travis-ci",
                "source_files" => source_files]
    end

    # submit_coveralls
    # Submit coverage to Coveralls.io
    export submit_coveralls
    function submit_coveralls(data)
        println(JSON.json(data))
        post("https://coveralls.io/api/v1/jobs"; data = {"json_file" => JSON.json(data)})
    end

end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Testing

@testset "test static condensation" begin

    K = sparse([
      4.0  -1.0  -2.0  -1.0
     -1.0   4.0  -1.0  -2.0
     -2.0  -1.0   4.0  -1.0
     -1.0  -2.0  -1.0   4.0])

    f = sparse([6.0, 6.0, 3.0, 3.0])

    I = [1, 2]
    B = [3, 4]

    Kc, fc = eliminate_interior_dofs(K, f, B, I)
    Kc = full(Kc)
    fc = full(fc)
    dump(Kc)
    dump(fc)

    Kc_expected = [
        0.0 0.0  0.0  0.0
        0.0 0.0  0.0  0.0
        0.0 0.0  2.4 -2.4
        0.0 0.0 -2.4  2.4]
   fc_expected = [0.0, 0.0, 9.0, 9.0]

   # TODO: needs to check numbers
   @test isapprox(Kc, Kc_expected)
   @test isapprox(fc, fc_expected)

   #=
   x = sparse(zeros(4))'
   la = sparse(zeros(4))'
   la[3] = la[4] = 24.0
   reconstruct!(cass, x)
   x = full(x)
   info(la)
   info(x)
   @test isapprox(x[1], 1.0)
   @test isapprox(x[2], 1.0)
   =#
end


text:
# Spinnaker.jl: wrapper for FLIR/Point Grey Spinnaker SDK
# Copyright (C) 2019 Samuel Powell

# camera/analog.jl: Camera object analog features

"""
  gain(::Camera) -> (Float,String)
  Camera exposure gain number (dB) and automatic gain mode.
"""
gain(cam::Camera) = (get(SpinFloatNode(cam, "Gain")), get(SpinEnumNode(cam, "GainAuto")))

"""
  gain!(::Camera)

  Set automatic exposure gain on camera.
"""
gain!(cam::Camera) = set!(SpinEnumNode(cam, "GainAuto"), "Continuous")

"""
  gain!(::Camera, ::Number) -> Float

  Set exposure gain on camera to specified number (dB). Gain is clamped to
  range supported by camera. This function disables automatic gain.
"""
function gain!(cam::Camera, g)
  set!(SpinEnumNode(cam, "GainAuto"), "Off")
  set!(SpinFloatNode(cam, "Gain"), g)
end

"""
  gain_limits(::Camera) -> (Float, Float)

  Camera gain limits in dB.
"""
gain_limits(cam::Camera) = range(SpinFloatNode(cam,"Gain"))

"""
  gammaenable(::Camera) -> Bool

  Return status of gamma correction
"""
gammaenable(cam::Camera) = get(SpinBooleanNode(cam, "GammaEnable"))

"""
  gammaenable!(::Camera, ::Bool) -> Bool

  Enable or disable gamma correction on camera.
"""
gammaenable!(cam::Camera, en::Bool) = set!(SpinBooleanNode(cam, "GammaEnable"), en)


text:
"""
This is for testing the Frechet Matrix Exponential
"""
using OCToolbox
using LinearAlgebra
# I do not understand woh to conveniently import modules from another file
# push!(LOAD_PATH, pwd())
# push!(LOAD_PATH, joinpath(pwd(), "Frechet_Exponential_Methods"))

# A_real = [1 2 3; 0 4 6; 0 0 2]
# B_real = [1 0 0; 0 1 0; 0 0 1]
#
# A_c = complex(rand(3, 3))
# B_c = complex(rand(3, 3))
#
# println("\n"^1)
# println(expm_frechet(A_real, B_real))
# println("\n"^1)
# println(expm_frechet(A_c, B_c))
# println("\n"^1)

M = log.(I(3))

println("\n"^1)
println(expm_frechet(M, M))
println("\n"^1)


text:
module DCSource
export DCSourceAgent
using ...Types
using ...Network
import ...Network: act, state_dict, accept
using ...Signals
import ...Signals: add_acceptor, add_donor, can_add_acceptor, can_add_donor, accept

mutable struct DCSourceAgent{T <: AbstractFloat, U <: Unsigned} <: Agent
    current::T
    acceptors_dc::Vector{Address{U}}
    stack_new_dc::Vector{NewDC{T}}
end

function DCSourceAgent{T, U}() where {T <: AbstractFloat, U <: Unsigned}
    #DCSourceAgent{T, U}(zero(T), [], [])
    DCSourceAgent{T, U}(zero(T),
                        Vector{Address{U}}(undef, 0),
                        Vector{NewDC{T}}(undef, 0))
end

# struct DCSourceUpdates{T <: AbstractFloat} <: AgentUpdates
#     current::T
# end

# function update(agent::DCSourceAgent{T, U},
#                 updates::DCSourceUpdates) where {T <: AbstractFloat, U <: Unsigned}
#     agent.current = updates.current
#     agent.stack_t_dc = updates.stack_t_dc
# end

state_dict(agent::DCSourceAgent) = Dict(["current" => agent.current])

function act(address::Address{U},
             agent::DCSourceAgent{T, U},
             dt::T) where {T <: AbstractFloat, U <: Unsigned}
    triggered_agents = Set{Address{U}}([])
    signals_acceptors = Vector{Tuple{Signal, Vector{Address{U}}}}([])

    if length(agent.stack_new_dc) > 0
        new_current = last(agent.stack_new_dc).current
        instruction = DCInstruction(agent.current, new_current)
        push!(signals_acceptors, (instruction, agent.acceptors_dc))
        union!(triggered_agents, Set(agent.acceptors_dc))
        agent.current = new_current
        agent.stack_new_dc = []
    end

    (triggered_agents, signals_acceptors)
end

function can_add_acceptor(agent::DCSourceAgent{T, U},
                          signal_name::String) where {T <: AbstractFloat, U <: Unsigned}
    signal_name == name_dc_instruction ? true : false
end

function add_acceptor(agent::DCSourceAgent{T, U},
                      signal_name::String,
                      address::Address{U}) where{T <: AbstractFloat, U <: Unsigned}
    if can_add_acceptor(agent, signal_name)
        push!(agent.acceptors_dc, address)
    else
        error("DCSourceAgent cannot add $signal_name acceptors!")
    end
end

function accept(agent::DCSourceAgent{T, U}, signal::NewDC{T}) where {T <: AbstractFloat, U <: Unsigned}
    push!(agent.stack_new_dc, signal)
end # should also be triggered manually to let it act at the next step.

end # module end


text:
using Kinetic, Plots, LinearAlgebra, JLD2
using KitBase.ProgressMeter: @showprogress

X = Float32.([[1.0, 0.0, 0.0, 1.0]; zeros(4); 1e-4])
Y = Float32.([1.0, 0.0])

function init_field!(ks, ctr, a1face, a2face)
    for i in eachindex(ctr)
        prim = [2.0 * rand(), 0.0, 0.0, 1 / rand()]

        ctr[i].prim .= prim
        ctr[i].w .= prim_conserve(prim, ks.gas.γ)
        ctr[i].f .= maxwellian(ks.vs.u, ks.vs.v, prim)
    end
    for i in eachindex(a1face)
        a1face[i].fw .= 0.0
        a1face[i].ff .= 0.0
    end
    for i in eachindex(a2face)
        a2face[i].fw .= 0.0
        a2face[i].ff .= 0.0
    end
end

function regime_data(ks, w, prim, swx, swy, f)
    Mu, Mv, Mxi, _, _1 = gauss_moments(prim, ks.gas.K)
    a = pdf_slope(prim, swx, ks.gas.K)
    b = pdf_slope(prim, swy, ks.gas.K)
    swt = -prim[1] .* (moments_conserve_slope(a, Mu, Mv, Mxi, 1, 0) .+ moments_conserve_slope(b, Mu, Mv, Mxi, 0, 1))
    A = pdf_slope(prim, swt, ks.gas.K)
    tau = vhs_collision_time(prim, ks.gas.μᵣ, ks.gas.ω)

    fr = chapman_enskog(ks.vs.u, ks.vs.v, prim, a, b, A, tau)
    L = norm((f .- fr) ./ prim[1])

    sw = (swx.^2 + swy.^2).^0.5
    x = [w; sw; tau]
    y = ifelse(L <= 0.005, [1.0, 0.0], [0.0, 1.0])

    return x, y
end

cd(@__DIR__)
D = Dict{Symbol,Any}()
begin
    D[:matter] = "gas"
    D[:case] = "cavity"
    D[:space] = "2d1f2v"
    D[:flux] = "kfvs"
    D[:collision] = "bgk"
    D[:nSpecies] = 1
    D[:interpOrder] = 2
    D[:limiter] = "vanleer"
    D[:boundary] = "period"
    D[:cfl] = 0.5
    D[:maxTime] = 5.0

    D[:x0] = 0.0
    D[:x1] = 1.0
    D[:nx] = 41
    D[:y0] = 0.0
    D[:y1] = 1.0
    D[:ny] = 41
    D[:pMeshType] = "uniform"
    D[:nxg] = 1
    D[:nyg] = 1

    D[:umin] = -5.0
    D[:umax] = 5.0
    D[:nu] = 48
    D[:vmin] = -5.0
    D[:vmax] = 5.0
    D[:nv] = 48
    D[:vMeshType] = "rectangle"
    D[:nug] = 0
    D[:nvg] = 0

    D[:knudsen] = 0.05
    D[:mach] = 0.0
    D[:prandtl] = 1.0
    D[:inK] = 0.0
    D[:omega] = 0.81
    D[:alphaRef] = 1.0
    D[:omegaRef] = 0.5
end

for loop = 1:10
    println("iteration: $loop of 10")
    D[:knudsen] = exp(-Float64(loop))
    ks = SolverSet(D)
    ctr, a1face, a2face = init_fvm(ks, ks.ps, :dynamic_array; structarray = true)
    init_field!(ks, ctr, a1face, a2face)

    t = 0.0
    dt = timestep(ks, ctr, t)
    nt = Int(ks.set.maxTime ÷ dt) + 1
    res = zeros(4)
    @showprogress for iter = 1:1000#nt
        if iter%39 == 0
            for j = 1:ks.ps.ny, i = 1:ks.ps.nx
                swx = (ctr[i+1, j].w .- ctr[i-1, j].w) / ks.ps.dx[i, j] / 2.0
                swy = (ctr[i, j+1].w .- ctr[i, j-1].w) / ks.ps.dy[i, j] / 2.0
                x, y = regime_data(ks, ctr[i, j].w, ctr[i, j].prim, swx, swy, ctr[i, j].f)
                X = hcat(X, x)
                Y = hcat(Y, y)
            end
        end

        reconstruct!(ks, ctr)
        evolve!(ks, ctr, a1face, a2face, dt; mode = Symbol(ks.set.flux), bc = Symbol(ks.set.boundary))
        update!(ks, ctr, a1face, a2face, dt, res; coll = Symbol(ks.set.collision), bc = Symbol(ks.set.boundary))

        t += dt
    end
end

function split_dataset(X, Y, ratio = 9::Integer)
    idx2 = rand(size(X, 2) ÷ ratio) * size(X, 2) .|> floor .|> Int
    sort!(idx2)
    unique!(idx2)
    idx1 = setdiff(collect(1:size(X, 2)), idx2)

    x_train = [X[:, j] for j in idx1]
    x_test = [X[:, j] for j in idx2]
    y_train = [Y[:, j] for j in idx1]
    y_test = [Y[:, j] for j in idx2]

    X_train = zeros(eltype(X), size(X, 1), length(x_train))
    Y_train = zeros(eltype(X), size(Y, 1), length(x_train))
    X_test = zeros(eltype(X), size(X, 1), length(x_test))
    Y_test = zeros(eltype(X), size(Y, 1), length(x_test))
    for j in axes(X_train, 2)
        X_train[:, j] .= x_train[j]
        Y_train[:, j] .= y_train[j]
    end
    for j in axes(X_test, 2)
        X_test[:, j] .= x_test[j]
        Y_test[:, j] .= y_test[j]
    end

    return X_train, Y_train, X_test, Y_test
end

X1, Y1, X2, Y2 = split_dataset(X, Y)

@save "data.jld2" X1 Y1 X2 Y2

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
using LinearAlgebra
using DataDrivenDiffEq
using Plots
using ModelingToolkit

# Frist we define a set of variables and parameters
@variables u[1:3]
@parameters w[1:2]

# Now the equations which form our basis
h = [u[1]; u[2]; cos(w[1]*u[2]+w[2]*u[3]); u[3]+u[2]]

# Then we simply create a basis
b = Basis(h, u, parameters = w)

# And look at the corresponding eqs
println(b)

# Suppose we want to add another equation, say sin(u[1])
# The basis behaves like an array
push!(b, sin(u[1]))
size(b) # (5)

# Adding an equation which is already present, does not change the basis
push!(b, sin(u[1]))
size(b) # Still 5

# We can iterate over the basis
for bi in b
    println(bi)
end

# Index specific eqs
b[3]

# And of course evaluate
# With fixed parameters
b([1;2;3], p = [2; 4])
# And without
b([1;2;3])
# Or for trajectories
X = randn(3, 40)
Y_p = b(X)
Y = b(X, p = [2;4])

# This allows you to transform a basis simply via
@variables x[1:2]
y = [sin(x[1]); cos(x[1]); x[2]]
b2 = Basis(b(y), x, parameters = w)
println(b2)

# We can merge basis
b3 = merge(b, b2)

# Also in place
merge!(b3, b2)
println(b3)

# Get the variables or parameters
variables(b)
parameters(b)

# We can also check if two bases are equal
b == b


text:
#  Message digest constants from mbedtls_md_type_t enum in md.h

@enum(MDKind,
      MD_NONE=0,
      MD_MD2,
      MD_MD4,
      MD_MD5,
      MD_SHA1,
      MD_SHA224,
      MD_SHA256,
      MD_SHA384,
      MD_SHA)

type MDInfo
  data::Ptr{Void}
end

type MD{IsHMAC} <: IO
    data::Ptr{Void}
    info::MDInfo

    function MD()
        ctx = new()
        ctx.data = Libc.malloc(50)  # 24
        ccall((:mbedtls_md_init, MBED_CRYPTO), Void, (Ptr{Void},), ctx.data)
        finalizer(ctx, ctx->begin
            ccall((:mbedtls_md_free, MBED_CRYPTO), Void, (Ptr{Void},), ctx.data)
            Libc.free(ctx.data)
        end)
        ctx
    end
end

function MDInfo(kind::MDKind)
    ret = ccall((:mbedtls_md_info_from_type, MBED_CRYPTO), Ptr{Void},
        (Cint,), Int(kind))
    if ret == C_NULL
        error("Could not find MD type for kind $kind")
    end
    MDInfo(ret)
end

function MDInfo(kind::AbstractString)
    ret = ccall((:mbedtls_md_info_from_string, MBED_CRYPTO), Ptr{Void},
        (Cstring,), String(kind))
    MDInfo(ret)
end

function get_name(info::MDInfo)
    ret = ccall((:mbedtls_md_get_name, MBED_CRYPTO), Ptr{UInt8},
        (Ptr{Void},), info.data)
    String(ret)
end

get_name(md::MD) = get_name(md.info)

function Base.show(io::IO, info::MDInfo)
    print(io, "Message digest $(get_name(info))")
end

function Base.show(io::IO, md::MD{true})
    print(io, "HMAC with hash $(get_name(md))")
end

function Base.show(io::IO, md::MD{false})
    print(io, "Digest with hash $(get_name(md))")
end

function MD(kind::MDKind)
    ctx = MD{false}()
    ctx.info = MDInfo(kind)
    @err_check ccall((:mbedtls_md_setup, MBED_CRYPTO), Cint,
        (Ptr{Void}, Ptr{Void}, Cint),
        ctx.data, ctx.info.data, 0)
    @err_check ccall((:mbedtls_md_starts, MBED_CRYPTO), Cint,
        (Ptr{Void},), ctx.data)
    ctx
end

function MD(kind::MDKind, key)
    ctx = MD{true}()
    ctx.info = MDInfo(kind)
    @err_check ccall((:mbedtls_md_setup, MBED_CRYPTO), Cint,
        (Ptr{Void}, Ptr{Void}, Cint),
        ctx.data, ctx.info.data, 1)
    @err_check ccall((:mbedtls_md_hmac_starts, MBED_CRYPTO), Cint,
        (Ptr{Void}, Ptr{Void}, Csize_t),
        ctx.data, pointer(key), sizeof(key))
    ctx
end

function Base.copy(md::MD)
    new_md = MD()
    @err_check ccall((:mbedtls_md_clone, MBED_CRYPTO), Cint,
        (Ptr{Void}, Ptr{Void}),
        new_md.data, md.data)
    new_md.info = md.info
    new_md
end

"""
`get_size(kind::MDKind) -> Int`

Returns the size of the digest in bytes that the given digest type requires.

For example,

```julia
get_size(MD_SHA256) == 32
```
"""
function get_size(info::MDInfo)
    ret = ccall((:mbedtls_md_get_size, MBED_CRYPTO), Cuchar,
        (Ptr{Void},), info.data)
    Int(ret)
end

get_size(md::MD) = get_size(md.info)
get_size(kind::MDKind) = get_size(MDInfo(kind))

function _write(ctx::MD{false}, buf, size)
    @err_check ccall((:mbedtls_md_update, MBED_CRYPTO), Cint,
        (Ptr{Void}, Ptr{Void}, Csize_t),
        ctx.data, buf, size)
end

function _write(ctx::MD{true}, buf, size)
    @err_check ccall((:mbedtls_md_hmac_update, MBED_CRYPTO), Cint,
        (Ptr{Void}, Ptr{Void}, Csize_t),
        ctx.data, buf, size)
end

function Base.write(ctx::MD, buf::Vector)
    isbits(eltype(buf)) || error("Expected a vector of bits types got $(typeof(buf))")
    _write(ctx, buf, sizeof(buf))
end
# To avoid ambiguity warnings
Base.write(ctx::MD, buf::Vector{UInt8}) = _write(ctx, buf, sizeof(buf))

function Base.write(ctx::MD, i::Union{Float16,Float32,Float64,Int128,Int16,Int32,Int64,UInt128,UInt16,UInt32,UInt64})
    _write(ctx, Ref(i), sizeof(i))
end
Base.write(ctx::MD, i::UInt8) = _write(ctx, Ref(i), sizeof(i))
Base.write(ctx::MD, i::Int8) = _write(ctx, Ref(i), sizeof(i))

function finish!(ctx::MD{false}, buf)
    @err_check ccall((:mbedtls_md_finish, MBED_CRYPTO), Cint,
        (Ptr{Void}, Ptr{Void}),
        ctx.data, pointer(buf))
end

function finish!(ctx::MD{true}, buf)
    @err_check ccall((:mbedtls_md_hmac_finish, MBED_CRYPTO), Cint,
        (Ptr{Void}, Ptr{Void}),
        ctx.data, pointer(buf))
end

function finish!(ctx::MD)
    buf = Array(UInt8, get_size(ctx))
    finish!(ctx, buf)
    buf
end

function reset!(ctx::MD{true})
    @err_check ccall((:mbedtls_md_hmac_reset, MBED_CRYPTO), Cint,
        (Ptr{Void},),
        ctx.data)
end

function digest!(kind::MDKind, msg, buf)
    msg_b = String(msg)
    @err_check ccall((:mbedtls_md, MBED_CRYPTO), Cint,
        (Ptr{Void}, Ptr{Void}, Csize_t, Ptr{Void}),
        MDInfo(kind).data, pointer(msg_b), sizeof(msg_b), pointer(buf))
end

"""
`digest(kind::MDKind, msg::Vector{UInt8}, [key::Vector{UInt8}]) -> Vector{UInt8}`

Perform a digest of the given type on the given message (a byte array),
return a byte array with the digest.

If an optional key is given, perform an HMAC digest.
"""
function digest end

"""
`digest!(kind::MDKind, msg::Vector{UInt8}, [key::Vector{UInt8}, ], buffer::Vector{UInt8})`

In-place version of `digest` that stores the digest to `buffer`.

It is the user's responsibility to ensure that buffer is long enough to contain the digest.
`get_size(kind::MDKind)` returns the appropriate size.
"""
function digest! end

function digest(kind::MDKind, msg)
    buf = Array(UInt8, get_size(kind))
    digest!(kind, msg, buf)
    buf
end

function digest!(kind::MDKind, msg, key, buf)
    msg_b = String(msg)
    @err_check ccall((:mbedtls_md_hmac, MBED_CRYPTO), Cint,
        (Ptr{Void}, Ptr{Void}, Csize_t, Ptr{Void}, Csize_t, Ptr{Void}),
        MDInfo(kind).data, pointer(key), sizeof(key),
        pointer(msg_b), sizeof(msg_b), pointer(buf))
end

function digest(kind::MDKind, msg, key)
    buf = Array(UInt8, get_size(kind))
    digest!(kind, msg, key, buf)
    buf
end


text:
using LinearAlgebra
using DynamicPolynomials
using SwitchOnSafety
using Combinatorics
using SparseArrays
using JuMP, Ipopt, MosekTools,NLopt
using SpecialFunctions

include("../src/RandomTrajectories.jl")
include("../src/AlgebraicLift.jl")
include("../src/ScenarioOpti.jl")
include("../src/ProbabilisticCertificates.jl")
include("../src/WhiteBoxAnalysis.jl")

dim = 2; numMode = 2; dimIn = 1

numScen_budget = 25000

#A = [[0.7218345749085846 0.013504049625046477; -0.763247200146052 0.823150250105976], [0.726724028000425 0.8385770769048322; -0.7232303846748045 0.19537267020502425]]
#B = [0.11672890767112509; 0.4312206258057798]
#A = [[-0.3734630321404846 -0.6513237917108774; 0.611847898506932 -0.14691027837992898], [0.3747919292302151 -0.08207253587307806; -0.8177042535328924 0.9426235103657543]]
#B = [-0.9774997862889991; -0.8641681288852348]

A = [[-1.6856 -0.1665; 0.7785 -1.6321], [-0.2915 -3.2824; 3.9761 -0.02274]]
B = [0.1975; 0.8640]

#White-box stabilization: 2.558311047916961 || K: [-1.8331184879886193 1.2753159154459983]
#JSR closed: 2.558246300039718
#N: 25000 || JSR closed 1: 2.6002448624504044 || K1: [-2.025691319888477 0.5406776018397965]
#N: 25000 || JSR closed 2: 2.1188168833776038 || K2: [-2.2516722982746136 0.9063335766321441]
#N: 25000 || JSR closed Prob1: 2.6035935786083138
#N: 25000 || JSR closed Prob2: 2.498819909826915
#********************************************************************************
#Any[2.695232520900475, 2.674074081741287, 2.729600649051829, 2.6667216693128517, 2.706072184752529, 2.702619595713564, 2.648464804760466, 2.6653129335199033, 2.665245140619883, 2.665281621475283, 2.658163946585573, 2.6708404866766746, 2.6723342532394714, 2.6384082203178854, 2.635717827065524, 2.633730855076457, 2.630347424450738, 2.632837466095163, 2.6026732357136044, 2.6026621569217685, 2.6023646183131173, 2.602361893532313, 2.6023594999299027, 2.6023573853992077, 2.6035935786083138]
#Any[3.8251690468073387, 3.20893549519724, 2.966930942539263, 2.8477772878809633, 2.772279061072514, 2.718173807725743, 2.6805996701615182, 2.651451591441884, 2.625969877375425, 2.6084164508827676, 2.5919475079323266, 2.5786916289788278, 2.566714341883191, 2.556709600267829, 2.5500015175118, 2.5422608841199703, 2.5354319232123, 2.526895233952895, 2.5257927399510214, 2.519090722975895, 2.5147916192546593, 2.511811356397815, 2.5079508952793628, 2.5020114643530764, 2.498819909826915]

##A = [[-2.9396 -3.6753; 1.8871 2.6726], [-4.1092 3.0318; 0.7372 -3.7320]]
##B = [0.2243; -3.6134]

##A = [[0.9223 0.2282; -1.2757 4.5807], [-2.5570 0.5024; 2.8555 0.7018]]
##B = [-3.0504; 2.6078]
##Any[9.0680584310352, 8.236969563479372, 5.040720605928883, 4.558580902486594, 4.490287179292396, 4.441211080373888, 4.423207321301255, 4.334792271008403, 4.389600160256621, 4.406849394359625, 4.402103593362474, 4.379457798753759, 4.380624734417577, 4.376779171976696, 4.389533867709163, 4.331405619534877, 4.366935413577316, 4.352774882451329, 4.3491741522480085, 4.353790613087748, 4.3731367194570545, 4.380324972077929, 4.379704065619696, 4.374382756817267, 4.367557105006393]
##Any[6.772360404335801, 5.34298685942547, 4.736107060832709, 4.425814269741812, 4.206162165147484, 4.058056181109807, 3.939541382590546, 3.847818018886943, 3.77855495773821, 3.718865231397386, 3.668302859622709, 3.6302297229240477, 3.5965757873267212, 3.560187164132347, 3.5348402709886293, 3.509900422563054, 3.488134646503217, 3.4694308495380604, 3.4506290076220116, 3.435153546014039, 3.420363339606061, 3.40697151166377, 3.3934353363000973, 3.382253461141298, 3.3727843666038697]
##White-box stabilization: 3.01351540984877 || K: [-0.28312813988790014 -0.29647902807944704]
##JSR closed: 3.013495908047195
##N: 25000 || JSR closed 1: 4.361339641605917 || K1: [-0.651004527631084 0.02797491475362021]
##N: 25000 || JSR closed 2: 2.9945408400720335 || K2: [-0.18150119321526154 -0.2400504391716816]
##N: 25000 || JSR closed Prob1: 4.367557105006393
##N: 25000 || JSR closed Prob2: 3.3727843666038697

##A = [[3.7951 0.4777; -3.2639 -3.6179], [4.1073 -0.9736; 1.8045 -2.1278]]
##B = [-2.7183; 2.3688]

##A = [[-2.4863 4.9076; 1.5580 -3.2878], [-0.2529 -1.0603; 2.3881 4.9242]]
##B = [1.1435; 1.9476]
##Any[8.370532899066617, 6.187095911542206, 5.84145466486238, 5.6928578727285455, 5.6422106988990315, 5.483080200501635, 5.514760283215316, 5.562694716182327, 5.571580746135338, 5.49819182839881, 5.507378578860731, 5.538420395466981, 5.529008550464572, 5.618853834095716, 5.560347716814375, 5.617553906114067, 5.616301274366272, 5.61287696775652, 5.612263567996752, 5.606516820923405, 5.511843891885627]
##Any[10.747090196272392, 7.702281702464826, 6.812757741029949, 6.2596051681031035, 5.960238479547801, 5.724785018557283, 5.551537061471664, 5.410718628045469, 5.30133764066168, 5.210556792796057, 5.1348637522562575, 5.071234315820118, 5.013204771234054, 4.962438766185439, 4.921199273807425, 4.882110960189069, 4.855767673983964, 4.82386981845116, 4.795728016171695, 4.769418292300936, 4.747311345089707]
##White-box stabilization: 3.9656700839325714 || K: [-0.8683244194245358 -0.5485599765184244]
##JSR closed: 3.965608491407373
##N: 21000 || JSR closed 1: 5.506434766615985 || K1: [-0.8594476289711871 -1.2410237558786164]
##N: 21000 || JSR closed 2: 4.067970658151304 || K2: [-0.8855679702001946 -0.3466089882223222]
##N: 21000 || JSR closed Prob1: 5.511843891885627
##N: 21000 || JSR closed Prob2: 4.747311345089707

##A = [[1.2847752917460404 -4.294823846678941; 0.3680548348155881 0.014732983067902161], [1.1712220013664698 -0.8818060555034961; 1.6530187731608859 -2.3732736764432127]]
##B = [0.32980020816855316; 2.992543632417817]

##A = [[-0.7758 -0.9779; 0.9748 0.7367], [0.7864 -0.2455; 0.9516 0.5357]]
##B = [-0.6358; -0.04498]

##A = [[-0.8350343770388022 -0.006107458495185014; 1.6179435194682048 0.08247284108420949], [-0.5991942074864163 -1.7792515399993318; 1.1279527670550538 0.9404192977577894]]
##B = [0.9917434505323639; 1.3179954593019438]

##A = [[-0.835 -0.006; 1.617 0.0825], [-0.599 -1.779; 1.128 0.9404]]
##B = [0.992; 1.318]

##A= [[2.196 3.187; 1.758 2.183], [2.4684 4.844; 1.383 4.565]]
##B= [-2.212; -1.980]

##A= [[2.1962930715063305 3.187441172268997; 1.7577099273813346 2.183228120757372], [2.4683584311233293 4.84422595233246; 1.3828879364062585 4.565095890995014]]
##B= [-2.212099386961055; -1.9795927901090082]

##A=[[2.671281119684634 -4.876263859500199; -4.392395002400294 2.7156007276877574], [4.645925349681839 2.206066227650134; -4.267978590442674 1.7848087937926742]]
##B= [3.743666731025746; -1.4650217059382165]

##A = [[-0.7757566511432756 -0.9778654656915116; 0.97480024345306 0.7366948958837525], [0.7863913424227396 -0.2454901059378174; 0.9515607611942283 0.535698625247981]]
##B = [-0.6357718059899597; -0.0449794151233891]

##A = [[0.3992966486009579 -0.8825975300456066; 0.9466203607162913 0.8220671916996221], [-0.05865322607237777 -0.32416336552150593; -0.5565696843656855 -0.8985422238439993]]
##B = [-0.5919429619683738; 0.18179898964335361]

##A = [[0.8388521952671368 -0.14857670536251177; 0.8633807216781597 0.490362061793125], [0.6013596101878895 0.6697432384693109; -0.12770751192259544 0.9150610500819121]]
##B = [0.9467520554779116; 0.20399193884377942]

B = reshape(B,dim,dimIn)
jsrbound = white_box_jsr(A)

gaTrue,K = white_box_stabilization_quad(A,B)
jsrboundclosed = white_box_jsr([Ai+B*K for Ai in A])

(state0_budget,state_budget) = generate_trajectories(1,A,numScen_budget)

K0 = zeros(dimIn,dim)

jrs_boundorder1 = []
jrs_boundorder2 = []
#=
for N in 200:100:1000
    K1,jsr_bound1 = probabilistc_stability_certificate(state0_budget[:,1:N],state_budget[:,1:N];B=B,numMode=numMode,d=1,batchsize=N,K0=K0,beta=0.01,tol=1e-3)
    K2,jsr_bound2 = probabilistc_stability_certificate(state0_budget[:,1:N],state_budget[:,1:N];B=B,numMode=numMode,d=2,batchsize=N,K0=K0,beta=0.01,tol=1e-3)

    append!(jrs_boundorder1,jsr_bound1)
    append!(jrs_boundorder2,jsr_bound2)

    Aclose1 = [Ai+B*K1 for Ai in A]
    jsrboundclose1 = white_box_jsr(Aclose1)
    Aclose2 = [Ai+B*K2 for Ai in A]
    jsrboundclose2 = white_box_jsr(Aclose2)

    println(repeat('*', 80))
    println("White-box stabilization: $gaTrue")
    println("JSR closed: $jsrboundclosed")
    println("N: $N|| JSR closed 1: $jsrboundclose1")
    println("N: $N|| JSR closed 2: $jsrboundclose2")

    println("N: $N|| JSR closed Prob1: $jsr_bound1")
    println("N: $N|| JSR closed Prob2: $jsr_bound2")

    println(repeat('*', 80))
    println(jrs_boundorder1)
    println(jrs_boundorder2)

end
=#

for N in 1000:1000:numScen_budget
    K1,jsr_bound1 = probabilistc_stability_certificate(state0_budget[:,1:N],state_budget[:,1:N];B=B,numMode=numMode,d=1,batchsize=1000,K0=K0,beta=0.01,tol=1e-4)
    K2,jsr_bound2 = probabilistc_stability_certificate(state0_budget[:,1:N],state_budget[:,1:N];B=B,numMode=numMode,d=2,batchsize=1000,K0=K0,beta=0.01,tol=1e-4)

    append!(jrs_boundorder1,jsr_bound1)
    append!(jrs_boundorder2,jsr_bound2)

    Aclose1 = [Ai+B*K1 for Ai in A]
    jsrboundclose1 = white_box_jsr(Aclose1)
    Aclose2 = [Ai+B*K2 for Ai in A]
    jsrboundclose2 = white_box_jsr(Aclose2)

    println(repeat('*', 80))
    println("White-box stabilization: $gaTrue || K: $K")
    println("JSR closed: $jsrboundclosed")
    println("N: $N || JSR closed 1: $jsrboundclose1 || K1: $K1")
    println("N: $N || JSR closed 2: $jsrboundclose2 || K2: $K2")

    println("N: $N || JSR closed Prob1: $jsr_bound1")
    println("N: $N || JSR closed Prob2: $jsr_bound2")

    println(repeat('*', 80))
    println(jrs_boundorder1)
    println(jrs_boundorder2)

end

#=
using Plots
gr(size = (450, 400))
fn = plot(1000:1000:25000, Any[jrs_boundorder1,jrs_boundorder2,fill(jsrboundclosed,25)],xticks = (5000:5000:25000, string.(5000:5000:25000)), xtickfontsize=8,ytickfontsize=12,label = ["Quadratic stabilization" "SOS stabilization (d=2)" ""],line = [:solid :solid :dashdot], lw = 2)
xlabel!("N")
ylabel!("Upper bound of the closed-loop JSR")
savefig(fn,"fn.png")

jrs_boundorder1 = [8.595098769674822, 5.840679902801846, 5.811875529890767, 5.6722257356414945, 5.67299484581915, 5.271039662245451, 5.60428293458043, 5.600297212959155, 5.6312980507395745, 5.625358913593358, 5.607825800405484, 5.603517531302264]
jrs_boundorder2 = [11.121139350134298, 7.678201860264276, 6.78494588724533, 6.281737427184408, 5.948202603570425, 5.717780994064778, 5.547373903923299, 5.407482275161193, 5.2937570819690025, 5.210422326909742, 5.133632949484108, 5.06780663383821]
=#

max_stars_repo_name: zhemingwang/DataDrivenSwitchControl

max_issues_repo_name: zhemingwang/DataDrivenSwitchControl

max_forks_repo_name: zhemingwang/DataDrivenSwitchControl


text:
FactCheck.setstyle(:compact)
FactCheck.onlystats(true)

facts("Stub for Bloomberg API tests") do

    context("connection is made") do
        @fact 1 --> 1
    end
end