text:
export densitymatrix, ket

"""
    densitymatrix(net, prob, norm=true)

Returns the Density matrix encoded by the neural network `net`, and normalizes
it if `norm==true`.
"""
function densitymatrix(net::NeuralNetwork, hilb::AbstractSuperOpBasis, norm=true)
    if !(net isa CachedNet)
        net = cached(net)
    end

    v = state(hilb, net)

    p_hilb = physical(hilb)
    ρ = zeros(out_type(net), spacedimension(p_hilb), spacedimension(p_hilb))
    #if v isa DiagonalStateWrapper
    #    v = v.parent
    #end
    for i=1:spacedimension(p_hilb)
        set_index!(row(v), p_hilb, i)
        for j=1:spacedimension(p_hilb)
            set_index!(col(v), p_hilb, j)
            ρ[i,j] = exp(net(v))
        end
    end
    if norm
        ρ ./= tr(ρ)
    end
    ρ
end

densitymatrix(net::NeuralNetwork, hilb::AbstractHilbert, norm=true) =
    densitymatrix(net, SuperOpSpace(hilb), norm)

Base.Matrix(net::NeuralNetwork, hilb, norm=true) = densitymatrix(net, hilb, norm)
Base.Vector(net::MatrixNet, hilb, norm=true) = vec(densitymatrix(net, hilb, norm))

"""
    ket(net, prob, norm=true)

Returns the state (ket) encoded by the neural network `net`, and normalizes
it if `norm==true`.
"""
function ket(net, hilb::AbstractHilbert, norm=true)
    if !(net isa CachedNet)
        net = cached(net)
    end

    psi = zeros(out_type(net), spacedimension(hilb))

    v = state(hilb, net)
    for i=1:spacedimension(hilb)
        set_index!(v, hilb, i)
        psi[i] = exp(net(v))
    end
    if norm
        normalize!(psi)
    end
    return psi
end

ket(net::MatrixNet, hilb::AbstractSuperOpBasis, norm=true) =
    ket(net, hilb, norm)

ket(net::MatrixNet, hilb::AbstractHilbert, norm=true) =
    ket(net, SuperOpSpace(hilb), norm)

Base.Vector(net::NeuralNetwork, hilb, norm=true) = ket(net, hilb, norm)

QuantumOpticsBase.DenseOperator(net::NeuralNetwork, hilb::AbstractHilbert, norm=true) =
    DenseOperator(convert(CompositeBasis, hilb), Matrix(net, hilb, norm))

QuantumOpticsBase.Ket(hilb::AbstractHilbert, net::KetNeuralNetwork, norm=true) =
    Ket(convert(CompositeBasis, hilb), ket(net, hilb, norm))

QuantumOpticsBase.expect(Ô::AbsLinearOperator, state) =
    expect(SparseOperator(Ô), state)
##


text:
module TestIO
    using Base.Test
    using DataFrames, Compat

    #test_group("We can read various file types.")

    data = joinpath(dirname(@__FILE__), "data")

    filenames = ["$data/blanklines/blanklines.csv",
                 "$data/compressed/movies.csv.gz",
                 "$data/newlines/embedded_os9.csv",
                 "$data/newlines/embedded_osx.csv",
                 "$data/newlines/embedded_windows.csv",
                 "$data/padding/space_after_delimiter.csv",
                 "$data/padding/space_around_delimiter.csv",
                 "$data/padding/space_before_delimiter.csv",
                 "$data/quoting/empty.csv",
                 "$data/quoting/escaping.csv",
                 "$data/quoting/quotedcommas.csv",
                 "$data/scaling/10000rows.csv",
                 "$data/utf8/corrupt_utf8.csv",
                 "$data/utf8/short_corrupt_utf8.csv",
                 "$data/utf8/utf8.csv"]

    for filename in filenames
        try
            df = readtable(filename)
        catch
            error(@sprintf "Failed to read %s\n" filename)
        end
    end

    #test_group("We get the right size, types, values for a basic csv.")

    filename = "$data/scaling/movies.csv"
    df = readtable(filename)

    @test size(df) == (58788, 25)

    @test df[1, 1] === 1
    @test df[1, 2] == "\$"
    @test df[1, 3] === 1971
    @test df[1, 4] === 121
    @test df[1, 5] === NA
    @test df[1, 6] === 6.4
    @test df[1, 7] === 348
    @test df[1, 8] === 4.5
    @test df[1, 9] === 4.5
    @test df[1, 10] === 4.5
    @test df[1, 11] === 4.5
    @test df[1, 12] === 14.5
    @test df[1, 13] === 24.5
    @test df[1, 14] === 24.5
    @test df[1, 15] === 14.5
    @test df[1, 16] === 4.5
    @test df[1, 17] === 4.5
    @test df[1, 18] == ""
    @test df[1, 19] === 0
    @test df[1, 20] === 0
    @test df[1, 21] === 1
    @test df[1, 22] === 1
    @test df[1, 23] === 0
    @test df[1, 24] === 0
    @test df[1, 25] === 0

    @test df[end, 1] === 58788
    @test df[end, 2] == "xXx: State of the Union"
    @test df[end, 3] === 2005
    @test df[end, 4] === 101
    @test df[end, 5] === 87000000
    @test df[end, 6] === 3.9
    @test df[end, 7] === 1584
    @test df[end, 8] === 24.5
    @test df[end, 9] === 4.5
    @test df[end, 10] === 4.5
    @test df[end, 11] === 4.5
    @test df[end, 12] === 4.5
    @test df[end, 13] === 14.5
    @test df[end, 14] === 4.5
    @test df[end, 15] === 4.5
    @test df[end, 16] === 4.5
    @test df[end, 17] === 14.5
    @test df[end, 18] == "PG-13"
    @test df[end, 19] === 1
    @test df[end, 20] === 0
    @test df[end, 21] === 0
    @test df[end, 22] === 0
    @test df[end, 23] === 0
    @test df[end, 24] === 0
    @test df[end, 25] === 0

    #test_group("readtable handles common separators and infers them from extensions.")

    df1 = readtable("$data/separators/sample_data.csv")
    df2 = readtable("$data/separators/sample_data.tsv")
    df3 = readtable("$data/separators/sample_data.wsv")
    df4 = readtable("$data/separators/sample_data_white.txt", separator = ' ')

    @test df1 == df2
    @test df2 == df3
    @test df3 == df4

    readtable("$data/quoting/quotedwhitespace.txt", separator = ' ')

    #test_group("readtable handles common newlines.")

    df = readtable("$data/newlines/os9.csv")
    @test isequal(readtable("$data/newlines/osx.csv"), df)
    @test isequal(readtable("$data/newlines/windows.csv"), df)

    @test isequal(df, readtable("$data/newlines/os9.csv", skipblanks = false))
    @test isequal(df, readtable("$data/newlines/osx.csv", skipblanks = false))
    @test isequal(df, readtable("$data/newlines/windows.csv", skipblanks = false))

    #test_group("readtable treats rows as specified.")

    df1 = readtable("$data/comments/before_after_data.csv", allowcomments = true)
    df2 = readtable("$data/comments/middata.csv", allowcomments = true)
    df3 = readtable("$data/skiplines/skipfront.csv", skipstart = 3)
    df4 = readtable("$data/skiplines/skipfront.csv", skipstart = 4, header = false)
    names!(df4, names(df1))
    df5 = readtable("$data/comments/before_after_data_windows.csv", allowcomments = true)
    df6 = readtable("$data/comments/middata_windows.csv", allowcomments = true)
    df7 = readtable("$data/skiplines/skipfront_windows.csv", skipstart = 3)
    df8 = readtable("$data/skiplines/skipfront_windows.csv", skipstart = 4, header = false)
    names!(df8, names(df1))
    # df9 = readtable("$data/skiplines/skipfront.csv", skipstart = 3, skiprows = 5:6)
    # df10 = readtable("$data/skiplines/skipfront.csv", skipstart = 3, header = false, skiprows = [4, 6])
    # names!(df10, names(df1))

    @test df2 == df1
    @test df3 == df1
    @test df4 == df1

    # Windows EOLS
    @test df5 == df1
    @test df6 == df1
    @test df7 == df1
    @test df8 == df1

    # @test df9 == df1[3:end]
    # @test df10 == df1[[1, 3:end]]

    function normalize_eol!(df)
        for (name, col) in eachcol(df)
            if eltype(col) <: String
                df[name] = map(s -> replace(s, "\r\n", "\n"), col)
            end
        end
        df
    end

    osxpath = "$data/skiplines/complex_osx.csv"
    winpath = "$data/skiplines/complex_windows.csv"

    opts1 = @compat Dict{Any,Any}(:allowcomments => true)
    opts2 = @compat Dict{Any,Any}(:skipstart => 4, :skiprows => [6, 7, 12, 14, 17], :skipblanks => false)

    df1 = readtable(osxpath; opts1...)
    # df2 = readtable(osxpath; opts2...)
    df1w = readtable(winpath; opts1...)
    # df2w = readtable(winpath; opts2...)

    # @test df2 == df1
    @test normalize_eol!(df1w) == df1
    # @test normalize_eol!(df2w) == df1

    opts1[:nrows] = 3
    opts2[:nrows] = 3

    @test readtable(osxpath; opts1...) == df1[1:3, :]
    # @test readtable(osxpath; opts2...) == df1[1:3, :]
    @test normalize_eol!(readtable(winpath; opts1...)) == df1[1:3, :]
    # @test readtable(winpath; opts2...) == df1[1:3, :]

    # opts2[:header] = false
    # opts2[:skipstart] = 5

    # df2b = readtable(path; opts2...)
    # names!(df2b, names(df1))

    # @test df2b == df1[1:3]

    #test_group("readtable handles custom delimiters.")

    readtable("$data/skiplines/skipfront.csv", allowcomments = true, commentmark = '%')

    readtable("$data/separators/sample_data.csv", quotemark = Char[])
    @test_throws ErrorException readtable("$data/newlines/embedded_osx.csv", quotemark = Char[])
    df = readtable("$data/quoting/single.csv", quotemark = ['\''])
    @test df == readtable("$data/quoting/mixed.csv", quotemark = ['\'', '"'])

    # df = readtable("$data/decimal/period.csv")
    # @test df[2, :A] == 0.3
    # @test df[2, :B] == 4.0

    # @test df == readtable("$data/decimal/comma.tsv", decimal = ',')

    #test_group("readtable column names.")

    ns = [:Var1, :Var2, :Var3, :Var4, :Var5]
    df = readtable("$data/typeinference/mixedtypes.csv")
    names!(df, ns)
    @test df == readtable("$data/typeinference/mixedtypes.csv", names = ns)

    df = readtable("$data/separators/sample_data.csv", header = false, names = ns[1:3])
    @test df[1, :Var1] == 0
    df = readtable("$data/separators/sample_data.csv", names = ns[1:3])
    @test df[1, :Var1] == 1

    #test_group("Properties of data frames returned by readtable method.")

    # Readtable ignorepadding
    io = IOBuffer("A , \tB  , C\n1 , \t2, 3\n")
    @test readtable(io, ignorepadding = true) == DataFrame(A = 1, B = 2, C = 3)

    # Readtable c-style escape options

    df = readtable("$data/escapes/escapes.csv", allowescapes = true)
    @test df[1, :V] == "\t\r\n"
    @test df[2, :V] == "\\\\t"
    @test df[3, :V] == "\\\""

    df = readtable("$data/escapes/escapes.csv")
    @test df[1, :V] == "\\t\\r\\n"
    @test df[2, :V] == "\\\\t"
    @test df[3, :V] == "\\\""

    # df = readtable("$data/escapes/escapes.csv", escapechars = ['"'], nrows = 2)
    # @test df[1, :V] == "\\t\\r\\n"
    # @test df[2, :V] == "\\\\\\\\t"

    # Readtable with makefactors active should only make factors from columns
    # of strings.
    filename = "$data/factors/mixedvartypes.csv"
    df = readtable(filename, makefactors = true)

    @test typeof(df[:factorvar]) == PooledDataArray{UTF8String,UInt32,1}
    @test typeof(df[:floatvar]) == DataArray{Float64,1}

    # Readtable shouldn't silently drop data when reading highly compressed gz.
    df = readtable("$data/compressed/1000x2.csv.gz")
    @test size(df) == (1000, 2)

    # Readtable type inference
    filename = "$data/typeinference/bool.csv"
    df = readtable(filename)
    @test typeof(df[:Name]) == DataArray{UTF8String,1}
    @test typeof(df[:IsMale]) == DataArray{Bool,1}
    @test df[:IsMale][1] == true
    @test df[:IsMale][4] == false

    filename = "$data/typeinference/standardtypes.csv"
    df = readtable(filename)
    @test typeof(df[:IntColumn]) == DataArray{Int,1}
    @test typeof(df[:IntlikeColumn]) == DataArray{Float64,1}
    @test typeof(df[:FloatColumn]) == DataArray{Float64,1}
    @test typeof(df[:BoolColumn]) == DataArray{Bool,1}
    @test typeof(df[:StringColumn]) == DataArray{UTF8String,1}

    filename = "$data/typeinference/mixedtypes.csv"
    df = readtable(filename)
    @test typeof(df[:c1]) == DataArray{UTF8String,1}
    @test df[:c1][1] == "1"
    @test df[:c1][2] == "2.0"
    @test df[:c1][3] == "true"
    @test typeof(df[:c2]) == DataArray{Float64,1}
    @test df[:c2][1] == 1.0
    @test df[:c2][2] == 3.0
    @test df[:c2][3] == 4.5
    @test typeof(df[:c3]) == DataArray{UTF8String,1}
    @test df[:c3][1] == "0"
    @test df[:c3][2] == "1"
    @test df[:c3][3] == "f"
    @test typeof(df[:c4]) == DataArray{Bool,1}
    @test df[:c4][1] == true
    @test df[:c4][2] == false
    @test df[:c4][3] == true
    @test typeof(df[:c5]) == DataArray{UTF8String,1}
    @test df[:c5][1] == "False"
    @test df[:c5][2] == "true"
    @test df[:c5][3] == "true"

    # Readtable defining column types
    filename = "$data/definedtypes/mixedvartypes.csv"

    df = readtable(filename)
    @test typeof(df[:n]) == DataArray{Int,1}
    @test df[:n][1] == 1
    @test typeof(df[:s]) == DataArray{UTF8String,1}
    @test df[:s][1] == "text"
    @test typeof(df[:f]) == DataArray{Float64,1}
    @test df[:f][1] == 2.3
    @test typeof(df[:b]) == DataArray{Bool,1}
    @test df[:b][1] == true

    df = readtable(filename, eltypes = [Int64, UTF8String, Float64, Bool])
    @test typeof(df[:n]) == DataArray{Int64,1}
    @test df[:n][1] == 1
    @test typeof(df[:s]) == DataArray{UTF8String,1}
    @test df[:s][1] == "text"
    @test df[:s][4] == "text ole"
    @test typeof(df[:f]) == DataArray{Float64,1}
    @test df[:f][1] == 2.3
    @test typeof(df[:b]) == DataArray{Bool,1}
    @test df[:b][1] == true
    @test df[:b][2] == false

    df = readtable(filename, eltypes = [Int64, UTF8String, Float64, UTF8String])
    @test typeof(df[:n]) == DataArray{Int64,1}
    @test df[:n][1] == 1.0
    @test isna(df[:s][3])
    @test typeof(df[:f]) == DataArray{Float64,1}
    # Float are not converted to int
    @test df[:f][1] == 2.3
    @test df[:f][2] == 0.2
    @test df[:f][3] == 5.7
    @test typeof(df[:b]) == DataArray{UTF8String,1}
    @test df[:b][1] == "T"
    @test df[:b][2] == "FALSE"

    # Readtable name normalization
    abnormal = "\u212b"
    ns = [:Å, :_B_C_, :_end]
    @test !in(symbol(abnormal), ns)

    io = IOBuffer(abnormal*",%_B*\tC*,end\n1,2,3\n")
    @test names(readtable(io)) == ns

    # Test writetable with NA and compare to the results
    tf = tempname()
    isfile(tf) && rm(tf)
    df = DataFrame(A = @data([1,NA]), B = @data(["b", NA]))
    writetable(tf, df)
    @test readcsv(tf) == ["A" "B"; 1 "b"; "NA" "NA"]

    # Test writetable with nastring set and compare to the results
    isfile(tf) && rm(tf)
    writetable(tf, df, nastring="none")
    @test readcsv(tf) == ["A" "B"; 1 "b"; "none" "none"]

    # Test writetable with append
    df1 = DataFrame(a = @data([1, 2, 3]), b = @data([4, 5, 6]))
    df2 = DataFrame(a = @data([1, 2, 3]), b = @data([4, 5, 6]))
    df3 = DataFrame(a = @data([1, 2, 3]), c = @data([4, 5, 6])) # 2nd column mismatch
    df3b = DataFrame(a = @data([1, 2, 3]), b = @data([4, 5, 6]), c = @data([4, 5, 6])) # number of columns mismatch

    # Would use joinpath(tempdir(), randstring()) to get around tempname
    # creating a file on Windows, but Julia 0.3 has no srand() to unset the
    # seed set in test/data.jl -- annoying for local testing.
    tf = tempname()
    isfile(tf) && rm(tf)

    # Written as normal if file doesn't exist
    writetable(tf, df1, append = true)
    @test readtable(tf) == df1

    # Written as normal if file is empty
    open(io -> print(io, ""), tf, "w")
    writetable(tf, df1, append = true)
    @test readtable(tf) == df1

    # Appends to existing file if append == true
    writetable(tf, df1)
    writetable(tf, df2, header = false, append = true)
    @test readtable(tf) == vcat(df1, df2)

    # Overwrites file if append == false
    writetable(tf, df1)
    writetable(tf, df2)
    @test readtable(tf) == df2

    # Enforces matching column names iff append == true && header == true
    writetable(tf, df2)
    @test_throws KeyError writetable(tf, df3, append = true)
    writetable(tf, df3, header = false, append = true)

    # Enforces matching column count if append == true
    writetable(tf, df3)
    @test_throws DimensionMismatch writetable(tf, df3b, header = false, append = true)
