text:
using ChaosTools
using Test, StaticArrays, OrdinaryDiffEq, LinearAlgebra

println("\nTesting poincare sections...")

@testset "Poincare SOS" begin
    @testset "Henon Heiles" begin
        ds = Systems.henonheiles([0, .483000, .278980390, 0] )
        psos = poincaresos(ds, (2, 0.0), 1000.0, direction = +1)
        xcross = psos[:, 2]
        @test length(xcross) > 1
        for x in xcross
            @test abs(x) < 1e-3
        end

        # @inline Vhh(q1, q2) = 1//2 * (q1^2 + q2^2 + 2q1^2 * q2 - 2//3 * q2^3)
        # @inline Thh(p1, p2) = 1//2 * (p1^2 + p2^2)
        # @inline Hhh(q1, q2, p1, p2) = Thh(p1, p2) + Vhh(q1, q2)
        # @inline Hhh(u::AbstractVector) = Hhh(u...)
        #
        # E = [Hhh(p) for p in psos]
        #
        # @test std(E) < 1e-10
        # @test maximum(@. E - E[1]) < 1e-10
    end
    @testset "Gissinger crazy plane" begin
        gis = Systems.gissinger([2.32865, 2.02514, 1.98312])
        # Define appropriate hyperplane for gissinger system
        ν = 0.1
        Γ = 0.9 # default parameters of the system

        # I want hyperperplane defined by these two points:
        Np(μ) = SVector{3}(sqrt(ν + Γ*sqrt(ν/μ)), -sqrt(μ + Γ*sqrt(μ/ν)), -sqrt(μ*ν))
        Nm(μ) = SVector{3}(-sqrt(ν + Γ*sqrt(ν/μ)), sqrt(μ + Γ*sqrt(μ/ν)), -sqrt(μ*ν))

        # Create hyperplane using normal vector to vector connecting points:
        gis_plane(μ) = [cross(Np(μ), Nm(μ))..., 0]

        μ = 0.12
        set_parameter!(gis, 1, μ)
        psos = poincaresos(gis, gis_plane(μ), 5000.0, Ttr = 200.0, direction = +1)
        @test length(psos) > 1
        ε = estimate_boxsizes(psos; w = 0, z = -1)
        @test generalized_dim(psos, ε; q = 2) < 1
    end
    @testset "beginning on the plane" begin

        @inbounds @inline function ż(z, p, t)
            A, B, D = p
            p₀, p₂ = z[1:2]
            q₀, q₂ = z[3:4]

            return SVector{4}(
                -A * q₀ - 3 * B / √2 * (q₂^2 - q₀^2) - D * q₀ * (q₀^2 + q₂^2),
                -q₂ * (A + 3 * √2 * B * q₀ + D * (q₀^2 + q₂^2)),
                A * p₀,
                A * p₂
            )
        end
        z0 = SVector{4}(2.3499921651423565, -9.801351029039825, 0.0, -2.7230316965872268)

        ds = ContinuousDynamicalSystem(ż, z0, (A=1, B=0.55, D=0.4))
        psos = poincaresos(ds, (3, 0), 10., direction=+1, idxs=[2,4])
        @test length(psos) > 0
        @test size(psos)[2] == 2
    end
    @testset "Henon Heiles iterated Poincaré map" begin
        ds = Systems.henonheiles([0, 0.483000, 0.278980390, 0] )
        pmap = poincaremap(ds, (2, 0.0); direction=+1)
        next = step!(pmap)
        @test abs(next[2]) < 1e-6
        # Test that we can re-init
        reinit!(pmap, ds.u0)
        next2 = step!(pmap)
        @test !isnothing(next2)
        @test abs(next2[2]) < 1e-6
        @test next[1] ≈ next2[1] atol = 1e-3
        # Test an ill defined hyperplane
        pmap = poincaremap(ds, (2, 5.0), 100.0; direction=+1)
        @test isnothing(step!(pmap))
    end
end

@testset "Produce OD" begin
  @testset "Shinriki" begin
    ds = Systems.shinriki([-2, 0, 0.2])

    pvalues = range(19,stop=22,length=11)
    parameter = 1
    i = 1
    j = 2
    tf = 200.0

    de = (abstol=1e-9, reltol = 1e-9)
    output = produce_orbitdiagram(ds, (j, 0.0), i, parameter, pvalues; tfinal = tf,
    Ttr = 100.0, printparams = false, direction = +1, de...)
    @test length(output) == length(pvalues)

    v = round.(output[1], digits = 4)
    s = collect(Set(v))
    @test length(s) == 1
    @test s[1] == -0.856

    v = round.(output[4], digits = 4)
    s = Set(v)
    @test length(s) == 2
    @test s == Set([-0.376, -1.2887])

  end
end

@testset "PoincareSOS of trajectory" begin
    ds = Systems.lorenz()
    tr = trajectory(ds, 1000.0; Ttr = 10.0)
    psos = poincaresos(tr, (2, 0.0))
    @test all(abs.(psos[:, 2]) .≤ 1e-15)
    psos = poincaresos(tr, [0, 1.0, 0, 0])
    @test all(abs.(psos[:, 2]) .≤ 1e-15)
    psos = poincaresos(tr, [1.0, 1.0, 0, 0])
    g = generalized_dim(psos; q = 0)
    @test g ≤ 1
end


text:
###test a Data object with one variable, and generating Float32 values
d0 = data(:array,0.0f0:0.1f0:1.0f0,zeros(Float32,11))
@test numvalues(d0) == 11
@test numvars(d0) == 1
@test generate!(d0).values == zeros(Float32,11)
@test dataindex(d0) == collect(0.0f0:0.1f0:1.0f0)
@test datavalues(d0) == zeros(Float32,11)
@test eltype(d0) == Float32

###measurment time
t1 = 0.0:0.1:1.0

###data generating functions
f1(t,a) = hcat(sin(a*t),cos(a*t))
f1!(r,t,a) = (for i=1:length(t) r[i,:] = [sin(a*t[i]),cos(a*t[i])] end ; r)
f2!(r,μ::Float64,Σ::Float64) = (d = Distributions.Normal(μ,Σ) ; rand!(d,r) ; r)

###data array
a1 = hcat(sin(2*pi*t1),cos(2*pi*t1))

###create the Data objects
d1 = data(:array,t1,a1)
d2 = data(:function,t1,f1,t1,2π)
d2! = data(:function!,t1,zeros(length(t1),2),f1!,t1,2π)
d3! = data(:function!,t1,zeros(length(t1),3),f2!,1.0,0.1)
d4! = data(:function!,t1,zeros(length(t1),3),rand!)

@test numvalues(d1) == numvalues(d2) == numvalues(d2!) == numvalues(d3!) == numvalues(d4!) == 11
@test numvars(d1) == numvars(d2) == numvars(d2!) == 2 && numvars(d3!) == numvars(d4!) == 3
@test datavalues(d1) == datavalues(d2) == datavalues(d2!) == a1
@test datavalues(generate!(d1)) == datavalues(generate!(d2)) == datavalues(generate!(d2!)) == a1
@test dataindex(d1) == dataindex(d2) == dataindex(d2!) == dataindex(d3!) == dataindex(d4!) == collect(t1)
@test eltype(d1) == eltype(d2) == eltype(d2!) == eltype(d3!) == eltype(d4!) == Float64
@test (srand(5643) ; datavalues(generate!(d3!))) == (srand(5643) ; rand(Distributions.Normal(1.0,0.1),11,3)) == datavalues(d3!)
@test (srand(5643) ; datavalues(generate!(d4!))) == (srand(5643) ; rand(11,3)) == datavalues(d4!)

###test obtaining a new copy of the datavector
e1 = datavalues(d1) ; e1[1] = 10.0
e2 = datavalues(d2) ; e2[1] = 10.0
@test e1 == e2 && e1 == datavalues(d1) && e2 == datavalues(d2)
@test e1 == datavalues(generate!(d1)) && e2 == datavalues(generate!(d2))

###noise model
s1 = [0.1,0.2]
n1 = noise(:gaussian,s1)

###apply a noise model to model data
c1 = datavalues(d1)
c2 = datavalues(d2)
srand(48032) ; m1 = copy(datavalues(d1)) ; for i=1:numvars(d1) m1[:,i] += rand(n1.distributions[i],numvalues(d1)) end ; m1
srand(48032) ; m2 = applynoise!(n1,copy(datavalues(d1)))
@test_approx_eq c1 datavalues(d1)
@test_approx_eq c2 datavalues(d2)
@test_approx_eq m1 m2

