text:
using CoverageTools
using Test

if Base.VERSION < v"1.1"
    isnothing(x) = false
    isnothing(x::Nothing) = true
end

withenv("DISABLE_AMEND_COVERAGE_FROM_SRC" => nothing) do

@testset "iscovfile" begin
    # test our filename matching. These aren't exported functions but it's probably
    # a good idea to have explicit tests for them, as they're used to match files
    # that get deleted
    @test CoverageTools.iscovfile("test.jl.cov")
    @test CoverageTools.iscovfile("test.jl.2934.cov")
    @test CoverageTools.iscovfile("/home/somebody/test.jl.2934.cov")
    @test !CoverageTools.iscovfile("test.ji.2934.cov")
    @test !CoverageTools.iscovfile("test.jl.2934.cove")
    @test !CoverageTools.iscovfile("test.jicov")
    @test !CoverageTools.iscovfile("test.c.cov")
    @test CoverageTools.iscovfile("test.jl.cov", "test.jl")
    @test !CoverageTools.iscovfile("test.jl.cov", "other.jl")
    @test CoverageTools.iscovfile("test.jl.8392.cov", "test.jl")
    @test CoverageTools.iscovfile("/somedir/test.jl.8392.cov", "/somedir/test.jl")
    @test !CoverageTools.iscovfile("/otherdir/test.jl.cov", "/somedir/test.jl")
end # testset

@testset "isfuncexpr" begin
    @test CoverageTools.isfuncexpr(:(f() = x))
    @test CoverageTools.isfuncexpr(:(function() end))
    @test CoverageTools.isfuncexpr(:(function g() end))
    @test CoverageTools.isfuncexpr(:(function g() where {T} end))
    @test !CoverageTools.isfuncexpr("2")
    @test !CoverageTools.isfuncexpr(:(f = x))
    @test CoverageTools.isfuncexpr(:(() -> x))
    @test CoverageTools.isfuncexpr(:(x -> x))
    @test CoverageTools.isfuncexpr(:(f() where A = x))
    @test CoverageTools.isfuncexpr(:(f() where A where B = x))
end # testset

@testset "Processing coverage" begin
    cd(dirname(@__DIR__)) do
        datadir = joinpath("test", "data")
        # Process a saved set of coverage data...
        r = process_file(joinpath(datadir, "CoverageTools.jl"))

        # ... and memory data
        malloc_results = analyze_malloc(datadir)
        filename = joinpath(datadir, "testparser.jl.9172.mem")
        @test malloc_results == [CoverageTools.MallocInfo(96669, filename, 2)]

        lcov = IOBuffer()
        # we only have a single file, but we want to test on the Vector of file results
        LCOV.write(lcov, FileCoverage[r])
        expected = read(joinpath(datadir, "tracefiles", "expected.info"), String)
        if Sys.iswindows()
            expected = replace(expected, "\r\n" => "\n")
            expected = replace(expected, "SF:test/data/CoverageTools.jl\n" => "SF:test\\data\\CoverageTools.jl\n")
        end
        @test String(take!(lcov)) == expected

        # LCOV.writefile is a short-hand for writing to a file
        lcov = joinpath(datadir, "lcov_output_temp.info")
        LCOV.writefile(lcov, FileCoverage[r])
        expected = read(joinpath(datadir, "tracefiles", "expected.info"), String)
        if Sys.iswindows()
            expected = replace(expected, "\r\n" => "\n")
            expected = replace(expected, "SF:test/data/CoverageTools.jl\n" => "SF:test\\data\\CoverageTools.jl\n")
        end
        @test String(read(lcov)) == expected
        # tear down test file
        rm(lcov)

        # test that reading the LCOV file gives the same data
        lcov = LCOV.readfolder(datadir)
        @test length(lcov) == 1
        r2 = lcov[1]
        r2_filename = r2.filename
        if Sys.iswindows()
            r2_filename = replace(r2_filename, '/' => '\\')
        end
        @test r2_filename == r.filename
        @test r2.source == ""
        @test r2.coverage == r.coverage[1:length(r2.coverage)]
        @test all(isnothing, r.coverage[(length(r2.coverage) + 1):end])
        lcov2 = [FileCoverage(r2.filename, "sourcecode", CoverageTools.CovCount[nothing, 1, 0, nothing, 3]),
                 FileCoverage("file2.jl", "moresource2", CoverageTools.CovCount[1, nothing, 0, nothing, 2]),]
        lcov = merge_coverage_counts(lcov, lcov2, lcov)
        @test length(lcov) == 2
        r3 = lcov[1]
        @test r3.filename == r2.filename
        @test r3.source == "sourcecode"
        r3cov = CoverageTools.CovCount[x === nothing ? nothing : x * 2 for x in r2.coverage]
        r3cov[2] += 1
        r3cov[3] = 0
        r3cov[5] = 3
        @test r3.coverage == r3cov
        r4 = lcov[2]
        @test r4.filename == "file2.jl"
        @test r4.source == "moresource2"
        @test r4.coverage == lcov2[2].coverage

        # Test a file from scratch
        srcname = joinpath("test", "data", "testparser.jl")
        covname = srcname*".cov"
        # clean out any previous coverage files. Don't use clean_folder because we
        # need to preserve the pre-baked coverage file CoverageTools.jl.cov
        clean_file(srcname)
        cmdstr = "include($(repr(srcname))); using Test; @test f2(2) == 4"
        run(`$(Base.julia_cmd()) --startup-file=no --code-coverage=user -e $cmdstr`)
        r = process_file(srcname, datadir)

        target = if VERSION >= v"1.5-"
            CoverageTools.CovCount[nothing, 1, nothing, 0, nothing, 0, nothing, nothing, nothing, nothing, 0, nothing, nothing, nothing, nothing, 0, 0, nothing, nothing, 0, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing]
        else
            CoverageTools.CovCount[nothing, 2, nothing, 0, nothing, 0, nothing, nothing, nothing, nothing, 0, nothing, nothing, nothing, nothing, nothing, 0, nothing, nothing, 0, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing]
        end
        target_disabled = map(x -> (x !== nothing && x > 0) ? x : nothing, target)
        @test r.coverage == target

        covtarget = (sum(x -> x !== nothing && x > 0, target), sum(!isnothing, target))
        @test get_summary(r) == covtarget
        if VERSION >= v"1.5-"
            @test get_summary(process_folder(datadir)) == (98, 107)
        else
            @test get_summary(process_folder(datadir)) == (98, 106)
        end

        r_disabled = withenv("DISABLE_AMEND_COVERAGE_FROM_SRC" => "yes") do
            process_file(srcname, datadir)
        end

        @test r_disabled.coverage == target_disabled
        amend_coverage_from_src!(r_disabled.coverage, r_disabled.filename)
        @test r_disabled.coverage == target

        # Handle an empty coverage vector
        emptycov = FileCoverage("", "", [])
        @test get_summary(emptycov) == (0, 0)

        @test isempty(CoverageTools.process_cov(joinpath("test", "fakefile"), datadir))

        # test clean_folder
        # set up the test folder
        datadir_temp = joinpath("test", "data_temp")
        cp(datadir, datadir_temp)
        # run clean_folder
        clean_folder(datadir_temp)
        # .cov files should be deleted
        @test !isfile(joinpath(datadir_temp, "CoverageTools.jl.cov"))
        # other files should remain untouched
        @test isfile(joinpath(datadir_temp, "CoverageTools.jl"))
        # tear down test data
        rm(datadir_temp; recursive=true)

        # test exclusion markers
        srcname = joinpath("test", "exclusions.jl")
        covname = srcname*".cov"
        clean_file(srcname)
        cmdstr = "include($(repr(srcname))); main()"
        run(`$(Base.julia_cmd()) --startup-file=no --code-coverage=user -e $cmdstr`)
        r = withenv("DISABLE_AMEND_COVERAGE_FROM_SRC" => "yes") do
            process_file(srcname, "test")
        end
        # FIXME: coverage information for this function is useless with Julia 1.0
        if VERSION >= v"1.1-"
            if VERSION >= v"1.6-"
                @test r.coverage == [1, 1, 1, 1, nothing, 1, nothing, 1, nothing]
            elseif VERSION >= v"1.5-"
                @test r.coverage == [nothing, 1, 1, 1, nothing, 1, nothing, 1, nothing]
            else
                @test r.coverage == [nothing, 2, 1, 1, nothing, 1, nothing, 1, nothing, nothing]
            end
            amend_coverage_from_src!(r.coverage, r.filename)
            if VERSION >= v"1.6-"
                @test r.coverage == [1, 1, nothing, 1, nothing, nothing, nothing, 1, nothing]
            elseif VERSION >= v"1.5-"
                @test r.coverage == [nothing, 1, nothing, 1, nothing, nothing, nothing, 1, nothing]
            else
                @test r.coverage == [nothing, 2, nothing, 1, nothing, nothing, nothing, 1, nothing, nothing]
            end
        end
    end

    srcdir = joinpath(@__DIR__, "BustedPackage", "src")
    clean_folder(srcdir)
    cmdstr = "using BustedPackage; BustedPackage.greet()"
    cd("BustedPackage") do
        run(`$(Base.julia_cmd()) --startup-file=no --code-coverage=user --project -e $cmdstr`)
    end
    bustedfile = joinpath(srcdir, "parseerr.jl")
    if VERSION.major == 1 && VERSION.minor < 5
        msg = "parsing error in $bustedfile:7: space before \"[\" not allowed in \"i [\""
    elseif VERSION.major == 1 && VERSION.minor == 5
        msg = "parsing error in $bustedfile:7: space before \"[\" not allowed in \"i [\" at none:4"
    else
        msg = "parsing error in $bustedfile:7: invalid iteration specification"
    end
    @test_throws Base.Meta.ParseError(msg) process_file(bustedfile, srcdir)
end # testset

@testset "malloc.jl" begin
    @testset "types" begin
        a = CoverageTools.MallocInfo(1, "", 1)
        b = CoverageTools.MallocInfo(2, "", 1)
        @test CoverageTools.sortbybytes(a, b)
    end
end # testset

@testset "parser.jl" begin
    @testset "function_body_lines!" begin
        @testset "ast.head == :module" begin
            code = """
                module Foo
                eval() = "foo"
                eval(x) = "bar"
                f() = "baz"
                end # module
                """
            ast = Meta.parse(code)
            # remove the top-level line number nodes
            filter!(x -> !(x isa LineNumberNode), ast.args[end].args)
            coverage = CoverageTools.CovCount[]
            lineoffset = 1
            infunction = false
            flines = []
            CoverageTools.function_body_lines!(flines, ast, coverage, lineoffset, infunction)
            @test flines == [4]
        end

        @testset "ast.head == :line" begin
            ast = Expr(:line, 12345, nothing)
            coverage = CoverageTools.CovCount[]
            lineoffset = 1
            infunction = true
            flines = []
            CoverageTools.function_body_lines!(flines, ast, coverage, lineoffset, infunction)
            @test flines == [12345]
        end
    end
end

@testset "CoverageTools.jl" begin
    @testset "clean_file" begin
        mktempdir() do tmp_dir
            cd(tmp_dir) do
                source_file = joinpath(tmp_dir, "foo.jl")
                cov_file = joinpath(tmp_dir, "foo.jl.12345.cov")
                touch(source_file)
                touch(cov_file)
                @test isfile(source_file)
                @test isfile(cov_file)
                CoverageTools.clean_file(source_file)
                @test isfile(source_file)
                @test !isfile(cov_file)
            end
        end
    end
end

end # withenv


text:
module multi
using BenchmarkTools
using PyPlot
PyPlot.pygui(true)