end

ext: jl


text:
n, p = 500, 5
((X, y, θ), (X1, y1, θ1)) = generate_binary(n, p; seed=52661)

@testset "Logreg" begin
    # No intercept
    λ = 5.0
    lr = LogisticRegression(λ; fit_intercept=false)
    J  = objective(lr, X, y)
    o  = LogisticLoss() + λ * L2Penalty()
    @test J(θ) == o(y, X*θ, θ)
    θ_newton   = fit(lr, X, y, solver=Newton())
    θ_newtoncg = fit(lr, X, y, solver=NewtonCG())
    θ_lbfgs    = fit(lr, X, y, solver=LBFGS())
    @test isapprox(J(θ),          191.949109, rtol=1e-5)
    @test isapprox(J(θ_newton),   184.704499, rtol=1e-5)
    @test isapprox(J(θ_newtoncg), 184.704499, rtol=1e-5)
    @test isapprox(J(θ_lbfgs),    184.704499, rtol=1e-5)

    # With intercept
    lr1 = LogisticRegression(λ; penalize_intercept=true)
    J   = objective(lr1, X, y1)
    θ_newton   = fit(lr1, X, y1, solver=Newton())
    θ_newtoncg = fit(lr1, X, y1, solver=NewtonCG())
    θ_lbfgs    = fit(lr1, X, y1, solver=NewtonCG())
    @test isapprox(J(θ1),         236.95486, rtol=1e-5)
    @test isapprox(J(θ_newton),   230.47826, rtol=1e-5)
    @test isapprox(J(θ_newtoncg), 230.47826, rtol=1e-5)
    @test isapprox(J(θ_lbfgs),    230.47826, rtol=1e-5)

    # with intercept and not penalizing it
    lr1 = LogisticRegression(λ)
    J   = objective(lr1, X, y1)
    θ_newton   = fit(lr1, X, y1, solver=Newton())
    θ_newtoncg = fit(lr1, X, y1, solver=NewtonCG())
    θ_lbfgs    = fit(lr1, X, y1, solver=NewtonCG())
    @test isapprox(J(θ1),         231.44778, rtol=1e-5)
    @test isapprox(J(θ_newton),   226.44523, rtol=1e-5)
    @test isapprox(J(θ_newtoncg), 226.44523, rtol=1e-5)
    @test isapprox(J(θ_lbfgs),    226.44523, rtol=1e-5)

    if DO_COMPARISONS
        # This checks that the parameters recovered using Sklearn lead
        # to a similar loss than the one given by our code to verify the
        # correctness of the code.
        lr_sk_ncg = SKLEARN_LM.LogisticRegression(
                        C=1.0/λ, solver="newton-cg", random_state=155151)
        lr_sk_ncg.fit(X, y1)
        θ1_sk_ncg = vcat(lr_sk_ncg.coef_[:], lr_sk_ncg.intercept_)
        lr_sk_lbfgs = SKLEARN_LM.LogisticRegression(
                        C=1.0/λ, solver="lbfgs", random_state=155151)
        lr_sk_lbfgs.fit(X, y1)
        θ1_sk_lbfgs = vcat(lr_sk_lbfgs.coef_[:], lr_sk_lbfgs.intercept_)
        # Comparing with ours
        @test isapprox(J(θ1_sk_ncg),   226.4, rtol=1e-2)
        @test isapprox(J(θ1_sk_lbfgs), 226.4, rtol=1e-2)
    end
end

n, p, c = 500, 5, 4
((X, y, θ), (X1, y1, θ1)) = generate_multiclass(n, p, c; seed=525)

@testset "Multinomial" begin
    # No intercept
    λ = 5.0
    mnr = MultinomialRegression(λ; fit_intercept=false)
    J   = objective(mnr, X, y; c=c)
    θ_newtoncg = fit(mnr, X, y, solver=NewtonCG())
    θ_lbfgs    = fit(mnr, X, y, solver=R.LBFGS())
    @test isapprox(J(θ),          428.25787, rtol=1e-5)
    @test isapprox(J(θ_newtoncg), 408.90295, rtol=1e-5)
    @test isapprox(J(θ_lbfgs),    408.90295, rtol=1e-5)

    #  With intercept
    mnr = MultinomialRegression(λ; penalize_intercept=true)
    J   = objective(mnr, X, y1; c=c)
    θ_newtoncg = fit(mnr, X, y1, solver=NewtonCG())
    θ_lbfgs    = fit(mnr, X, y1, solver=R.LBFGS())
    @test isapprox(J(θ),          1415.66534, rtol=1e-5)
    @test isapprox(J(θ_newtoncg),  375.95659, rtol=1e-5)
    @test isapprox(J(θ_lbfgs),     375.95659, rtol=1e-5)

    mnr = MultinomialRegression(λ)
    J   = objective(mnr, X, y1; c=c)
    θ_newtoncg = fit(mnr, X, y1, solver=NewtonCG())
    θ_lbfgs    = fit(mnr, X, y1, solver=R.LBFGS())
    @test isapprox(J(θ),          1414.01629, rtol=1e-5)
    @test isapprox(J(θ_newtoncg),  372.42948, rtol=1e-5)
    @test isapprox(J(θ_lbfgs),     372.42948, rtol=1e-5)

    if DO_COMPARISONS
        lr_sk_ncg = SKLEARN_LM.LogisticRegression(
                        C=1.0/λ, solver="newton-cg", multi_class="multinomial",
                        random_state=551551)
        lr_sk_ncg.fit(X, y1)
        θ1_sk_ncg = vec(vcat(lr_sk_ncg.coef_', lr_sk_ncg.intercept_'))
        lr_sk_lbfgs = SKLEARN_LM.LogisticRegression(
                        C=1.0/λ, solver="lbfgs", multi_class="multinomial",
                        random_state=551551)
        lr_sk_lbfgs.fit(X, y1)
        θ1_sk_lbfgs = vec(vcat(lr_sk_ncg.coef_', lr_sk_ncg.intercept_'))
        # Comparing with ours
        @test isapprox(J(θ1_sk_ncg),   374.697, rtol=1e-3)
        @test isapprox(J(θ1_sk_lbfgs), 374.697, rtol=1e-3)
    end
end

n, p = 500, 100
(_, (X, y, θ)) = generate_binary(n, p;  seed=52551, sparse=0.1)

@testset "Logreg/EN" begin
    ρ = 0.8
    α = 0.03
    λ = α * (1 - ρ) * n
    γ = α * ρ * n
    enlr = LogisticRegression(λ, γ; penalize_intercept=true)
    J    = objective(enlr, X, y)
    θ_fista = fit(enlr, X, y)
    θ_ista  = fit(enlr, X, y, solver=ISTA())
    @test isapprox(J(θ),       273.58196, rtol=1e-5)
    @test isapprox(J(θ_fista), 226.92630, rtol=1e-5)
    @test isapprox(J(θ_ista),  226.92627, rtol=1e-5)
    # sparsity
    @test nnz(θ)       == 8  # generating θ
    @test nnz(θ_fista) == 15
    @test nnz(θ_ista)  == 15

    # pure l1 regularization (unclear how sklearn does the mixing of l1/l2 for logreg)
    enlr = LogisticRegression(γ; penalty=:l1)
    J    = objective(enlr, X, y)
    θ_fista = fit(enlr, X, y)
    @test isapprox(J(θ),       255.83313, rtol=1e-5)
    @test isapprox(J(θ_fista), 218.79232, rtol=1e-5)
    @test nnz(θ_fista) == 14

    if DO_COMPARISONS
        enlr_sk = SKLEARN_LM.LogisticRegression(
                    penalty="elasticnet", C=1.0/γ, l1_ratio=1, solver="saga",
                    random_state=6618)
        enlr_sk.fit(X, y)
        θ_sk = vcat(enlr_sk.coef_[:], enlr_sk.intercept_)
        @test isapprox(J(θ_sk), 218.79, rtol=1e-3)
        @test nnz(θ_sk) == 14
    end
end

n, p, c = 1000, 100, 3
(_, (X, y, θ)) = generate_multiclass(n, p, c;  seed=33, sparse=0.1)

@testset "Multin/EN" begin
    λ = 10
    γ = 50
    enmnr = MultinomialRegression(λ, γ; penalize_intercept=true)
    J     = objective(enmnr, X, y; c=c)
    θ_fista = fit(enmnr, X, y)
    θ_ista  = fit(enmnr, X, y)
    @test isapprox(J(θ),       1614.90802, rtol=1e-5)
    @test isapprox(J(θ_fista),  931.85109, rtol=1e-5)
    @test isapprox(J(θ_ista),   931.85109, rtol=1e-5)
    @test nnz(θ_fista) == 14
    @test nnz(θ_ista)  == 14

    # pure l1 regularization for sklearn comp
    enmnr = MultinomialRegression(γ; penalty=:l1)
    J     = objective(enmnr, X, y; c=c)
    θ_fista = fit(enmnr, X, y)
    @test isapprox(J(θ),       1492.82898, rtol=1e-5)
    @test isapprox(J(θ_fista),  917.72513, rtol=1e-5)
    @test nnz(θ_fista) == 15

    if DO_COMPARISONS
        enmnr_sk = SKLEARN_LM.LogisticRegression(
                    penalty="elasticnet", C=1.0/γ, l1_ratio=1, solver="saga", multi_class="multinomial", random_state=1616)
        enmnr_sk.fit(X, y)
        θ_sk = enmnr_sk.coef_'[:]
        @test isapprox(J(θ_sk), 921.54, rtol=1e-3)
        @test nnz(θ_sk) == 14
    end
end


text:
function l0EM_orig(X, y, lambda = 2;
  epsilon = 1e-4, deltaTreshold = 1e-5,
  nonnegative = false, maxIterations = 10000)

  N, M = size(X)
  lambda_eye = lambda * eye(N)
  xt = X'
  theta = xt * inv(X * xt .+ lambda_eye) * y

  if nonnegative
    set_negative_to_zero!(theta)
  end

  iterations = 0

  while iterations < maxIterations
    iterations += 1

    # E-step:
    eta = theta

    # M-step:
    eta_squared = eta .^ 2
    xt_eta = broadcast(*, eta_squared, xt)
    theta = xt_eta * inv(X * xt_eta .+ lambda_eye) * y

    # Ensure non-negativity if required
    if nonnegative
      set_negative_to_zero!(theta)
    end

    # We have converged if change is too small => break.
    if norm(theta - eta) < deltaTreshold
      break
    end
  end

  theta[abs(theta) .< epsilon] = 0.0

  theta
end