###test the loglikelihood of the data relative to the noise model
@test_approx_eq loglikelihood(n1,m1,datavalues(d1)) sum(broadcast((μ,σ,x)->Distributions.logpdf(Distributions.Normal(μ,σ),x),datavalues(d1),sqrt(s1)',m1))

println()
println("====================")
println("Test show() function")
show(d1)
show(d2)
show(d3)
show(d4)
show(n1)
println("End  show() function")
println("====================")
println()

nothing

ext: jl


text:
function AbstractSelfFeedingLoader(arg0::JString, arg1::DataProvidersManager)
    return AbstractSelfFeedingLoader((JString, DataProvidersManager), arg0, arg1)
end


text:
# draw histogram of sensor data 600

module DrawHist600
    using Plots, DataFrames, CSV
    pyplot()

    function main(is_test=false)
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_600.txt")
        df_600_mm = CSV.read(data_path, DataFrame,
                             header=["date", "time", "ir", "lidar"],
                             delim=' ')

        bin_min_max = maximum(df_600_mm.lidar) - minimum(df_600_mm.lidar)
        histogram(df_600_mm.lidar, bins=bin_min_max, label="histogram")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/complex_dist/histogram_600mm/histogram_600_mm.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
@testset "interpolate.interp1d" begin
    x = collect(0:10)
    y = exp.(-x./3.0)
    f = interpolate.interp1d(x, y)
    @test f([0.5, 0.6]) == [0.8582656552868946, 0.8299187863442736]
    @test f(0.5)[1] == 0.8582656552868946
end


text:
#Lowpass and Bandpass Demonstration Script
include("dsexample1_fn.jl")
#Reminder: nlev=M+1

function wait_for_user(dsm)
	println("\nSimulation complete for modulator:\n\t$dsm")
	println("\n*** Press ENTER to continue ***")
	readline(stdin)
end

#2nd-order lowpass
dsm = RealDSM(order=2, OSR=16, M=8, Hinf=2.0, opt=0)
dsexample1(dsm, LiveDemo=true)
wait_for_user(dsm)

#5th-order lowpass
dsm = RealDSM(order=5, OSR=16, M=8, Hinf=2.0, opt=0)
dsexample1(dsm, LiveDemo=true)
wait_for_user(dsm)

#5th-order lowpass with optimized zeros
dsm = RealDSM(order=5, OSR=16, M=8, Hinf=2.0, opt=1)
dsexample1(dsm, LiveDemo=true)
wait_for_user(dsm)

#5th-order lowpass with optimized zeros and larger Hinf
dsm = RealDSM(order=5, OSR=16, M=8, Hinf=3.0, opt=1)
dsexample1(dsm, LiveDemo=true)
wait_for_user(dsm)

#7th-order lowpass; Hinf=2
dsm = RealDSM(order=7, OSR=8, M=16, Hinf=2.0, opt=1)
dsexample1(dsm, ampdB=-6.0, LiveDemo=true)
wait_for_user(dsm)

#7th-order lowpass; Hinf=8
dsm = RealDSM(order=7, OSR=8, M=16, Hinf=8.0, opt=1)
dsexample1(dsm, ampdB=-6.0, LiveDemo=true)
wait_for_user(dsm)

#6th-order bandpass
f0 = 1/6 #Normalized center frequency
dsm = RealDSM(order=6, OSR=16, M=8, f0=f0, Hinf=2.0, opt=1)
dsexample1(dsm, ampdB=-6.0, LiveDemo=true)
wait_for_user(dsm)

:END_OF_DEMO


text:
for (Func,Len) in ((:(Base.sum),:complexlength),(:linesum,:arclength))
    @eval begin
        function $Func{CC<:Chebyshev,DD}(G::Function,u::Fun{JacobiWeight{CC,DD}},z)
            d,α,β,n=domain(u),u.space.α,u.space.β,2ncoefficients(u)
            vals,t = ichebyshevtransform(pad(u.coefficients,n)),points(d,n)
            if α == β == -0.5
                return 0.5*$Len(d)*map(z->mean(G(z,t).*vals),z)*π
            else
                sp,p = space(u),plan_chebyshevtransform(complex(vals))
                map(z->$Func(Fun(chebyshevtransform(G(z,t).*vals,p),sp)),z)
            end
        end
    end
end

function logkernel{CC<:Chebyshev,DD}(G::Function,u::Fun{JacobiWeight{CC,DD}},z)
    sp,n=space(u),2ncoefficients(u)
    vals,t = ichebyshevtransform(pad(u.coefficients,n)),points(sp,n)
    p = plan_chebyshevtransform(complex(vals))
    return map(z->logkernel(Fun(chebyshevtransform(G(z,t).*vals,p),sp),z),z)
end

function cauchy{CC<:Chebyshev,DD}(G::Function,u::Fun{JacobiWeight{CC,DD}},z)
    sp,n=space(u),2ncoefficients(u)
    vals,t = ichebyshevtransform(pad(u.coefficients,n)),points(sp,n)
    p = plan_chebyshevtransform(vals)#complex(vals))
    return map(z->cauchy(Fun(chebyshevtransform(G(z,t).*vals,p),sp),z),z)
end

for TYP in (:Fourier,:Laurent)
    @eval begin
        function Base.sum{DD}(G::Function,u::Fun{$TYP{DD}},z)
            d,n=domain(u),2ncoefficients(u)
            vals,t = values(pad(u,n)),points(d,n)
            if isa(d,Circle)
              return map(z->mean(G(z,t).*vals.*t),z)*2π*im
            else
              return map(z->mean(G(z,t).*vals),z)*arclength(d)
            end
        end

        function linesum{DD}(G::Function,u::Fun{$TYP{DD}},z)
            d,n=domain(u),2ncoefficients(u)
            vals,t = values(pad(u,n)),points(d,n)
            map(z->mean(G(z,t).*vals),z)*arclength(d)
        end

#        function logkernel{DD}(G::Function,u::Fun{$TYP{DD}},z)
#            sp,n=space(u),2ncoefficients(u)
#            vals,t = values(pad(u,n)),points(sp,n)
#            p = plan_transform(sp,vals)
#            return map(z->logkernel(Fun(transform(sp,G(z,t).*vals,p),sp),z),z)
#        end
    end
end

function logkernel{DD}(G::Function,sp::Fourier{DD},u,z)
    n=2length(u)
    vals,t = values(pad(Fun(u,sp),n)),points(sp,n)
    p = plan_transform(sp,vals)
    return map(z->logkernel(sp,transform(sp,G(z,t).*vals,p),z),z)
end
logkernel{DD}(G::Function,sp::Laurent{DD},u,z)=logkernel(G,Fun(Fun(u,sp),Fourier),z)

for Func in (:(Base.sum),:linesum,:logkernel,:cauchy)
    @eval begin
        $Func{F<:Fun}(G::Function,u::Vector{F},z)=mapreduce(u->$Func(G,u,z),+,u)
        $Func(G::Function,sp::PiecewiseSpace,u,z)=$Func(G,vec(Fun(u,sp)),z)
        $Func(G::Function,u::Fun,z)=$Func(G,space(u),coefficients(u),z)
        $Func{PS<:PolynomialSpace}(G::Function,sp::JacobiWeight{PS},f,z)=$Func(G,Fun(Fun(f,sp),JacobiWeight(sp.α,sp.β,Chebyshev(domain(sp)))),z)
    end
end

function logkernel{LS,RR<:Arc}(G::Function,sp::Space{LS,RR},u,z)
    n=2length(u)
    vals,t = itransform(sp,pad(u,n)),points(sp,n)
    p = plan_transform(sp,complex(vals))
    return map(z->logkernel(sp,transform(sp,G(z,t).*vals,p),z),z)
end

function linesum{LS,RR<:Arc}(G::Function,sp::Space{LS,RR},u,z)
    n=2length(u)
    vals,t = itransform(sp,pad(u,n)),points(sp,n)
    p = plan_transform(sp,complex(vals))
    return map(z->linesum(sp,transform(sp,G(z,t).*vals,p),z),z)
end


text:
using DelayEmbeddings
using Test, StaticArrays
println("Testing custom QR-decomposition...")

@testset "QR-decomposition" begin
    tol = 1e-10
    @testset "Known matrix" begin
        A = [1. 1 0;
             1  0 1;
             0  1 1]

        QA, RA = DynamicalSystemsBase.qr_sq(A)
        QtA =  [1/√2   1/√6  -1/√3;
                1/√2  -1/√6   1/√3;
                0      2/√6   1/√3]
        RtA =  [2/√2,  3/√6,  2/√3]
        for i in length(QA)
            @test isapprox(abs(QA[i]), abs(QtA[i]), rtol = tol)
        end
        for i in 1:3
            @test abs(RA[i,i]) ≈ abs(RtA[i])
        end
    end
    @testset "Random matrix" begin
        for i in 1:10

            A = rand(5,5)
            QA, RA = DynamicalSystemsBase.qr_sq(A)
            QtA, RtA = qr(A)
            for i in length(QA)
                @test isapprox(abs(QA[i]), abs(QtA[i]), rtol = tol)
            end
            for i in 1:3
                @test isapprox(abs(RA[i,i]), abs(RtA[i,i]), rtol = tol)
            end
        end
    end
end


text:
@testset "constructor" begin
    p = Properties()
    @test typeof(p.properties) == Dict{String,String}
end

@testset "put string" begin
    p = Properties()
    setprop(p, "k1", "key one")
    @test getprop(p, "k1") == "key one"

    setprop(p, " k 2 ", raw"C:\Users\igor\elefant\.jl")
    @test getprop(p, " k 2 ") == "C:\\Users\\igor\\elefant\\.jl"

    setprop(p, "k1", "another key")
    @test getprop(p, "k1") == "another key"
end

@testset "put !string" begin
    p = Properties()
    setprop(p, "int", 1)
    setprop(p, "float64", 1.0)
    setprop(p, "float32", 1.0f0)

    @test getprop(p, "int") == "1"
    @test getprop(p, "float64") == "1.0"
    @test getprop(p, "float32") == "1.0"
end

@testset "get" begin
    p = Properties()
    setprop(p, "", 15)
    setprop(p, "key1", "")
    #TODO: docs: can't put nothing, because it can't be printed
    @test getprop(p, "") == "15"
    @test getprop(p, "", "0") == "15"
    @test getprop(p, "", nothing) == "15"
    @test getprop(p, "15", "0") == "0"

    @test getprop(p, "nokey") == nothing
    @test getprop(p, "nokey", nothing) == nothing
    @test getprop(p, "nokey", "default") == "default"
end

@testset "get!" begin
    p = Properties()
    setprop(p, "key1", "val1")

    @test getprop!(p, "key1", "def1") == "val1"
    @test getprop(p, "key1") == "val1"
    @test getprop(p, "key2") == nothing

    @test getprop!(p, "key2", "def2") == "def2"
    @test getprop(p, "key2") == "def2"
end


text:
"""
    x = spqr_rank_deflation(...)

SPQR_RANK_DEFLATION constructs pseudoinverse or basic solution using deflation.

 Called from spqr_basic and spqr_cod after these routines call spqr,
 spqr_rank_ssi and spqr_rank_form_basis.  The input parameters are as used
 in spqr_basic and spqr_cod.
 The parameter call_from indicates the type of call:
      call_from = 1 indicates a call from spqr_basic
      call_from = 2 indicates a call from spqr_cod
 Output:
   x -- for call_from = 1 x is a basic solution to
                 min || b - A * x ||.               (1)
        The basic solution has n - (rank returned by spqr) zero components.
        For call_from = 2 x is an approximate pseudoinverse solution to (1).
 Not user-callable.

 Algorithm:   R * wh = ch or R' * wh = ch is solved where ch is described
              in the code and R comes from the QR factorizations in
              spqr_basic or spqr_cod. R is triangular and potentially
              numerically singular with left and right singular vectors for
              small singular values stored in U and V.  When R is numerically
              singular deflation (see SIAM SISC, 11:519-530, 1990) to
              calculate an approximate truncated singular value solution to
              the triangular system.  Orthogonal transformations
              are applied to wh to obtain the solutions x to (1).

"""
function spqr_rank_deflation(call_from, R::AbstractMatrix{T}, U, V, C, m, n, rank_spqr,
                    numerical_rank, prow, pcol, Q::AbstractMatrix{T}; opts...) where T

# disable nearly-singular matrix warnings, and save the current state
## TODO: how? warning_state = warning ('off', 'MATLAB:nearlySingularMatrix') ;

    opts, nsvals_large, start_with_A_transpose =
        get_opts(opts, :nsvals_large, :start_with_A_transpose)

    if isempty(C)
        return zeros(T, m)
    end
    # restrict to columns of U, V belonging to small singular Values of R (< tol)
    if numerical_rank < rank_spqr
        k = size(U, 2)
        U = view(U, 1:m, nsvals_large+1:k)
        V = view(V, 1:m, nsvals_large+1:k)
    end

   if start_with_A_transpose || call_from == 1

       ch = view(C, 1:rank_spqr, 1:size(C,2))
       if numerical_rank == rank_spqr
           wh = R \ ch
       else
           # use deflation (see SIAM SISC, 11:519-530, 1990) to calculate an
           # approximate truncated singular value solution to R * wh = ch
           wh = ch - U * (U' * ch)
           wh = R \ wh
           wh = wh - V * (V' * wh)
        end
        if call_from == 2
            wh[prow,:] = wh;
        end
        wh = [wh ; zeros(T, n - rank_spqr,size(C,2)) ];
        if call_from == 1
            x = wh
            x[pcol,:] = x
        else
            x = spqr_qmult(Q, wh, 1)
        end

    else

        ch = view(C, prow, 1:size(C,2))
        if numerical_rank == rank_spqr
            wh = R' \ ch
        else
            # use deflation (see SIAM SISC, 11:519-530, 1990) to calculate an
            # approximate truncated singular value solution to R' * wh = ch
            wh = ch - V * (V' * ch)
            wh = R' \ wh
            wh = wh - U * (U' * wh)
        end
        wh = [wh ; zeros(T, n - rank_spqr, size(C,2)) ];
        wh = spqr_qmult(Q, wh, 1)
        x[pcol,:] = wh

    end

    # restore the warning back to what it was
    ## TODO: warning (warning_state) ;

    x
end


text:
module DataFile

using Avro.Schemas

export VALID_CODECS,
       OBJECT_CONTAINER_FILE_MAGIC,
       META_CODEC_KEY,
       META_SCHEMA_KEY,
       MAGIC_SCHEMA,
       SYNC_SCHEMA,
       METADATA_SCHEMA

const VALID_CODECS = Set(["null", "deflate"])

"""
Magic constant that is written to the top of every object container file. This
is the hexadecimal form of 'Obj1'.
"""
const OBJECT_CONTAINER_FILE_MAGIC = [0x4f, 0x62, 0x6a, 0x01]

const META_CODEC_KEY = "avro.codec"
const META_SCHEMA_KEY = "avro.schema"

const MAGIC_SCHEMA = Schemas.FixedSchema(Schemas.FullName("Magic"), 4)
const META_SCHEMA = Schemas.MapSchema(Schemas.STRING)
const SYNC_SCHEMA = Schemas.FixedSchema(Schemas.FullName("Sync"), 16)
const METADATA_SCHEMA =
    Schemas.RecordSchema(
        Schemas.FullName("Header", "org.apache.avro.file"),
        [
            Schemas.Field("magic", 0, MAGIC_SCHEMA),
            Schemas.Field("meta", 1, META_SCHEMA),
            Schemas.Field("sync", 2, SYNC_SCHEMA)
        ]
    )

include("codec.jl")
include("data_file_writer.jl")
include("data_file_reader.jl")

using Avro.DataFile.Writers
using Avro.DataFile.Readers

end


text:
# A comparison between ℓ_0 and ℓ_1 trend filtering

using Plots
using Convex
using Mosek
using EllZeroTrendFiltering

g = snp500_data()[1:300]

N = length(g)

## ℓ_1 trend filtering, see "ℓ_1 Trend Filtering", Kim et al. (2009)
H =  spdiagm(0=>ones(N-2), 1=>-2*ones(N-2), 2=>ones(N-2))

x = Variable(N)

λ = 2.0
problem = minimize(sumsquares(x - g) + λ * norm(H*x, 1))

solve!(problem, Mosek.Optimizer)
#solve!(problem, SCS.Optimizer)
println(problem.status)

x_l1 = evaluate(x)[:]
f_l1 = sum((x_l1 - g).^2)

## ℓ_0 trend filtering
ζ = 0.02
I_l0, Y_l0, f_l0 = fit_pwl_regularized(g, ζ)

plot(g, label="SP500", size=(1000,550));
plot!(x_l1, lw=2.5, label="l1 reg. lambda = $λ");
plot!(I_l0, Y_l0, lw=2.5, label="l0 reg. zeta = $ζ")

# plot(H*x.value);
# plot!(I_l0[2:end], 0.1*diff(Y_l0) ./ diff(I_l0), m=:o, lw=0)
#


text:
using SHA, JSON
using Restful, Logging
using UUIDs
using HTTP, URIParser

import Restful: json

difficulty = Int8(4)

struct Transaction
    sender::String
    recipient::String
    amount::Float64
end

mutable struct Block
    index::Int64
    timestamp::Int64
    transaction_list::Array{Transaction}
    proof::Int64
    previous_hash::String
end

mutable struct Blockchain
    chain::Array{Block}
    current_transaction::Array{Transaction}
    nodes::Array{String}
end

function init(blockchain::Blockchain)
    blockchain.chain = []
    blockchain.current_transaction = []
    blockchain.nodes = []
    new_block(blockchain=blockchain, proof=100, previous_hash="genesis_block_hash")
end

function new_block(;blockchain::Blockchain, proof::Int64, previous_hash = nothing)
    if previous_hash == nothing
        previous_hash = blockhash(blockchain.chain[end])
    end
    block = Block(length(blockchain.chain)+1, round(Int64, time()), blockchain.current_transaction, proof, previous_hash)

    push!(blockchain.chain, block)
    blockchain.current_transaction = []
    return block
end

function new_transaction(;blockchain::Blockchain, sender::String, recipient::String, amount::Float64)
    new_tx = Transaction(sender, recipient, amount)
    push!(blockchain.current_transaction, new_tx)
    return length(blockchain.chain) + 1
end

function blockhash(block::Block)
    return bytes2hex(sha256(JSON.json(block)))
end

function proof_of_work(last_proof)
    proof = 0
    while valid_proof(last_proof, proof) == false
        proof += 1
    end
    return proof
end

function valid_proof(last_proof, proof)
    header = "0"^difficulty
    guess = string(last_proof) * string(proof)
    guess_hash = bytes2hex(sha256(guess))
    return guess_hash[1:difficulty] == header
end

function register_node(blockchain::Blockchain, address::String)
    push!(blockchain.nodes, address)
end

function valid_chain(chain::Array{Block})
    last_block = chain[end]
    current_index = 1
    while current_index < length(chain)
        block = chain[current_index]
        if block.previous_hash != blockhash(last_block)
            return false
        end

        if valid_proof(last_block.proof, block.proof) == false
            return false
        end

        last_block = block
        current_index += 1
    end

    return true
end

function resolve_conflict(blockchain::Blockchain)
    neighbours = blockchain.nodes
    new_chain = nothing

    my_height = length(blockchain.chain)

    for node in neighbours
        response = HTTP.request("GET", "http://" * node * "/chain")

        if response.status == 200
            response_body = JSON.Parser.parse(String(response.body))
            height = response_body[1]["length"]
            chain = response_body[2]["chain"]

            if height > my_height && valid_chain(chain)
                my_height = height
                new_chain = chain
            end
        end
    end

    if new_chain != nothing
        blockchain.chain = new_chain
        return true
    end

    return false
end

const app = Restful.app()

blockchain = Blockchain([], [], [])
node_identifier = replace(string(uuid4()), "-" => "")
init(blockchain)

app.get("/mine", json) do req, res, route
    last_block = blockchain.chain[end]
    last_proof = last_block.proof
    proof = proof_of_work(last_proof)

    new_transaction(blockchain=blockchain, sender="0", recipient=node_identifier, amount=1.0)

    previous_hash = blockhash(last_block)
    block = new_block(blockchain=blockchain, proof=proof, previous_hash=previous_hash)
    res.json(Dict("message" => "New Block Forged",
    "index" => block.index,
    "transaction" => block.transaction_list,
    "proof" => block.proof,
    "previous_hash" => block.previous_hash) |> collect)
end

app.post("/transactions/new", json) do req, res, route
    required = ["sender", "recipient", "amount"]
    data = JSON.parse(req.body)

    if all(i->(i in keys(data)), required) == false
        res.code(400)
    else
        block_id = new_transaction(blockchain=blockchain, sender=data["sender"], recipient=data["recipient"], amount=data["amount"])
        res.json(Dict("message" => "Transaction will be added to Block " * string(block_id)) |> collect)
    end
end

app.get("/chain", json) do req, res, route
    res.json(Dict("chain"=>blockchain.chain, "length"=>length(blockchain.chain)) |> collect)
end

app.post("/nodes/register", json) do req, res, route
    data = JSON.parse(req.body)
    nodes = data["nodes"]

    if nodes == nothing
        res.code(400)
    end

    for node in nodes
        register_node(blockchain, node)
    end

    res.json(Dict("message" => "Nodes will be added to Blockchain") |> collect)
end

app.get("/nodes/resolve", json) do req, res, route
    replaced = resolve_conflict(blockchain)
    if replaced == true
        res.json(Dict("message" => "Our chain was replaced") |> collect)
    else
        res.json(Dict("message" => "Our chain is main chain") |> collect)
    end
end

@async with_logger(SimpleLogger(stderr, Logging.Warn)) do
    app.listen("127.0.0.1", 3001)
end

max_stars_repo_name: HuLaTown/blockchain

max_issues_repo_name: HuLaTown/blockchain

max_forks_repo_name: HuLaTown/blockchain


text:
using Test
using Mods

p = 23
@test Mod(2,p) == Mod(25,p)
@test Mod(2,p)+Mod(21,p) == Mod(0,p)
@test Mod(10,p)*Mod(5,p) == Mod(4,p)
@test Mod(10,p)/Mod(5,p) == Mod(2,p)
@test Mod(17,p) * inv(Mod(17,p)) == Mod(1,p)
@test Mod(17^6,p)==Mod(17,p)^6
@test Mod(17,p)^(-2) == inv(Mod(17,p))^2
@test Mod(17,p) == 17

q = 91
a = Mod(17,p)
b = Mod(32,q)

x = CRT(a,b)
@test Mod(x.val,p)==a
@test Mod(x.val,q)==b
@test inv(a)*a == 1

p = 9223372036854775783   # This is a big prime
x = Mod(-2,p)
@test x*x == 4
@test x+x == -4
@test x/x == 1
@test x^4 == 16
@test x^(p-1) == 1   # Fermat Little Theorem test
@test 2x == x+x
@test x-x == 0
y = inv(x)
@test x*y == 1
@test x+p == x
@test x*p == 0
@test p-x == x-2x


text:
using LinearAlgebra
println("testing transforms")

# function test_affine
moving_pts = [1.0 1.0;
			2.0 3.0;
			3.0 3.0]
fixed_pts = [1.0 1.0;
			2.0 3.0;
			3.0 3.0]
tform = calculate_affine(moving_pts, fixed_pts)
@test tform ≈ Matrix(1.0I,2,2)

tform = calculate_rigid(moving_pts, fixed_pts)
@test tform ≈ Matrix(1.0I,2,2)

tform = calculate_translation(moving_pts, fixed_pts)
@test tform ≈ Matrix(1.0I,3,3)

moving_pts = [1.0 1.0;
			2.0 3.0;
			3.0 3.0]
fixed_pts = [5.0 1.0;
			6.0 3.0;
			7.0 3.0]
shift = [1 0 0; 0 1 0; 4 0 1]

tform = calculate_translation(moving_pts, fixed_pts)
@test tform ≈ shift


text:
using Polynomials, Roots
export get_polynomial
export get_polynomial_coefficients
export compute_fwd_theta
export compute_bnd_rel_bwd_err
export compute_bwd_theta

"""
    p = get_polynomial(graph::Compgraph)

Return the polynomial underlying the computational `graph`. The `graph` is
assumed to involve only linear combinations and multiplications.

`p` is polynomial of type `Polynomial` from the package
[`Polynomials.jl`](https://juliamath.github.io/Polynomials.jl/stable/).

See also [`get_polynomial_coefficients`](@ref).
"""
function get_polynomial(graph::Compgraph)
    graph = Compgraph(Any, graph) # Allow for symbolic values
    p = eval_graph(graph, Polynomial("x"))
    return p
end

"""
    c = get_polynomial_coefficients(graph::Compgraph)

Return the coefficients of the polynomial underlying the computational `graph`.
The `graph` is assumed to involve only linear combinations and multiplications.

The coefficients `c` is a vector containing the coefficients as expressed in
the monomial basis and sorted from the leading coefficient to the constant term.

See also [`get_polynomial`](@ref).
"""
function get_polynomial_coefficients(graph::Compgraph)
    return get_polynomial(graph).coeffs
end

"""
    (e_fwd,theta)=compute_fwd_theta(
        graph::Compgraph{T},
        f;
        coefftype = T,
        tolerance = eps(coefftype) / 2,
        theta_init = big"0.1",
    )

Return relative forward error and corresponding theta for approximation to `f`.

The first output ``e_{fwd}`` is the function that computes the relative forward
error``e_{fwd}(z)=|f(z) - p(z)| / |z|``, where ``p`` is the polynomial
underlying the computational graph `graph`. This is meaningful only if ``p``
approximates ``f`` in some sense.

The second output ``θ`` is the largest positive real number such that
``e_{fwd}(θ) ≤ tolerance``, where ``tolerance`` is by default the unit roundoff
of the data type `coefftype`, which in turn defaults to the type of the
coefficients of `graph`.

The value of ``θ`` is approximated by using the built-in function `fzero` with
starting value set to `theta_init`. By default, this root-finding procedure
uses high-precision arithmetic.
"""
function compute_fwd_theta(
    graph::Compgraph{T},
    f;
    coefftype = T,
    tolerance = eps(coefftype) / 2,
    theta_init = big"0.1",
) where {T}
    # Obtain coefficients of polynomial.
    coeff = get_polynomial_coefficients(graph)

    # Convert coefficients to required type.
    coeff = convert.(coefftype, coeff)
    if isreal(coeff)
        coeff = real(coeff)
    end

    # Find point where the relative forward error equals tolerance.
    p = Polynomial(abs.(coeff), :x)
    e_fwd(z) = abs.(f.(z) - p.(z)) ./ abs.(z)
    g(z) = e_fwd(z) - tolerance
    theta_fwd = fzero(g, theta_init, verbose = true)
    return e_fwd, convert(coefftype, theta_fwd)
end

"""
    e_bwd=compute_bnd_rel_bwd_err(
        f,
        graph::Compgraph{T};
        coefftype = T,
        numterms = 100,
        numdigits = 100,
    )

Compute a bound on the relative backward error of the function `f`.

The returned function ``e_{bwd}`` is a bound on the relative backward error of
the polynomial underlying the computational graph `graph`, seen as an
approximant to `f`. For the underlying polynomial ``p``, the
function returns a bound on `|δ|` where `δ` is such that ``f(z+δ) = p(z)``.

The first argument `f` is a symbol. Currently supported values include:

  - `:exp` The bound is computed by means of the identity δ = log((exp(-z)
    p(z)-1)+1). The code computes the series expansion of the right-hand side of the
    equation, truncates it to the first `numterms` coefficients, bounds each
    coefficient of the ensuing polynomial with its absolute value. The computation
    uses `numdigits` digits of precision.
"""
function compute_bnd_rel_bwd_err(
    f,
    graph::Compgraph{T};
    coefftype = T,
    numterms = 100,
    numdigits = 100,
) where {T}
    if f != :exp
        error("Currently :exp is the only supported function.")
    end
    # Set precision to numdigits decimal digits.
    setprecision(Integer(ceil(log2(10.0^numdigits))))

    # Obtain coefficients of polynomial p.
    coeff = get_polynomial_coefficients(graph)

    # Convert coefficients of p to required type.
    coeff = convert.(coefftype, coeff)
    if isreal(coeff)
        coeff = real(coeff)
    end

    # Series expansion of polynomial approximant p.
    papproximant = Polynomial(coeff)

    # (Truncated) series expansions of exp(-z) and log(1+z).
    expminusz =
        (-1) .^ (0:1:numterms) ./
        factorial.(collect(big(0.0):big(1.0):numterms))
    pexpminusz = Polynomial(expminusz)
    logzplusone = [0; (-1) .^ (0:1:numterms-1) ./ (1:1:numterms)]
    plogzplusone = Polynomial(logzplusone)

    # From exp(z+δ) ≈ p(z), approximate δ ≈ log((exp(-z) p(z)-1)+1).
    # The variables are as follows:
    #    * pexpzpz: the coefficients of exp(-z) p(z)
    #    * presult: the coefficients of log((exp(-z) p(z)-1)+1)
    pexpzpz = Polynomial((pexpminusz*papproximant).coeffs[1:numterms])
    presult = Polynomial((plogzplusone(pexpzpz - 1)).coeffs[1:numterms])

    # Compute bound on |δ|.
    bnd_bwd_err = Polynomial(abs.(presult.coeffs))
    e_bwd(z) = abs.(bnd_bwd_err(z)) ./ abs.(z)

    return e_bwd
end

"""
    (e_bwd,theta)=compute_bwd_theta(;
        bnd_rel_err = compute_bnd_rel_bwd_err(:exp, graph),
        tolerance = eps(coefftype) / 2,
        theta_init = big"0.2",
        use_log = false,
    )

    (e_bwd,theta)=compute_bwd_theta(
        graph::Compgraph{T};
        coefftype = T,
        numterms = 100,
        numdigits = 100,
        tolerance = eps(coefftype) / 2,
        theta_init = big"0.2",
        use_log = false,
    )

Compute bound on relative backward error with corresponding theta.

The first output ``e_{bwd}`` is a function that returns a bound on the relative
backward error of the polynomial underlying the computational graph `graph`,
seen as an approximant to the exponential. For the underlying polynomial ``p``,
the function returns a bound on `|δ|` where `δ` is such that exp(z+δ) = p(z).
This function is computed using `compute_bnd_rel_bwd_err(:graph)`.

The second output ``θ`` is the largest positive real number such that
``e_{bwd}(θ) ≤ tolerance``, where ``tolerance`` is by default the unit roundoff
of the data type `coefftype`, which in turn defaults to the type of the
coefficients of `graph`.

The value of ``theta`` is estimated by approximately solving the equation
``e_{bwd}(z) = tolerance``, using the built-in function `fzero` with starting
value set to `theta_init`. By default this root-finding procedure uses
high-precision arithmetic.

If the kwarg `use_log` is set to `true`, then the value of ``theta`` is computed
by approximating a solution to the equation ``log(e_{bwd}(z)) = log(tolerance)``
instead.

The alternative form accepts a function that returns a bound on the relative
backward error. By default, the function is constructed with
`compute_bnd_rel_bwd_err(:exp, ...)` and the default values for the kwargs in
the first form.
"""
function compute_bwd_theta(;
    bnd_rel_err = compute_bnd_rel_bwd_err(:exp, graph),
    tolerance = eps(coefftype) / 2,
    theta_init = big"0.2",
    use_log = false,
) where {T}
    # Find point where bound on relative backward error equals tolerance.
    e_bwd = bnd_rel_err
    h(z) = use_log ? log(e_bwd(z) / tolerance) : e_bwd(z) - tolerance
    theta_bwd = fzero(h, theta_init)

    return e_bwd, theta_bwd
end

function compute_bwd_theta(
    graph::Compgraph{T};
    coefftype = T,
    numterms = 100,
    numdigits = 100,
    tolerance = eps(coefftype) / 2,
    theta_init = big"0.2",
    use_log = false,
) where {T}
    return compute_bwd_theta(
        bnd_rel_err = compute_bnd_rel_bwd_err(
            :exp,
            graph;
            coefftype = coefftype,
            numterms = numterms,
            numdigits = numdigits,
        ),
        tolerance = tolerance,
        theta_init = theta_init,
        use_log = use_log,
    )
end


text:
using Test
using Dates: now
using ExtensibleScheduler
using ExtensibleScheduler: Job, get_job_id, JobConfig

@testset "Job" begin

    function return_args(id)
        "Returned from return_args $id"
    end

    id1 = get_job_id()
    action1 = Action(return_args, 1)
    @test run(action1) == "Returned from return_args 1"
    trigger1 = Trigger(DateTime(2010, 1, 1))
    name1 = "Textual description of job1"
    priority1 = 0
    dt_created1 = now()
    dt_updated1 = dt_created1
    dt_next_fire1 = DateTime(0)
    n_triggered1 = 0
    config1 = JobConfig()
    job1 = Job(id1, action1, trigger1, name1, priority1, dt_created1, dt_updated1, dt_next_fire1, n_triggered1, config1)

end


text:
# general interface functions
"""
    grid(D)

Return the grid associated to a derivative operator `D`.
"""
function grid end

"""
    accuracy_order(D)

Return the order of accuracy of a derivative operator `D`. For SBP finite difference
operators, this refers to the interior order of accuracy.
"""
function accuracy_order end

"""
    derivative_order(D)

Return the order of the derivative associated to the derivative operator `D`.
For example, it will return `1` for a first-derivative SBP operator.
"""
function derivative_order end

"""
    source_of_coefficients(D)

Return the source of coefficients of the derivative operator `D`. If you use the
operator `D` for your research, please cite this source in addition to
[`SummationByPartsOperators`](@ref).
"""
source_of_coefficients(D) = SummationByPartsOperators

"""
    left_boundary_weight(D)

Return the left-boundary weight of the (diagonal) mass matrix `M` associated to
the derivative operator `D`.
"""
function left_boundary_weight end

"""
    right_boundary_weight(D)

Return the left-boundary weight of the (diagonal) mass matrix `M` associated to
the derivative operator `D`.
"""
function right_boundary_weight end

function Base.summary(io::IO, D::AbstractDerivativeOperator)
    print(io, nameof(typeof(D)), "(derivative_order=", derivative_order(D),
              ", accuracy_order=", accuracy_order(D), ")")
end

abstract type AbstractExecutionMode end
"""
    SafeMode()

A safe execution mode relying only on basic functionality of Julia.
"""
struct SafeMode <: AbstractExecutionMode end
"""
    FastMode()

A (probably) faster execution mode that might depend on packages such as
[LoopVectorization.jl](https://github.com/JuliaSIMD/LoopVectorization.jl).
"""
struct FastMode <: AbstractExecutionMode end
"""
    ThreadedMode()

An execution mode using multiple threads and possibly further optimizations, cf.
[`FastMode`](@ref).
"""
struct ThreadedMode <: AbstractExecutionMode end

# TODO: deprecated in v0.5
_parallel_to_mode(::Val{:threads}) = ThreadedMode()
_parallel_to_mode(::Val{:serial}) = FastMode()

derivative_order(coefficients::AbstractDerivativeCoefficients) = coefficients.derivative_order
accuracy_order(coefficients::AbstractDerivativeCoefficients) = coefficients.accuracy_order
Base.eltype(coefficients::AbstractDerivativeCoefficients{T}) where {T} = T
LinearAlgebra.issymmetric(coefficients::AbstractDerivativeCoefficients) = coefficients.symmetric

derivative_order(D::AbstractDerivativeOperator) = derivative_order(D.coefficients)
accuracy_order(D::AbstractDerivativeOperator) = accuracy_order(D.coefficients)
Base.eltype(D::AbstractDerivativeOperator{T}) where {T} = T
Base.real(D::AbstractDerivativeOperator) = real(eltype(D))
LinearAlgebra.issymmetric(D::AbstractDerivativeOperator) = issymmetric(D.coefficients)
function Base.size(D::AbstractDerivativeOperator)
    N = length(grid(D))
    (N, N)
end
function Base.size(D::AbstractDerivativeOperator, i::Int)
    if i <= 0
        error("arraysize: dimension out of range")
    elseif i <= 2
        length(grid(D))
    else
        1
    end
end

@inline grid(D::AbstractDerivativeOperator) = D.grid
@inline grid(D::AbstractPeriodicDerivativeOperator) = D.grid_compute

Base.@propagate_inbounds function mul!(dest, D::AbstractDerivativeOperator, u)
    mul!(dest, D, u, one(recursive_bottom_eltype(dest)))
end

function Base.:*(D::AbstractDerivativeOperator, u)
    @boundscheck begin
        @argcheck size(D,1) == size(D,2) == length(u) DimensionMismatch
    end
    T = typeof(one(eltype(D)) * first(u))
    dest = similar(u, T); fill!(dest, zero(eltype(dest)))
    @inbounds mul!(dest, D, u)
    dest
end

function Base.Matrix(D::AbstractDerivativeOperator{T}) where {T}
    v = Array{T}(undef, size(D, 2)...)
    fill!(v, zero(eltype(v)))
    A = Array{T}(undef, size(D)...)
    for i in 1:size(D,2)
        v[i] = one(T)
        # Using a view here can cause problems with FFT based operators.
        # Since this part is not performance critical, we can also just use copies.
        # mul!(view(A,:,i), D, v)
        A[:,i] .= D * v
        v[i] = zero(T)
    end
    A
end

function SparseArrays.sparse(D::AbstractDerivativeOperator{T}) where {T}
    M, N = size(D)
    rowind = Vector{Int}()
    nzval = Vector{T}()
    colptr = Vector{Int}(undef, N+1)
    v = fill(zero(T), N)
    dest = Array{T}(undef, M)

    for i = 1:N
        v[i] = one(T)
        mul!(dest, D, v)
        js = findall(!iszero, dest)
        colptr[i] = length(nzval)+1
        if length(js) > 0
            append!(rowind, js)
            append!(nzval, dest[js])
        end
        v[i] = zero(T)
    end
    colptr[N+1] = length(nzval)+1

    return SparseMatrixCSC(M, N, colptr, rowind, nzval)
end

"""
    compute_coefficients(u, D::AbstractDerivativeOperator)

Compute the nodal values of the function `u` at the grid associated to the
derivative operator `D`.
"""
function compute_coefficients(u, D::AbstractDerivativeOperator)
    x = grid(D)
    xmin = first(x)
    xmax = last(x)
    uval = Array{typeof(u((xmin+xmax)/2))}(undef, size(x)...)
    compute_coefficients!(uval, u, D)
end

"""
    compute_coefficients!(uval, u, D::AbstractDerivativeOperator)

Compute the nodal values of the function `u` at the grid associated to the
derivative operator `D` and stores the result in `uval`.
"""
function compute_coefficients!(uval, u, D::AbstractDerivativeOperator)
    uval .= u.(grid(D))
end

"""
    evaluate_coefficients(u, D::AbstractDerivativeOperator)

Evaluates the nodal coefficients `u` at a grid associated to the derivative
operator `D`.
Returns `xplot, uplot`, where `xplot` contains the nodes and `uplot` the
corresponding values of `u`.
"""
function evaluate_coefficients(u, D::AbstractDerivativeOperator)
    x = grid(D)
    xplot = Array{eltype(x)}(undef, size(x)...)
    uplot = Array{eltype(u)}(undef, size(x)...)

    evaluate_coefficients!(xplot, uplot, u, D)
end

"""
    evaluate_coefficients!(xplot, uplot, u, D::AbstractDerivativeOperator)

Evaluates the nodal coefficients `u` at a grid associated to the derivative
operator `D` and stores the result in `xplot, uplot`.
Returns `xplot, uplot`, where `xplot` contains the nodes and `uplot` the
corresponding values of `u`.
"""
function evaluate_coefficients!(xplot, uplot, u, D::AbstractDerivativeOperator)
    @argcheck size(uplot) == size(xplot)
    @argcheck size(uplot) == size(grid(D))

    xplot .= grid(D)
    uplot .= u

    xplot, uplot
end

"""
    integrate(func, u, D::AbstractPeriodicDerivativeOperator)

Map the function `func` to the coefficients `u` and integrate with respect to
the quadrature rule associated with the derivative operator `D`.
"""
function integrate(func, u::AbstractVector, D::AbstractPeriodicDerivativeOperator)
    @boundscheck begin
        length(u) == length(grid(D))
    end
    @unpack Δx = D

    res = sum(func, u)

    Δx * res
end

"""
    SumOfDerivativeOperators

Sum several derivative operators lazily.
"""
struct SumOfDerivativeOperators{T,N,Operators<:Tuple{Vararg{AbstractDerivativeOperator{T},N}}} <: AbstractDerivativeOperator{T}
    operators::Operators

    function SumOfDerivativeOperators(operators::Operators) where {T,N,Operators<:Tuple{Vararg{AbstractDerivativeOperator{T},N}}}
        @argcheck all(i->size(operators[i]) == size(first(operators)), eachindex(operators)) DimensionMismatch
        @argcheck all(i->grid(operators[i]) ≈ grid(first(operators)), eachindex(operators)) ArgumentError
        new{T,N,Operators}(operators)
    end
end

SumOfDerivativeOperators(ops...) = SumOfDerivativeOperators(ops)

Base.size(sum::SumOfDerivativeOperators) = size(first(sum.operators))
Base.size(sum::SumOfDerivativeOperators, i::Int) = size(first(sum.operators), i)
function Base.length(::Type{SumOfDerivativeOperators{T,N,Operators}}) where {T,N,Operators}
    N
end
grid(sum::SumOfDerivativeOperators) = grid(first(sum.operators))

function Base.:+(D1::AbstractDerivativeOperator, D2::AbstractDerivativeOperator)
    SumOfDerivativeOperators(D1, D2)
end

function Base.:+(sum::SumOfDerivativeOperators, D::AbstractDerivativeOperator)
    SumOfDerivativeOperators(sum.operators..., D)
end

function Base.:+(D::AbstractDerivativeOperator, sum::SumOfDerivativeOperators)
    SumOfDerivativeOperators(D, sum.operators...)
end

function Base.:+(sum1::SumOfDerivativeOperators, sum2::SumOfDerivativeOperators)
    SumOfDerivativeOperators(sum1.operators..., sum2.operators...)
end

function Base.show(io::IO, sum::SumOfDerivativeOperators)
    print(io, "Sum of ", length(sum.operators), " operators")
    if get(io, :compact, false) == false
        print(io, ":")
        for D in sum.operators
            print(io, "\n", D)
        end
    end
end

@unroll function mul!(dest::AbstractVector, sum::SumOfDerivativeOperators, u::AbstractVector,
                      α, β)
    @unpack operators = sum
    @boundscheck begin
        @argcheck size(first(operators), 2) == length(u) DimensionMismatch
        @argcheck size(first(operators), 1) == length(dest) DimensionMismatch
    end

    @inbounds mul!(dest, operators[1], u, α, β)
    @unroll for i in 1:length(sum)
        if i != 1
            @inbounds mul!(dest, operators[i], u, α, one(β))
        end
    end

    nothing
end

@unroll function mul!(dest::AbstractVector, sum::SumOfDerivativeOperators, u::AbstractVector,
                      α)
    @unpack operators = sum
    @boundscheck begin
        @argcheck size(first(operators), 2) == length(u) DimensionMismatch
        @argcheck size(first(operators), 1) == length(dest) DimensionMismatch
    end

    @inbounds mul!(dest, operators[1], u, α)
    @unroll for i in 1:length(sum)
        if i != 1
            @inbounds mul!(dest, operators[i], u, α, one(α))
        end
    end

    nothing
end


text:
abstract type AbstractLSODAHandle end

mutable struct ContextHandle <: AbstractLSODAHandle
    ctx::lsoda_context_t
    freed::Bool
    function (::Type{ContextHandle})(ctx::lsoda_context_t)
        h = new(ctx,false)
        finalizer(release_handle,h)
        return h
    end
end
release_handle(ch::ContextHandle) = lsoda_free(ch)

function lsoda_free(ch::ContextHandle)
    if !ch.freed
        lsoda_free(ch.ctx)
        ch.freed = true
    end
    nothing
end

# Now wrap the rest of the APIs for convenience
lsoda_reset(ch::ContextHandle) = lsoda_reset(ch.ctx)
lsoda_prepare(ch::ContextHandle,opt::lsoda_opt_t) = lsoda_prepare(ch.ctx,opt)
lsoda(ch::ContextHandle,y::Vector,t::Vector{Float64},tout) = lsoda(ch.ctx,y,t,tout)


text:
using PMI
using Test

PMI.init()

nranks = PMI.get_size()
rank = PMI.get_rank()

kvs = PMI.KVS()
if rank == 0
    PMI.put!(kvs, "primary", string(rank))
    PMI.commit!(kvs)
end
PMI.barrier()

@test PMI.get(kvs, "primary") == "0"

PMI.finalize()


text:
@testset "Message dispatcher" begin

    if Sys.iswindows()
        global_socket_name1 = "\\\\.\\pipe\\jsonrpc-testrun1"
    elseif Sys.isunix()
        global_socket_name1 = joinpath(tempdir(), "jsonrpc-testrun1")
    else
        error("Unknown operating system.")
    end

    request1_type = JSONRPC.RequestType("request1", Foo, String)
    request2_type = JSONRPC.RequestType("request2", Nothing, String)
    notify1_type = JSONRPC.NotificationType("notify1", String)

    global g_var = ""

    server_is_up = Base.Condition()

    server_task = @async try
        server = listen(global_socket_name1)
        notify(server_is_up)
        sock = accept(server)
        global conn = JSONRPC.JSONRPCEndpoint(sock, sock)
        global msg_dispatcher = JSONRPC.MsgDispatcher()

        msg_dispatcher[request1_type] = (conn, params) -> begin
            @test JSONRPC.is_currently_handling_msg(msg_dispatcher)
            params.fieldA == 1 ? "YES" : "NO"
        end
        msg_dispatcher[request2_type] = (conn, params) -> JSONRPC.JSONRPCError(-32600, "Our message", nothing)
        msg_dispatcher[notify1_type] = (conn, params) -> g_var = params

        run(conn)

        for msg in conn
            JSONRPC.dispatch_msg(conn, msg_dispatcher, msg)
        end
    catch err
        Base.display_error(stderr, err, catch_backtrace())
    end

    wait(server_is_up)

    sock2 = connect(global_socket_name1)
    conn2 = JSONRPCEndpoint(sock2, sock2)

    run(conn2)

    JSONRPC.send(conn2, notify1_type, "TEST")

    res = JSONRPC.send(conn2, request1_type, Foo(fieldA=1, fieldB="FOO"))

    @test res == "YES"
    @test g_var == "TEST"

    @test_throws JSONRPC.JSONRPCError(-32600, "Our message", nothing) JSONRPC.send(conn2, request2_type, nothing)

    close(conn2)
    close(sock2)
    close(conn)

    fetch(server_task)

    # Now we test a faulty server

    if Sys.iswindows()
        global_socket_name2 = "\\\\.\\pipe\\jsonrpc-testrun2"
    elseif Sys.isunix()
        global_socket_name2 = joinpath(tempdir(), "jsonrpc-testrun2")
    else
        error("Unknown operating system.")
    end

    server_is_up = Base.Condition()

    server_task2 = @async try
        server = listen(global_socket_name2)
        notify(server_is_up)
        sock = accept(server)
        global conn = JSONRPC.JSONRPCEndpoint(sock, sock)
        global msg_dispatcher = JSONRPC.MsgDispatcher()

        msg_dispatcher[request2_type] = (conn, params)->34 # The request type requires a `String` return, so this tests whether we get an error.

        run(conn)

        for msg in conn
            @test_throws ErrorException("The handler for the 'request2' request returned a value of type $Int, which is not a valid return type according to the request definition.") JSONRPC.dispatch_msg(conn, msg_dispatcher, msg)
        end
    catch err
        Base.display_error(stderr, err, catch_backtrace())
    end

    wait(server_is_up)

    sock2 = connect(global_socket_name2)
    conn2 = JSONRPCEndpoint(sock2, sock2)

    run(conn2)

    @test_throws JSONRPC.JSONRPCError(-32603, "The handler for the 'request2' request returned a value of type $Int, which is not a valid return type according to the request definition.", nothing) JSONRPC.send(conn2, request2_type, nothing)

    close(conn2)
    close(sock2)
    close(conn)

    fetch(server_task)

end


text:
using Test, Random, LinearAlgebra, SparseArrays

using Yao
using QuAlgorithmZoo

@testset "QFT" begin
    include("QFT.jl")
end

@testset "CircuitBuild" begin
    include("CircuitBuild.jl")
end

@testset "RotBasis" begin
    include("RotBasis.jl")
end

@testset "Grover" begin
    include("Grover.jl")
end

@testset "PhaseEstimation" begin
    include("PhaseEstimation.jl")
end

@testset "HHL" begin
    include("HHL.jl")
end
@testset "diff Eq" begin
    include("lin_diffEq_test.jl")
end

@testset "QCOptProblem" begin
    include("QCOptProblem.jl")
end

@testset "hamiltonian solvers" begin
    include("hamiltonian_solvers.jl")
end

@testset "hadamard test" begin
    include("HadamardTest.jl")
end

@testset "Sequence" begin
    include("Sequence.jl")
end

@testset "Diff" begin
    include("Diff.jl")
end


text:
using PyPlot

nt = 2048
t = collect(linspace(0,10,nt))
s0 = 1.0
period = 0.3
om0 = 2pi/period
q = 1000.0
alpha = s0*om0*(1 + im/sqrt(4.*q^2-1.))
beta = om0/2/q*(1+im*sqrt(4.*q^2-1.))

# Now build matrix & Cholesky:
mat = zeros(nt,nt)
for i=1:nt
  for j=1:nt
    dt = abs.(t[i]-t[j])
    mat[i,j]=exp(-real(beta)*dt)*(real(alpha)*cos(imag(beta)*dt)+imag(alpha)*sin(imag(beta)*dt))
  end
end
chmat=transpose(chol(mat))

nfft = nt
ndev = randn(nt)
sim=*(chmat,ndev)
#plot(t,sim,color="blue")

numuniquepts = ceil(Int64,(nfft+1)/2)

spec = fft(sim)
spec = spec[1:numuniquepts]
spec = abs.(spec).^2/nfft
if rem(nfft,2) == 1
  spec[2:end]=spec[2:end]*2
else
  spec[2:end-1]=spec[2:end-1]*2
end

f = (0:numuniquepts-1)/(t[2]-t[1])/nfft

loglog(f,spec)

# Now compute noisy spectrum:
noise = 50.0
simn = sim + noise*randn(nt)

specn = fft(simn)
specn = specn[1:numuniquepts]
specn = abs.(specn).^2/nfft
if rem(nfft,2) == 1
  specn[2:end]=specn[2:end]*2
else
  specn[2:end-1]=specn[2:end-1]*2
end

loglog(f,specn)
om = 2pi*f
spect = s0*om0^4./((om.^2-om0^2).^2+om0^2.*om.^2./q^2)
loglog(f,spect)


text:
#===============================================================================

    AmpInput.jl

    defines struct of data and all analysis parameters
    to be passed to amp_analysis() in amplification.jl

    the constructor is intended as the only interface
    to the amplification analysis and the only place
    where argument defaults are applied

    Author: Tom Price
    Date:   July 2019

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

import DataStructures.OrderedDict
import StaticArrays: SMatrix, SVector
import Ipopt: IpoptSolver #, NLoptSolver

#===============================================================================
    defaults >>
===============================================================================#

## defaults for solver
DEFAULT_AMP_SOLVER                      = IpoptSolver(print_level = 0, max_iter = 35)
const DEFAULT_AMP_SOLVER_PRINT_PREFIX   = ""

## default values for baseline model
const DEFAULT_AMP_MODEL_DEF_FUNC        = m -> SFC_MODEL_DEFS[m]
const DEFAULT_AMP_MODEL                 = SFCModel
const DEFAULT_AMP_BL_METHOD             = l4_enl
const DEFAULT_AMP_FALLBACK_FUNC         = median
const DEFAULT_AMP_MIN_RELIABLE_CYC      = Int_T(5)         ## >= 1
const DEFAULT_AMP_BL_CYC_BOUNDS         = Vector{Int_T}()
const DEFAULT_AMP_CQ_FLUO_METHOD        = cp_dr1

## defaults for quantification model
const DEFAULT_AMP_QUANT_METHOD          = l4_enl
const DEFAULT_AMP_DENSER_FACTOR         = Int_T(3)         ## must be an integer
const DEFAULT_AMP_CQ_METHOD             = Cy0

## default for set_qt_fluos!()
const DEFAULT_AMP_QT_PROB               = Float_T(0.9)

## defaults for report_cq!()
## note for default scaled_max_dr1_lb:
## 'look like real amplification, scaled_max_dr1 0.00894855, ip223, exp. 75, well A7, channel 2`
const DEFAULT_AMP_BEFORE_128X           = false
const DEFAULT_AMP_MAX_BSF_LB            = Int_T(4356)      ## ⎫
const DEFAULT_AMP_MAX_DR1_LB            = Int_T(472)       ## ⎬ must be integers
const DEFAULT_AMP_MAX_DR2_LB            = Int_T(41)        ## ⎭
const DEFAULT_AMP_SCALED_MAX_BSF_LB     = Float_T(0.086)
const DEFAULT_AMP_SCALED_MAX_DR1_LB     = Float_T(0.0089)
const DEFAULT_AMP_SCALED_MAX_DR2_LB     = Float_T(0.000689)

## defaults for process_ad()
## see AssignGenosResult.jl

## defaults for output
const DEFAULT_AMP_OUTPUT_FORMAT         = pre_json_output
const DEFAULT_AMP_OUTPUT_OPTION         = short
const DEFAULT_AMP_MODEL_OUTPUT_STRUCT   = AmpShortModelResults
DEFAULT_AMP_REPORTER                    = roundoff(JSON_DIGITS)

#===============================================================================
    Field definitions >>
===============================================================================#

## name, DataType, default value
const AMP_FIELD_DEFS = [
    ## data
    Field(:raw,                 RawData{<: Union{Int_T,Float_T}}),
    Field(:num_cycles,          Int_T),
    Field(:num_wells,           Int_T),
    Field(:num_channels,        Int_T),
    Field(:cycles,              SVector{N,Int_T} where {N}),
    Field(:wells,               SVector{W,Symbol} where {W}),
    Field(:channels,            SVector{C,Int_T} where {C}),
    Field(:calibration_data,    CalibrationData{<: NumberOfChannels, <: Union{Int_T,Float_T}}),

    ## calibration parameters
    Field(:calibration_args,    CalibrationParameters,      DEFAULT_CAL_ARGS),

    ## solver parameters
    Field(:solver,              IpoptSolver,                DEFAULT_AMP_SOLVER),
    Field(:ipopt_print_prefix,  String,                     DEFAULT_AMP_SOLVER_PRINT_PREFIX),

    ## amplification model parameters
    Field(:amp_model,           Type{<: AmpModel},          DEFAULT_AMP_MODEL),
    Field(:SFC_model_def_func,  Function,                   DEFAULT_AMP_MODEL_DEF_FUNC),
    Field(:bl_method,           SFCModelName,               DEFAULT_AMP_BL_METHOD),
    Field(:bl_fallback_func,    Function,                   DEFAULT_AMP_FALLBACK_FUNC),
    Field(:min_reliable_cyc,    Int_T,                        DEFAULT_AMP_MIN_RELIABLE_CYC,   "min_reliable_cyc"),
    Field(:baseline_cyc_bounds, Union{Vector,Array{Vector,2}},
                                                            DEFAULT_AMP_BL_CYC_BOUNDS,      "baseline_cyc_bounds"),
    Field(:quant_method,        SFCModelName,               DEFAULT_AMP_QUANT_METHOD),
    Field(:denser_factor,       Int_T,                        DEFAULT_AMP_DENSER_FACTOR),
    Field(:cq_method,           CqMethod,                   DEFAULT_AMP_CQ_METHOD),
    Field(:qt_prob,             Float_T,                    DEFAULT_AMP_QT_PROB),
    Field(:before_128x,         Bool,                       DEFAULT_AMP_BEFORE_128X),
    Field(:max_bsf_lb,          Int_T,                        DEFAULT_AMP_MAX_BSF_LB,         "min_fluomax"),
    Field(:max_dr1_lb,          Int_T,                        DEFAULT_AMP_MAX_DR1_LB,         "min_D1max"),
    Field(:max_dr2_lb,          Int_T,                        DEFAULT_AMP_MAX_DR2_LB,         "min_D2max"),
    Field(:scaled_max_bsf_lb,   Float_T,                    DEFAULT_AMP_SCALED_MAX_BSF_LB),
    Field(:scaled_max_dr1_lb,   Float_T,                    DEFAULT_AMP_SCALED_MAX_DR1_LB),
    Field(:scaled_max_dr2_lb,   Float_T,                    DEFAULT_AMP_SCALED_MAX_DR2_LB),

    ## allelic discrimination parameters
    Field(:categ_well_vec,      Vector{Pair{Symbol,Any}},   DEFAULT_AMP_CATEG_WELL_VEC),
    Field(:ctrl_well_dict,      OrderedDict{Vector{Int_T},Vector{Int_T}},
                                                            DEFAULT_AMP_CTRL_WELL_DICT),
    ## output format parameters
    Field(:out_format,          OutputFormat,               DEFAULT_AMP_OUTPUT_FORMAT),
    Field(:amp_output,          AmpOutputOption,            DEFAULT_AMP_OUTPUT_OPTION),
    Field(:amp_model_results,   Type{<: AmpModelResults},   DEFAULT_AMP_MODEL_OUTPUT_STRUCT),
    Field(:reporting,           Function,                   DEFAULT_AMP_REPORTER)]

#===============================================================================
    macro calls >>
===============================================================================#

## generate struct and constructor
SCHEMA = AMP_FIELD_DEFS
@make_struct_from_SCHEMA AmpInput Input
@make_constructor_from_SCHEMA AmpInput

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
# # [Bending of an Initially Curved Beam](@id curved)
#
# This example is a common benchmark problem for the geometrically exact bending of
# nonlinear beams.
#
# ![](../assets/curved-drawing.svg)
#
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`curved.ipynb`](@__NBVIEWER_ROOT_URL__/examples/curved.ipynb).
#-

using GXBeam, LinearAlgebra

## problem constants
R = 100
L = R*pi/4 # inches
h = w = 1 # inches
E = 1e7 # psi Young's Modulus
ν = 0.0
G = E/(2*(1+ν))

## beam starting point, frame, and curvature
start = [0, 0, 0]
frame = [0 -1 0; 1 0 0; 0 0 1]
curvature = [0, 0, -1/R]

## cross section properties
A = h*w
Ay = A
Az = A
Iyy = w*h^3/12
Izz = w^3*h/12
J = Iyy + Izz

## discretize the beam
nelem = 16
ΔL, xp, xm, Cab = discretize_beam(L, start, nelem;
    frame = frame,
    curvature = curvature)

## force
P = 600 # lbs

## index of left and right endpoints of each beam element
pt1 = 1:nelem
pt2 = 2:nelem+1

## compliance matrix for each beam element
compliance = fill(Diagonal([1/(E*A), 1/(G*Ay), 1/(G*Az), 1/(G*J), 1/(E*Iyy),
    1/(E*Izz)]), nelem)

## create assembly of interconnected nonlinear beams
assembly = Assembly(xp, pt1, pt2, compliance=compliance, frames=Cab,
    lengths=ΔL, midpoints=xm)

## create dictionary of prescribed conditions
prescribed_conditions = Dict(
    ## fixed left endpoint
    1 => PrescribedConditions(ux=0, uy=0, uz=0, theta_x=0, theta_y=0, theta_z=0),
    ## force on right endpoint
    nelem+1 => PrescribedConditions(Fz = P)
)

## perform static analysis
system, converged = static_analysis(assembly;
    prescribed_conditions = prescribed_conditions)

## post-process results
state = AssemblyState(system, assembly;
    prescribed_conditions = prescribed_conditions)

println("Tip Displacement: ", state.points[end].u)
println("Tip Displacement (Bathe and Bolourch): [-13.4, -23.5, 53.4]")

# The calculated tip displacements match those reported by Bathe and Bolourch in "Large
# Displacement Analysis of Three-Dimensional Beam Structures" closely, thus verifying our
# implementation of geometrically exact beam theory.
#
# We can visualize the deformed geometry and inspect the associated point and element data
# using ParaView.

write_vtk("curved-visualization", assembly, state)

# ![](../assets/curved-visualization.png)


text:
using Base.Iterators: repeated

# Abstract Interface

"""
    abstract type IncrementalCycleTracker

The supertype for incremental cycle detection problems. The abstract type
constructor IncrementalCycleTracker(G) may be used to automatically select
a specific incremental cycle detection algorithm. See [`add_edge_checked!`](@ref)
for a usage example.
"""
abstract type IncrementalCycleTracker{I} <: AbstractGraph{I} end

function (::Type{IncrementalCycleTracker})(s::AbstractGraph{I}; dir=:out) where {I}
    # TODO: Once we have more algorithms, the poly-algorithm decision goes here.
    # For now, we only have Algorithm N.
    return DenseGraphICT_BFGT_N{something(dir == :in, false)}(s)
end

# Cycle Detection Interface
"""
    add_edge_checked!([f!,], ict::IncrementalCycleTracker, v, w)

Using the incremental cycle tracker, ict, check whether adding the edge `v=>w`.
Would introduce a cycle in the underlying graph. If so, return false and leave
the ict intact. If not, update the underlying graph and return true.

# Optional `f!` Argument

By default the `add_edge!` function is used to update the underlying graph.
However, for more complicated graphs, users may wish to manually specify the
graph update operation. This may be accomplished by passing the optional `f!`
callback arhgument. This callback is called on the underlying graph when no
cycle is detected and is required to modify the underlying graph in order to
effectuate the proposed edge addition.

# Batched edge additions

Optionally, either `v` or `w` (depending on the `in_out_reverse` flag) may be a
collection of vertices representing a batched addition of vertices sharing a
common source or target more efficiently than individual updates.

## Example

```jldoctest
julia> G = SimpleDiGraph(3)

julia> ict = IncrementalCycleTracker(G)
BFGT_N cycle tracker on {3, 0} directed simple Int64 graph

julia> add_edge_checked!(ict, 1, 2)
true

julia> collect(edges(G))
1-element Vector{Graphs.SimpleGraphs.SimpleEdge{Int64}}:
 Edge 1 => 2

julia> add_edge_checked!(ict, 2, 3)
true

julia> collect(edges(G))
2-element Vector{Graphs.SimpleGraphs.SimpleEdge{Int64}}:
 Edge 1 => 2
 Edge 2 => 3

julia> add_edge_checked!(ict, 3, 1) # Would add a cycle
false

julia> collect(edges(G))
2-element Vector{Graphs.SimpleGraphs.SimpleEdge{Int64}}:
Edge 1 => 2
Edge 2 => 3
```
"""
function add_edge_checked! end

to_edges(v::Integer, w::Integer) = (v=>w,)
to_edges(v::Integer, ws) = zip(repeated(v), ws)
to_edges(vs, w::Integer) = zip(vs, repeated(w))

add_edge_checked!(ict::IncrementalCycleTracker, vs, ws) = add_edge_checked!(ict, vs, ws) do g
    foreach(((v, w),)->add_edge!(g, v, w), to_edges(vs, ws))
end

# Utilities
"""
    struct TransactionalVector

A vector with one checkpoint that may be reverted to by calling `revert!`. The setpoint itself
is set by calling `commit!`.
"""
struct TransactionalVector{T} <: AbstractVector{T}
    v::Vector{T}
    log::Vector{Pair{Int, T}}
    TransactionalVector(v::Vector{T}) where {T} =
        new{T}(v, Vector{Pair{Int, T}}())
end

function commit!(v::TransactionalVector)
    empty!(v.log)
    return nothing
end

function revert!(vec::TransactionalVector)
    for (idx, val) in reverse(vec.log)
        vec.v[idx] = val
    end
    return nothing
end

function Base.setindex!(vec::TransactionalVector, val, idx)
    oldval = vec.v[idx]
    vec.v[idx] = val
    push!(vec.log, idx=>oldval)
    return nothing
end
Base.getindex(vec::TransactionalVector, idx) = vec.v[idx]
Base.size(vec::TransactionalVector) = size(vec.v)

# Specific Algorithms

const bibliography = """
## References

[BFGT15] Michael A. Bender, Jeremy T. Fineman, Seth Gilbert, and Robert E. Tarjan. 2015
    A New Approach to Incremental Cycle Detection and Related Problems.
    ACM Trans. Algorithms 12, 2, Article 14 (December 2015), 22 pages.
    DOI: http://dx.doi.org/10.1145/2756553
"""

## Bender, Algorithm N

"""
    struct DenseGraphICT_BFGT_N

Implements the "Naive" (Algorithm N) Bender-Fineman-Gilbert-Tarjan one-way line search incremental cycle detector
for dense graphs from [BFGT15] (Section 3).

$bibliography
"""
struct DenseGraphICT_BFGT_N{InOutReverse, I, G<:AbstractGraph{I}} <: IncrementalCycleTracker{I}
    graph::G
    levels::TransactionalVector{Int}
    DenseGraphICT_BFGT_N{InOutReverse}(g::G) where {InOutReverse, I, G<:AbstractGraph{I}} =
        new{InOutReverse, I, G}(g, TransactionalVector(fill(0, nv(g))))
end
function Base.show(io::IO, ict::DenseGraphICT_BFGT_N)
    print(io, "BFGT_N cycle tracker on ")
    show(io, ict.graph)
end

function topological_sort(ict::DenseGraphICT_BFGT_N{InOutReverse}) where {InOutReverse}
    # The ICT levels are a weak topological ordering, so a sort of the levels
    # will give a topological sort of the vertices.
    perm = sortperm(ict.levels)
    InOutReverse && (perm = reverse(perm))
    return perm
end

# Even when both `v` and `w` are integer, we know that `v` would come first, so
# we prefer to check for `v` as the cycle vertex in this case.
add_edge_checked!(f!, ict::DenseGraphICT_BFGT_N{false}, v::Integer, ws) =
    _check_cycle_add!(f!, ict, to_edges(v, ws), v)
add_edge_checked!(f!, ict::DenseGraphICT_BFGT_N{true}, vs, w::Integer) =
    _check_cycle_add!(f!, ict, to_edges(vs, w), w)

### [BFGT15] Algorithm N
#
# Implementation Notes
#
# This is Algorithm N from [BFGT15] (Section 3), plus limited patching support and
# a number of standard tricks. Namely:
#
# 1. Batching is supported as long as there is only a single source or destination
#    vertex. General batching is left as an open problem. The reason that the
#    single source/dest batching is easy to add is that we know that either the
#    source or the destination vertex is guaranteed to be a part of any cycle
#    that we may have added. Thus we're guaranteed to encounter one of the two
#    verticies in our cycle validation and the rest of the algorithm goes through
#    as usual.
# 2. We opportunistically traverse each edge when we see it and only add it
#    to the worklist if we know that traversal will recurse further.
# 3. We add some early out checks to detect we're about to do redundant work.
function _check_cycle_add!(f!, ict::DenseGraphICT_BFGT_N{InOutReverse}, edges, v) where {InOutReverse}
    g = ict.graph
    worklist = Pair{Int, Int}[]
    # TODO: In the case where there's a single target vertex, we could saturate
    # the level first before we assign it to the tracked vector to save some
    # log space.
    for (v, w) in edges
        InOutReverse && ((v, w) = (w, v))
        if ict.levels[v] < ict.levels[w]
            continue
        end
        v == w && return false
        ict.levels[w] = ict.levels[v] + 1
        push!(worklist, v=>w)
    end
    while !isempty(worklist)
        (x, y) = popfirst!(worklist)
        xlevel = ict.levels[x]
        ylevel = ict.levels[y]
        if xlevel >= ylevel
            # The xlevel may have been incremented further since we added this
            # edge to the worklist.
            ict.levels[y] = ylevel = xlevel + 1
        elseif ylevel > xlevel + 1
            # Some edge traversal scheduled for later already incremented this
            # level past where we would have been. Delay processing until then.
            continue
        end
        for z in (InOutReverse ? inneighbors(g, y) : outneighbors(g, y))
            if z == v
                revert!(ict.levels)
                return false
            end
            if ylevel >= ict.levels[z]
                ict.levels[z] = ylevel + 1
                push!(worklist, y=>z)
            end
        end
    end
    commit!(ict.levels)
    f!(g)
    return true
end


text:
a = 1

if a < 10 a = a + 1

println(a)


text:
#############################################################
#
# Suggested Users Functions
#
#############################################################

function get_user_suggestions_slug(slug::AbstractString; options=Dict{AbstractString, AbstractString}())

    r = get_oauth("https://api.twitter.com/1.1/users/suggestions/$(slug).json", options)

    #return to_USERS(r["users"])

    return r.status == 200 ? to_USERS(Requests.json(r)["users"]) : r

end

function get_user_suggestions_slug_members(slug::AbstractString; options=Dict{AbstractString, AbstractString}())

    r = get_oauth("https://api.twitter.com/1.1/users/suggestions/$(slug)/members.json", options)

    #return to_USERS(r)

    return r.status == 200 ? to_USERS(Requests.json(r)) : r

end


text:
using DynamicHMCModels
#=
using LinearAlgebra
using StaticArrays
using TransformVariables
using LogDensityProblems
using DynamicHMC
using Parameters
using Random
using MCMCChains
=#

function LogLikelihoodMin(ω,dyDep)

 dt = 0.1
 cD = [0 -1im ; 1im 0]
 H = (ω/2.) * cD
 rho=[0.05 0 ; 0 0.95]
 M = I - ((cD'*cD)/2) * dt + (cD * dyDep) - 1im * H * dt
 newRho = M * rho * M'
 lklhood = real(tr(newRho))- (ω* dt/2)^2
 return log(lklhood)

end

struct Experiment
      dyDep::Float64
end

function (problem::Experiment)((ω,)::NamedTuple{(:ω,)})
      @unpack dyDep  = problem        # extract the data
      LogLikelihoodMin(ω,dyDep)
end

dyDepObs=0.690691
p1 = Experiment(dyDepObs)

println(p1((ω=.4,)))

trans_single = as((ω=as(Real, 2, 4),))
P1 = TransformedLogDensity(trans_single, p1)
∇P1 = ADgradient(:ForwardDiff, P1)

# Sample 4 chains

a3d = Array{Float64, 3}(undef, 1000, 1, 4);
for j in 1:4
  global results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P1, 1000;
    reporter = NoProgressReport())
  global posterior = P1.transformation.(results.chain)

  for i in 1:1000
    a3d[i, 1, j] = values(posterior[i].ω)
  end
end

# Create MCMCChains object

parameter_names = ["ω"]
sections =   Dict(
  :parameters => parameter_names,
)
chns = create_mcmcchains(a3d, parameter_names, sections, start=1)
show(chns)

println()
DynamicHMC.Diagnostics.EBFMI(results.tree_statistics) |> display

println()
DynamicHMC.Diagnostics.summarize_tree_statistics(results.tree_statistics)


text:
module JungleHelperJungleSinkingPlatform

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Platform (Sinking, Jungle) (Jungle Helper)" => Ahorn.EntityPlacement(
        Maple.SinkingPlatform,
        "rectangle",
        Dict{String, Any}(
          "texture" => "JungleHelper/jungle"
        )

    )
)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
import LinearAlgebra

function close_to_zero(value::T) where {T<:AbstractFloat}
    return abs(value) < 1e-2
end

function se3_inv(T::Matrix{<:AbstractFloat})
    C = T[1:3, 1:3]
    d = T[1:3, 4]
    return [C' -C'*d; 0 0 0 1]
end

# function skew(x::Vector{<:AbstractFloat})
#     return [
#         0 -x[3] x[2]
#         x[3] 0 -x[1]
#         -x[2] x[1] 0
#     ]
# end

function so3_wedge(ϕ::Vector{<:AbstractFloat})
    return [
        0 -ϕ[3] ϕ[2]
        ϕ[3] 0 -ϕ[1]
        -ϕ[2] ϕ[1] 0
    ]
end

function so3_exp(psi::Vector{<:AbstractFloat})
    angle = norm(psi, 2)
    axis = psi / angle
    s = sin(angle)
    c = cos(angle)
    return c * I(3) + (1 - c) * axis .* axis' - s * so3_wedge(axis)
end

function so3_vee(C::Matrix{<:AbstractFloat})
    return [C[3, 2], C[1, 3], C[2, 1]]
end

function so3_log(C::Matrix{<:AbstractFloat})

    cos_angle = 0.5 * tr(C) - 0.5
    cos_angle = clamp(cos_angle, -1.0, 1.0)
    angle = acos(cos_angle)

    if close_to_zero(angle)
        tmp = C - I
    else
        tmp = (angle / (2 * sin(angle))) * (C - C')
    end

    return so3_vee(tmp)
end

function so3_jacobian(psi::Vector{<:AbstractFloat})
    angle = norm(psi)

    if close_to_zero(angle)
        return I + 0.5 * so3_wedge(psi)
    end

    axis = psi / angle
    s = sin(angle)
    c = cos(angle)

    return (s / angle) * I(3) + (1 - s / angle) * axis * axis' + (1 - c) / angle * so3_wedge(axis)
end

function so3_inv_left_jacobian(phi::Vector{<:AbstractFloat})
    dof = 3
    angle = norm(phi)

    # Near phi==0, use first order Taylor expansion
    if close_to_zero(angle)
        return I(3) - 0.5 * so3_wedge(phi)
    end

    axis = phi / angle
    half_angle = 0.5 * angle
    cot_half_angle = 1.0 / tan(half_angle)

    return half_angle * cot_half_angle * I(dof) + (1 - half_angle * cot_half_angle) * (axis * axis') - half_angle * so3_wedge(axis)

end

# Might have to update this
function se3_exp(d::Vector{<:AbstractFloat}, psi::Vector{<:AbstractFloat})
    J = so3_jacobian(psi)
    C = so3_exp(psi)
    return [C J*d; 0 0 0 1]
end

function se3_log(T::Matrix{<:AbstractFloat})
    psi = so3_log(T[1:3, 1:3])

    J = so3_inv_left_jacobian(psi)
    r = inv(J) * T[1:3, 4]

    return [r; psi]
end

# Functions taken from the python liegroups packages

function se3_wedge(ξ::Vector{<:AbstractFloat})
    dim = 4
    Ξ = zeors(size(ξ, 1), dim, dim)
    Ξ[:, 1:3, 1:3] = so3_wedge(ξ[:, 4:6])
    Ξ[:, 0:3, 3] = ξ[:, 1:3]
    return Ξ
end

function se3_odot(p::Vector{<:AbstractFloat})
    dim = 4
    dof = 6

    result = zeros(size(p, 1), dof)

    if size(p, 1) == dim - 1
        result[1:3, 1:3] = I(3)
        result[1:3, 4:6] = so3_wedge(-p)
    else
        result[1:3, 1:3] = I(3)
        result[1:3, 4:6] = so3_wedge(-p)
    end

    return result
end

function se3_adjoint(T::Matrix{<:AbstractFloat})
    C = T[1:3, 1:3]
    r = T[1:3, 4]
    return [
        C so3_wedge(r)*C
        zeros(3, 3) C
    ]
end

max_stars_repo_name: jandaa/state-estimation

max_issues_repo_name: jandaa/state-estimation

max_forks_repo_name: jandaa/state-estimation


text:
using Literate

function preprocess(str)
    str = replace(str, "# PREAMBLE" => "")
    str = replace(str, "# PKG_SETUP" => "")
    return str
end

exampledir = @__DIR__
Literate.notebook(joinpath(exampledir, "acrobot.jl"), exampledir, execute=false, preprocess=preprocess)
Literate.notebook(joinpath(exampledir, "cartpole.jl"), exampledir, execute=false, preprocess=preprocess)
Literate.notebook(joinpath(exampledir, "hopper.jl"), exampledir, execute=false, preprocess=preprocess)
Literate.notebook(joinpath(exampledir, "planar_push.jl"), exampledir, execute=false, preprocess=preprocess)
Literate.notebook(joinpath(exampledir, "rocket.jl"), exampledir, execute=false, preprocess=preprocess)


text:
module helloIBox
using JSON
using HTTP
using Base64
using StructArrays
using Dates
using Sockets

export  start_server, getData, getTag, getDataTree,getChildsCount,
        getAttr, getStructData,
        dg_new

include("io/pack.jl")
include("helloServer.jl")
include("IBoxAPI.jl")
include("timeseries.jl")
include("segments.jl")
include("features.jl")
include("datasets.jl")
include("undoredo.jl")
include("datagroup.jl")
include("manualmarks.jl")

end # module


text:
using Colors, DelimitedFiles, Test, Statistics
const G = GPlot

isin(s, str) = @test occursin(str, s)
notisin(s, str) = @test !occursin(str, s)

# disable continuous preview in tests
_bk_cont_preview = G.GP_ENV["CONT_PREVIEW"]
_bk_palette = G.GP_ENV["PALETTE"]
continuous_preview(false)

if VERSION < v"1.1"
    eachrow(A::AbstractArray) = (view(A, j, :) for j ∈ axes(A, 1))
end

function m2b(a,b)
    ismissing(a) && ismissing(b) && return true
    isinf(a) && isinf(b) && return true
    isnan(a) && isnan(b) && return true
    return (a == b)
end

function checkzip(z, v::VecOrMat)
    for (zi, vi) ∈ zip(z, eachrow(v))
      all(m2b(zi[j], vi[j]) for j ∈ eachindex(vi)) || (@show zi; @show vi; return false)
    end
    return true
end


text:
using Revise
using VolterraGP
using BenchmarkTools
using Flux
using Zygote
using LinearAlgebra

D = 3
C = 2
P = 2

train, test = generate_toy_data(10)
gp = GaussianProcess(scaledEQs, D, C, P, train)

# @btime negloglikelihood(gp)
# @btime gradient(negloglikelihood, gp)

# @btime VolterraGP.kernel(0.1, 0.2, 1, 1, gp)
# _, back = Zygote.pullback(VolterraGP.kernel, 0.1, 0.2, 1, 1, gp)
# @btime back(I)

# @btime VolterraGP.full_cov(0.1, 0.2, 1, 1, gp)
# _, back = Zygote.pullback(VolterraGP.full_cov, 0.1, 0.2, 1, 1, gp)
# @btime back(I)

# @btime sum(VolterraGP.get_phi_cov(0.1, 0.2, 2, 2, 1, 1, gp))
# @btime gradient(pars -> sum(VolterraGP.get_phi_cov(pars...)), (0.1, 0.2, 2, 2, 1, 1, gp))

# @btime sum(VolterraGP.fill_phi(gp.base_kernel, ones(4), ones(4, 2), [1., 1.]))
# @btime gradient(pars -> sum(VolterraGP.fill_phi(pars...)), (gp.base_kernel, ones(4), ones(4, 2), [0.1]))

# phi = VolterraGP.get_phi_cov(0.1, 0.2, 2, 2, 1, 1, gp)

# @btime VolterraGP.kan_rv_prod(phi)
# @btime gradient(VolterraGP.kan_rv_prod, phi)

# @btime VolterraGP.scaledEQs(0.1, 0.2, ones(2), ones(2), [1.])
# @btime gradient(VolterraGP.scaledEQs, 0.1, 0.2, ones(2), ones(2), [1.])
# println(VolterraGP.scaledEQs_adj(2., 0.1, 0.2, ones(2), ones(2), [1.]))

# # a = Zygote.forward_jacobian(negloglikelihood, gp)
# function sig_like(sig, gp)
#     gp.dpars.σ = sig
#     negloglikelihood(gp)
# end
# a = Zygote.forward_jacobian(x -> sig_like(x, gp), ones(3))
# f = x -> x' * x


text:
module Cerner

import HTTP
import JSON3
import Base64

include("system_accounts.jl")

end # module


text:
"""
$SIGNATURES

Find footnotes refs and defs and eliminate the ones that don't verify the appropriate regex.
For a footnote ref: `\\[\\^[a-zA-Z0-0]+\\]` and `\\[\\^[a-zA-Z0-0]+\\]:` for the def.
"""
function validate_footnotes!(tokens::Vector{Token})
    fn_refs = Vector{Token}()
    rm      = Int[]
    for (i, τ) in enumerate(tokens)
        τ.name == :FOOTNOTE_REF || continue
        # footnote ref [^1]:
        m = match(r"^\[\^[a-zA-Z0-9]+\](:)?$", τ.ss)
        if !isnothing(m)
            if !isnothing(m.captures[1])
                # it's a def
                tokens[i] = Token(:FOOTNOTE_DEF, τ.ss)
            end
            # otherwise it's a ref, leave as is
        else
            # delete
            push!(rm, i)
        end
    end
    deleteat!(tokens, rm)
    return nothing
end

"""
$SIGNATURES

Verify that a given string corresponds to a well formed html entity.
"""
function validate_html_entity(ss::AS)
    !isnothing(match(r"&(?:[a-z0-9]+|#[0-9]{1,6}|#x[0-9a-f]{1,6});", ss))
end

"""
$(SIGNATURES)

Given a candidate header block, check that the opening `#` is at the start of a line, otherwise
ignore the block.
"""
function validate_headers!(tokens::Vector{Token})::Nothing
    isempty(tokens) && return
    s = str(tokens[1].ss) # does not allocate
    rm = Int[]
    for (i, τ) in enumerate(tokens)
        τ.name in MD_HEADER_OPEN || continue
        # check if it overlaps with the first character
        fromτ = from(τ)
        fromτ == 1 && continue
        # otherwise check if the previous character is a linereturn
        prevc = s[prevind(s, fromτ)]
        prevc == '\n' && continue
        push!(rm, i)
    end
    deleteat!(tokens, rm)
    return
end

postprocess_link(s::String) = replace(r"")

"""
$(SIGNATURES)

Keep track of link defs.
"""
function validate_and_store_link_defs!(blocks::Vector{OCBlock})::Nothing
    isempty(blocks) && return
    rm = Int[]
    parent = str(blocks[1])
    for (i, β) in enumerate(blocks)
        if β.name == :LINK_DEF
            # incremental backward look until we find a `[` or a `\n` if `\n` first, discard
            ini  = prevind(parent, from(β))
            k    = ini
            char = '\n'
            while k ≥ 1
                char = parent[k]
                char ∈ ('[','\n') && break
                k = prevind(parent, k)
            end
            if char == '['
                # redefine the full block
                ftk = Token(:FOO,subs(""))
                # we have a [id]: lk add it to PAGE_LINK_DEFS
                id = subs(parent, nextind(parent, k), ini)
                # issue #266 in case there's formatting in the link
                id = jd2html(id, internal=true)
                id = replace(id, r"^<p>"=>"")
                id = replace(id, r"<\/p>\n$"=>"")
                lk = β |> content |> strip |> string
                PAGE_LINK_DEFS[id] = lk
                # replace the block by a full one so that it can be fully
                # discarded in the process of md blocks
                blocks[i] = OCBlock(:LINK_DEF, ftk=>ftk, subs(parent, k, to(β)))
            else
                # discard
                push!(rm, i)
            end
        end
    end
    deleteat!(blocks, rm)
    return nothing
end