# println(Threads.nthreads())

function write_tid()
    n_dim = 5_000_000
    vec = zeros(Float64, n_dim)
    Threads.@threads for kk = 1:n_dim
        vec[kk] = Threads.threadid()
    end
    # for kk = 1:n_dim
    #     vec[kk] = Threads.threadid()
    # end
end

function mandelbrot(x::T, y::T; max_iter::Int = 100_000) where {T <: Real}
    # x and y are interpreted as complex numbers
    # z = x + i*y
    # z * z = (x^2 - y^2, 2*x*y)

    z0_x = zero(T)
    z0_y = zero(T)

    critical_norm = 2 * one(T)

    x_old = x
    y_old = y
    n_old = x^2 + y^2

    x_new = zero(T)
    y_new = zero(T)

    iter = 0
    for outer iter = 1:max_iter
        x_new = x_old^2 - y_old^2 + x
        y_new = 2 * x_old * y_old + y
        n_new = x_new^2 + y_new^2
        if n_new > critical_norm
            break
        end
        # if abs(n_new - n_old) <= 10 * eps(n_old)
        #     break
        # end
        n_old = n_new
        x_old = x_new
        y_old = y_new
    end

    return iter

end

function mandelbrot_matrix(;
    max_iter = 100,
    x_min = -2.0,
    x_max = 1.0,
    y_min = -1.0,
    y_max = -y_min,
    n_y_pixels = 1000,
    n_x_pixels = ceil(Int, (x_max - x_min) / (y_max - y_min) * n_y_pixels)
    )

    x_range = range(x_min, x_max, n_x_pixels)
    y_range = range(y_min, y_max, n_y_pixels)

    pixels = zeros(Float64, n_y_pixels, n_x_pixels)
    Threads.@sync for ii = 1:n_x_pixels # Threads.@threads
        Threads.@spawn for jj = 1:n_y_pixels
            @inbounds pixels[jj, ii] = (max_iter - mandelbrot(x_range[ii], y_range[jj]; max_iter = max_iter) ) / max_iter
        end
    end

    return pixels
end

function test_tid()

    b = @benchmark write_tid()
    show(stdout, MIME("text/plain"), b)
    println("\n\n")
    @time write_tid()
end

function plot_mandelbrot()

    @benchmark pixels = mandelbrot_matrix(n_y_pixels = 3000, x_min = -1.5, x_max = -1.2, y_min = -0.1)

    # PyPlot.figure()
    # PyPlot.imshow(pixels, cmap = PyPlot.cm.cubehelix )

end

end

ext: jl

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using JuMP, BARON, Base.Test

m = Model(solver=BaronSolver())
ub = [2, 2, 1]
@variable(m, 12 ≤ x[i=1:4] ≤ 60, start = 24)

@NLconstraints(m, begin
    x[3] ≤ x[4]
    x[2] ≤ x[1]
end)

@NLobjective(m, Min, 6.931 - x[1]*x[2]/(x[3]*x[4])^2 + 1)

solve(m)


text:
# Hinge loss function

export SqrHingeLoss

"""
**Squared Hinge loss**

  SqrHingeLoss(y, μ=1.0)

Returns the function
```math
f(x) = μ⋅∑_i \\max\\{0, 1 - y_i ⋅ x_i\\}^2,
```
where `y` is an array and `μ` is a positive parameter.
"""
struct SqrHingeLoss{R <: Real, T <: AbstractArray{R}} <: ProximableFunction
  y::T
  mu::R
  function SqrHingeLoss{R, T}(y::T, mu::R) where {R <: Real, T <: AbstractArray{R}}
    if mu <= 0
      error("parameter mu must be positive")
    else
      new(y, mu)
    end
  end
end

is_separable(f::SqrHingeLoss) = true
is_convex(f::SqrHingeLoss) = true
is_smooth(f::SqrHingeLoss) = true

SqrHingeLoss(b::T, mu::R=1.0) where {R <: Real, T <: AbstractArray{R}} = SqrHingeLoss{R, T}(b, mu)

(f::SqrHingeLoss)(x::T) where {R <: Real, T <: AbstractArray{R}} = f.mu*sum(max.(zero(R), (one(R) .- f.y.*x)).^2)

function gradient!(y::AbstractArray{R}, f::SqrHingeLoss{R, T}, x::AbstractArray{R}) where {R <: Real, T}
	sum = zero(R)
	for i in eachindex(x)
		zz = 1-f.y[i]*x[i]
		z = max(zero(R), zz)
		y[i] = z .> 0 ? -2*f.mu*f.y[i]*zz : 0
		sum += z^2
	end
	return f.mu*sum
end

function prox!(z::AbstractArray{R}, f::SqrHingeLoss{R, T}, x::AbstractArray{R}, gamma::R=one(R)) where {R, T}
    v = zero(R)
    for k in eachindex(x)
        if f.y[k]*x[k] >= 1
            z[k] = x[k]
        else
            z[k] = (x[k] + 2*f.mu*gamma*f.y[k])/(1+2*f.mu*gamma*f.y[k]^2)
            v += (1-f.y[k]*z[k])^2
        end
    end
    return f.mu*v
end

fun_name(f::SqrHingeLoss) = "squared hinge loss"
fun_dom(f::SqrHingeLoss) = "AbstractArray{Real}"
fun_expr(f::SqrHingeLoss) = "x ↦ μ * sum( max(0,1-b*x_i)^2, i=1,...,n )"
fun_params(f::SqrHingeLoss) = "b = $(typeof(f.y)), μ = $(f.mu)"

function prox_naive(f::SqrHingeLoss{R, T}, x::AbstractArray{R}, gamma::R=one(R)) where {R, T}
    flag = f.y.*x .<= 1
    z = copy(x)
    z[flag] = (x[flag] .+ 2 .* f.mu.*gamma.*f.y[flag])./(1. + 2 .* f.mu.*gamma.*f.y[flag].^2)
    return z, f.mu*sum(max.(0.0, 1 .- f.y.*z).^2)
end


text:
using Femshop
using Test

@testset "Femshop.jl" begin
    init_femshop("poisson1d");

    useLog("poisson1dlog")

    domain(1)                      # dimension
    functionSpace(order=3)         # basis function polynomial order

    mesh(LINEMESH, elsperdim=20)   # build uniform LINEMESH with 20 elements

    u = variable("u")              # make a scalar variable with symbol u
    testSymbol("v")                # sets the symbol for a test function

    boundary(u, 1, DIRICHLET, 0)  # boundary condition for BID 1 is Dirichlet with value 0

    coefficient("f", "-100*pi*pi*sin(10*pi*x)*sin(pi*x) - pi*pi*sin(10*pi*x)*sin(pi*x) + 20*pi*pi*cos(10*pi*x)*cos(pi*x)")
    weakForm(u, "-grad(u)*grad(v) - f*v")

    solve(u);

    allerr = zeros(size(Femshop.grid_data.allnodes,2));
    exact(x) = sin(10*pi*x)*sin(pi*x);

    for i=1:size(Femshop.grid_data.allnodes,2)
        x = Femshop.grid_data.allnodes[1,i];
        err = abs(u.values[i] - exact(x));
        allerr[i] = err;
    end
    maxerr = maximum(abs, allerr);
    println("max error = "*string(maxerr));

    @test(maxerr < 0.01)

    # Write your tests here.
end


text:
abstract type AbstractSymTensor{S, T, N} end
const AbstractSymMatrix{S, T} = AbstractSymTensor{S, T, 2}

"""
A symmetric tensor is represented as a tensor map over the space of
V_row (codomain) with tensor product structure of (((V1 ⊗ V2) ⊗ V3) ⊗
...) and the space of V_col (domain) with tensor product structure of
(...⊗ (V_{n-3} ⊗ (V_{n-1} ⊗ V_n))). The blocks is also represented as a
matrix
"""

mutable struct SymTensor{S, T, N} <: AbstractSymTensor{S, T, N}
    charge   :: S
    space    :: NTuple{N, VectorSpace{S}}
    blocks   :: SortedDict{Sector{S, N}, Array{T, N}}

    function SymTensor(charge :: S,
                       space  :: NTuple{N, VectorSpace{S}},
                       blocks :: SortedDict{Sector{S, N}, Array{T, N}}) where{S,T,N}
        sects, sizes = _allsectorsandsizes(charge, space)
        i = 1
        for (s,d) in blocks
            s == sects[i] || throw("$(keys(blocks)) SectorMismatch(), $s vs $(sects[i])")
            size(d) == sizes[i] || throw(SizeMismatch())
            i += 1
        end
        new{S,T,N}(charge, space, blocks)
    end
    SymTensor{S,T,N}(c,s,d::SortedDict{Sector{S, N}, Array{T, N}}) where{S,T,N} = new{S,T,N}(c,s,d)
end

# function SymTensor{S,T1,N}(A::SymTensor{S,T2,N}) where {S,T1, T2, N}
#     SymTensor{S,T1,N}(A.charge, A.space, SortedDict{Sector{S,N}, T1}(A.blocks))
# end

const SymMatrix{S, T} = SymTensor{S, T, 2}
SymMatrix(c::S,
          s::NTuple{2, VectorSpace{S}},
          d::SortedDict{Sector{S, 2}, Array{T, 2}}) where{S, T} = SymMatrix{S,T}(c,s,d)

const SymVector{S, T} = SymTensor{S, T, 1}

const U1Tensor{T, N} = SymTensor{U1Charge, T, N}
const U1Matrix{T} = SymMatrix{U1Charge, T}
const U1Vector{T} = SymVector{U1Charge, T}

mutable struct SymDiagonal{S, T<:Number} <: AbstractSymMatrix{S, T}
    charge :: S
    space  :: NTuple{2, VectorSpace{S}}
    blocks   :: SortedDict{Sector{S, 2}, Diagonal{T,Vector{T}}}

    function SymDiagonal(charge :: S,
                         space  :: NTuple{2, VectorSpace{S}},
                         blocks   :: SortedDict{Sector{S,2}, Diagonal{T,Vector{T}}}) where{S,T}
        sects, sizes = _allsectorsandsizes(charge, space)
        i = 1
        for (s,d) in blocks
            s == sects[i] || throw("$blocks SectorMismatch(), $s vs $(sects[i])")
            size(d) == sizes[i] || throw(SizeMismatch())
            i += 1
        end
        # dims(first(space)) == dims(last(space)) ||
        #     throw("SymDiagonal has to be square $space")
        new{S,T}(charge, space, blocks)
    end
end

SymDiagonal{S, T}(c,s,b) where{S,T} = SymDiagonal(c,s,b)
const U1Diagonal{T} = SymDiagonal{Int, T}

@inline diag(A::SymDiagonal) = sort(vcat([diag(blk) for (c,blk) in A.blocks]...), rev=true)

function SymVector(charge::S, v::Vector{T}, isdual::Bool=false) where {S<:AbstractCharge, T}
    V = VectorSpace{S}([charge => length(v)], isdual)
    SymVector{S, T}(charge, (V,), SortedDict([Sector(charge) => v]))
end

function convert(::Type{SymTensor{S, T1, N}},
                 A::SymTensor{S, T2, N}) where {S, T1, T2, N}
    SymTensor{S, T1, N}(A.charge, A.space, A.blocks)
end

@inline space(A::AbstractSymTensor) = A.space
@inline space(A::AbstractSymTensor, l::Int) = A.space[l]
@inline LinearAlgebra.rank(::AbstractSymTensor{S, T, N}) where {S, T, N} = N
@inline charge(A::AbstractSymTensor) = A.charge
@inline sectors(A::AbstractSymTensor) = collect(keys(A.blocks))
@inline blocks(A::AbstractSymTensor) = A.blocks