function l0EM_us(X, y, lambda = 2;
  epsilon = 1e-4, deltaTreshold = 1e-5,
  nonnegative = false, maxIterations = 10000)

  N, M = size(X)
  lambda_eye = UniformScaling(lambda)
  xt = X'
  theta = xt * inv(X * xt .+ lambda_eye) * y

  if nonnegative
    set_negative_to_zero!(theta)
  end

  iterations = 0

  while iterations < maxIterations
    iterations += 1

    # E-step:
    eta = theta

    # M-step:
    eta_squared = eta .^ 2
    xt_eta = broadcast(*, eta_squared, xt)
    theta = xt_eta * inv(X * xt_eta .+ lambda_eye) * y

    # Ensure non-negativity if required
    if nonnegative
      set_negative_to_zero!(theta)
    end

    # We have converged if change is too small => break.
    if norm(theta - eta) < deltaTreshold
      break
    end
  end

  theta[abs(theta) .< epsilon] = 0.0

  theta
end

function l0EM_solve(X, y, lambda = 2;
  epsilon = 1e-4, deltaTreshold = 1e-5,
  nonnegative = false, maxIterations = 10000)

  N, M = size(X)
  lambda_eye = lambda * eye(N)
  xt = X'
  theta = xt * ((X * xt .+ lambda_eye) \ y)

  if nonnegative
    set_negative_to_zero!(theta)
  end

  iterations = 0

  while iterations < maxIterations
    iterations += 1

    # E-step:
    eta = theta

    # M-step:
    eta_squared = eta .^ 2
    xt_eta = broadcast(*, eta_squared, xt)
    theta = xt_eta * ((X * xt_eta .+ lambda_eye) \ y)

    # Ensure non-negativity if required
    if nonnegative
      set_negative_to_zero!(theta)
    end

    # We have converged if change is too small => break.
    if norm(theta - eta) < deltaTreshold
      break
    end
  end

  theta[abs(theta) .< epsilon] = 0.0

  theta
end

function l0EM_us_solve(X, y, lambda = 2;
  epsilon = 1e-4, deltaTreshold = 1e-5,
  nonnegative = false, maxIterations = 10000)

  N, M = size(X)
  lambda_eye = UniformScaling(lambda)
  xt = X'
  theta = xt * ((X * xt .+ lambda_eye) \ y)

  if nonnegative
    set_negative_to_zero!(theta)
  end

  iterations = 0

  while iterations < maxIterations
    iterations += 1

    # E-step:
    eta = theta

    # M-step:
    eta_squared = eta .^ 2
    xt_eta = broadcast(*, eta_squared, xt)
    theta = xt_eta * ((X * xt_eta .+ lambda_eye) \ y)

    # Ensure non-negativity if required
    if nonnegative
      set_negative_to_zero!(theta)
    end

    # We have converged if change is too small => break.
    if norm(theta - eta) < deltaTreshold
      break
    end
  end

  theta[abs(theta) .< epsilon] = 0.0

  theta
end

function l0EM_adaptus_solve(X, y, lambda = 2;
  epsilon = 1e-4, deltaTreshold = 1e-5,
  nonnegative = false, maxIterations = 10000)

  N, M = size(X)
  if N >= M
    lambda_eye = UniformScaling(lambda)
  else
    lambda_eye = lambda * eye(N)
  end
  xt = X'
  theta = xt * ((X * xt .+ lambda_eye) \ y)

  if nonnegative
    set_negative_to_zero!(theta)
  end

  iterations = 0

  while iterations < maxIterations
    iterations += 1

    # E-step:
    eta = theta

    # M-step:
    eta_squared = eta .^ 2
    xt_eta = broadcast(*, eta_squared, xt)
    theta = xt_eta * ((X * xt_eta .+ lambda_eye) \ y)

    # Ensure non-negativity if required
    if nonnegative
      set_negative_to_zero!(theta)
    end

    # We have converged if change is too small => break.
    if norm(theta - eta) < deltaTreshold
      break
    end
  end

  theta[abs(theta) .< epsilon] = 0.0

  theta
end

function set_negative_to_zero!(v)
  v[v .< 0.0] = 0.0
  v
end

function format_time(t)
  if t < 1e-11
    return @sprintf("%.2e nsec", t/1e-9)
  end
  sec_steps = [(1e-7, 1e-9, "nsec"), (1e-4, 1e-6, "usec"),
    (1e-1, 1e-3, "msec"), (60, 1, "sec"), (3600, 60, "mins")]
  for i in 1:length(sec_steps)
    limit, divisor, name = sec_steps[i]
    if t < limit
      return @sprintf("%.2f %s", t/divisor, name)
    end
  end
  @sprintf("%.2f hours", t/3600)
end

using DataFrames
using HypothesisTests

coefficient_of_variation(v) = std(v) / mean(v)
mad(v) = median(abs(v .- median(v)))
rcv(v) = mad(v) / median(v)
iqr(v) = quantile(v, 0.75) - quantile(v, 0.25)
riqr(v) = iqr(v) / median(v)

SkipIndicator = -100
WP = symbol("p-value")
PM = symbol("+/- %")

summarize(df::AbstractDataFrame) = DataFrame(
  Avg = mean(df[:Time]),
  PM = 100.0 * coefficient_of_variation(df[:Time]),
  #RCV = rcv(df[:Time]),
  #RIQR = riqr(df[:Time]),
  Max = maximum(df[:Time]),
  Min = minimum(df[:Time]),
  Reps = length(df[:Time])
)

#significance_indicator(p) = (p < 0.001) ? "***" : ((p < 0.01) ? "**" : ((p < 0.05) ? "*" : ((p == 1.0) ? "N/A" : "")))
significance_indicator(p) = (p < 0.001) ? "Yes (p < 0.001)" : ((p < 0.01) ? "Yes (p < 0.01)" : ((p < 0.05) ? "Yes (p < 0.05)" : ((p == 1.0) ? "N/A" : "No")))
speed_indication(t) = (t < 1) ? "." : ((t < 10) ? "," : ((t < 300) ? ":" : "|"))

calc_pvalue(fastest, times) = begin
  len = min(length(fastest), length(times))
  # We use a right sided p-value since we know times is at least as large as fastest, on average.
  pvalue(SignTest( (times .- fastest[1:len]), 0 ), tail = :right)
  # Wilcoxon/MannWhitneyUTest not really right to use since executions are paired...
  #pvalue(MannWhitneyUTest(times, fastest), tail = :both)
end

mean_speed_diff(fastest, times) = mean(times ./ fastest[1:length(times)])

summarize_measurements(times, indices, n, rep, descs, sortByAvg = false) = begin
  df = DataFrame({times, indices}, [:Time, :Function])

  # Skip rows that we have not yet filled in / measured:
  df = df[df[:Function] .> 0, :]

  # Skip rows that we have indicated as skipped
  df = df[!(df[:Time] .== SkipIndicator), :]

  # Summarize and calc relative speed
  dfs = by(df, :Function, summarize)
  dfs[:Relative] = dfs[:Avg] / minimum(dfs[:Avg])

  # Find fastest function
  fastest = dfs[:Function][findfirst(dfs[:Avg], minimum(dfs[:Avg]))]

  # Compare all functions to the fastest with a wilcoxon test
  times_fastest = df[:Time][df[:Function] .== fastest]
  #dfs[WP] = map((j) -> pvalue(MannWhitneyUTest(times_fastest, df[:Time][df[:Function] .== j])), 1:n)
  dfs[WP] = map((j) -> calc_pvalue(times_fastest, df[:Time][df[:Function] .== j]), 1:n)

  # Calc relative for each paired rep and take mean of that:
  #dfs[:Rel2] = map((j) -> mean_speed_diff(times_fastest, df[:Time][df[:Function] .== j]), 1:n)

  # Make information more readable
  dfs[:Function] = descs
  if sortByAvg
    dfs = sort(dfs, cols = :Avg)
  end
  dfs[:Avg] = map(format_time, dfs[:Avg])
  #dfs[:Median] = map(format_time, dfs[:Median])
  dfs[symbol("Rel.")] = map((v) -> signif(v, 3), dfs[:Relative])
  map( (s) -> dfs[s] = map((v) -> signif(v, 2), dfs[s]), [WP] )
  dfs[:Max] = map(format_time, dfs[:Max])
  dfs[:Min] = map(format_time, dfs[:Min])
  ss = symbol("Signif slower?")
  dfs[ss] = map(significance_indicator, dfs[WP])
  dfs[PM] = map((v) -> signif(v, 3), dfs[:PM])

  dfs[:,[:Function, :Avg, symbol("Rel."), :Reps, PM, :Max, :Min, WP, ss]]
end

spaces(num) = join( map((i) -> " ", 1:num), "" )

# Compare the execution speed of a set of functions on the same randomly generated input.
# Run only as many replications of each function as is needed to show that it is slower than
# the fastest of the functions.
function speed_compare(funcs::Dict{ASCIIString, Function}, inputsGen::Function;
  MaxReplications = 15,
  MinReplications = 3,
  alpha = 0.05)

  # Create description strings of equal length
  descs = collect(keys(funcs))
  max_desc_len = maximum(map(length, descs))
  formatted_descs = map((s) -> s * spaces(max_desc_len - length(s)), descs)

  fs = map((d) -> funcs[d], descs)
  n = length(fs)

  # Force JIT compilation
  inputs = inputsGen()
  map((i) -> fs[i](inputs...), 1:n)

  times = zeros(Float64, n * MaxReplications)
  indices = zeros(Int, n * MaxReplications)
  still_run = map((i) -> true, 1:n)

  local rep # since we need it in return statement below

  for rep in 1:MaxReplications
    inputs = inputsGen()
    for i in shuffle(collect(1:n))
      index = (i - 1) * MaxReplications + rep
      indices[index] = i
      if still_run[i]
        f = fs[i]
        times[index] = @elapsed f(inputs...)
        # println("  $(i). $(descs[i]): $(format_time(times[index]))")
      else
        times[index] = SkipIndicator
        # println("  $(i). $(descs[i]): skipped")
      end
      print(speed_indication(times[index]))
    end

    if rep >= MinReplications
      dfs = summarize_measurements(times, indices, n, rep, descs)
      # Only run those which do not have p-values below alpha compared
      # to the fastest running functions:
      still_run = still_run & (dfs[WP] .> alpha)
      if sum(still_run) < 2
        break
      end
    end
  end

  return summarize_measurements(times, indices, n, rep, formatted_descs, true)
end

alts = {
  "original" => l0EM_orig,
  "solve" => l0EM_solve,
  "unisc & solve" => l0EM_us_solve,
}
alts = convert(Dict{ASCIIString, Function}, alts)

input_gen_func(N, M, k, epsilon = 0.10) = () -> begin
  X = randn(N, M)
  #actual_theta = shuffle([linspace(1,k,k), zeros(M-k)][:])
  actual_theta = shuffle([10.0 * rand(k,1), zeros(M-k)][:])
  y = X * actual_theta + epsilon * randn(N, 1)
  (X, y)
end

alts2 = {
  "original" => l0EM_orig,
  "solve" => l0EM_solve,
  "unisc & solve" => l0EM_us_solve,
  # "ad. unisc & solve" => l0EM_adaptus_solve,
}
alts2 = convert(Dict{ASCIIString, Function}, alts2)

for N in [1e2, 1e3]
  for M in [1e2, 1e3, 1e4]
    for k in [4, 16]
      println("\nN = $(int(N)), M = $(int(M)), k = $(int(k))")
      showall( speed_compare(alts2, input_gen_func(int(N), int(M), k)) )
      println("")
    end
  end
end


text:
using RNLA, Test, LinearAlgebra

@testset "rpca_test1" begin
    m = 10
    k = 10
    X = rand(m, k)
    U, S, W1 = svd(X)
    Z1 = U * Diagonal(S)
    Λ1 = (S .^ 2) ./ (m - 1)

    W2, Λ2, Z2 = rpca(X, k)

    @test norm(abs.(W1) - abs.(W2)) < sqrt(eps())
    @test norm(abs.(Λ1) - abs.(Λ2)) < sqrt(eps())
    @test norm(abs.(Z1) - abs.(Z2)) < sqrt(eps())
end

@testset "rpca_test2" begin
    m = 100
    n = 1000
    k = 50
    X = rand(m, n)
    U, S, W1 = svd(X)
    W1 = W1[:, 1:k]
    Z1 = U[:, 1:k] * Diagonal(S[1:k])
    Λ1 = (S[1:k] .^ 2) ./ (m - 1)

    W2, Λ2, Z2 = rpca(X, k)

    @test norm(abs.(W1) - abs.(W2)) < sqrt(eps())
    @test norm(abs.(Λ1) - abs.(Λ2)) < sqrt(eps())
    @test norm(abs.(Z1) - abs.(Z2)) < sqrt(eps())
end

@testset "rpca_test3" begin
    m = 1000
    n = 100
    k = 50
    X = rand(m, n)
    U, S, W1 = svd(X)
    W1 = W1[:, 1:k]
    Z1 = U[:, 1:k] * Diagonal(S[1:k])
    Λ1 = (S[1:k] .^ 2) ./ (m - 1)

    W2, Λ2, Z2 = rpca(X, k)

    @test norm(abs.(W1) - abs.(W2)) < sqrt(eps())
    @test norm(abs.(Λ1) - abs.(Λ2)) < sqrt(eps())
    @test norm(abs.(Z1) - abs.(Z2)) < sqrt(eps())
end


text:
module PartedArrays
export
    PartedArray,
    PartedVector,
    PartedMatrix,
    create_partition,
    create_partition2

include("partitioning.jl")
include("partedarray.jl")

end # module


text:
using StaticArrays, NLSolvers, LinearAlgebra
OPT_PROBS = Dict()

#### Himmelblau

OPT_PROBS["himmelblau"] = Dict()
OPT_PROBS["himmelblau"]["array"] = Dict()

function himmelblau!(x)
    fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
    return fx
end
function himmelblau_batched_f!(X)
    F = map(himmelblau!, X)
    return F
end
function himmelblau_batched_f!(F, X)
    map!(himmelblau!, F, X)
    return F
end
function himmelblau_g!(∇f, x)
    ∇f[1] =
        4.0 * x[1]^3 + 4.0 * x[1] * x[2] - 44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
    ∇f[2] =
        2.0 * x[1]^2 + 2.0 * x[2] - 22.0 + 4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
    ∇f
end

function himmelblau_fg!(∇f, x)
    ∇f = himmelblau_g!(∇f, x)
    fx = himmelblau!(x)
    return fx, ∇f
end
function himmelblau_fgh!(∇f, ∇²f, x)
    ∇²f = himmelblau_h!(∇²f, x)
    fx, ∇f = himmelblau_fg!(∇f, x)
    return fx, ∇f, ∇²f
end
function himmelblau_h!(∇²f, x)
    ∇²f[1, 1] = 12.0 * x[1]^2 + 4.0 * x[2] - 44.0 + 2.0
    ∇²f[1, 2] = 4.0 * x[1] + 4.0 * x[2]
    ∇²f[2, 1] = ∇²f[1, 2]
    ∇²f[2, 2] = 2.0 + 4.0 * x[1] + 12.0 * x[2]^2 - 28.0
    return ∇²f
end
function himmelblau_hv!(hv, x, v)
    hv[1] =
        (12.0 * x[1]^2 + 4.0 * x[2] - 44.0 + 2.0) * v[1] + (4.0 * x[1] + 4.0 * x[2]) * v[2]
    hv[2] =
        (4.0 * x[1] + 4.0 * x[2]) * v[1] + (2.0 + 4.0 * x[1] + 12.0 * x[2]^2 - 28.0) * v[2]
    return hv
end

OPT_PROBS["himmelblau"]["array"]["x0"] = [3.0, 1.0]
OPT_PROBS["himmelblau"]["array"]["mutating"] = ScalarObjective(
    himmelblau!,
    himmelblau_g!,
    himmelblau_fg!,
    himmelblau_fgh!,
    himmelblau_h!,
    himmelblau_hv!,
    himmelblau_batched_f!,
    nothing,
)

OPT_PROBS["himmelblau"]["staticarray"] = Dict()