@inline size(A::AbstractSymTensor) = Tuple(dim(V) for V in A.space)
@inline size(A::AbstractSymTensor, l::Int) = dim(A.space[l])

@inline issimilar(A::T, B::T) where {T<:AbstractSymTensor} =
    A.charge == B.charge && A.space == B.space

@inline isequal(A::T, B::T) where {T<:AbstractSymTensor} =
    issimilar(A, B) && isequal(A.blocks, B.blocks)

@inline isapprox(A::AbstractSymTensor, B::AbstractSymTensor) =
    issimilar(A, B) && all(isapprox(A.blocks[c], B.blocks[c]) for c in keys(A.blocks))

@inline ==(A::T, B::T) where {T<:AbstractSymTensor} = isequal(A, B)

@inline eltype(::AbstractSymTensor{S, T}) where {S, T} = T
@inline eltype(::Type{<:AbstractSymTensor{S, T}}) where {S, T} = T

@inline vtype(::AbstractSymTensor{S, T}) where {S, T} = S
@inline vtype(::Type{<:AbstractSymTensor{S, T}}) where {S, T} = S
@inline vtype(::Type{<:Array{T}}) where {T} = Trivial

#@inline size(A::SymVector) = size(nzblks[1])

# function index_sector(A::AbstractSymTensor{T, N},
#                       s::NTuple{N, Int}) where{T<:Number, N}
#     index = searchsortedfirst(A.sects, s, lt=_sectorlessthan)
#     (index > length(A.sects) || s != A.sects[index]) &&
#         error("sector not found!")
#     index
# end

@inline Base.getindex(A::AbstractSymTensor{S,T,N},
                      s::Sector{S, N}) where {S,T,N} = getindex(A.blocks, s)

function Base.setindex!(A::AbstractSymTensor{S,T,N},
                        blk::Array{T, N},
                        s::Sector{S, N}) where {S,T,N}
    haskey(A.blocks, s) || throw("Block doesn't exist")
    A.blocks[s] = blk
    A
end

SymTensor(f::Function, charge::S, space::NTuple{N,VectorSpace{S}}) where {S,N} =
    SymTensor(f, Float64, charge, space)

function SymTensor(f::Function, ::Type{T}, charge::S,
                   space::NTuple{N,VectorSpace{S}}) where {S,T,N}
    sects, sizes = _allsectorsandsizes(charge, space)
    blocks = SortedDict([sects[i] => f(T,sizes[i]) for i in eachindex(sects)])
    SymTensor(charge, space, blocks)
end

function rand(::Type{T}, charge::S, space::NTuple{N, VectorSpace{S}};
              seed::Int=1911) where {S,T,N}
    sects, sizes = _allsectorsandsizes(charge, space)
    rng = MersenneTwister(seed)
    blocks = SortedDict([sects[i] => rand(rng, T, sizes[i]) for i in eachindex(sects)])
    SymTensor(charge, space, blocks)
end

rand(::Type{T}, space::NTuple{N, VectorSpace{S}}) where {S,T,N} =
    rand(T, zero(S), space)
rand(charge::Int, space::NTuple{N, VectorSpace{S}}) where {S,N} =
    rand(Float64, charge, space)
rand(space::NTuple{N, VectorSpace{S}}) where {S,N} =
    rand(Float64, zero(0), space)

fill(x::T, space::NTuple{N, VectorSpace{S}}) where {S,T,N} =
    fill(x, zero(S), space)
function fill(x::T, charge::S, space::NTuple{N, VectorSpace{S}}) where {S,T,N}
    sects, sizes = _allsectorsandsizes(charge, space)
    blocks = SortedDict([sects[i] => fill(x, sizes[i]) for i in eachindex(sects)])
    SymTensor(charge, space, blocks)
end
fill(x::T, space::NTuple{N, TrivialVectorSpace}) where{T, N} = fill(x, dim.(space))

function fill!(A::SymTensor{S,T,N}, x::T) where {S,T,N}
    for b in values(A.blocks)
        fill!(b, x)
    end
    A
end

function fill_linearindex(charge::S, space::NTuple{N, VectorSpace{S}}) where {S,N}
    sects, sizes = _allsectorsandsizes(charge, space)
    blocks = SortedDict{Sector{S, N}, Array{Int, N}}()
    p = 0
    for i in eachindex(sects)
        sz = sizes[i]
        blocks[sects[i]] = reshape(collect(p+1:p+prod(sz)), sz...)
        p += prod(sz)
    end
    SymTensor(charge, space, blocks)
end
fill_linearindex(space::NTuple{N, VectorSpace{S}}) where{S,N} =
    fill_linearindex(zero(S), space)

##TODO: Mix this with the usual UnitScaling
@inline eye(V::VectorSpace{S}) where{S} = eye(Float64, V)
@inline eye(::Type{T}, V::VectorSpace{S}) where {S, T} =
    SymMatrix(zero(S), (V, dual(V)),
              SortedDict([Sector(c, c)=>Matrix{T}(I,d,d) for (c,d) in V]))

"""
    dual(A)

Make the dual of tensor. This operation negates the: tensor charge,
leg signs, sectors as well! One can set whether or not it also
conjugates the blocks with `conjugate` which is true by default.

"""
function dual(A::AbstractSymTensor; conjugate::Bool=true)
    if !conjugate
        return typeof(A)(inv(A.charge),
                         dual.(A.space),
                         A.blocks)
        #SortedDict([inv(s)=>b for (s,b) in A.blocks]))
    end
    typeof(A)(inv(A.charge),
              dual.(A.space),
              A.blocks)
    #SortedDict([inv(s)=>conj(b) for (s,b) in A.blocks]))
end
dual(A::AbstractArray) = conj(A)

function mapcharges(f::Function, A::AbstractSymTensor)
    SymTensor(f(A.charge),
              mapcharges.(f, A.space),
              typeof(A.blocks)(mapcharge(f, s),d for (s,d) in A.blocks))
end

function mapcharges(f::NTuple{N, Function},
                    A::AbstractSymTensor{S,T,N}) where{S,T<:Number,N}
    space = Tuple(mapcharges(f[i], A.space[i]) for i in 1:N)

    sects = [Tuple(f[i](s[i]) for i in 1:N) for s in sectors(A)]
    charges = [sum(s, isdual.(space)) for s in sectors(A)]
    !all([charges[i] == charges[1] for i in eachindex(charges)]) &&
        error("mapcharges functions are inconsistent!")
    SymTensor(charges[1], space,
              typeof(A.blocks)(Sector(Tuple(convert(S,f[i](s[i])) for i in 1:N))=>d for (s,d) in A.blocks))
end

conj(A::AbstractSymTensor) =
    typeof(A)(A.charge, A.space, conj(A.blocks))

function array(A::AbstractSymTensor{S,T,N}) where {S,T,N}
    arrep = zeros(eltype(A), dim.(space(A))...)
    for (s,b) in A.blocks
        arrep[layout(A.space, s)...] = b
    end
    arrep
end

@inline *(A::AbstractSymTensor, a::T) where {T<:Number} =
    typeof(A)(A.charge, A.space,
              SortedDict([s=>a.*b for (s,b) in A.blocks]))

@inline *(a::T, A::AbstractSymTensor) where {T<:Number} = *(A, a)

"construct an additional SymTensor similar to A, possibly with a
different scalar type T."
function similar(A::AbstractSymTensor, T::Type=eltype(A))
    typeof(A)(A.charge, A.space,
              SortedDict([s=>similar(d, T) for (s,d) in A.blocks]))
end

"copy the contents of A to a preallocated SymTensor B"
function copyto!(B::T, A::T) where {T<:AbstractSymTensor}
    issimilar(B, A) && error("Can not copy to a  non-similar SymTensor!")
    for st in onlysemitokens(B.blocks)
        copyto!(B.blocks[st], A.blocks[st])
    end
    B
end

"out of place scalar multiplication; multiply SymTensor A with scalar α
and store the result in B"
function mul!(B::T, A::T, α) where {T<:AbstractSymTensor}
    issimilar(A, B) || error("oops")
    semitsA = collect(onlysemitokens(A.blocks))
    semitsB = collect(onlysemitokens(B.blocks))
    for st in onlysemitokens(B.blocks)
        mul!(B.blocks[st], A.blocks[st], α)
    end
    B
end

"in-place scalar multiplication of A with α; in particular with α =
false, A is initialized with all zeros"
function rmul!(A::T, α) where {T<:AbstractSymTensor}
    for st in onlysemitokens(A.blocks)
        rmul!(A.blocks[st],  α)
    end
    A
end

"Stores in B the result of α*A + B"
function axpy!(α,
               A::AbstractSymTensor{S,T1,N},
               B::AbstractSymTensor{S,T2,N}) where {S,T1<:Number, T2<:Number, N}
    issimilar(A, B) || error("axpy! The two matrices are not simliar!")
    for st in onlysemitokens(B.blocks)
        B.blocks[st] = α .* A.blocks[st] + B.blocks[st]
    end
    B
end

"Stores in B the result of α*A + β*B"
function axpby!(α,
                A::AbstractSymTensor{S,T1,N},
                β,
                B::AbstractSymTensor{S,T2,N}) where {S,T1<:Number, T2<:Number, N}
    issimilar(A, B) || error("axpy! The two matrices are not simliar!")
    for st in onlysemitokens(B.blocks)
        B.blocks[st] = α .* A.blocks[st] + β .* B.blocks[st]
    end
    B
end

"""
compute the inner product of two similar SymTensors which conjugates
the second one
"""
function dot(A::AbstractSymTensor{S,T1,N},
             B::AbstractSymTensor{S,T2,N}) where{S,T1<:Number, T2<:Number, N}
    issimilar(A, B) || error("The two SymTensors have to be similar to dot!")
    sum([dot(A.blocks[st], B.blocks[st]) for st in onlysemitokens(A.blocks)])
end

" compute the 2-norm of a  AbstractSymTensor"
function norm(A::AbstractSymTensor)
    sqrt(Float64(dot(A, A)))
end

normalize!(S::SymDiagonal) = rmul!(S, 1/norm(S))
normalize!(S::LinearAlgebra.Diagonal) = rmul!(S, 1/norm(S))

# function dropdummyleg(A::AbstractSymTensor, l::Int)
#     isdummy(A.space[l]) || error("space is not dummy!")

#     N = rank(A)
#     SymTensor(A.charge,
#               (A.space[1:l-1]..., A.space[l+1:N]...),
#               typeof(A.blocks)((s[1:l-1]...,s[l+1:N]...),
#     reshape(d, size(d)[1:l-1]...,size(d)[l+1:N]...) for (s,d) in A.blocks)
#               end
#               end

function Base.show(io::IO, A::AbstractSymTensor)
    show(io, typeof(A))
    num_ch = length(A.space)
    println(io)
    #println("$(size(A)) $(typeof(A)) with $(length(A.sects)) blocks: ")
    for (s,b) in A.blocks
        print(io, "$s with size: ")
        show(io, "text/plain", b)
        println(io)
    end
    nothing
end

function trimspace!(A::AbstractSymTensor, l::Int)
    N = rank(A)
    1 <= l <= N || error("not a leg index!")
    leg = A.legs[l]
    s = leg.sign
    legcharges = BitSet()
    for sect in A.sects
        push!(legcharges, s*sect[l])
    end
    indexes = _findindexes_sorted(leg.chrs, collect(legcharges))
    leg = U1Space(leg.sign, leg.chrs[indexes], leg.dims[indexes])
    A.legs = (A.legs[1:l-1]..., leg, A.legs[l+1:N]...)
    A
end

# function negateleg(A::AbstractSymTensor, l::Int)
#     N = rank(A)
#     0 < l <= N || error("leg $l doesn't exist!")
#     legs = (A.legs[1:l-1]..., negate(A.legs[l]), A.legs[l+1:N]...)
#     SymTensor{eltype(A),N}(A.charge, legs, sects, A.nzblks)
# end

function _findindexes_sorted(A::Vector{T}, elements::Vector{T}) where{T}
    #println(A, elements)
    p = 1
    indexes = Int[]
    for i in eachindex(A)
        if A[i] == elements[p]
            push!(indexes, i)
            if p==length(elements) break end
            p += 1
        elseif A[i] > elements[p]
            error("_findindexes_sorted ", A, elements)
        end
    end
    indexes
end

# """
#     SymTensor(signs, sector, nzblock)

# Make a SymTensor with the given possible charges and only one nonzero
# sector. Find the dimensions from that.

# """
# function SymTensor(signs::NTuple{N, Int},
#                    sector::NTuple{N, Int},
#                    nzblock::Array{T, N}) where {T<:Number, N}

#     chrs = [0, 1]
#     charge = sum(signs .* sector)
#     _sector_is_allowed(charge, signs, sector) ||
#         error("sector is not allowed ", sector)
#     legs = U1Space[]
#     legs = Tuple([U1Space(signs[n], [sector[n]], [size(nzblock, n)])
#                   for n in eachindex(sector)])
#     SymTensor(charge, Tuple(legs), [sector], [nzblock])
# end

# """
#     findindexes_sorted(A, elements

# Find the index of a set of sorted elements in a sorted vector `A`.
# """
# function _trimlegs(legs::NTuple{N, U1Space}, sects::Vector{NTuple{N, Int}}) where{N}
#     legcharges = Tuple(BitSet() for i=1:N)
#     for sect in sects
#         for i in 1:N
#             push!(legcharges[i], sect[i])
#         end
#     end
#     #print(legcharges)
#     trimmedlegs = U1Space[]
#     for lidx in 1:N
#         leg = legs[lidx]
#         indexes = findindexes_sorted(leg.chrs, collect(legcharges[lidx]))
#         push!(trimmedlegs, U1Space(leg.sign, leg.chrs[indexes], leg.dims[indexes]))
#     end
#     # println("trimming :")
#     # println(legs)
#     # println(trimmedlegs)
#     Tuple(trimmedlegs)
# end

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
abstract type CorrelationStructure end

# not mutable because it will never change
struct Independence <: CorrelationStructure
    N::Number
    α
end

Independence(N, α = 0) = Independence(N, α)

mutable struct Exchangeable <: CorrelationStructure
    N::Real
    α̂::Real
end

mutable struct Unstructered <: CorrelationStructure
    N::Real
    α̂
end

(X::Independence)() = Matrix{Int}(I, X.N, X.N)

(X::Unstructered)() = reshape(X.α̂, (X.N, X.N))

function (X::Exchangeable)()
    R = zeros(X.N, X.N)
    for j in 1:X.N
        for k in 1:X.N
            R[j, k] = j == k ? 1 : X.α̂
        end
    end
    return R
end

function updateα(::Type{Independence}, n, m, e, ϕ)
    Independence(m)()
end

function updateα(::Type{Exchangeable}, n, m, e, ϕ)
    R = updateα(Unstructered, n, m, e, ϕ)
    mean([R[r] for r in CartesianIndices(R) if r[1] ≠ r[2]])
end

# n is the number of clusters
# m is the number of observations within clusters
# e is the pearsons residuals within a cluster
# ϕ is well... ϕ
# returns the unstructured correlation matrix
function updateα(::Type{Unstructered}, n, m, e, ϕ)
    R = zeros(m, m)
    for i in 1:n
        for j in 1:m
            for k in 1:m
                R[j, k] += e[i][j]*e[i][k]
            end
        end
    end
    R = R / (n * ϕ)
    R[diagind(R)] .= 1
    return R
end


text:
include("../src/ConicHulls.jl")

include("test_det.jl")
include("test_isevenperm.jl")

include("test_basic.jl")
include("test_cases.jl")
include("test_random_adds.jl")


text:
# Julia wrapper for header: /usr/local/include/libavutil/display.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_display_rotation_get,
    av_display_rotation_set,
    av_display_matrix_flip

function av_display_rotation_get(matrix::NTuple{9,Int32})
    ccall((:av_display_rotation_get,libavutil),Cdouble,(NTuple{9,Int32},),matrix)
end

function av_display_rotation_set(matrix::NTuple{9,Int32},angle::Cdouble)
    ccall((:av_display_rotation_set,libavutil),Void,(NTuple{9,Int32},Cdouble),matrix,angle)
end

function av_display_matrix_flip(matrix::NTuple{9,Int32},hflip::Integer,vflip::Integer)
    ccall((:av_display_matrix_flip,libavutil),Void,(NTuple{9,Int32},Cint,Cint),matrix,hflip,vflip)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
#####################
###### Solving ######
#####################

const sosdebug = true

# NOTE: doesn't return the result of ex.
# This allows ex to contain return statements.
macro timeifdebug(str,ex)
    if sosdebug
        quote
            @printf("%s -- ", $(esc(str)))
            Libc.flush_cstdio()
            tic()
            $(esc(ex))
            toc()
            Libc.flush_cstdio()
        end
    else
        quote $(esc(ex)) end
    end
end

# The row-redundancy test works as follows: we choose a random 1-d projection,
# and keep a dictionary of existing rows sorted by their value in that 1-d
# projection. Each time we get a new row, we take its value in that projection
# and fully compare that row with the projection-nearby rows. This should result in
# very few superfluous full-row comparisons.
# In order to get the nearby rows as indexed by projection value, we use the SortedDict
# type from DataStructures.jl.
function row_redundancy_object(wid)
    randproj = rand(1,wid)
    rowsd = SortedDict(Dict{Float64,SparseMatrixCSC{Float64,Int64}}())
    (randproj, rowsd)
end

function row_redundancy_check(redobj, row)
    (randproj, rowsd) = redobj
    rval = (randproj * row)[1,1]

    st = searchsortedfirst(rowsd,rval)
    while status((rowsd,st)) == 1
        (k,v) = deref((rowsd,st))
        if k - rval > 1e-7
            break
        end
        if norm(v - row) < 1e-7
            return false
        end
        st = advance((rowsd,st))
    end

    st = searchsortedlast(rowsd,rval)
    while status((rowsd,st)) == 1
        (k,v) = deref((rowsd,st))
        if rval - k > 1e-7
            break
        end
        if norm(v - row) < 1e-7
            return false
        end
        st = regress((rowsd,st))
    end

    rowsd[rval] = row
    true
end

#  Make the magic happen.
function sossolve(prog :: Program, d :: Int64; solver="csdp", call="csdp")

    # sanity check
    if(d <= 0)
        throw(ArgumentError("Degree must be positive"))
    elseif(d < deg(prog.objective))
        throw(ArgumentError("Degree must be at least that of the objective"))
    elseif(d % 2 == 1)
        throw(ArgumentError("Degree must be even"))
    elseif(length(prog.vars) == 0)
        throw(ArgumentError("Program is empty"))
    end

    @timeifdebug "Enumerating monomials" begin
        # choose an order for indexing non-1 monomials
        monall = [ monoms(prog,i) for i in 0:d ]
        mon0d2 = vcat( monall[ 1:(div(d,2)+1) ]... ) # monomials of degree from 0 to d/2
    end

    @timeifdebug "Symmetrizing monomials" begin
        omap, revomap = monom_orbits(vcat(monall...), prog.symmetries)
    end

    @timeifdebug "Writing constraint-by-orbit matrix" begin
        # write the (sparse) constraint-by-orbit occurence matrix C, removing redundant constraints
        CI = Array{Int64,1}()
        CJ = Array{Int64,1}()
        CV = Array{Float64,1}()
        redund = row_redundancy_object(length(revomap))
        i = 1
        for poly in prog.constraints
            pd = deg(poly)
            if pd > d
                println("Warning: a constraint has degree $pd, larger than relaxation degree $d")
            end

            # promote to higher degree in all possible ways
            for promdeg in 0:(d-pd)
                for monom in monall[promdeg+1] # +1 because monall is an array indexed from 1

                    row = spzeros(length(revomap),1)
                    for (k,v) in poly
                        orbit = omap[k*monom]
                        row[orbit,1] += v
                    end

                    # add this row, but only if it's actually new
                    if row_redundancy_check(redund, row)
                        rv = rowvals(row)
                        append!(CI,fill(i,length(rv)))
                        append!(CJ,rv)
                        append!(CV,nonzeros(row))
                        i += 1
                    end

                end
            end
        end

        # we have to specify width `length(revomap)` in case higher moments /
        # orbits don't show up in the promoted constraints.
        # XXX In principle, though, it'd make the linear algebra faster to
        # exclude these columns, and manually re-include them after the linear
        # algebra. At some point we should do this; we can also reduce columns
        # by joining orbits by equality constraints (e.g. boolean constraints)
        C = sparse(CI,CJ,CV,i-1,length(revomap))
    end
    @printf("Size of C is %s\n", size(C))

    @timeifdebug "Manipulating data structures" begin
        # write our objective in moments
        O = zeros(length(revomap))
        for (k,v) in prog.objective
            O[omap[k]] += v
        end

        # separate the constant column from the rest
        Cconst = C[:,1]
        C = C[:,2:end]
        Oconst = O[1]
        O = O[2:end]
        revomap = revomap[2:end]

        # decomposition-to-orbit map
        domap = Dict{Tuple{Int64,Int64},Int64}()
        for a in 1:length(mon0d2)
            for b in max(a,2):length(mon0d2)
                domap[(a,b)] = omap[mon0d2[a]*mon0d2[b]]-1
            end
        end

        omap = 0 # done with this
    end

    # find an initial solution for moments -- satisfying affine constraints, but not necessarily PSD
    @timeifdebug "Computing initial value" if size(C,1) > 0
        initial = vec(C \ full(Cconst))
        # test that this actually satisfies constraints. if not, declare infeasibility
        if norm(Cconst - C * initial) > 1e-6
            return sos_sol_infeasible(prog,d)
        end
    else
        initial = zeros(size(C,2))
    end

    # compute the nullspace
    # XXX really want to do this in a sparse way, from the QR decomposition of B
    @timeifdebug "Computing ideal basis" begin
        B = nullspace(full(C))
    end

    # we will need to adjust our objective by the following constant term, which is missing from the SDP
    adjust_obj = Oconst - dot(O, initial)

    # set up problem & solver instance
    sdp = SDPSession(prog.maximize, size(B,2), length(mon0d2), adjust_obj)

    # write objective value
    @timeifdebug "Writing objective -- " begin
        for i in 1:size(B,2)
            sdp_obj!(sdp, dot(O, B[:,i]))
        end
    end

    # dual objective = primal constant-term
    @timeifdebug "Writing initial value -- " begin
        sdp_con!(sdp, 0, 1, 1, -1.0)
        for ((a,b),orbit) in domap
            val = initial[orbit]
            if(abs(val) > 1e-8)
                sdp_con!(sdp, 0, a, b, val)
            end
        end
    end

    # dual constraints = primal building-blocks
    @timeifdebug "Writing ideal basis -- " begin
        for i in 1:size(B,2)
            # rewrite B[:,i] as a matrix in d/2 x d/2 moment decompositions
            for ((a,b),orbit) in domap
                val = B[orbit,i]
                if(abs(val) > 1e-8)
                    sdp_con!(sdp, i, a, b, val)
                end
            end
        end
    end

    # free some memory
    C = Cconst = O = Oconst = B = initial = monall = mond = mondrev = revomap = domap = 0

    # solve the SDP
    println("Solving SDP with $(sdp.nmoments)^2 entries and $(sdp.nconstraints) constraints")
    @time sdpsol = sdp_solve(sdp,call=call)

    # build & return a solution object
    @timeifdebug "Building solution object" begin
        if sdpsol.status == :Infeasible
            return sos_sol_unbounded(prog,d)
        elseif sdpsol.status == :Unbounded
            return sos_sol_infeasible(prog,d)
        end
        moments = Dict{SoSMonom,Float64}()
        for i in 1:length(mon0d2)
            for j in i:length(mon0d2)
                moments[mon0d2[i]*mon0d2[j]] = sdpsol.primalmatrix[i,j]
            end
        end
        status = sdpsol.status == :Infeasible ? :Unbounded : sdpsol.status == :Unbounded ? :Infeasible : sdpsol.status
        return SoSSolution(prog, d, sdpsol.primalobj, sdpsol.dualobj, moments, sdpsol.primalmatrix, sdpsol.dualmatrix, status)
    end
end

rowdist(r1,r2) = halfrowdist(r1,r2) + halfrowdist(r2,r1)
function halfrowdist(r1 :: Dict{Int64,Float64}, r2 :: Dict{Int64,Float64})
    dist = 0.0
    for (k,v) in r1
        dist += (v - get(r2,k,0.0))^2
    end
    return dist
end

max_stars_repo_name: willperry/SumOfSquaresOptimization


text:
# # Using with FluxTraining.jl
#
# *DLPipelines.jl* works great together with
# [*FluxTraining.jl*](https://github.com/lorenzoh/FluxTraining.jl).
#
# We'll train an image classifier to show how.
#
# If you want to follow along, you'll have to install some other libraries

using Pkg
Pkg.add("DataLoaders")
Pkg.add("Flux")
Pkg.add("LearnBase")
Pkg.add("MLDataPattern")
Pkg.add(url="https://github.com/lorenzoh/DLDatasets.jl")
Pkg.add(url="https://github.com/lorenzoh/FluxTraining.jl")
Pkg.add(url="https://github.com/lorenzoh/ModelUtils.jl")
Pkg.add(url="https://github.com/lorenzoh/FluxModels.jl")

# and import the following:
# {cell=main}

using DataLoaders: DataLoader
using DLPipelines
using DLDatasets: ImageNette, loaddataset, metadata
using Flux: ADAM, logitcrossentropy
using FluxModels: xresnet18
using FluxTraining: Learner, fit!
using LearnBase: getobs
using MLDataPattern: splitobs

# Recall the 4 things needed to start training with *FluxTraining.jl*:
#
# - a model to train,
# - a loss function,
# - an optimizer; and
# - training and validation data iterators
#
# For the first three we use:
#
# {cell=main, result=false}

lossfn = logitcrossentropy
optimizer = ADAM()
model = xresnet18()

# The last part is where *DLPipelines.jl* comes in.
#
# As in the [`ImageClassification` overview](../methods/imageclassification.jl),
# we first need to load an image classification dataset:
#
# {cell=main, result = false}

dataset = loaddataset(ImageNette, "v2_160px")
categories = metadata(ImageNette).labels

# Now we construct an [`ImageClassification`](#) object with
# the proper configuration:
#
# {cell=main}

method = ImageClassification(categories, sz = (128, 128))

# After splitting the dataset into training and validation split, we can use
# [`MethodDataset`](#) to get data containers that will encode the
# observations.
#
# {cell=main}

traindataset_, validdataset_ = splitobs(dataset, 0.8)
traindataset = MethodDataset(traindataset_, method, Training())
validdataset = MethodDataset(validdataset_, method, Validation())

# By passing either [`Training`](#) or [`Validation`](#) as [`Context`](#)
# to the datasets, the correct transformations will be automatically applied.
# For example, a random crop will be used for the `traindataset`, while
# images in `validdataset` will always be cropped from the center.
#
# Loading an observation from either will now correctly return a normalized
# 3D-array with the image data and a one-hot encoded category:
#
# {cell=main}

x, y = getobs(traindataset, 1)
summary.((x, y))

# To finally get an iterator over batches, we can use `DataLoader` from
# [DataLoaders.jl](https://lorenzoh.github.io/DataLoaders.jl/dev/README.html).
# It will also make sure to prefetch the data on background threads so that
# the training loop isn't slowed down waiting for the next batch of data.
#
# {cell=main, result=false}

batchsize = 16
traindl = DataLoader(traindataset, batchsize)
validdl = DataLoader(validdataset, batchsize)

# Now we can create a `Learner`:
#
# {cell=main}

learner = Learner(model, (traindl, validdl), optimizer, lossfn)

# You can pass any number of
# [callbacks](https://lorenzoh.github.io/FluxTraining.jl/dev/docs/callbacks/reference.html)
# to `Learner`, for example `ToGPU()` to utilize your GPU when training.
#
# And that's it, you can start training with a call to `fit!`:

fit!(learner, 10)

# To find out more about what you can do with *FluxTraining.jl*, see
# [its documentation](https://lorenzoh.github.io/FluxTraining.jl/dev/README.html).


text:
# TODO: set up our own tfuncs

function CC.builtin_tfunction(interp::JETInterpreter, @nospecialize(f), argtypes::Array{Any,1},
                              sv::InferenceState) # `JETInterpreter` isn't overloaded on `return_type`
    ret = @invoke builtin_tfunction(interp::AbstractInterpreter, f, argtypes::Array{Any,1},
                                    sv::Union{InferenceState,Nothing})

    if f === throw
        # here we only report a selection of "serious" exceptions, i.e. those should be
        # reported even if they may be caught in actual execution;
        # other general `throw` calls will be reported within `_typeinf(interp::JETInterpreter, frame::InferenceState)`
        # only when they are not caught by control flow, which is judged by whether if the
        # final return type of `sv` is annotated as `Bottom` or not
        if length(argtypes) ≥ 1
            a = first(argtypes)
            if isa(a, Const)
                v = a.val
                if isa(v, UndefKeywordError)
                    report!(interp, UndefKeywordErrorReport(interp, sv, v, get_cur_linfo(sv)))
                end
            end
        end

        return ret
    elseif f === getfield
        # getfield is so common, let's special case it
        if 2 ≤ length(argtypes) ≤ 3
            obj, fld = argtypes
            if isa(fld, Const)
                name = fld.val
                if isa(name, Symbol)
                    if isa(obj, Const) && (mod = obj.val; isa(mod, Module))
                        if isdefined(mod, name)
                            # TODO: add report pass here (for performance linting)
                        else
                            # report access to undefined global variable
                            report!(interp, GlobalUndefVarErrorReport(interp, sv, mod, name))
                            # return Bottom
                        end
                    elseif ret === Bottom
                        # general case when an error is detected by the native `getfield_tfunc`
                        typ = widenconst(obj)
                        report!(interp, NoFieldErrorReport(interp, sv, typ, name))
                        return ret
                    end
                end
            end
        end
    elseif f === fieldtype
        # the valid widest possible return type of `fieldtype_tfunc` is `Union{Type,TypeVar}`
        # because fields of unwrapped `DataType`s can legally be `TypeVar`s,
        # but this will cause lots of false positive `NoMethodErrorReport` reports for
        # inference with accessing to fields of abstract types since most methods don't
        # expect `TypeVar` (e.g. `@report_call readuntil(stdin, 'c')`);
        # JET.jl further widens this case to `Any` and give up further analysis rather than
        # noisy and sound analysis
        # xref: https://github.com/JuliaLang/julia/pull/38148
        if ret === Union{Type, TypeVar}
            return Any
        end
    end

    if isa(ret, VirtualGlobalVariable)
        # propagate virtual global variable

        add_backedge!(ret.li, sv)
        # this might be `Bottom`, but hopefully the error on this variable is already reported,
        # so we don't check `InvalidBuiltinCallErrorReport` for this pass
        return ret.t
    elseif ret === Bottom
        # XXX: for general case, JET just relies on the (maybe too persmissive) return type
        # from native tfuncs to report invalid builtin calls and probably there're lots of
        # false negatives
        report!(interp, InvalidBuiltinCallErrorReport(interp, sv, argtypes))
    end

    return ret
end

# `return_type_tfunc` internally uses `abstract_call` to model `return_type` function and
# here we shouldn't pass `JETInterpreter` to it; otherwise we may get false error reports from
# the  `abstract_call`, which isn't the abstraction of actual execution, thus here we just
# check if the call of `return_type` is valid or not
function CC.return_type_tfunc(interp::JETInterpreter, argtypes::Vector{Any}, sv::InferenceState)
    if length(argtypes) !== 3
        # invalid argument number, let's report and return error result (i.e. `Bottom`)
        report!(interp, NoMethodErrorReport(interp,
                                            sv,
                                            false,
                                            # this is not necessary to be computed correctly, though
                                            argtypes_to_type(argtypes),
                                            ))
        return Bottom
    else
        # don't recursively pass on `JETInterpreter` via `@invoke_native`
        return return_type_tfunc(interp.native, argtypes, sv)
    end
end


text:
using Test
using TensorNetworkAD
using TensorNetworkAD: diaglocalhamiltonian, energy, expectationvalue, optimiseipeps,
                       hamiltonian, indexperm_symmetrize, num_grad
using OMEinsum, Zygote, Random
using LinearAlgebra: svd, norm
using Optim, LineSearches

@testset "non-interacting" begin
    h = diaglocalhamiltonian([1,-1.0])
    as = (rand(3,3,3,3,2) for _ in 1:100)
    @test all(a -> -1 < energy(h,SquareIPEPS(a); χ=5, tol=0, maxit=10)/2 < 1, as)

    h = diaglocalhamiltonian([1,-1.0])
    a = zeros(2,2,2,2,2) .+ 1e-12 * randn(2,2,2,2,2)
    a[1,1,1,1,2] = randn()
    @test energy(h,SquareIPEPS(a); χ=4, tol=1e-12, maxit=100)/2 ≈ -1

    a = zeros(2,2,2,2,2) .+ 1e-12 * randn(2,2,2,2,2)
    a[1,1,1,1,1] = randn()
    @test energy(h,SquareIPEPS(a); χ=10, tol=0, maxit=300)/2 ≈ 1

    a = zeros(2,2,2,2,2) .+ 1e-12 * randn(2,2,2,2,2)
    a[1,1,1,1,2] = a[1,1,1,1,1] = randn()
    @test abs(energy(h,SquareIPEPS(a); χ=10, tol=0, maxit=300)) < 1e-9

    grad = let energy = x -> real(energy(h, SquareIPEPS(x); χ=8, tol=0, maxit=3))
        res = optimize(energy,
            Δ -> Zygote.gradient(energy,Δ)[1], a, LBFGS(m=20), inplace = false)
    end
    @test grad != Nothing

    hdiag = [0.3,0.1,-0.43]
    h = diaglocalhamiltonian(hdiag)
    a = SquareIPEPS(randn(2,2,2,2,3))
    res = optimiseipeps(a, h; χ=4, tol=0, maxit=100,
        optimargs = (Optim.Options(f_tol=1e-6, show_trace=false),))
    e = minimum(res)/2
    @test isapprox(e, minimum(hdiag), atol=1e-3)
end

@testset "ising" begin
    Random.seed!(3)
    h = zeros(2,2,2,2)
    h[1,1,2,2] = h[2,2,1,1] = 1
    h[2,2,2,2] = h[1,1,1,1] = -1
    ipeps = SquareIPEPS(randn(2,2,2,2,2))
    a = indexperm_symmetrize(ipeps)
    res = optimiseipeps(a, h; χ=4, tol=0, maxit=100,
        optimargs = (Optim.Options(f_tol=1e-6, show_trace=false),))
    e = minimum(res)
    @test isapprox(e,-1, atol=1e-3)

    h = zeros(2,2,2,2)
    h[1,1,2,2] = h[2,2,1,1] = 1
    h[2,2,2,2] = h[1,1,1,1] = -1
    randu, s,  = svd(randn(2,2))
    h = ein"(((abcd,ai),bj),ck),dl -> ijkl"(h,randu,randu',randu,randu')
    ipeps = SquareIPEPS(randn(2,2,2,2,2))
    a = indexperm_symmetrize(ipeps)
    res = optimiseipeps(a, h; χ=6, tol=0, maxit=200,
        optimargs = (Optim.Options(f_tol=1e-6, show_trace=false),))
    e = minimum(res)
    @test isapprox(e,-1, atol=1e-3)

    # comparison with results from https://github.com/wangleiphy/tensorgrad
    h = hamiltonian(TFIsing(1.0))
    ipeps = SquareIPEPS(randn(2,2,2,2,2))
    a = indexperm_symmetrize(ipeps)
    res = optimiseipeps(a, h; χ=5, tol=0, maxit=100,
        optimargs = (Optim.Options(f_tol=1e-6, show_trace=false),))
    e = minimum(res)
    @test isapprox(e, -2.12566, atol = 1e-3)

    h = hamiltonian(TFIsing(0.5))
    ipeps = SquareIPEPS(randn(2,2,2,2,2))
    a = indexperm_symmetrize(ipeps)
    res = optimiseipeps(a, h; χ=5, tol=0, maxit=100,
        optimargs = (Optim.Options(f_tol=1e-6, show_trace=false),))
    e = minimum(res)
    @test isapprox(e, -2.0312, atol = 1e-2)

    Random.seed!(1)
    h = hamiltonian(TFIsing(2.0))
    ipeps = SquareIPEPS(randn(2,2,2,2,2))
    a = indexperm_symmetrize(ipeps)
    res = optimiseipeps(a, h; χ=6, tol=1e-9, maxit=100,
        optimargs = (Optim.Options(f_tol=1e-8, show_trace=false),))
    e = minimum(res)
    @test isapprox(e, -2.5113, atol = 1e-3)
end

@testset "heisenberg" begin
    # comparison with results from https://github.com/wangleiphy/tensorgrad
    Random.seed!(2)
    h = hamiltonian(Heisenberg())
    ipeps = SquareIPEPS(randn(2,2,2,2,2))
    a = indexperm_symmetrize(ipeps)
    res = optimiseipeps(a, h; χ=5, tol=0, maxit=100,
        optimargs = (Optim.Options(f_tol=1e-6, show_trace=false),))
    e = minimum(res)
    @test isapprox(e, -0.66023, atol = 1e-3)

    h = hamiltonian(Heisenberg(2.0, 2.0, 1.0))
    ipeps = SquareIPEPS(randn(2,2,2,2,2))
    a = indexperm_symmetrize(ipeps)
    res = optimiseipeps(a, h; χ=6, tol=0, maxit=100, #optimmethod = Optim.LBFGS(),
        optimargs = (Optim.Options(f_tol = 1e-6, show_trace = false),))
    e = minimum(res)
    @test isapprox(e, -1.190, atol = 1e-2)

    h = hamiltonian(Heisenberg(0.5, 0.5, 2.0))
    ipeps = SquareIPEPS(randn(2,2,2,2,2))
    a = indexperm_symmetrize(ipeps)
    res = optimiseipeps(a, h; χ=5, tol=0, maxit=100,
        optimargs = (Optim.Options(f_tol = 1e-6, show_trace = false),))
    e = minimum(res)
    @test isapprox(e, -1.0208, atol = 1e-3)
end

@testset "gradient" begin
    Random.seed!(0)
    h = hamiltonian(Heisenberg())
    ipeps = SquareIPEPS(randn(2,2,2,2,2))
    a = indexperm_symmetrize(ipeps)
    gradzygote = first(Zygote.gradient(a) do x
        energy(h,x; χ=4, tol=0, maxit=100)
    end).bulk
    gradnum = num_grad(a.bulk, δ=1e-3) do x
        energy(h, SquareIPEPS(x); χ=4, tol=0, maxit=100)
    end

    @test isapprox(gradzygote, gradnum, atol=1e-3)
end

@testset "complex" begin
    Random.seed!(2)
    h = hamiltonian(Heisenberg())
    ipeps = SquareIPEPS(randn(2,2,2,2,2))
    a = indexperm_symmetrize(ipeps)
    ca = SquareIPEPS(a.bulk .+ 0im)
    @test energy(h,a; χ=4, tol=1e-12, maxit=100) ≈ energy(h,ca; χ=4, tol=1e-12, maxit=100)
    ϕ = exp(1im * rand()* 2π)
    ca = SquareIPEPS(a.bulk .* ϕ)
    @test energy(h,ca; χ=4, tol=1e-12, maxit=100) ≈ energy(h,a; χ=4, tol=1e-12, maxit=100)

    gradzygote = first(Zygote.gradient(a) do x
        real(energy(h,x; χ=4, tol=1e-12,maxit=100))
    end)

    ca = SquareIPEPS(a.bulk .+ 0im)
    @test gradzygote.bulk ≈ first(Zygote.gradient(ca) do x
        real(energy(h,x; χ=4, tol=1e-12,maxit=100))
    end).bulk

    Random.seed!(2)
    # real
    h = hamiltonian(Heisenberg())
    ipeps = SquareIPEPS(randn(2,2,2,2,2))
    a = indexperm_symmetrize(ipeps)
    res1 = optimiseipeps(a, h; χ=20, tol=1e-12, maxit=100,
        optimargs = (Optim.Options(f_tol=1e-6, store_trace = true, show_trace=false),));

    # complex
    ipeps = SquareIPEPS(randn(2,2,2,2,2) .+ randn(2,2,2,2,2) .* 1im)
    a = indexperm_symmetrize(ipeps)
    res2 = optimiseipeps(a, h; χ=20, tol=1e-12, maxit=100,
        optimargs = (Optim.Options(f_tol=1e-6,store_trace = true,  show_trace=false, allow_f_increases=true),),
        optimmethod = Optim.LBFGS(
            m = 10,
            alphaguess = LineSearches.InitialStatic(alpha=1, scaled=true),
            linesearch = LineSearches.Static())
        );
    @test isapprox(minimum(res1), minimum(res2), atol = 1e-3)
end


text:
using NNlib, Flux, LinearAlgebra, DiffEqBase

function createUser(req::HTTP.Request)
  data = JSON2.read(IOBuffer(HTTP.payload(req)))
  r = findall(x -> isequal(x.email, data.email), USERS)
  if !isempty(r)
    println("log: cannot create user $(data.email)")
    return HTTP.Response(400, "{ \"message\": \"email already used\" }")
  end
  user = User(getNextUserId(), data.name, data.password, data.email, 0)
  USERS[user.id] = user
  println("log: create user $(user.email)")
  return HTTP.Response(200, JSON2.write(user))
end

function getUser(req::HTTP.Request)
    data = JSON2.read(IOBuffer(HTTP.payload(req)))
    token = data.token
    r = findall(x -> (x.token == token), USERS)
    if isempty(r)
      return HTTP.Response(400, "{ \"message\": \"user not found\" }")
    else
      println("log: get user $(USERS[r[1]].email)")
      return HTTP.Response(200, JSON2.write(USERS[r[1]]))
    end
end

function login(req::HTTP.Request)
  input = JSON2.read(IOBuffer(HTTP.payload(req)))
  r = findall(x -> (x.email == input.email && x.password == input.password), USERS)
  if isempty(r)
    return HTTP.Response(400, "{ \"message\": \"email not found\" }")
  else
    token = rand(UInt32)
    USERS[r[1]].token = token
    return HTTP.Response(200, "{ \"token\": $(token) }")
  end
end

function get_day_offset(epoch)
  d = unix2datetime(epoch)
  return day(d) - 4 # 4 means 4 june
end

function add_schedule(req::HTTP.Request)
  data = JSON2.read(IOBuffer(HTTP.payload(req)))
  token, starttime, endtime, state, district = data.token, data.starttime, data.endtime, data.state, data.district
  r = findall(x -> (x.token == token), USERS)
  if isempty(r)
    return HTTP.Response(400, "{ \"message\": \"user not found\" }")
  else
    user  = USERS[r[1]]
    pred = prediction(state, district, get_day_offset(starttime))
    if pred isa Nothing
      sched = Schedule(getNextScheduleId(), user.id, starttime, endtime, state * "_" * district, 0, Int(round(time())))
    else
      sched = Schedule(getNextScheduleId(), user.id, starttime, endtime, state * "_" * district, pred[1][1], Int(round(time())))
    end
    SCHEDULES[sched.id] = sched
    return HTTP.Response(200, JSON2.write(sched))
  end
end

function get_shedules(req::HTTP.Request)
  data = JSON2.read(IOBuffer(HTTP.payload(req)))
  token = data.token
  r = findall(x -> (x.token == token), USERS)
  if isempty(r)
    return HTTP.Response(400, "{ \"message\": \"user not found\" }")
  else
    user  = USERS[r[1]]
    idxs = findall(x -> (x.uid == user.id), SCHEDULES)
    out = Array{Schedule}(undef, 0)
    for i in idxs
      push!(out, SCHEDULES[i])
    end
    return HTTP.Response(200, JSON2.write(out))
  end
end

function get_status(req::HTTP.Request)
  data = JSON2.read(IOBuffer(HTTP.payload(req)))
  district, state = data.district, data.state
  vector = Vector(collect(1:7))
  pred = Covid19Modelling.prediction(state, district, vector)
  if pred isa Nothing
    infected, rd = zeros(7), zeros(7)
  else
    infected, rd = pred
  end
  out = Dict("infected" => infected .+ rd, "rd" => rd)
  return HTTP.Response(200, [Pair("Access-Control-Allow-Origin", "*")]; body=JSON2.write(out))
end

function index(req::HTTP.Request)
  file = "Server/app" * HTTP.unescapeuri(req.target[1:end])
  return isfile(file) ? HTTP.Response(200, read(file)) : HTTP.Response(404)
end

HTTP.@register(ROUTER, "POST",   "/api/v1/register", createUser)
HTTP.@register(ROUTER, "POST",   "/api/v1/login", login)
HTTP.@register(ROUTER, "POST",   "/api/v1/user", getUser)
HTTP.@register(ROUTER, "POST",   "/api/v1/schedule/add", add_schedule)
HTTP.@register(ROUTER, "POST",   "/api/v1/schedule/list", get_shedules)
HTTP.@register(ROUTER, "POST",   "/api/v1/status", get_status)
HTTP.@register(ROUTER, "GET",    "/*", index)


text:
import StatsFuns: log2π

@node MvNormalMeanPrecision Stochastic [ out, (μ, aliases = [ mean ]), (Λ, aliases = [ invcov, precision ]) ]

conjugate_type(::Type{ <: MvNormalMeanPrecision }, ::Type{ Val{ :out } }) = MvNormalMeanPrecision
conjugate_type(::Type{ <: MvNormalMeanPrecision }, ::Type{ Val{ :μ } })   = MvNormalMeanPrecision
conjugate_type(::Type{ <: MvNormalMeanPrecision }, ::Type{ Val{ :Λ } })   = Wishart

# default method for mean-field assumption
@average_energy MvNormalMeanPrecision (q_out::Any, q_μ::Any, q_Λ::Any) = begin
    # naive: (ndims(q_out) * log2π - mean(logdet, q_Λ) + tr(mean(q_Λ)*(v_out + v_mean + (m_out - m_mean)*(m_out - m_mean)'))) / 2
    dim = ndims(q_out)

    m_mean, v_mean = mean_cov(q_μ)
    m_out, v_out   = mean_cov(q_out)
    m_Λ            = mean(q_Λ)

    result = zero(promote_type(eltype(q_out), eltype(q_μ), eltype(q_Λ)))
    result += dim * log2π
    result -= mean(logdet, q_Λ)
    @inbounds for k1 in 1:dim, k2 in 1:dim
        # optimize trace operation (indices can be interchanges because of symmetry)
        result += m_Λ[k1,k2] * (v_out[k1,k2] + v_mean[k1,k2] + (m_out[k2] - m_mean[k2]) * (m_out[k1] - m_mean[k1]))
    end
    result /= 2

    return result
end

# specialized method for mean-field assumption with q_Λ::Wishart
@average_energy MvNormalMeanPrecision (q_out::Any, q_μ::Any, q_Λ::Wishart) = begin
    # m_out, v_out = mean_cov(q_out)
    # m_mean, v_mean = mean_cov(q_μ)
    # return (ndims(q_out) * log2π - mean(logdet, q_Λ) + tr(mean(q_Λ)*(v_out + v_mean + (m_out - m_mean)*(m_out - m_mean)'))) / 2
    dim = ndims(q_out)

    m_mean, v_mean = mean_cov(q_μ)
    m_out, v_out   = mean_cov(q_out)
    df_Λ, S_Λ      = params(q_Λ)  # prevent allocation of mean matrix

    result = zero(promote_type(eltype(q_out), eltype(q_μ), eltype(q_Λ)))

    @inbounds for k1 in 1:dim, k2 in 1:dim
        # optimize trace operation (indices can be interchanges because of symmetry)
        result += S_Λ[k1,k2] * (v_out[k1,k2] + v_mean[k1,k2] + (m_out[k2] - m_mean[k2]) * (m_out[k1] - m_mean[k1]))
    end
    result *= df_Λ
    result += dim * log2π
    result -= mean(logdet, q_Λ)
    result /= 2

    return result
end

# default method for structured mean-field assumption
@average_energy MvNormalMeanPrecision (q_out_μ::Any, q_Λ::Any) = begin
    # naive: @views (d*log2π - mean(logdet, q_Λ) + tr(mean(q_Λ)*( V[1:d,1:d] - V[1:d,d+1:end] - V[d+1:end,1:d] + V[d+1:end,d+1:end] + (m[1:d] - m[d+1:end])*(m[1:d] - m[d+1:end])' ))) / 2
    dim = div(ndims(q_out_μ), 2)

    m, V = mean_cov(q_out_μ)
    m_Λ  = mean(q_Λ)

    result = zero(promote_type(eltype(q_out_μ), eltype(q_Λ)))
    result += dim * log2π
    result -= mean(logdet, q_Λ)
    @inbounds for k1 in 1:dim, k2 in 1:dim
        # optimize trace operation (indices can be interchanges because of symmetry)
        result += m_Λ[k1,k2] * (V[k1,k2] + V[dim+k1,dim+k2] - V[dim+k1,k2] - V[k1,dim+k2] + (m[k1] - m[dim+k1]) * (m[k2] - m[dim+k2]))
    end
    result /= 2

    return result
end

# specialized method for structured mean-field assumption with q_Λ::Wishart
@average_energy MvNormalMeanPrecision (q_out_μ::Any, q_Λ::Wishart) = begin
    # naive: (ndims(q_out) * log2π - mean(logdet, q_Λ) + tr(mean(q_Λ)*(v_out + v_mean + (m_out - m_mean)*(m_out - m_mean)'))) / 2
    dim = div(ndims(q_out_μ), 2)

    m, V       = mean_cov(q_out_μ)
    df_Λ, S_Λ  = params(q_Λ)     # prevent allocation of mean matrix

    result = zero(promote_type(eltype(q_out_μ), eltype(q_Λ)))

    @inbounds for k1 in 1:dim, k2 in 1:dim
        # optimize trace operation (indices can be interchanges because of symmetry)
        result += S_Λ[k1,k2] * (V[k1,k2] + V[dim+k1,dim+k2] - V[dim+k1,k2] - V[k1,dim+k2] + (m[k1] - m[dim+k1]) * (m[k2] - m[dim+k2]))
    end
    result *= df_Λ
    result += dim * log2π
    result -= mean(logdet, q_Λ)
    result /= 2

    return result
end


text:
abstract type AbstractIFParameter end
@snn_kw struct IFParameter{FT=Float32} <: AbstractIFParameter
    τm::FT = 20ms
    τe::FT = 5ms
    τi::FT = 10ms
    Vt::FT = -50mV
    Vr::FT = -60mV
    El::FT = Vr
end

abstract type AbstractIF end

@snn_kw mutable struct IF{VFT=Vector{Float32},VBT=Vector{Bool}} <: AbstractIF
    param::IFParameter = IFParameter()
    N::Int32 = 100
    v::VFT = param.Vr .+ rand(N) .* (param.Vt - param.Vr)
    ge::VFT = zeros(N)
    gi::VFT = zeros(N)
    fire::VBT = zeros(Bool, N)
    I::VFT = zeros(N)
    records::Dict = Dict()
end

"""
    [Integrate-And-Fire Neuron](https://neuronaldynamics.epfl.ch/online/Ch1.S3.html)
"""
IF

function integrate!(p::IF, param::IFParameter, dt::Float32)
    @unpack N, v, ge, gi, fire, I = p
    @unpack τm, τe, τi, Vt, Vr, El = param
    @inbounds for i = 1:N
        v[i] += dt * (ge[i] + gi[i] - (v[i] - El) + I[i]) / τm
        ge[i] += dt * -ge[i] / τe
        gi[i] += dt * -gi[i] / τi
    end
    @inbounds for i = 1:N
        fire[i] = v[i] > Vt
        v[i] = ifelse(fire[i], Vr, v[i])
    end
end


text:
n1 = BayesianNetwork([], [])
n1 = BayesianNetwork()
pd1 = ProbabilityDistribution(["head","tails"],[0.5,0.5])
pd2 = ProbabilityDistribution(["head","tails"],[0.3,0.7])
a1 = DBayesianNode(:hulu, pd1)
a2 = DBayesianNode(:bulu, pd2)

@test a1.label == :hulu

@test length(n1.binclist) == 0
@test length(n1.finclist) == 0

ie1 = in_edges(a1, n1)
@test length(ie1) == 0

@test in_degree(a1, n1) == 0

oe1 = out_edges(a1, n1)
@test length(oe1) == 0

@test out_degree(a1, n1) == 0

add_node!(n1, a1)

@test has_cached_result(n1,P(:hulu)) == true
@test length(n1.binclist) == 1
@test length(n1.finclist) == 1

add_node!(n1, a2)

@test length(n1.binclist) == 2
@test length(n1.finclist) == 2

@test a1.index == 1
@test a2.index == 2
@test length(n1.nodes) == 2

e1 = add_edge!(n1, a1, a2)

@test source(e1).label == :hulu
@test edge_index(e1) == 1
@test length(n1.edges) == 1
@test n1.edges[1] == e1
@test length(in_neighbors(a2,n1)) == 1
@test length(out_neighbors(a2,n1)) == 0
@test length(out_neighbors(a1,n1)) == 1
@test out_neighbors(a1,n1)[1].label == :bulu
@test in_degree(a2,n1) != out_degree(a2,n1)

@test symbols_in_network(n1, [:hulu, :bulu]) == true
@test symbols_in_network(n1, [:hulu, :bjarke]) == false

@test add_probability!(n1, P(:hulu|:bulu), ProbabilityDensityDistribution(states(a1), [x -> x+1, x -> x+2])) == true
@test add_probability!(n1, P(:hulu|:bjarke), ProbabilityDensityDistribution(states(a1), [x -> x+1, x -> x+2])) == false
@test add_probability!(n1, P(:bjarke|:hulu), ProbabilityDensityDistribution(states(a1), [x -> x+1, x -> x+2])) == false

pd3 = ProbabilityDistribution(["head","tails"],[0.25,0.75])
pd4 = ProbabilityDistribution(["head","tails"],[0.75,0.25])
b1 = DBayesianNode(:bjarke, pd3)
b2 = DBayesianNode(:hesthaven, pd4)
b3 = CBayesianNode(:esben, x->x^0.5)

n2 = BayesianNetwork([b1,b2,b3],[])

@test b1.index == 1

pd5 = ProbabilityDistribution(["head","tails"],[0.46,0.54])
pd6 = ProbabilityDistribution(["head","tails"],[0.99,0.01])
c1 = DBayesianNode(:thomas, pd5)
c2 = DBayesianNode(:panum, pd6)
c3 = CBayesianNode(:moller, x->x^0.5)

@test_approx_eq_eps probability(c3,5) 2.236 0.0001

@test_throws CBayesianNode(:moller, x-> "I'm not entirely sure what I'm doing")

n3 = BayesianNetwork([c1,c2],[])

add_node!(n3, c3)
@test_throws add_edge!(n3,c3,c2)

e1 = add_edge!(n3,c2,c3)
@test typeof(e1) == BayesianEdge
@test length(edges(n3)) == 1
@test source(e1) == c2

@test length(nodes(n3)) == 3

@test_throws add_node!(n3,b1)
@test node_index(b1) == 1

###################################
#Tests that should cause errors
###################################
@test_throws add_edge!(n3,b1,b2)
@test_throws add_edge!(n3,c3,c1)
@test_throws q1 = DBayesianNode("hej", pd3)
@test_throws q1 = CBayesianNode("Hej", x -> 1)

###################################

pd7 = ProbabilityDistribution(["head","tails"],[0.46,0.54])
pd8 = ProbabilityDistribution(["head","tails"],[0.99,0.01])
pd9 = ProbabilityDistribution(["head","tails"],[0.99,0.01])
pd10 = ProbabilityDistribution(["head","tails"],[0.99,0.01])
pd11 = ProbabilityDistribution(["head","tails"],[0.99,0.01])
pd12 = ProbabilityDistribution(["head","tails"],[0.99,0.01])
pd13 = ProbabilityDistribution(["head","tails"],[0.99,0.01])
pd14 = ProbabilityDistribution(["head","tails"],[0.99,0.01])
d1 = DBayesianNode(:d1, pd7)
d2 = DBayesianNode(:d2, pd8)
d3 = DBayesianNode(:d3, pd9)
d4 = DBayesianNode(:d4, pd10)
d5 = DBayesianNode(:d5, pd11)
d6 = DBayesianNode(:d6, pd12)
d7 = DBayesianNode(:d7, pd13)
d8 = DBayesianNode(:d8, pd14)
n4 = BayesianNetwork([d1,d2,d3,d4,d5,d6,d7,d8],[])

add_edge!(n4,d2,d1)
add_edge!(n4,d3,d1)
add_edge!(n4,d6,d7)
add_edge!(n4,d7,d2)
add_edge!(n4,d8,d5)
add_edge!(n4,d1,d5)

@test legal_configuration(n4, CPD([:d1], [:d2,:d3])) == true
@test legal_configuration(n4, CPD([:d2,:d3], [:d1])) == true
@test legal_configuration(n4, CPD([:d1], [:d2,:d3,:d5,:d6,:d7,:d8])) == true
@test legal_configuration(n4, CPD([:d1], [:d2,:d3,:d8])) == true
@test legal_configuration(n4, CPD([:d2,:d3], [:d4])) == false
@test legal_configuration(n4, CPD([:d2], [:d4, :d1])) == false

@test has_cached_result(n4, P(:d9)) == false
@test has_cached_result(n4, P(:d1|:d2)) == false
@test has_cached_result(n4, P(:d4|:d2)) == false

####################################

pdbhs1 = ProbabilityDistribution(["head", "tails"], [0.5,0.5])
pdbhs2 = ProbabilityDistribution(["head", "tails"], [0.75,0.25])
pdbhs3 = ProbabilityDistribution(["head", "tails"], [0.33,0.67])
abhs1 = DBayesianNode(:R, pdbhs1)
abhs2 = DBayesianNode(:S, pdbhs2)
abhs3 = DBayesianNode(:D, pdbhs3)

nbhs1 = BayesianNetwork([abhs1, abhs2, abhs3])
add_edge!(nbhs1, abhs1, abhs2)
add_edge!(nbhs1, abhs1, abhs3)
add_edge!(nbhs1, abhs2, abhs3)

@test check_requirements(nbhs1, P(:R|:S)) == false
@test check_requirements(nbhs1, P(:S|:R)) == false
add_probability!(nbhs1, P(:R|:S), ProbabilityDistribution(["head", "tails"], [0.33, 0.67]))
@test check_requirements(nbhs1, P(:S|:R)) == true

@test check_requirements(nbhs1, P(:D|[:R,:S])) == false
add_probability!(nbhs1, P([:R,:S]|:D), ProbabilityDistribution(["head", "tails"], [0.20, 0.80]))
@test check_requirements(nbhs1, P(:D|[:R,:S])) == true
add_probability!(nbhs1, P([:R,:S]), ProbabilityDistribution(["head", "tails"], [0.40, 0.60]))
@test check_requirements(nbhs1, P(:D|[:R,:S])) == true

####################################
## Test query
####################################

DB=DBayesianNode(:D, ["d1","d2","d_water"])

an1=CBayesianNode(:a1)
an2=CBayesianNode(:a2)
an3=CBayesianNode(:a3)

bn1=BayesianNetwork([DB,an1,an2,an3])

add_edge!(bn1,DB,an1)
add_edge!(bn1,DB,an2)
add_edge!(bn1,DB,an3)

analysis=[an1,an2,an3]

# Distributions
using Distributions

for a in analysis
    a_d1_dist = Normal(1, 0.2)
    a_d2_dist = Normal(10, 0.2)

    samples = [ rand(dist) for i=1:100, dist in [a_d1_dist, a_d2_dist] ]

    a_d_water_pdf = x -> 1/(maximum(samples) - minimum(samples))

    add_probability!(bn1, P(a.label|:D), ProbabilityDensityDistribution(states(DB), [x -> 1, x -> 1, a_d_water_pdf]))
end

add_probability!(bn1, P(:D), ProbabilityDistribution(["d1","d2","d_water"], [0.4,0.4,0.2]))

testCPD1 = P(:D|Dict{Symbol,Any}({:a1 => 30, :a2 => 250}))
diagProbs1 = query(bn1,testCPD1)

@test_approx_eq diagProbs1["d1"] 0.4
@test_approx_eq diagProbs1["d2"] 0.4


text:
#=

  Euler #46 in Julia.

  Problem 46:
  """
  It was proposed by Christian Goldbach that every odd composite number can be
  written as the sum of a prime and twice a square.

  9 = 7 + 2×1^2
  15 = 7 + 2×2^2
  21 = 3 + 2×3^2
  25 = 7 + 2×3^2
  27 = 19 + 2×2^2
  33 = 31 + 2×1^2

  It turns out that the conjecture was false.

  What is the smallest odd composite that cannot be written as the
  sum of a prime and twice a square?
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# 0.00057412s
function euler46a()

    for i in 3:2:10000
        if !isPrime(i)
            s = round(Int,sqrt(i/2))
            found = 0
            for j in 1:s
                if found === 0
                    ts = j*j*2
                    if isPrime(abs(i-ts))
                        found = 1
                    end
                end
            end
            if found == 0
                return i
            end
        end
    end
    # return res
    return nothing
end

run_euler(euler46a)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
@doc """Checks whether a solver is applicable to a test case""" ->
function isapplicable(s::Solver, tc::TestCaseExplicit)
    out = _isapplicable(s,tc)
    if hasmass(tc)
        out = out && s.eq_type>=explicit_mass_eq
    end
    return out
end
function _isapplicable(s::Solver, tc::TestCase)
    out = true
    # stiffness
    out = out && s.stiff >= tc.stiff
    # dae
    out = out && s.dae>=tc.dae
    return out
end
isapplicable(tc::TestCaseExplicit, s::Solver) = isapplicable(s,tc)

# test whether a self consistent test-run:
isconsistent(tr::TestRunAdapt) = isapplicable(tr.solver, tr.tc) && isadaptive(tr)
isconsistent(tr::TestRunFixedStep) = isapplicable(tr.solver, tr.tc) && !isadaptive(tr)
isconsistent(ts::TestSuite) = isapplicable(ts.solver, ts.tc)

# ###
# # saving to files

# function save(res::TestResults)...

# function save(tc::


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

type Text
    atlas
    offsets
    uvs
    positions
    colors
    scales
end
length(rt::Text) = length(rt.text)
function sizehint!(t:Text, size::Int)
    for elem in (:text, :positions, :colors, :scales, :uvs, :offsets)
        sizehint!(getfield(t, elem), size)
    end
    nothing
end
const TextDefault = Dict(
    :atlas => GLVisualize.get_texture_atlas()
    :color => RGBA{Float32}(0,0,0,1)
    :scale => Vec2f0(50)
    :font => GLVisualize.DEFAULT_FONT_FACE
)

type RichText
    defaults::FixedKeyDict
    text::Text
    inlined_objects::Vector
    links::Vector
    cursor
end
function sizehint!(t:RichText, size::Int)
    sizehint!(t.text, size)
    nothing
end
length(rt::RichtText) = length(rt.text)
make_iter(x) = repeat(x)
make_iter(x::AbstractArray) = x

function get_iter(defaultfunc, dictlike, key)
    make_iter(get(defaultfunc, dictlike, key))
end

"""
Inserts text at [`position`] with `style`
"""
function insert!(text::RichText, text2::AbstractString, position=text.cursor, style=EMPTY_STYLE)
    positions = get_iter(()->positions(text, text2), style, :positions)
    offsets   = get_iter(()->offsets(text, text2),   style, :offsets)
    colors    = get_iter(()->font_color(text),       style, :colors)
    fonts     = get_iter(()->font(text),             style, :fonts)
    sizes     = get_iter(()->font_size(text),        style, :fontsizes)

    for elem in zip(text2, positions, colors, fonts, sizes)
        push!(text.text, elem)
    end
end

"""
Inserts code with highlighting
"""
function insert!(text::RichText, code::Code,
        position=text.cursor,
        color_scheme=color_scheme(text, code)
    )

end

"""
Inlines any `object` in `area` at `position`.
"""
function insert!(text::RichText, object, position=text.cursor, area=nothing)

end

function push!(text::RichText, values::RichtTextValues)
    for props in zip(text2, positions, colors, fonts, sizes)
        push!(text, props)
    end
    nothing
end


text:
# Autogenerated wrapper script for AprilTags_jll for x86_64-apple-darwin
export libapriltag

JLLWrappers.@generate_wrapper_header("AprilTags")
JLLWrappers.@declare_library_product(libapriltag, "@rpath/libapriltag.3.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libapriltag,
        "lib/libapriltag.3.1.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using BinDeps

@BinDeps.setup

libgit2 = library_dependency("libgit2")
version = "0.21.0"
provides(Sources,
         URI("https://github.com/libgit2/libgit2/archive/v$version.tar.gz"),
         libgit2,
         unpacked_dir="libgit2-$version")

# Windows binaries built in MSYS2 via:
# /d/code/CMake-2.8/bin/cmake .. -DCMAKE_INSTALL_PREFIX=$PWD/usr -G"MSYS Makefiles"
# /d/code/CMake-2.8/bin/cmake --build . --target install
provides(Binaries, URI("http://sourceforge.net/projects/juliadeps-win/files/libgit2-v$version.7z"),
         libgit2, unpacked_dir="usr$WORD_SIZE/bin", os = :Windows)

prefix = joinpath(BinDeps.depsdir(libgit2),"usr")
srcdir = joinpath(BinDeps.depsdir(libgit2),"src","libgit2-$version")

provides(SimpleBuild,
    (@build_steps begin
        GetSources(libgit2)
        @build_steps begin
            ChangeDirectory(srcdir)
            `cmake . -DCMAKE_INSTALL_PREFIX=$prefix
                     -DCMAKE_BUILD_TYPE=Release
                     -DTHREADSAFE=ON
                     -DBUILD_CLAR=OFF`
            `make install`
        end
    end), libgit2, os = :Unix)

@BinDeps.install [:libgit2 => :libgit2]


text:
using Distributions
using JSON, ForwardDiff, Calculus, PDMats # test dependencies
using Test
using Distributed
using Random
using StatsBase

tests = [
    "mvnormal",
    "mvlognormal",
    "types",
    "utils",
    "samplers",
    "categorical",
    "univariates",
    "continuous",
    "fit",
    "multinomial",
    "binomial",
    "poissonbinomial",
    "dirichlet",
    "dirichletmultinomial",
    "mvtdist",
    "kolmogorov",
    "edgeworth",
    "matrix",
    "vonmisesfisher",
    "conversion",
    "mixture",
    "gradlogpdf",
    "truncate",
    "noncentralt",
    "locationscale",
    "quantile_newton",
    "semicircle",
    "qq",
    "product",
    "truncnormal",
    "discretenonparametric",
    "functionals"
]

printstyled("Running tests:\n", color=:blue)

using Random
Random.seed!(345679)

res = map(tests) do t
    @eval module $(Symbol("Test_", t))
    using Distributions
    using JSON, ForwardDiff, Calculus, PDMats # test dependencies
    using Test
    using Random
    Random.seed!(345679)
    using LinearAlgebra
    using StatsBase
    include($t * ".jl")
    end
    return
end

# print method ambiguities
println("Potentially stale exports: ")
display(Test.detect_ambiguities(Distributions))
println()


text:
include("../utils/one-NN.jl")
include("../utils/distancias.jl")

export LearnerEuclideanOneNN
export WeightedLearnerEuclideanOneNN

"""
    LearnerEuclideanOneNN(data, labels)

Devuelve un clasificador 1-NN basado en la distancia Euclídea
"""
function LearnerEuclideanOneNN(data, labels)
    return f(x)= OneNN(x, EuclideanDistance, data, labels)
end

"""
    WeightedLearnerEuclideanOneNN(w,data, labels)

Devuelve un clasificador 1-NN basado en la distancia Euclídea
ponderada a `w`

"""
function WeightedLearnerEuclideanOneNN(w,data, labels)
    distance(x,y)=WeightedEuclideanDistance(x,y,w)
    return f(x)= OneNN(x, distance, data, labels)
end


text:
lines = readlines("input.txt")
numbers = parse.(Int, lines)

sums = sum.(zip(numbers[1:end-2], numbers[2:end-1], numbers[3:end]))
println(sum(sums[2:end] .> sums[1:end-1]))