function himmelblau_f(x)
    fx = (x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2
    return fx
end

function himmelblau_fg(∇f, x)
    ∇f = himmelblau_g(∇f, x)
    fx = himmelblau_f(x)
    return fx, ∇f
end

function himmelblau_g(∇f, x)
    ∇f1 = 4.0 * x[1]^3 + 4.0 * x[1] * x[2] - 44.0 * x[1] + 2.0 * x[1] + 2.0 * x[2]^2 - 14.0
    ∇f2 = 2.0 * x[1]^2 + 2.0 * x[2] - 22.0 + 4.0 * x[1] * x[2] + 4.0 * x[2]^3 - 28.0 * x[2]
    ∇f = @SVector([∇f1, ∇f2])
    return ∇f
end

function himmelblau_fgh(∇f, ∇²f, x)
    ∇²f11 = 12.0 * x[1]^2 + 4.0 * x[2] - 42.0
    ∇²f12 = 4.0 * x[1] + 4.0 * x[2]
    ∇²f21 = 4.0 * x[1] + 4.0 * x[2]
    ∇²f22 = 12.0 * x[2]^2 + 4.0 * x[1] - 26.0
    ∇²f = @SMatrix([∇²f11 ∇²f12; ∇²f21 ∇²f22])
    fx, ∇f = himmelblau_fg(∇f, x)

    return fx, ∇f, ∇²f
end

sl = @SVector([3.0, 1.0])
state0newton = (@SVector([2.0, 2.0]), I + sl * sl')
OPT_PROBS["himmelblau"]["staticarray"]["x0"] = sl
OPT_PROBS["himmelblau"]["staticarray"]["state0"] = state0newton
OPT_PROBS["himmelblau"]["staticarray"]["static"] = ScalarObjective(
    himmelblau_f,
    himmelblau_g,
    himmelblau_fg,
    himmelblau_fgh,
    nothing,
    nothing,
    nothing,
    nothing,
)

### Exponential

exponential!(x) = exp((2.0 - x[1])^2) + exp((3.0 - x[2])^2)
function exponential_g!(g, x)
    g[1] = -2.0 * (2.0 - x[1]) * exp((2.0 - x[1])^2)
    g[2] = -2.0 * (3.0 - x[2]) * exp((3.0 - x[2])^2)
    return g
end
function exponential_h!(H, x)
    H[1, 1] = 2.0 * exp((2.0 - x[1])^2) * (2.0 * x[1]^2 - 8.0 * x[1] + 9)
    H[1, 2] = 0.0
    H[2, 1] = 0.0
    H[2, 2] = 2.0 * exp((3.0 - x[2])^2) * (2.0 * x[2]^2 - 12.0 * x[2] + 19)
    return H
end
function exponential_hv!(Hv, x, v)
    Hv[1, 1] = (2.0 * exp((2.0 - x[1])^2) * (2.0 * x[1]^2 - 8.0 * x[1] + 9)) * v[1]
    Hv[2, 2] = (2.0 * exp((3.0 - x[2])^2) * (2.0 * x[2]^2 - 12.0 * x[2] + 19)) * v[2]
    return Hv
end
function exponential_fg!(g, x)
    fx = exponential!(x)
    g = exponential_g!(g, x)
    return fx, g
end
function exponential_fgh!(g, H, x)
    fx, g = exponential_fg!(g, x)
    H = exponential_h!(H, x)
    return fx, g, H
end

OPT_PROBS["exponential"] = Dict()
OPT_PROBS["exponential"]["array"] = Dict()
# Byttet om på x og H
OPT_PROBS["exponential"]["array"]["x0"] = [0.0, 0.0]
OPT_PROBS["exponential"]["array"]["mutating"] = ScalarObjective(
    exponential!,
    exponential_g!,
    exponential_fg!,
    exponential_fgh!,
    exponential_h!,
    exponential_hv!,
    nothing,
    nothing,
)

OPT_PROBS["laplacian"] = Dict()
OPT_PROBS["laplacian"]["array"] = Dict()
# Byttet om på x og H
OPT_PROBS["laplacian"]["array"]["x0(n)"] = n -> zeros(n)

plap(U; n = length(U)) = (n - 1) * sum((0.1 .+ diff(U) .^ 2) .^ 2) - sum(U) / (n - 1)
plap1(U; n = length(U), dU = diff(U), dW = 4 .* (0.1 .+ dU .^ 2) .* dU) =
    (n - 1) .* ([0.0; dW] .- [dW; 0.0]) .- ones(n) / (n - 1)
precond(x::Vector) = precond(length(x))
precond(n::Number) =
    spdiagm(-1 => -ones(n - 1), 0 => 2 * ones(n), 1 => -ones(n - 1)) * (n + 1)
_f(x) = plap([0; x; 0])
function _fg(∇f, x)
    fx = _f(x)
    copyto!(∇f, (plap1([0; x; 0]))[2:end-1])
    fx, ∇f
end

OPT_PROBS["laplacian"]["array"]["mutating"] = OptimizationProblem(
    ScalarObjective(_f, nothing, _fg, nothing, nothing, nothing, nothing, nothing),
)

function theta(x)
    if x[1] > 0
        return atan(x[2] / x[1]) / (2.0 * pi)
    else
        return (pi + atan(x[2] / x[1])) / (2.0 * pi)
    end
end

function fletcher_powell_f(x)
    theta_x = theta(x)
    fx = 100.0 * ((x[3] - 10.0 * theta_x)^2 + (sqrt(x[1]^2 + x[2]^2) - 1.0)^2) + x[3]^2
    return fx
end
function fletcher_powell_g(∇f, x)
    T = eltype(x)
    theta_x = theta(x)

    if x[1]^2 + x[2]^2 == 0
        dtdx1 = T(0)
        dtdx2 = T(0)
    else
        dtdx1 = -x[2] / (T(2) * pi * (x[1]^2 + x[2]^2))
        dtdx2 = x[1] / (T(2) * pi * (x[1]^2 + x[2]^2))
    end
    ∇f1 =
        -2000.0 * (x[3] - 10.0 * theta_x) * dtdx1 +
        200.0 * (sqrt(x[1]^2 + x[2]^2) - 1) * x[1] / sqrt(x[1]^2 + x[2]^2)
    ∇f2 =
        -2000.0 * (x[3] - 10.0 * theta_x) * dtdx2 +
        200.0 * (sqrt(x[1]^2 + x[2]^2) - 1) * x[2] / sqrt(x[1]^2 + x[2]^2)
    ∇f3 = 200.0 * (x[3] - 10.0 * theta_x) + 2.0 * x[3]
    ∇f = @SVector[∇f1, ∇f2, ∇f3]

    return ∇f
end
function fletcher_powell_fg(∇f, x)
    T = eltype(x)
    theta_x = theta(x)

    if x[1]^2 + x[2]^2 == 0
        dtdx1 = T(0)
        dtdx2 = T(0)
    else
        dtdx1 = -x[2] / (T(2) * pi * (x[1]^2 + x[2]^2))
        dtdx2 = x[1] / (T(2) * pi * (x[1]^2 + x[2]^2))
    end
    ∇f1 =
        -2000.0 * (x[3] - 10.0 * theta_x) * dtdx1 +
        200.0 * (sqrt(x[1]^2 + x[2]^2) - 1) * x[1] / sqrt(x[1]^2 + x[2]^2)
    ∇f2 =
        -2000.0 * (x[3] - 10.0 * theta_x) * dtdx2 +
        200.0 * (sqrt(x[1]^2 + x[2]^2) - 1) * x[2] / sqrt(x[1]^2 + x[2]^2)
    ∇f3 = 200.0 * (x[3] - 10.0 * theta_x) + 2.0 * x[3]
    ∇f = @SVector[∇f1, ∇f2, ∇f3]

    fx = 100.0 * ((x[3] - 10.0 * theta_x)^2 + (sqrt(x[1]^2 + x[2]^2) - 1.0)^2) + x[3]^2

    return fx, ∇f
end
const sv3 = @SVector[0.1, 0.0, 0.0]

const fp_static = ScalarObjective(
    fletcher_powell_f,
    fletcher_powell_g,
    fletcher_powell_fg,
    nothing,
    nothing,
    nothing,
    nothing,
    nothing,
)

const state0 = (@SVector[-0.5, 0.0, 0.0], I + sv3 * sv3')
OPT_PROBS["fletcher_powell"] = Dict()
OPT_PROBS["fletcher_powell"]["staticarray"] = Dict()
# Byttet om på x og H
OPT_PROBS["fletcher_powell"]["staticarray"]["x0"] = state0
OPT_PROBS["fletcher_powell"]["staticarray"]["static"] =
    OptimizationProblem(fp_static; inplace = false)

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
using DataFrames
using RDatasets

function col_normalize(X::Array{Float64})
    xmin = minimum(X, dims=1)
    xmax = maximum(X, dims=1)
    for i in 1:size(X, 2)
        X[:, i] = (X[:, i] .- xmin[i]) ./ (xmax[i] - xmin[i])
    end
    Array{Float64}(X')
end

function make_x_y(data::DataFrame, seed::Int64)
    x_data = convert(Array{Float64}, data[1:(size(data, 2) - 1)])
    X = col_normalize(x_data)

    labels = unique(data[end])
    Y = convert(Array{Int64}, indexin(data[end], labels))

    rng = MersenneTwister(seed)
    inds = randperm(rng, length(Y))
    X[:, inds], Y[inds]
end

function get_iris(seed::Int64)
    make_x_y(dataset("datasets", "iris"), seed)
end

function get_diabetes(seed::Int64)
    make_x_y(vcat(dataset("MASS", "Pima.te"), dataset("MASS", "Pima.tr")),
             seed)
end

function get_glass(seed::Int64)
    make_x_y(dataset("MASS", "fgl"), seed)
end


text:
using LinearAlgebra
using Random
using Statistics
using ProgressMeter

PauliMatrices = [[0 1; 1 0], [0 -1im; 1im 0], [1 0; 0 -1]]
Idmatrix = [1 0; 0 1] .+ 0.0im

# kron(PauliMatrices[1], PauliMatrices[1])
⊗(x, y) = kron(x, y)

############################### hamiltonian ###########################

function Ham_Jz_ops(Jz_arr::Array{Float64, 1}, N::Int64)
    term_ops_sum = zeros(ComplexF64, 2^N, 2^N)

    index = 1
    for i in 1:N
        a = 3
        term_ops = fill(Idmatrix, N)
        term_ops[i] = PauliMatrices[a]
        term_ops[mod(i, N) + 1] = PauliMatrices[a]
        term_ops_sum[:,:] += Jz_arr[i] * foldl(⊗, term_ops)

        index += 1
    end

    return term_ops_sum
end

function Ham_hz_ops(hz_arr::Array{Float64, 1}, N::Int64)
    term_ops_sum = zeros(ComplexF64, 2^N, 2^N)

    index = 1
    for i in 1:N
        a = 3
        term_ops = fill(Idmatrix, N)
        term_ops[i] = PauliMatrices[a]
        term_ops_sum[:,:] += hz_arr[i] * foldl(⊗, term_ops)

        index += 1
    end

    return term_ops_sum
end

function Ham_Jx_ops(Jx_arr::Array{Float64, 1}, N::Int64)
    term_ops_sum = zeros(ComplexF64, 2^N, 2^N)

    index = 1
    for i in 1:N
        a = 1
        term_ops = fill(Idmatrix, N)
        term_ops[i] = PauliMatrices[a]
        term_ops_sum[:,:] += Jx_arr[i] * foldl(⊗, term_ops)

        index += 1
    end

    return term_ops_sum
end


text:
struct DDPSimulation{NS,NC,T}
	prob::DDP
	sol::AbstractDDPSolution{NS,NC}
	value::Union{Array{T,3}, Nothing}
	states::Array{T,3}
	policy::Array{T,3}
end
value(sim::DDPSimulation) = sim.value
states(sim::DDPSimulation) = sim.states
policy(sim::DDPSimulation) = sim.policy

"""Construct a DataFrame from the simulation."""
function DataFrame(sim::DDPSimulation)
	(nStateVars, nPeriods, nFirms) = size(states(sim))
	nChoiceVars = size(policy(sim), 1)

	mFirm = repeat(collect(1:nFirms)', nPeriods, 1)
	mPeriod = repeat(collect(0:nPeriods-1), 1, nFirms)

	df = DataFrame()
	df.firm = mFirm[:]
	df.period = mPeriod[:]

	[df[Symbol("state_" * string(i))] = states(sim)[i,:,:][:] for i = 1 : nStateVars]
	[df[Symbol("choice_" * string(i))] = policy(sim)[i,:,:][:] for i = 1 : nChoiceVars]
	if value(sim) != nothing
		df[:value] = value(sim)[1,:,:][:]
	end
	return df
end

"""Construct an `variable x periods x firms` array from the simulation.
The order of the variables is states, choices and value."""
function Array(sim::DDPSimulation)
	if value(sim) == nothing
		return cat(states(sim), policy(sim), dims=1)
	else
		return cat(states(sim), policy(sim), value(sim), dims=1)
	end
end


text:
module day12

include("inputs.jl")

ex = inputs.ex_d12()
input = inputs.input_d12()

mutable struct Ship
    pos::Array{Int}
    heading::Array{Int}
    waypoint::Array{Int}
end

parse_instruction(line) = line[1], parse(Int, line[2:end])
rot(θ) = [cosd(θ) -sind(θ); sind(θ) cosd(θ)]

function heading!(ship, act, val)
    if     act == 'N' ship.pos += [0;1] * val
    elseif act == 'E' ship.pos += [1;0] * val
    elseif act == 'S' ship.pos += [0;-1] * val
    elseif act == 'W' ship.pos += [-1;0] * val
    elseif act == 'R' ship.heading = rot(-val) * ship.heading
    elseif act == 'L' ship.heading = rot(val) * ship.heading
    elseif act == 'F' ship.pos += ship.heading * val
    end
end

function waypoint!(ship, act, val)
    if     act == 'N' ship.waypoint += [0;1] * val
    elseif act == 'E' ship.waypoint += [1;0] * val
    elseif act == 'S' ship.waypoint += [0;-1] * val
    elseif act == 'W' ship.waypoint += [-1;0] * val
    elseif act == 'R' ship.waypoint = rot(-val) * ship.waypoint
    elseif act == 'L' ship.waypoint = rot(val) * ship.waypoint
    elseif act == 'F' ship.pos += ship.waypoint * val
    end
end

function navigate(input, move)
    insts = parse_instruction.(input)
    ship = Ship([0;0], [1; 0], [10,1])
    [move(ship, inst...) for inst in insts]
    return sum(abs.(ship.pos))
end

part1(input) = navigate(input, heading!)
part2(input) = navigate(input, waypoint!)

@assert part1(ex) == 25
@assert part1(input) == 415
@info "Day 12, Part 1 answer: $(part1(input))"
@assert part2(ex) == 286
@assert part2(input) == 29401
@info "Day 12, Part 2 answer: $(part2(input))"

end


text:
module JuAFEM
using Reexport
@reexport using Tensors
@reexport using WriteVTK

using LinearAlgebra
using SparseArrays
using Base: @propagate_inbounds

include("exports.jl")

"""
Represents a reference shape which quadrature rules and interpolations are defined on.
Currently, the only concrete types that subtype this type are `RefCube` in 1, 2 and 3 dimensions,
and `RefTetrahedron` in 2 and 3 dimensions.
"""
abstract type AbstractRefShape end

struct RefTetrahedron <: AbstractRefShape end
struct RefCube <: AbstractRefShape end

"""
Abstract type which has `CellValues` and `FaceValues` as subtypes
"""
abstract type Values{dim,T,refshape} end
abstract type CellValues{dim,T,refshape} <: Values{dim,T,refshape} end
abstract type FaceValues{dim,T,refshape} <: Values{dim,T,refshape} end

include("utils.jl")

# Interpolations
include("interpolations.jl")

# Quadrature
include("Quadrature/quadrature.jl")

# FEValues
include("FEValues/cell_values.jl")
include("FEValues/face_values.jl")
include("FEValues/common_values.jl")
include("FEValues/face_integrals.jl")

# Grid
include(joinpath("Grid", "grid.jl"))
include(joinpath("Grid", "grid_generators.jl"))
include(joinpath("Grid", "coloring.jl"))

# Dofs
include(joinpath("Dofs", "DofHandler.jl"))
include(joinpath("Dofs", "MixedDofHandler.jl"))
include(joinpath("Dofs", "ConstraintHandler.jl"))

include("iterators.jl")

# Assembly
include("assembler.jl")

# Projection
include("L2_projection.jl")

# Export
include("Export/VTK.jl")

# Point Evaluation
include(joinpath("PointEval", "PointEvalHandler.jl"))

# Other
include("deprecations.jl")

end # module


text:
module TrackError

using ..Util, ..Tagging, ..Conditioning, ..Traits
import ..Conditioning
export applytrackerr, condvar, condvarapply, tagerror, initerror

# Options.
#> condition on logerr, explicitly

# Options
# 1. Always initialize with correct type -- best performance
# 2. Have nothing, small union ⊥ ∪ Float64
# 3. Vector

"""
`applytrackerr(x, ω, initerr)`

Is `ω` in the domain of `x`?
Returns `(x(ω), v)` where `v = ω ∈ dom(x)`
"""
function applytrackerr(x, ω, initerr)
  ω_ = tagerror(ω, initerr) # zt: surprisingly expensive
  fx = x(ω_)
  (fx = fx, err = ω_.tags.err.val)
end

"Equivalent to `condvar(x)(ω)`, but more efficient"
condvarapply(x, ω, initerr) = applytrackerr(x, ω, initerr).err
condvarapply(x, ω, errtype::Type{T}) where T =
  condvarapply(x, ω, initerror(T))

logcondvarapply(x, ω) = condvarapply(x, ω)

function initerror end
initerror(::Type{Bool}) = true
"""
`condvar(x)`

Random variable that `x` is conditioned on.
I.e. if `x` can be expressed as `x | y`, then `y = condvar(x)`
"""
condvar(x, initerr) = ω -> condvarapply(x, ω, initerr)

"Update `err` by conjoining current error with `b`"
conjoinerror!(err::Box, b) = err.val &= b

# Track error for this ω
dotrack(ω) = !haskey(ω.tags, :donttrack) && haskey(ω.tags, :err)

function Conditioning.condf(::trait(Err), ω, x, y)
  dotrack(ω) && conjoinerror!(ω.tags.err, y(ω))
  x(ω)
end

function Conditioning.cond!(::trait(Err), ω, bool)
  dotrack(ω) && conjoinerror!(ω.tags.err, bool)
  bool
end

# conjoinerror!(ω.tags.err, y(ω))
tagerror(::trait(Err), ω, initerr) = tag(ω,  (err = Box(initerr),))
tagerror(traits, ω, initerr) = tag(ω, (err = Box(initerr),))

"Tag `ω` with `err`"
tagerror(ω::Ω, initerr) where Ω = tagerror(traits(Ω), ω, initerr)

function combinetags(::Type{Val{:err}}, a, b)
  # @show a
  # @show b
  # @assert false
  @show a[]
  @show a[] & b[]
end

end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
n = 2
μ0 = zeros(n)
P0 = Diagonal(ones(n))

x0 = resample(μ0, P0, 1.0)

β = 10.0
μ1 = sample_mean(x0)
P1 = sample_covariance(x0, β)

A = Diagonal(ones(n))
B = zeros(n)

x1 = []
u0 = rand(1)
for j = 1:2n
	push!(x1, A * x0[j] + B * u0[1])
end

A' * P0 * A
β = 10.0
μ1 = sample_mean(x1)
P1 = sample_covariance(x1, β)


text:
@testset "Geo interface" begin
    pt = LibGEOS.Point(1.0, 2.0)
    @test GeoInterface.coordinates(pt) ≈ [1, 2] atol = 1e-5
    @test GeoInterface.geotype(pt) == :Point

    pt = LibGEOS.Point(1, 2)
    @test GeoInterface.coordinates(pt) ≈ [1, 2] atol = 1e-5
    @test GeoInterface.geotype(pt) == :Point

    pt = LibGEOS.readgeom("POINT EMPTY")
    @test GeoInterface.coordinates(pt) ≈ Float64[] atol = 1e-5
    @test GeoInterface.geotype(pt) == :Point

    mpt = LibGEOS.readgeom("MULTIPOINT(0 0, 10 0, 10 10, 11 10)")
    @test GeoInterface.coordinates(mpt) ==
          Vector{Float64}[[0, 0], [10, 0], [10, 10], [11, 10]]
    @test GeoInterface.geotype(mpt) == :MultiPoint

    coords = Vector{Float64}[[8, 1], [9, 1], [9, 2], [8, 2]]
    ls = LibGEOS.LineString(coords)
    @test GeoInterface.coordinates(ls) == coords
    @test GeoInterface.geotype(ls) == :LineString

    ls = LibGEOS.readgeom("LINESTRING EMPTY")
    @test GeoInterface.coordinates(ls) == []
    @test GeoInterface.geotype(ls) == :LineString

    mls = LibGEOS.readgeom("MULTILINESTRING ((5 0, 10 0), (0 0, 5 0))")
    @test GeoInterface.coordinates(mls) == [[[5, 0], [10, 0]], [[0, 0], [5, 0]]]
    @test GeoInterface.geotype(mls) == :MultiLineString

    coords = Vector{Float64}[[8, 1], [9, 1], [9, 2], [8, 2], [8, 1]]
    lr = LibGEOS.LinearRing(coords)
    @test GeoInterface.coordinates(lr) == coords
    @test GeoInterface.geotype(lr) == :LineString

    coords = Vector{Vector{Float64}}[
        Vector{Float64}[[0, 0], [10, 0], [10, 10], [0, 10], [0, 0]],
        Vector{Float64}[[1, 8], [2, 8], [2, 9], [1, 9], [1, 8]],
        Vector{Float64}[[8, 1], [9, 1], [9, 2], [8, 2], [8, 1]],
    ]
    polygon = LibGEOS.Polygon(coords)
    @test GeoInterface.coordinates(polygon) == coords
    @test GeoInterface.geotype(polygon) == :Polygon

    polygon = LibGEOS.readgeom("POLYGON EMPTY")
    @test GeoInterface.coordinates(polygon) == [[]]
    @test GeoInterface.geotype(polygon) == :Polygon

    multipolygon = LibGEOS.readgeom("MULTIPOLYGON(((0 0,0 10,10 10,10 0,0 0)))")
    @test GeoInterface.coordinates(multipolygon) ==
          Vector{Vector{Vector{Float64}}}[Vector{Vector{Float64}}[Vector{Float64}[
        [0, 0],
        [0, 10],
        [10, 10],
        [10, 0],
        [0, 0],
    ]]]
    @test GeoInterface.geotype(multipolygon) == :MultiPolygon

    geomcollection = LibGEOS.readgeom(
        "GEOMETRYCOLLECTION (POLYGON ((8 2, 10 10, 8.5 1, 8 2)), POLYGON ((7 8, 10 10, 8 2, 7 8)), POLYGON ((3 8, 10 10, 7 8, 3 8)), POLYGON ((2 2, 8 2, 8.5 1, 2 2)), POLYGON ((2 2, 7 8, 8 2, 2 2)), POLYGON ((2 2, 3 8, 7 8, 2 2)), POLYGON ((0.5 9, 10 10, 3 8, 0.5 9)), POLYGON ((0.5 9, 3 8, 2 2, 0.5 9)), POLYGON ((0 0, 2 2, 8.5 1, 0 0)), POLYGON ((0 0, 0.5 9, 2 2, 0 0)))",
    )
    collection = GeoInterface.geometries(geomcollection)
    coords = Vector{Vector{Vector{Float64}}}[
        Vector{Vector{Float64}}[Vector{Float64}[
            [8.0, 2.0],
            [10.0, 10.0],
            [8.5, 1.0],
            [8.0, 2.0],
        ]],
        Vector{Vector{Float64}}[Vector{Float64}[
            [7.0, 8.0],
            [10.0, 10.0],
            [8.0, 2.0],
            [7.0, 8.0],
        ]],
        Vector{Vector{Float64}}[Vector{Float64}[
            [3.0, 8.0],
            [10.0, 10.0],
            [7.0, 8.0],
            [3.0, 8.0],
        ]],
        Vector{Vector{Float64}}[Vector{Float64}[
            [2.0, 2.0],
            [8.0, 2.0],
            [8.5, 1.0],
            [2.0, 2.0],
        ]],
        Vector{Vector{Float64}}[Vector{Float64}[
            [2.0, 2.0],
            [7.0, 8.0],
            [8.0, 2.0],
            [2.0, 2.0],
        ]],
        Vector{Vector{Float64}}[Vector{Float64}[
            [2.0, 2.0],
            [3.0, 8.0],
            [7.0, 8.0],
            [2.0, 2.0],
        ]],
        Vector{Vector{Float64}}[Vector{Float64}[
            [0.5, 9.0],
            [10.0, 10.0],
            [3.0, 8.0],
            [0.5, 9.0],
        ]],
        Vector{Vector{Float64}}[Vector{Float64}[
            [0.5, 9.0],
            [3.0, 8.0],
            [2.0, 2.0],
            [0.5, 9.0],
        ]],
        Vector{Vector{Float64}}[Vector{Float64}[
            [0.0, 0.0],
            [2.0, 2.0],
            [8.5, 1.0],
            [0.0, 0.0],
        ]],
        Vector{Vector{Float64}}[Vector{Float64}[
            [0.0, 0.0],
            [0.5, 9.0],
            [2.0, 2.0],
            [0.0, 0.0],
        ]],
    ]
    @test map(GeoInterface.coordinates, collection) == coords
    @test GeoInterface.geotype(geomcollection) == :GeometryCollection

    geomcollection = LibGEOS.readgeom(
        "GEOMETRYCOLLECTION(MULTIPOINT(0 0, 0 0, 1 1),LINESTRING(1 1, 2 2, 2 2, 0 0),POLYGON((5 5, 0 0, 0 2, 2 2, 5 5)))",
    )
    collection = GeoInterface.geometries(geomcollection)
    coords = Vector[
        Vector{Float64}[[0.0, 0.0], [0.0, 0.0], [1.0, 1.0]],
        Vector{Float64}[[1.0, 1.0], [2.0, 2.0], [2.0, 2.0], [0.0, 0.0]],
        Vector{Vector{Float64}}[Vector{Float64}[
            [5.0, 5.0],
            [0.0, 0.0],
            [0.0, 2.0],
            [2.0, 2.0],
            [5.0, 5.0],
        ]],
    ]
    geotypes = [:MultiPoint, :LineString, :Polygon]
    for (i, item) in enumerate(collection)
        @test GeoInterface.coordinates(item) == coords[i]
        @test GeoInterface.geotype(item) == geotypes[i]
    end
    @test GeoInterface.geotype(geomcollection) == :GeometryCollection

    geomcollection = LibGEOS.readgeom("GEOMETRYCOLLECTION EMPTY")
    collection = GeoInterface.geometries(geomcollection)
    @test length(collection) == 0
    @test GeoInterface.geotype(geomcollection) == :GeometryCollection
end


text:
using Test
using DataStructures

using problems.linkedlist

@testset "dedupe 1 (immutable)" begin

    @test dedupe(list()) == list()

    @testset "case 1" begin
        actual = collect(dedupe(list(4, 6, 3, 1, 3, 2, 6, 1)))
        expected = [4, 6, 3, 1, 2]

        @test actual == expected
    end

    @testset "case 2" begin
        actual = collect(dedupe(list(4, 6, 3, 1, 3, 2, 6, 1, 7, 1, 2)))
        expected = [4, 6, 3, 1, 2, 7]

        @test actual == expected
    end

    @testset "case 3" begin
        actual = collect(dedupe(list(4, 6, 3, 1, 3, 2, 6, 1, 7, 1, 2)))
        expected = [4, 6, 3, 1, 2, 7]

        @test actual == expected
    end

    @testset "case 4" begin
        actual = collect(dedupe(list(1, 1, 1, 2, 3, 1, 1, 2, 3, 1, 4, 4)))
        expected = [1, 2, 3, 4]

        @test actual == expected
    end

    @testset "case 5" begin
        actual = collect(dedupe(list(1, 1, 'a', 2, 3, 1, 1, 2, 3, 'b', 4, 4)))
        expected = [1, 'a', 2, 3, 'b', 4]

        @test actual == expected
    end

end

@testset "dedupe 2 (mutable)" begin

    @test dedupe(MutableLinkedList()) == MutableLinkedList()

    @testset "case 1" begin

        lst = to_mutable_list([4, 6, 3, 1, 3, 2, 6, 1])

        actual = collect(dedupe(lst))
        expected = [4, 6, 3, 1, 2]

        @test actual == expected
    end

    @testset "case 2" begin

        lst = to_mutable_list([4, 6, 3, 1, 3, 2, 6, 1, 7, 1, 2])

        actual = collect(dedupe(lst))
        expected = [4, 6, 3, 1, 2, 7]

        @test actual == expected
    end

    @testset "case 3" begin

        lst = to_mutable_list([4, 6, 3, 1, 3, 2, 6, 1, 7, 1, 2])

        actual = collect(dedupe(lst))
        expected = [4, 6, 3, 1, 2, 7]

        @test actual == expected
    end

    @testset "case 4" begin
        lst = to_mutable_list([1, 1, 1, 2, 3, 1, 1, 2, 3, 1, 4, 4])

        actual = collect(dedupe(lst))
        expected = [1, 2, 3, 4]

        @test actual == expected
    end

    @testset "case 5" begin
        lst = to_mutable_list([1, 1, 'a', 2, 3, 1, 1, 2, 3, 'b', 4, 4])

        actual = collect(dedupe(lst))
        expected = [1, 'a', 2, 3, 'b', 4]

        @test actual == expected
    end
end

@testset "k-th last" begin

    @testset "case 1" begin

        lst = list(4, 6, 3, 1, 3, 2, 6, 1)

        actual = kth_last(lst, 0)
        expected = 1

        @test actual == expected
    end

    @testset "case 2" begin

        lst = list(4, 6, 3, 1, 3, 2, 6, 1)

        actual = kth_last(lst, 1)
        expected = 6

        @test actual == expected
    end

    @testset "case 3" begin

        lst = list(4, 6, 3, 1, 3, 2, 6, 1)

        actual = kth_last(lst, 2)
        expected = 2

        @test actual == expected
    end

    @testset "case 4" begin

        lst = list(4, 6, 3, 1, 3, 2, 6, 1)

        actual = kth_last(lst, 7)
        expected = 4

        @test actual == expected
    end

    @testset "case 5" begin

        lst = list(4, 6, 3, 1, 3, 2, 6, 1)

        actual = kth_last(lst, 6)
        expected = 6

        @test actual == expected
    end
end

@testset "delete middle" begin

    @testset "case 1" begin

        lst = list(4, 6, 3, 1, 3, 2, 6, 1)

        actual = delete_middle(lst, 4)
        expected = list(4, 6, 3, 3, 2, 6, 1)

        @test actual == expected
    end

    @testset "case 2" begin

        lst = list(4, 6, 2, 3, 1, 3, 2, 6, 1)

        actual = delete_middle(lst, 3)
        expected = list(4, 6, 3, 1, 3, 2, 6, 1)

        @test actual == expected
    end

    @testset "case 3" begin

        lst = list(4, 6, 2, 3, 1, 3, 2, 6, 1, 7)

        actual = delete_middle(lst, 6)
        expected = list(4, 6, 2, 3, 1, 2, 6, 1, 7)

        @test actual == expected
    end

    @testset "case 4" begin

        lst = list(4, 6, 2, 3, 1, 3, 2, 6, 1, 7, 8)

        actual = delete_middle(lst, 2)
        expected = list(4, 2, 3, 1, 3, 2, 6, 1, 7, 8)

        @test actual == expected
    end

    @testset "case 5" begin

        lst = list(4, 6, 2, 3, 1, 3, 2, 6, 1, 7, 8, 10)

        actual = delete_middle(lst, 10)
        expected = list(4, 6, 2, 3, 1, 3, 2, 6, 1, 8, 10)

        @test actual == expected
    end
end

@testset "partition" begin

    @testset "case 1" begin

        lst = list(4, 6, 3, 1, 3, 2, 6, 1)

        actual = partition(lst, 4)
        expected = list(1, 2, 3, 1, 3, 6, 6, 4)

        @test actual == expected
    end

    @testset "case 2" begin

        lst = list(4, 6, 2, 3, 1, 3, 2, 6, 1)

        actual = partition(lst, 3)
        expected = list(1, 2, 1, 2, 6, 3, 3, 6, 4)

        @test actual == expected
    end

    @testset "case 3" begin

        lst = list(4, 6, 2, 3, 1, 3, 2, 6, 1, 7)

        actual = partition(lst, 6)
        expected = list(1, 2, 3, 1, 3, 2, 4, 7, 6, 6)

        @test actual == expected
    end

    @testset "case 4" begin

        lst = list(4, 6, 2, 3, 1, 3, 2, 6, 1, 7, 8)

        actual = partition(lst, 2)
        expected = list(1, 1, 8, 7, 6, 2, 3, 3, 2, 6, 4)

        @test actual == expected
    end

    @testset "case 5" begin

        lst = list(4, 6, 2, 3, 1, 3, 2, 6, 1, 7, 8, 10)

        actual = partition(lst, 10)
        expected = list(8, 7, 1, 6, 2, 3, 1, 3, 2, 6, 4, 10)

        @test actual == expected
    end
end

max_stars_repo_name: bmoretz/Daily-Coding-Problem

max_issues_repo_name: bmoretz/Daily-Coding-Problem

max_forks_repo_name: bmoretz/Daily-Coding-Problem


text:
function work(job)
    endtime = time() + 0.1
    while time() < endtime # Blocking wait
        nothing
    end
    return (job, Threads.threadid())
end

mutable struct TestObj
    data :: Int
end


text:
function nessi_read_les(h5g::HDF5Group)
  nt = read(h5g, "nt")[]
  norb = read(h5g, "size1")[]

  les_ = read(h5g, "les")
  les = zeros(ComplexF64, norb, norb, nt+1, nt+1)

  indices = findall(triu(trues(nt+1,nt+1)))
  les[:,:,indices] .= les_

  for i in 1:norb
    for j in 1:norb
      les[i,j,:,:] .-= triu(les[i,j,:,:],1)'
    end
  end

  return les
end

function nessi_read_ret(h5g::HDF5Group)
  nt = read(h5g, "nt")[]
  norb = read(h5g, "size1")[]

  ret_ = read(h5g, "ret")
  ret = zeros(ComplexF64, norb, norb, nt+1, nt+1)

  indices = findall(triu(trues(nt+1,nt+1)))
  ret[:,:,indices] .= ret_

  for i in 1:norb
    for j in 1:norb
      ret[i,j,:,:] .-= triu(ret[i,j,:,:],1)'
    end
  end

  ret = permutedims(ret, (1,2,4,3))
  return ret
end

function nessi_read_tv(h5g::HDF5Group)
  nt = read(h5g, "nt")[]
  ntau = read(h5g, "ntau")[]
  norb = read(h5g, "size1")[]

  tv_ = read(h5g, "tv")
  tv = permutedims(reshape(tv_, (norb, norb, ntau+1, nt+1)),(1,2,4,3))
  return tv
end

function nessi_read_mat(h5g::HDF5Group)
  mat = read(h5g, "mat")
  return mat
end

function nessi_read_gf(h5g::HDF5Group)
  les = nessi_read_les(h5g)
  ret = nessi_read_ret(h5g)
  tv  = nessi_read_tv(h5g)
  mat = nessi_read_mat(h5g)
  return (;les, ret, tv, mat)
end

function TimeGF(les::AbstractMatrix, ret::AbstractMatrix, tv::AbstractMatrix, mat::AbstractVector, grid::TimeGrid)
  @assert grid.contour.domain == full_contour

  t = realtimes(grid)
  tau = imagtimes(grid)

  nt = length(t)
  ntau = length(tau)

  @assert size(les) == (nt,nt)
  @assert size(ret) == (nt,nt)
  @assert size(tv) == (nt,ntau)
  @assert size(mat) == (ntau,)

  G = TimeGF(grid) do t1, t2
    greater = heaviside(t1.val, t2.val)
    if ((t1.val.domain == forward_branch || t1.val.domain == backward_branch) &&
        (t2.val.domain == forward_branch || t2.val.domain == backward_branch))
      i = findfirst(isapprox(real(t1.val.val)), t)
      j = findfirst(isapprox(real(t2.val.val)), t)
      greater ? ret[i,j] + les[i,j] : les[i,j]
    elseif (t1.val.domain == imaginary_branch && (t2.val.domain == forward_branch || t2.val.domain == backward_branch))
      i = findfirst(isapprox(imag(t1.val.val)), -tau)
      j = findfirst(isapprox(real(t2.val.val)), t)
      conj(tv[j, ntau+1-i]) # akoi 19c
    elseif ((t1.val.domain == forward_branch || t1.val.domain == backward_branch) && t2.val.domain == imaginary_branch)
      i = findfirst(isapprox(real(t1.val.val)), t)
      j = findfirst(isapprox(imag(t2.val.val)), -tau)
      tv[i,j]
    else
      i = findfirst(isapprox(imag(t1.val.val)), -tau)
      j = findfirst(isapprox(imag(t2.val.val)), -tau)
      greater ? 1.0im * mat[i - j + 1] : -1.0im * mat[i - j + ntau]
    end
  end

  return G
end


text:
immutable InfSup{T<:Real}
    inf::T
    sup::T
end

# deliver all powers x^i for i = 1..n
function powers{T}(x::InfSup{T}, n::Integer)

	res = Array{InfSup{T}}(n)
    powers!(res, x, n)
end

function powers!{T}(res::Array{InfSup{T}}, x::InfSup{T}, n::Integer)
    if n <= 1
        res[n] = x
    else
        bc = count_ones(-one(n)) - leading_zeros(n-1)
        tp = one(n) << (bc - 1)
        powers!(res, x, tp)
        y = res[tp]
        for i = (tp+1)..n
            res[i] = y * res[i-tp]
        end
    end
    res
end

function *{T}(x::InfSup{T}, y::InfSup{T})


text:
function preprocess(PData::PDataST, H, g, n1, n2)
    return PData
end


text:
using Test
import sbp
using sbp.Collocated: boundary_matrix, build_operators_2d, build_covariant_basis

st=sbp.Strand1994

order=2
nx = 10
ny = 20
desc_D, desc_H = st.load_description(order=order)
x, h, Dx, Hx = st.build_operators(desc_D, desc_H, nx)
y, h, Dy, Hy = st.build_operators(desc_D, desc_H, ny)

@testset "Boundary matrix"  begin
        Bp = boundary_matrix(nx)
        @test size(Bp, 1) == nx
        @test size(Bp, 2) == nx
        @test isapprox(sum(Bp), 0)

        Bp = boundary_matrix(nx, false)
        @test size(Bp, 1) == nx
        @test size(Bp, 2) == nx
        @test isapprox(sum(Bp), 0)
end

@testset "2D operators" begin
        Bx = boundary_matrix(nx)
        By = boundary_matrix(ny)
        op = build_operators_2d(Dx, Dy, Hx, Hy, Bx, By)
        @test op.nx == nx
        @test op.ny == ny
        @test size(op.Dx, 1) == nx * ny
        @test size(op.Dx, 2) == nx * ny
        @test size(op.Dy, 1) == nx * ny
        @test size(op.Dy, 2) == nx * ny
end

@testset "Covariant basis" begin
        Bx = boundary_matrix(nx)
        By = boundary_matrix(ny)
        op = build_operators_2d(Dx, Dy, Hx, Hy, Bx, By)
        fx = ones(op.nx*op.ny)
        fy = ones(op.nx*op.ny)
        a = build_covariant_basis(fx, fy, op)

        x_r1 = zeros(op.nx * op.ny)
        x_r2 = zeros(op.nx * op.ny)
        y_r1 = zeros(op.nx * op.ny)
        y_r2 = zeros(op.nx * op.ny)

        @test size(a.x_r1) == size(a.x_r2)
        @test size(a.x_r1) == size(a.y_r1)
        @test size(a.x_r1) == size(a.y_r2)

        @test isapprox(a.x_r1, x_r1, atol=1.0)
        @test isapprox(a.x_r2, x_r2, atol=1.0)
        @test isapprox(a.y_r1, y_r1, atol=1.0)
        @test isapprox(a.y_r2, y_r2, atol=1.0)
end

max_stars_repo_name: ooreilly/sbpjl

max_issues_repo_name: ooreilly/sbpjl

max_forks_repo_name: ooreilly/sbpjl


text:
module MaxHelpingHandCustomSummitCheckpoint

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/CustomSummitCheckpoint" CustomSummitCheckpoint(x::Integer, y::Integer, firstDigit::String="zero", secondDigit::String="zero")

const placements = Ahorn.PlacementDict(
    "Custom Summit Checkpoint (max480's Helping Hand)" => Ahorn.EntityPlacement(
        CustomSummitCheckpoint
    )
)

const numberlist = String["zero", "one", "two", "three", "four", "five", "six", "seven", "eight", "nine", "minus", "questionmark"]

Ahorn.editingOptions(entity::CustomSummitCheckpoint) = Dict{String, Any}(
    "firstDigit" => numberlist,
    "secondDigit" => numberlist
)

baseSprite = "scenery/summitcheckpoints/base02.png"

function Ahorn.selection(entity::CustomSummitCheckpoint)
    x, y = Ahorn.position(entity)

    return Ahorn.getSpriteRectangle(baseSprite, x, y)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CustomSummitCheckpoint, room::Maple.Room)
    digit1 = get(entity.data, "firstDigit", 0)
    digit2 = get(entity.data, "secondDigit", 0)

    Ahorn.drawSprite(ctx, baseSprite, 0, 0)
    Ahorn.drawSprite(ctx, "MaxHelpingHand/summitcheckpoints/$digit1/numberbg.png", -2, 4)
    Ahorn.drawSprite(ctx, "MaxHelpingHand/summitcheckpoints/$digit1/number.png", -2, 4)
    Ahorn.drawSprite(ctx, "MaxHelpingHand/summitcheckpoints/$digit2/numberbg.png", 2, 4)
    Ahorn.drawSprite(ctx, "MaxHelpingHand/summitcheckpoints/$digit2/number.png", 2, 4)
end

end

max_stars_repo_name: swoolcock/MaxHelpingHand

max_issues_repo_name: swoolcock/MaxHelpingHand

max_forks_repo_name: swoolcock/MaxHelpingHand


text:
#+Documented_Item,Internal L:Documented_Item
# A *central* structure containing documented item
#
struct Documented_Item
    _filename::String
    _extracted_tag::Extract_Tag_Result
    _doc::Union{Nothing,
                Extracted_Comment}
    _code::Extracted_Item_Base
end

#+Documented_Item
filename(di::Documented_Item)::String = di._filename
#+Documented_Item
tokenized(di::Documented_Item)::Tokenized = di._extracted_tag._tokenized
#+Documented_Item
tags(di::Documented_Item)::Vector{String} = tags(di._extracted_tag)
#+Documented_Item
identifier(di::Documented_Item)::String = identifier(di._code)
#+Documented_Item
raw_string_doc(di::Documented_Item)::String = (di._doc==nothing) ? "" : raw_string(di._doc)
#+Documented_Item
raw_string_code(di::Documented_Item)::String = (di._code==nothing) ? "" : raw_string(di._code)
#+Documented_Item
raw_string_code_with_body(di::Documented_Item)::String = (di._code==nothing) ? "" : raw_string_with_body(di._code)
#+Documented_Item
contains_tag(di::Documented_Item,tag::String)::Bool = tag ∈ tags(di)
#contains_tag(di::Documented_Item,tag::String)::Bool = contains(==,tags(di),tag)
#+Documented_Item
# Checks if at least one tag of =tag_list= is contained in =di=
contains_tag(di::Documented_Item,tag_list::Array{String,1})::Bool = mapreduce(x->contains_tag(di,x),(x,y)->x||y,tag_list)
#+Documented_Item
line(di::Documented_Item)::Int = line(di._extracted_tag)
#+Documented_Item
skip(di::Documented_Item)::Int = skip(di._code)

#+Documented_Item
# Returns links as string if exists, "" otherwise
link(di::Documented_Item)::String = link(di._extracted_tag)

#+Documented_Item
is_documented_function(di::Documented_Item)::Bool = typeof(di._code) == Extracted_Function
#+Documented_Item
is_documented_structure(di::Documented_Item)::Bool = typeof(di._code) == Extracted_Struct
#+Documented_Item
is_documented_abstract_type(di::Documented_Item)::Bool = typeof(di._code) == Extracted_Abstract
#+Documented_Item
is_documented_enum_type(di::Documented_Item)::Bool = typeof(di._code) == Extracted_Enum
#+Documented_Item
is_documented_variable_type(di::Documented_Item)::Bool = typeof(di._code) == Extracted_Variable
#+Documented_Item
is_documented_macro_type(di::Documented_Item)::Bool = typeof(di._code) == Extracted_Macro

#+Links,Documented_Item
function create_file_org_link(di::Documented_Item)::String
    return create_file_org_link(filename(di),line(di))
end

#+Error,Documented_Item
# Use data from [[Documented_Item][]] to add a link after the message
function warning_message(di::Documented_Item,message::String)::Nothing
    warning_message("$(message) $(create_file_org_link(di))")
end

function create_documented_item(extracted_tag::Extract_Tag_Result;
                                filename::String="")::Documented_Item
    idx=skip(extracted_tag)
    tok=extracted_tag._tokenized
    doc=extract_comment(tok,idx)
    if doc!=nothing
        idx=skip(doc)
    else
        # for the moment do not overload logs
        # println("--- WARNING: undocumented stuff file: $(filename):$(first(startpos(tok[idx]))) $(untokenize(tok[idx]))")
    end
    code=extract_code(tok,idx)
    @assert code!=nothing

    return Documented_Item(filename,extracted_tag,doc,code)
end

function create_documented_item(tok::Tokenized,idx::Int;
                                filename::String="")::Documented_Item
    extracted_tag = extract_tag(tok,idx)
    @assert extracted_tag != nothing
    return create_documented_item(extracted_tag,filename=filename)
end

#+API L:create_documented_item_array_filename
#
# Reads a Julia code file and returns an array of documented items.
#
function create_documented_item_array(filename::String)::Array{Documented_Item,1}

    if VERSION < v"0.7"
        code = readstring(filename)
    else
        code = read(filename,String)
    end

    tok=tokenized(code)

    n=length(tok)
    idx=1
    toReturn=Array{Documented_Item,1}()

    while idx<=n
        extracted_tag = find_tag(tok,idx)

        # no more tag?
        #
        if extracted_tag==nothing
            break
        end

        try
            di=create_documented_item(extracted_tag,filename=filename)
            push!(toReturn,di)
        catch msg
            warning_message("cannot interpret $(filename):$(line(extracted_tag)) $(msg)")
        end

        # important to rescan after tag (and not after code):
        # tag1
        # struct foo
        #   tag2 <- catch me!
        #   constructor
        # end
        idx=skip(extracted_tag)
    end
    return toReturn
end

#+API
#
# Reads an array of Julia code files and returns an array of
# documented items.
#
# *Usage example:*
# #+BEGIN_SRC julia
# create_documented_item_array(["file1","file2",...])
# #+END_SRC
#
# *Note:* instead of a list of files you can also specify a directory,
# see [[create_documented_item_array_dir][]]
function create_documented_item_array(filename_list::Array{String,1})::Array{Documented_Item,1}
    docItem_array = Array{Documented_Item,1}()

    for file in filename_list
        docItem_array = vcat(docItem_array,create_documented_item_array(file))
    end

    return docItem_array
end

function scan_directory_return_file_names(directory_name::AbstractString)
    expanded_dirname=expanduser(directory_name)
    @assert isdir(expanded_dirname)

    file_names = Vector{String}()

    for (root, dirs, files) in walkdir(expanded_dirname)

        for file in files
            push!(file_names,joinpath(root, file))
        end

    end

    return file_names
end

function filter_julia_source_code_file_names(file_names::Array{<:AbstractString,1})
    file_name_pattern = r".jl$"
    return filter(x->occursin(file_name_pattern,x),file_names)
end

#+API L:create_documented_item_array_dir
#
# Reads all *.jl files in a directory and returns an array of
# documented items.
#
function create_documented_item_array_dir(dirname::AbstractString)
    # expanded_dirname=expanduser(dirname)
    # @assert isdir(expanded_dirname)
    # # tips from https://stackoverflow.com/questions/20484581/search-for-files-in-a-folder
    # files=filter(x->contains(x,r".jl$"), readdir(expanded_dirname))
    # map!(x->expanded_dirname*x,files,files)

    file_names=scan_directory_return_file_names(dirname)
    files=filter_julia_source_code_file_names(file_names)
    return create_documented_item_array(files)
end


text:
module FHist

export Hist1D, binedges, bincounts, bincenters, binerrors, nbins, integral
export sample, lookup, cumulative, normalize, restrict, rebin
export atomic_push!

export Hist2D, project, profile, transpose

using StatsBase, RecipesBase, UnicodePlots, Statistics
import LinearAlgebra: normalize, normalize!
using Base.Threads: SpinLock

const _default_overflow = false

struct Hist1D{T<:Real,E} <: AbstractHistogram{T,1,E}
    hist::Histogram{T,1,E}
    sumw2::Array{Float64, 1}
    hlock::SpinLock
    overflow::Bool
    function Hist1D(h::Histogram{T,1,E}, sw2 = copy(h.weights); overflow=_default_overflow) where {T,E}
        return new{T,E}(h, sw2, SpinLock(), overflow)
    end
end

struct Hist2D{T<:Real,E} <: AbstractHistogram{T,2,E}
    hist::Histogram{T,2,E}
    sumw2::Array{Float64, 2}
    hlock::SpinLock
    overflow::Bool
    function Hist2D(h::Histogram{T,2,E}, sw2 = copy(h.weights); overflow=_default_overflow) where {T,E}
        return new{T,E}(h, sw2, SpinLock(), overflow)
    end
end

include("./utils.jl")
include("./hist1d.jl")
include("./hist2d.jl")
include("./arithmatics.jl")
include("./plot.jl")
end


text:
C(x)=((x^2+1.)^0.5+x)^2
Shmin(Sv,Pp,μ)=((Sv-Pp)/C(μ))+Pp
Shmax(Sv,Pp,μ)=C(μ)*(Sv-Pp)+Pp

@userplot zobackogram

@recipe function f(h::zobackogram)

    Sv, Pp, μ = h.args

    Smin=Shmin(Sv,Pp,μ)
    Smax=Shmax(Sv,Pp,μ)

    #45° line. slope=1
    L1=[0,Smax*1.1]

    #Shmin line
    L2=[Smin Smin;Smin Sv]

    #Shmax Horizonal line
    L3=[Sv Smax;Smax Smax]

    #Shmax diagonal Line
    L4=[Smin Sv;Sv Smax]

    #Internal Lines
    L5=[Smin Sv; Sv Sv;Sv Smax]

    xlabel := "Shmin [psi/ft]"
    ylabel := "Shmax [psi/ft]"
    legend := :bottomright
    legendfontsize --> 6

    #45° line. slope=1
    @series begin
        seriestype := :path
        linewidth := 3
        linecolor := :black
        label := "45° Line"
        L1,L1
    end

    @series begin
        seriestype := :path
        linewidth := 2
        linecolor := :darkblue
        label := "Shmin Line"
        L2[:,1],L2[:,2]
    end

    @series begin
        seriestype := :path
        linewidth := 2
        linecolor := :darkred
        label := "Shmax Horizontal Line"
        L3[:,1],L3[:,2]
    end

    @series begin
        seriestype := :path
        linewidth := 2
        linecolor := :darkgreen
        label := "Shmax Strike Slip"
        L4[:,1],L4[:,2]
    end

    @series begin
        seriestype := :path
        linewidth := 1
        linestyle := :dash
        linecolor := :black
        label := "Limits"
        series_annotations := map(x->Plots.text(x ,:left,:top,6),["Shmin=$(round(Smin,digits=1))",
                                "Sv=$(round(Sv,digits=1))","Shmax=$(round(Smax,digits=1))"])
        L5[:,1],L5[:,2]
    end

end


text:
"""
    struct UnstructPSpace{A,B,C,D,E,F,G,H,I,J,K,L} <: AbstractPhysicalSpace
        cells::A # all information: cell, line, vertex
        points::B # locations of vertex points
        cellid::C # node indices of elements
        cellType::D # inner/boundary cell
        cellNeighbors::E # neighboring cells id
        cellFaces::F # cell edges id
        cellCenter::G # cell center location
        cellArea::H # cell size
        facePoints::I # ids of two points at edge
        faceCells::J # ids of two cells around edge
        faceCenter::K # edge center location
        faceType::L
    end

Physical space with unstructured mesh
"""
struct UnstructPSpace{
    A,
    B<:AbstractMatrix{<:AbstractFloat},
    C<:AbstractMatrix{<:Integer},
    D<:AbstractVector{<:Integer},
    E<:AbstractVector{<:AbstractFloat},
    F<:AbstractArray{<:AbstractFloat,3},
} <: AbstractPhysicalSpace
    cells::A # all information: cell, line, vertex
    points::B # locations of vertex points
    cellid::C # node indices of elements
    cellType::D # inner/boundary cell
    cellNeighbors::C # neighboring cells id
    cellFaces::C # cell edges id
    cellCenter::B # cell center location
    cellArea::E # cell size
    cellNormals::F # cell unit normal vectors
    facePoints::C # ids of two points at edge
    faceCells::C # ids of two cells around edge
    faceCenter::B # edge center location
    faceType::D # inner/boundary face
    faceArea::E # face area
end

function UnstructPSpace(file::T) where {T<:AbstractString}
    cells, points = read_mesh(file)
    p = mesh_connectivity_2D(cells, points)

    return UnstructPSpace((cells, points)..., p...)
end

function write_vtk(ks::T1, ctr) where {T1<:AbstractSolverSet}
    write_vtk(ks.ps, ctr)

    return nothing
end

function write_vtk(ps::UnstructPSpace, ctr)
    cdata = zeros(length(ctr), length(ctr[1].w))
    for i in eachindex(ctr)
        cdata[i, :] .= ctr[i].prim
        if size(cdata, 2) > 1
            cdata[i, end] = 1.0 / cdata[i, end]
        end
    end
    write_vtk(ps.points, ps.cellid, cdata)

    return nothing
end


text:
## Functions for centered moments of moments feature maps
#const m1s2 = [1,2,3,4,5];
const m2lim2 = 10;
function apply_mofm_full(X::Array{Float64}, rowcounts::Int64, mom::StatisticalMoments,options)
   lenresponse = length(mom.moments[1]);
   momnums = options.aux;
   s=size(momnums);
   count =1; i =1;
    if(length(s)==1)
        for i in momnums
            mm  = mom.moments[i];
            for j = 1:m2lim2
                X[rowcounts,count] = get_moment_vec(mm,j);
                count = count +1;
            end
        end
   else
        for mm in mom.moments
            for j = 1:m2lim2
                if(options.aux[i,j]==1)
                   X[rowcounts,count] = get_moment_vec(mm,j);
                   count = count +1;
                end
            end
            i = i +1;
        end
   end
end

function get_num_mofm__full_features(mom::StatisticalMoments, options)
    s=size(options.aux);
    if(length(s)==1)
        lenresponse =  length(options.aux)*m2lim2;
     else
        lenresponse = sum(options.aux);
    end
    return lenresponse;
end

# Defines the feature map
fmofm_full = FeatureMap(apply_mofm_full,get_num_mofm__full_features, "mofm_full_map");

#"returns the  mom_num-th moment of resps, where resps is a vector whose elements can be either numbers or vectors."
function get_moment_vec(resps, mom_num)
    mu = mean(resps);
    if(mom_num==1) return mu end;
    num_sim = length(resps);
    if(num_sim<2)
        println("Not enough simulations to calculate momments, return zeros")
        return 0.0*mu
    end
    moment = sum((mu - r).^mom_num for r in resps)./(num_sim-1.0)
    return sign(moment).*abs(moment).^(1.0/mom_num) # I think the sign does make a difference, as it changes with different labels
end


text:
# -*- coding: utf-8 -*-

# cd ~/.julia/v0.4/DirectX
# julia test/dx9adl_process.jl

include("dx9adl_demo_0001.jl")

println("test dx9adl with relocated dlls")
d9 = DirectX.connect(55 * 16, 55 * 9, "/private/dx9/dx9adl")
Demo_0001.demo_0001(d9)
DirectX.close(d9)
println("ok")


text:
module FinOptions

## Dependencies

using DataArrays
using DataFrames
using Base.Dates

export

        Black76,
        BlackScholes,
        BlackScholesDP,
        GBlackScholes,
        KirkSpreadOption,
        GeoAveRate,
        GDelta,
        Delta,
        Elasticity,
        Gamma,
        DGammaDVol,
        GammaP,
        DDeltaDVol,
        Vega,
        Vomma,
        VegaP,
        DVegaDVol,
        Theta,
        Rho,
        RhoFut,
        Rho2,
        Carry,
        Speed,
        StrikeDelta,
        SABRVolatility,
        RNDensity,
        PutCallParity,
        VannaVolga,
        SwingOptionPricer,
        StrikeLadder,
        VVParams,
        Aug16, Sep16, Oct16, Nov16, Dec16,
        Jan17, Feb17, Mar17, Apr17, May17, Jun17,
        Jul17, Aug17, Sep17, Oct17, Nov17, Dec17,
        Jan18, Feb18, Mar18, Apr18, May18, Jun18,
        Jul18, Aug18, Sep18, Oct18, Nov18, Dec18,
        Q416, Q117, Q217, Q317, Q417, Q118, Q218, Q318, Q418,
        W16, S17, W17, S18,
        CND

include("Black76.jl")
include("BlackScholes.jl")
include("BlackScholesDP.jl")
include("CND.jl")
#include("CNDEV.jl"),
include("GBlackScholes.jl")
include("GDelta.jl")
include("Greeks.jl")
include("PutCallParity.jl")
include("SABR.jl")
include("VannaVolga.jl")
include("SwingOptionPricer.jl")
include("StrikeLadder.jl")
include("Expiry.jl")

# package code goes here

type VVParams
        putStrike
        ATMStrike
        callStrike
        putVol
        ATMVol
        callVol
end

end # module


text:
# 2bit File Format
# ================
#
# Reader and writer of the 2bit file format.
#
# The 2bit file format stores multiple DNA sequences (up to 4 Gbp total) as a
# compact randomly-accessible format.
# See https://genome.ucsc.edu/FAQ/FAQformat.html#format7 for the details.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

module TwoBit

import BioCore
import BioSequences

include("record.jl")
include("reader.jl")
include("writer.jl")

end


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

# case 1: CS-FMU Simulation

myFMU = fmiLoad("SpringPendulum1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])

comp = fmiInstantiate!(myFMU; loggingOn=false)
@test comp != 0

# choose FMU or FMUComponent
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
@assert fmuStruct !== nothing "Unknown fmuStruct, environment variable `FMUSTRUCT` = `$envFMUSTRUCT`"

t_start = 0.0
t_stop = 8.0

# test without recording values (but why?)
solution = fmiSimulateCS(fmuStruct, t_start, t_stop; dt=1e-2)
@test solution.success

# test with recording values
solution = fmiSimulateCS(fmuStruct, t_start, t_stop; dt=1e-2, recordValues=["mass.s", "mass.v"])
@test solution.success
@test length(solution.values.saveval) == t_start:1e-2:t_stop |> length
@test length(solution.values.saveval[1]) == 2

t = solution.values.t
s = collect(d[1] for d in solution.values.saveval)
v = collect(d[2] for d in solution.values.saveval)
@test t[1] == t_start
@test t[end] == t_stop

# reference values from Simulation in Dymola2020x (Dassl)
@test s[1] == 0.5
@test v[1] == 0.0

if ENV["EXPORTINGTOOL"] == "Dymola/2020x" # ToDo: Linux FMU was corrupted
    @test s[end] ≈ 0.509219 atol=1e-1
    @test v[end] ≈ 0.314074 atol=1e-1
end

fmiUnload(myFMU)

# case 2: CS-FMU with input signal

function extForce(t)
    [sin(t)]
end

myFMU = fmiLoad("SpringPendulumExtForce1D", ENV["EXPORTINGTOOL"], ENV["EXPORTINGVERSION"])

comp = fmiInstantiate!(myFMU; loggingOn=false)
@test comp != 0

# choose FMU or FMUComponent
fmuStruct = nothing
envFMUSTRUCT = ENV["FMUSTRUCT"]
if envFMUSTRUCT == "FMU"
    fmuStruct = myFMU
elseif envFMUSTRUCT == "FMUCOMPONENT"
    fmuStruct = comp
end
@assert fmuStruct !== nothing "Unknown fmuStruct, environment variable `FMUSTRUCT` = `$envFMUSTRUCT`"

solution = fmiSimulateCS(fmuStruct, t_start, t_stop; dt=1e-2, recordValues=["mass.s", "mass.v"], inputValueReferences=["extForce"], inputFunction=extForce)
@test solution.success
@test length(solution.values.saveval) > 0
@test length(solution.values.t) > 0

@test t[1] == t_start
@test t[end] == t_stop

# reference values from Simulation in Dymola2020x (Dassl)
@test [solution.values.saveval[1]...] == [0.5, 0.0]
@test sum(abs.([solution.values.saveval[end]...] - [0.613371, 0.188633])) < 0.2
fmiUnload(myFMU)


text:
module ThArrays

using Libdl
using Requires

@static if Sys.islinux()
    using LibTorchCAPI_jll
elseif Sys.isapple()
    const libtorch_capi = :libtorch_capi
end

export TorchNumber, Tensor, Scalar, eltype_id,
    ThC, ThAD, TrackerAD, ThJIT,
    Device, CPU, CUDA, to, on

const PROJECT_DIR = (@__DIR__) |> dirname

function __init__()
    @static if Sys.isapple()
        push!(Libdl.DL_LOAD_PATH, joinpath(PROJECT_DIR, "deps/lib"))
        Libdl.dlopen(joinpath(PROJECT_DIR, "deps/lib/libtorch_capi"))
    else
        @static if VERSION < v"1.6.0"
            Libdl.dlopen(LibTorchCAPI_jll.libtorch_capi_path)
        end
    end
    @async handle_error_in_julia()
    @require Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" @eval include("compat/tracker.jl")
end

function handle_error_in_julia()
    err_handler = "jl_error"
    ccall((:set_error_handler, libtorch_capi),
          Cvoid, (Cstring, Csize_t),
          pointer(err_handler), length(err_handler))
end

const TYPE_MAP = Dict{Type, Int8}(
    ### float
    Float16 => 5,
    Float32 => 6,
    Float64 => 7,
    ### bool and char
    Bool => 11,
    # Char => 1, # Char in Julia is not single byte
    ### int
    Int8 => 1,
    # UInt8 => 1,
    Int16 => 2,
    # UInt16 => 2,
    Int32 => 3,
    # UInt32 => 3,
    Int64 => 4,
    # UInt64 => 4,
    # Int128 => ?,
    # UInt128 => ?,
)

const REVERSE_TYPE_MAP = Dict(reverse(p) for p in TYPE_MAP)

TorchNumber = Union{Float16, Float32, Float64,
                    Bool,
                    Int8, Int16, Int32, Int64}

include("tensor.jl")
include("scalar.jl")
include("thc/thc.jl")
include("common-methods.jl")
include("ad.jl")
include("th-jit.jl")

end


text:
#	GridWhale module
#	Copyright (c) 2020 Kronosaur Productions, LLC. All Rights Reserved.
#
#	This file provides functions for interacting with Google Drive.

export GoogleDriveFileResource
export input_resource
export load_resource_as

struct GoogleDriveFileResource
	resID::UInt32
end

function input_resource(res::Type{GoogleDriveFileResource})::GoogleDriveFileResource
	resid::UInt32 = ccall("InputGoogleDriveFile", Cuint, ())
	GoogleDriveFileResource(resid)
end

function load_resource_as(res::GoogleDriveFileResource, mimetype::AbstractString)
	if mimetype == "text/csv"
		rawtable = ccall("LoadResourceAs", Any, (Cuint, Cstring), res.resID, mimetype)

		colCount = size(rawtable, 2)
		rowCount = size(rawtable[2, 1])

		table = DataTable("", 0, 0, [], [])

		for i = 1:colCount
			table[rawtable[1, i]] = rawtable[2, i]
		end

		return table
	else
		return nothing
	end
end


text:
module QJuliaLanMR

using QJuliaInterface
using QJuliaEnums
using QJuliaBlas
using QJuliaReduce
using QJuliaSolvers

using LinearAlgebra
using MPI

norm2   = QJuliaReduce.gnorm2
rdot    = QJuliaReduce.reDotProduct
verbose = false

const ε = eps(Float64)

@inline function LanczosStep(Mat::Any, vk::AbstractArray, vkm1::AbstractArray, βold::Float64, σ::Float64, pk::AbstractArray)
    Mat(pk, vk)
    pk .=@. pk - βold*vkm1
    # Compute Lancsoz coefficient α
    α = rdot(vk, pk)
    pk .=@. pk - (α+σ)*vk
    # Compute Lancsoz coefficient β
    β = sqrt(rdot(pk, pk))
    if β  < ε ; error("Algorithm breakdown."); end
    vkm1 .=@. vk
    vk .=@. (1.0 / β)*pk
    return  α, β
end

function solver(x::AbstractArray, b::AbstractArray, Mat::Any, MatSloppy::Any, param::QJuliaSolvers.QJuliaSolverParam_qj)

    if verbose == true; println("Running Lanczos version of MR solver in maximum " , param.Nsteps, " step(s)."); end

    QJuliaReduce.set_blas_global_reduction(param.global_reduction)

    if MPI.Initialized() == false; error("MPI was not inititalized, copy source field to the solution."); end

    if (param.maxiter == 0) || (param.Nsteps == 0)
      if param.use_init_guess == false
        x .=@. 0.0
      end
      return
    end #if

    mixed = (param.dtype_sloppy != param.dtype)

    local r     = zero(typeof(x)(undef, length(x)))
    local Av    = zero(typeof(x)(undef, length(x)))
    # explicit Lanczos vectors
    local vk    = zero(typeof(x)(undef, length(x)))
    local vkm1  = zero(typeof(x)(undef, length(x)))
    # explicit Lanczos matrix
    # this is just (sub-)diagonal elements (but transposed when stored!):
    # #(0) β2      ... βk-1  βk
    # α1   α2   α3 ... αk-1  αk
    # β2   β3   β4 ... βk    βk+1
    #
    # R-matrix has the following structure
    # γ1 δ2 ϵ3  0  0  0  ...
    # 0  γ2 δ3 ϵ4  0  0  ...
    # 0  0  γ3 δ4  ϵ5 0  ...
    # 0 ................ γ_k

    local Tk    = zero(Matrix{Float64}(undef, (param.maxiter+1), 2))
    # some aux fields:
    local dk    = zero(typeof(x)(undef, length(x)))
    local dkm1  = zero(typeof(x)(undef, length(x)))

    bnorm = sqrt(norm2(b))  #Save norm of b
    rnorm = 0.0     #if zero source then we will exit immediately doing no work

    σ = param.shift

    if param.use_init_guess == true
      #r = b - Ax0
      Mat(r, x)
      r .=@. r - σ*x
      r .=@. b - r
      rnorm = sqrt(norm2(r))
    else
      rnorm = bnorm
      r .=@. b
    end

    if( rnorm == 0.0 ); return; end
    if verbose == true; println("MR: Initial residual = ", rnorm); end

    # Load init vectors
    vk   .=@. (1.0 / rnorm) * r
    # Givens parameters
    cs = -1.0; sn = 0.0
    # Matrix norm:
    Anorm = 0.0

    # Lancsoz coeffs
    α = 0.0; β = 0.0; βold = rnorm
    # aux parameters
    δk   = 0.0; ϵold = 0.0;
    #start iterations
    for k in 1:param.maxiter
      # Lanczos step:
      α, β = LanczosStep(Mat, vk, vkm1, βold, σ, Av); Tk[k,1] = α; Tk[k,2] = β
      # Givens QR
      # produce middle two enties in the last coloumn of Tk
      δ  = cs*δk + sn*α
      γk = sn*δk - cs*α
      # produce first two entries of T_{k+1}*e_{k+1}
      ϵ  =  sn*β
      δk = -cs*β
      # eliminate β
      γinv  = 1.0 / sqrt(γk*γk+β*β)
      sn    = β  * γinv
      cs    = γk * γinv
      # Solution update (using Av as tmp)
      Av   .=@. γinv*vkm1 - (γinv*δ)*dk - (γinv*ϵold)*dkm1
      dkm1 .=@. dk
      dk   .=@. Av
      #
      τk = rnorm*cs
      x .=@. x + τk*dk
      #residual norm and matrix norm:
      Arnorm = rnorm*sqrt(γk*γk+δk*δk)
      #
      rnorm  *= sn
      #
      Anorm = k == 1 ? sqrt(α*α+β*β) : max(Anorm, sqrt(βold*βold+α*α+β*β))

      if(verbose == true); println("Lanczos MR: ", k,"  residual norm = ", rnorm, "\t ( estimated matrix norm is ", Anorm,  ")"); end
      # Store current β and ϵ for the next cycle
      βold = β; ϵold = ϵ

    end #for itern

    Mat(r, x)

    r .=@. b - r
    rnorm = sqrt(norm2(r))

    param.true_res = rnorm / bnorm
    if verbose == true; println("MR: Relative residual: true = ", param.true_res); end

    QJuliaReduce.reset_blas_global_reduction()

    # Return computed matrix norm
    return Anorm

end #MR

end #QJuliaLanMR

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia