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

using Base.Test

srand(101)
debug = false #Turn on for more debugging info

#Pauli σ-matrices
for σ in map(Hermitian, Any[ eye(2), [0 1; 1 0], [0 -im; im 0], [1 0; 0 -1] ])
    @test ishermitian(σ)
end

# Hermitian matrix exponential/log
let A1 = randn(4,4) + im*randn(4,4)
    A2 = A1 + A1'
    @test expm(A2) ≈ expm(Hermitian(A2))
    @test logm(A2) ≈ logm(Hermitian(A2))
    A3 = A1 * A1' # posdef
    @test expm(A3) ≈ expm(Hermitian(A3))
    @test logm(A3) ≈ logm(Hermitian(A3))
end

let A1 = randn(4,4)
    A3 = A1 * A1'
    A4 = A1 + A1.'
    @test expm(A4) ≈ expm(Symmetric(A4))
    @test logm(A3) ≈ logm(Symmetric(A3))
    @test logm(A3) ≈ logm(Hermitian(A3))
end

let n=10
    areal = randn(n,n)/2
    aimg  = randn(n,n)/2
    debug && println("symmetric eigendecomposition")
    for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
        a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
        asym = a'+a                 # symmetric indefinite
        ε = εa = eps(abs(float(one(eltya))))

        x = randn(n)
        y = randn(n)
        b = randn(n,n)/2
        x = eltya == Int ? rand(1:7, n) : convert(Vector{eltya}, eltya <: Complex ? complex.(x, zeros(n)) : x)
        y = eltya == Int ? rand(1:7, n) : convert(Vector{eltya}, eltya <: Complex ? complex.(y, zeros(n)) : y)
        b = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex.(b, zeros(n,n)) : b)

        debug && println("\ntype of a: ", eltya, "\n")

        # constructor
        @test Symmetric(Symmetric(asym, :U))     === Symmetric(asym, :U)
        @test Hermitian(Hermitian(asym, :U))     === Hermitian(asym, :U)
        @test Symmetric(Symmetric(asym, :U), :U) === Symmetric(asym, :U)
        @test Hermitian(Hermitian(asym, :U), :U) === Hermitian(asym, :U)
        @test_throws ArgumentError Symmetric(Symmetric(asym, :U), :L)
        @test_throws ArgumentError Hermitian(Hermitian(asym, :U), :L)

        # similar
        @test isa(similar(Symmetric(asym)), Symmetric{eltya})
        @test isa(similar(Hermitian(asym)), Hermitian{eltya})
        @test isa(similar(Symmetric(asym), Int), Symmetric{Int})
        @test isa(similar(Hermitian(asym), Int), Hermitian{Int})
        @test isa(similar(Symmetric(asym), (3,2)), Matrix{eltya})
        @test isa(similar(Hermitian(asym), (3,2)), Matrix{eltya})
        @test isa(similar(Symmetric(asym), Int, (3,2)), Matrix{Int})
        @test isa(similar(Hermitian(asym), Int, (3,2)), Matrix{Int})

        # full
        @test asym == full(Hermitian(asym))

        # parent
        @test asym == parent(Hermitian(asym))

        # getindex
        @test asym[1,1] == Hermitian(asym)[1,1]
        @test asym[1,1] == Symmetric(asym)[1,1]

        #trace
        @test trace(asym) == trace(Hermitian(asym))

        # issymmetric, ishermitian
        if eltya <: Real
            @test issymmetric(Symmetric(asym))
            @test ishermitian(Symmetric(asym))
        end
        if eltya <: Complex
            @test ishermitian(Symmetric(b + b'))
        end

        #transpose, ctranspose
        if eltya <: Real
            @test transpose(Symmetric(asym)) == asym
        else
            @test transpose(Hermitian(asym)) == transpose(asym)
        end
        @test ctranspose(Symmetric(asym)) == Symmetric(conj(asym))
        @test ctranspose(Hermitian(asym)) == asym

        #tril/triu
        for di in -n:n
            @test triu(Symmetric(a+a.'),di) == triu(a+a.',di)
            @test tril(Symmetric(a+a.'),di) == tril(a+a.',di)
            @test triu(Hermitian(asym),di) == triu(asym,di)
            @test tril(Hermitian(asym),di) == tril(asym,di)
            @test triu(Symmetric(a+a.',:L),di) == triu(a+a.',di)
            @test tril(Symmetric(a+a.',:L),di) == tril(a+a.',di)
            @test triu(Hermitian(asym,:L),di) == triu(asym,di)
            @test tril(Hermitian(asym,:L),di) == tril(asym,di)
        end

        eltya == BigFloat && continue # Revisit when implemented in julia
        d, v = eig(asym)
        @test asym*v[:,1] ≈ d[1]*v[:,1]
        @test v*Diagonal(d)*v' ≈ asym
        @test isequal(eigvals(asym[1]), eigvals(asym[1:1,1:1]))
        @test abs.(eigfact(Hermitian(asym), 1:2)[:vectors]'v[:,1:2]) ≈ eye(eltya, 2)
        eig(Hermitian(asym), 1:2) # same result, but checks that method works
        @test abs.(eigfact(Hermitian(asym), d[1] - 1, (d[2] + d[3])/2)[:vectors]'v[:,1:2]) ≈ eye(eltya, 2)
        eig(Hermitian(asym), d[1] - 1, (d[2] + d[3])/2) # same result, but checks that method works
        @test eigvals(Hermitian(asym), 1:2) ≈ d[1:2]
        @test eigvals(Hermitian(asym), d[1] - 1, (d[2] + d[3])/2) ≈ d[1:2]
        @test full(eigfact(asym)) ≈ asym

        # relation to svdvals
        @test sum(sort(abs.(eigvals(Hermitian(asym))))) == sum(sort(svdvals(Hermitian(asym))))

        # cond
        @test cond(Hermitian(asym)) ≈ cond(asym)

        # det
        @test det(asym) ≈ det(Hermitian(asym, :U))
        @test det(asym) ≈ det(Hermitian(asym, :L))
        if eltya <: Real
            @test det(asym) ≈ det(Symmetric(asym, :U))
            @test det(asym) ≈ det(Symmetric(asym, :L))
        end
        @test det(a + a.') ≈ det(Symmetric(a + a.', :U))
        @test det(a + a.') ≈ det(Symmetric(a + a.', :L))

        # isposdef[!]
        @test isposdef(Symmetric(asym)) == isposdef(full(Symmetric(asym)))
        @test isposdef(Hermitian(asym)) == isposdef(full(Hermitian(asym)))
        if eltya != Int
            @test isposdef!(Symmetric(copy(asym))) == isposdef(full(Symmetric(asym)))
            @test isposdef!(Hermitian(copy(asym))) == isposdef(full(Hermitian(asym)))
        end

        # rank
        let A = a[:,1:5]*a[:,1:5]'
            # Make sure A is Hermitian even in the present of rounding error
            # xianyi/OpenBLAS#729
            A = (A' + A) / 2
            @test rank(A) == rank(Hermitian(A))
        end

        # mat * vec
        if eltya <: Complex
            @test Hermitian(asym)*x+y ≈ asym*x+y
        end
        if eltya <: Real && eltya != Int
            @test Symmetric(asym)*x+y ≈ asym*x+y
        end

        C = zeros(eltya,n,n)
        # mat * mat
        if eltya <: Complex
            @test Hermitian(asym) * a ≈ asym * a
            @test a * Hermitian(asym) ≈ a * asym
            @test Hermitian(asym) * Hermitian(asym) ≈ asym*asym
            @test_throws DimensionMismatch Hermitian(asym) * ones(eltya,n+1)
            Base.LinAlg.A_mul_B!(C,a,Hermitian(asym))
            @test C ≈ a*asym
        end
        if eltya <: Real && eltya != Int
            @test Symmetric(asym) * Symmetric(asym) ≈ asym*asym
            @test Symmetric(asym) * a ≈ asym * a
            @test a * Symmetric(asym) ≈ a * asym
            @test_throws DimensionMismatch Symmetric(asym) * ones(eltya,n+1)
            Base.LinAlg.A_mul_B!(C,a,Symmetric(asym))
            @test C ≈ a*asym
        end

        # solver
        @test Hermitian(asym)\x ≈ asym\x
        if eltya <: Real
            @test Symmetric(asym)\x ≈ asym\x
        end

        #inversion
        @test inv(Hermitian(asym)) ≈ inv(asym)
        if eltya <: Real && eltya != Int
            @test inv(Symmetric(asym)) ≈ inv(asym)
            @test inv(Hermitian(a)) ≈ inv(full(Hermitian(a)))
            @test inv(Symmetric(a)) ≈ inv(full(Symmetric(a)))
        end

        # conversion
        @test Symmetric(asym) == convert(Symmetric,Symmetric(asym))
        if eltya <: Real && eltya != Int
            typs = [Float16,Float32,Float64]
            for typ in typs
                @test Symmetric(convert(Matrix{typ},asym)) == convert(Symmetric{typ,Matrix{typ}},Symmetric(asym))
            end
        end
        if eltya <: Complex && eltya != Int
            typs = [Complex64,Complex128]
            for typ in typs
                @test Hermitian(convert(Matrix{typ},asym)) == convert(Hermitian{typ,Matrix{typ}},Hermitian(asym))
            end
        end

        #unsafe_getindex
        if eltya <: Real
            @test Symmetric(asym)[1:2,1:2] == asym[1:2,1:2]
        end
        @test Hermitian(asym)[1:2,1:2] == asym[1:2,1:2]
    end
end

#Issue #7647: test xsyevr, xheevr, xstevr drivers
for Mi7647 in (Symmetric(diagm(1.0:3.0)),
               Hermitian(diagm(1.0:3.0)),
               Hermitian(diagm(complex(1.0:3.0))),
               SymTridiagonal([1.0:3.0;], zeros(2)))
    debug && println("Eigenvalues in interval for $(typeof(Mi7647))")
    @test eigmin(Mi7647)  == eigvals(Mi7647, 0.5, 1.5)[1] == 1.0
    @test eigmax(Mi7647)  == eigvals(Mi7647, 2.5, 3.5)[1] == 3.0
    @test eigvals(Mi7647) == eigvals(Mi7647, 0.5, 3.5) == [1.0:3.0;]
end

#Issue #7933
let A7933 = [1 2; 3 4]
    B7933 = copy(A7933)
    C7933 = full(Symmetric(A7933))
    @test A7933 == B7933
end

# Issues #8057 and #8058
for f in (eigfact, eigvals, eig)
    for A in (Symmetric([0 1; 1 0]), Hermitian([0 im; -im 0]))
        @test_throws ArgumentError f(A, 3, 2)
        @test_throws ArgumentError f(A, 1:4)
    end
end

#Issue 10671
let A = [1.0+im 2.0; 2.0 0.0]
    @test !ishermitian(A)
    @test_throws ArgumentError Hermitian(A)
end

# Unary minus for Symmetric matrices
let A = Symmetric(randn(5,5))
    B = -A
    @test A + B ≈ zeros(5,5)
end

# 17780
let a = randn(2,2)
    a = a'a
    b = complex.(a,a)
    c = Symmetric(b)
    @test conj(c) == conj(Array(c))
    cc = copy(c)
    @test conj!(c) == conj(Array(cc))
    c = Hermitian(b + b')
    @test conj(c) == conj(Array(c))
    cc = copy(c)
    @test conj!(c) == conj(Array(c))
end

# 19225
let X = [1 -1; -1 1]
    for T in (Symmetric, Hermitian)
        Y = T(copy(X))
        _Y = similar(Y)
        copy!(_Y, Y)
        @test _Y == Y

        W = T(copy(X), :L)
        copy!(W, Y)
        @test W.data == Y.data
        @test W.uplo != Y.uplo

        W[1,1] = 4
        @test W == T([4 -1; -1 1])
        @test_throws ArgumentError (W[1,2] = 2)

        @test Y + I == T([2 -1; -1 2])
        @test Y - I == T([0 -1; -1 0])
        @test Y * I == Y

        @test Y + 1 == T([2 0; 0 2])
        @test Y - 1 == T([0 -2; -2 0])
        @test Y * 2 == T([2 -2; -2 2])
        @test Y / 1 == Y

        @test T([true false; false true]) + true == T([2 1; 1 2])
    end

    @test_throws ArgumentError Hermitian(X) + 2im*I
    @test_throws ArgumentError Hermitian(X) - 2im*I
end

max_stars_repo_name: odow/julia

max_issues_repo_name: odow/julia

max_forks_repo_name: odow/julia


text:
function parse_fun_expr(ex)
    if !(ex.head == :call)
        throw(ArgumentError("Invalid macro function definition"))
    end
    f_name = ex.args[1]
    vars = Symbol[]
    types = Union{Symbol,Expr,DataType}[]
    for var_ex in ex.args[2:end]
        if isa(var_ex, Expr) && var_ex.head == :(::)
            sym = var_ex.args[1]::Symbol
            typ = var_ex.args[2]::Union{Symbol,Expr,DataType}
        elseif isa(var_ex, Symbol)
            sym = var_ex
            typ = :Any
        else
            throw(ArgumentError("Can't parse variable list in " *
                                "macro function defintion"))
        end
        push!(vars, sym)
        push!(types, typ)
    end
    return f_name, vars, types
end

"""
Macro for calling CLBLAS C API functions. Usage example:

    @api2.blasfun clblasSaxpy(n::Csize_t, alpha::CL_float2,
                        X::cl.CL_mem, offx::Csize_t, incx::Cint,
                        Y::cl.CL_mem, offy::Csize_t, incy::Cint,
                        n_queues::cl.CL_uint,
                        queues::Ptr{cl.CL_command_queue},
                        n_events_in_wait_list::cl.CL_uint,
                        event_wait_list::Ptr{cl.CL_event},
                        events::Ptr{cl.CL_event})

Note: All parameter types are resolved in a macro definition environment,
      so only OpenCL types (e.g. cl.CL_float), CL_float2 and CL_double2
      are currently supported (which is essentially enough for wrapping
      CLBLAS C API)
"""
macro blasfun(expr)
    f, args, types = parse_fun_expr(expr)
    ex = Expr(:function)
    ex_fun = Expr(:(call))
    push!(ex_fun.args, f)
    for (arg, typ) in zip(args, types)
        var_ex = Expr(:(::))
        push!(var_ex.args, arg)
        push!(var_ex.args, typ)
        push!(ex_fun.args, var_ex)
    end
    push!(ex.args, ex_fun)

    ex_body = quote
        local err::cl.CL_int = ccall(($(string(f)), libCLBLAS),
                                     cl.CL_int,
                                     ($(map(eval, types)...),),
                                     $(args...))
        if err != cl.CL_SUCCESS
            throw(cl.CLError(err))
        end
        return err
    end

    push!(ex.args, ex_body)
    return esc(ex)
end

function kw_param_expr(name, dtype, value)
    Expr(:parameters, Expr(:kw, Expr(:(::), name, dtype), value))
end

"""
Shorter version of @blasfun with reduced number of parameters
"""
macro blasfun2(expr)
    f, args, types = parse_fun_expr(expr)
    ex = Expr(:function)
    ex_fun = Expr(:(call))
    push!(ex_fun.args, f)
    # optional list of events (kw args should go first in AST)
    push!(ex_fun.args,
          kw_param_expr(:inEvents, Vector{cl.CL_event}, cl.CL_event[]))
    # passed parameters
    for (arg, typ) in zip(args, types)
        var_ex = Expr(:(::))
        push!(var_ex.args, arg)
        push!(var_ex.args, typ)
        push!(ex_fun.args, var_ex)
    end
    # additional argument - list of queues
    push!(ex_fun.args, Expr(:(::), :inQueues, Vector{cl.CmdQueue}))
    push!(ex.args, ex_fun)

    ex_body = quote
        local ctx = cl.info(inQueues[1], :context)
        local num_queues = cl.CL_uint(length(inQueues))
        local queues = cl.CL_command_queue[queue.id for queue in inQueues]
        local num_events = cl.CL_uint(length(inEvents))
        local events = num_events == 0 ? C_NULL : Ref(inEvents)
        local ret_event = Ref{cl.CL_event}()

        # local err = $f($(args...), num_queues, pointer(queues), num_events,
        #                pointer(events), pointer(ret_event))
        local err::cl.CL_int = ccall(($(string(f)), libCLBLAS),
                                     cl.CL_int,
                                     ($(map(eval, types)...),
                                      # common types
                                      Cuint, Ptr{cl.CL_command_queue},
                                      Cuint, Ptr{cl.CL_event},
                                      Ptr{cl.CL_event},),
                                     $(args...),
                                     # common params
                                     num_queues, Ref(queues),
                                     num_events, events,
                                     ret_event)
        if err != cl.CL_SUCCESS
            throw(cl.CLError(err))
        end
        return cl.Event(ret_event[], retain=false)
    end

    push!(ex.args, ex_body)
    return esc(ex)
end

ext: jl


text:
using Test

# Include tests


text:
using Test, PolynomialBases
using InteractiveUtils

# regression test
for basis_type in subtypes(PolynomialBases.NodalBasis{PolynomialBases.Line})
    println(devnull, "  ", basis_type(1))
    for p in 0:20, T in (Float32, Float64)
        basis = try
            basis_type(p, T)
        catch m
            isa(m, DomainError) && continue
            throw(m)
        end
        xmin, xmax = 0.25, 0.9
        x = similar(basis.nodes)
        ξ = similar(x)
        map_from_canonical!(x, basis.nodes, xmin, xmax, basis)
        map_to_canonical!(ξ, x, xmin, xmax, basis)
        @test ξ ≈ basis.nodes atol=10*eps(T)
    end
end


text:
function readcigarmd2seqref_sam2pairwise(read, cigar, md)
    seq_str = convert(String, read)
    command = `./readcigarmd2seqref $seq_str $cigar $md`
    output = readchomp(command)
    seq, ref = split(output)
    return seq, ref
end

max_stars_repo_name: ChristianMichelsen/reference-free

max_issues_repo_name: ChristianMichelsen/reference-free

max_forks_repo_name: ChristianMichelsen/reference-free


text:
function pcGA(objfun::Function,
             vals::Vector{Int},
             actual::Bool
             )
    # Harbored Parameters
    N = length(vals)
    iterations = 50*N
    initPopulation = randperm
    populationSize = 50
    crossoverRate = 0.5
    mutationRate = 0.25
    selection = tournament(2)
    crossover = ox1
    mutation = swap2!
    ɛ = 0.1

    # Setup parameters
    elite = isa(ɛ, Int) ? ɛ : round(Int, ɛ * populationSize)
    fitFunc = inverseFunc(objfun)
    # Initialize population
    fitness = zeros(Float64,populationSize)
    population = Array(Vector{Int}, populationSize)
    offspring = similar(population)
    # Generate population
    for i in 1:populationSize
      population[i] = initPopulation(N)
      fitness[i] = fitFunc(vals[population[i]])
    end
    fitidx = sortperm(fitness, rev = true)
    # Generate and evaluate offspring
    itr = 1
    bestFitness = 0.0
    bestIndividual = 0
    while true
        # Select offspring
        selected = selection(fitness, populationSize)
        # Perform mating
        offidx = randperm(populationSize)
        for i in 1:2:populationSize
            j = (i == populationSize) ? i-1 : i+1
            if rand() < crossoverRate
                offspring[i], offspring[j] = crossover(population[selected[offidx[i]]], population[selected[offidx[j]]])
            else
                offspring[i], offspring[j] = copy(population[selected[i]]), copy(population[selected[j]])
            end
        end
        # Perform mutation
        for i in 1:populationSize
            if rand() < mutationRate
                mutation(offspring[i])
            end
        end
        # Elitism
        if elite > 0
          subs = sample(1:populationSize,elite,replace=false)
            for i in 1:elite
                offspring[subs[i]] = population[fitidx[i]]
            end
        end
        # New generation
        for i in 1:populationSize
            population[i] = copy(offspring[i])
            fitness[i] = fitFunc(vals[offspring[i]])
        end
        fitidx = sortperm(fitness, rev = true)
        bestIndividual = fitidx[1]
        curGenFitness = objfun(vals[population[fitidx[1]]])
        fittol = abs(bestFitness - curGenFitness)
        bestFitness = curGenFitness
        # if number of iterations more then specified
        if itr >= iterations
            break
        end
        itr += 1
    end
    if actual
      out = vals[population[fitidx[1]]]
    else
      out = population[fitidx[1]]
    end
    return out
end


text:
module debspkg

# Write your package code here.
using ForwardDiff

include("extra_file.jl")

export f,derivative_of_f1

end


text:
function bval_u(Te, Ti, trgt_sgn, H3)

    cs = trgt_sgn .* sqrt.(Te + Ti)

    return H3 * cs
end

function bval_phi(phi, Te, lcfs_avg, H1, H2, H3)

    lambda = 2.695
    phi_in = dot(lcfs_avg, phi)*ones(length(phi))
    phi_out = 2.695 * Te

    return H1*phi_in + (H2 + H3)*phi_out
end

function solve_pde(A, x0, b, xb, gc::GhostConditions)

    An, B = require_boundary_conditions(A, gc)
    bn = gc.Proj*b - B*xb
    x = jacobi_preconditioned_gmres(An, gc.Proj*x0, bn, 1, 1)
    # x = An \ bn
    return gc.Mirror*transpose(gc.Proj)*x + (I - gc.Mirror)*xb
end

function solve_diffusion(x, xb, k, Dxx, Dyy, gc::GhostConditions)

    M = I - k*(Dxx + Dyy)
    return solve_pde(M, x, x, xb, gc)
end

function solve_vorticity_eqn(phi, phi_b, w, n, lnn_x, lnn_y, Pi_xx, Pi_yy, ad, Dx, Dy, Dxx, Dyy, gc::GhostConditions)

    M = Dxx + Dyy + Diagonal(lnn_x) * Dx + Diagonal(lnn_y) * Dy
    b = (w - ad*(Pi_xx + Pi_yy)) ./ n

    return solve_pde(M, phi, b, phi_b, gc)
end

function solve_helmholtz_eqn(psi, A, de2, Dxx, Dyy, gc::GhostConditions)

    M = I - de2*(Dxx + Dyy)
    return solve_pde(M, psi, A, zeros(length(A)), gc)
end

function forward_euler(x, x_t, xb, K1, K2, dt, gc::GhostConditions)

    return (x + dt*(K1*x_t + K2*gc.Mirror*xb)) ./ (1 .+ dt*diag(K2))
end

function leapfrog!(lnn, lnTe, lnTi, u, w, A, phi, psi, n, Te, Ti, Pe, Pi, j, jn, Sn, STe, STi, Dx, Dy, Dxy, Dxx, Dyy, Dxxx, Dyyy,
                   Dxxy, Dxyy, Ds, Dss, dt, N_subcycle, K1, K2, H1, H2, H3, trgt_sgn, lcfs_avg, GC_nmann, GC_dchlt, GC_u, kdiff_lnn,
                   kdiff_lnTe, kdiff_lnTi, kdiff_u, kdiff_w, kdiff_A, am, ad, ki, ke, er, eg, ev, de2, eta)

    # homogeneous boundary values
    bval_hom = zeros(length(Te))

    # take derivatives of psi
    psi_x = Dx*psi
    psi_y = Dy*psi
    psi_xy = Dxy*psi
    psi_xx = Dxx*psi
    psi_yy = Dyy*psi

    # subcycle electron thermal conduction terms
    lnTe0 = lnTe[:,2]
    for _=1:N_subcycle

        # derivatives
        lnTe_x = Dx*lnTe[:,2]
        lnTe_y = Dy*lnTe[:,2]
        lnTe_xy = Dxy*lnTe[:,2]
        lnTe_xx = Dxx*lnTe[:,2]
        lnTe_yy = Dyy*lnTe[:,2]
        lnTe_s = partial_s(lnTe[:,2], lnTe_x, lnTe_y, psi_x, psi_y, Ds, am)
        lnTe_ss = partial_ss(lnTe[:,2], lnTe_s, lnTe_x, lnTe_y, lnTe_xy, lnTe_xx, lnTe_yy, psi_x, psi_y, psi_xy, psi_xx, psi_yy, Dss, Ds, Dx, Dy, am)

        lnTe_t = lnT_partial_thermal(lnTe_s, lnTe_ss, Te, n, ke) + STe
        lnTe[:,3] = forward_euler(0.5*(lnTe[:,1] + lnTe[:,2]), lnTe_t, bval_hom, K1, K2, dt/N_subcycle, GC_nmann)

        # derivatives again
        lnTe_x = Dx*lnTe[:,2]
        lnTe_y = Dy*lnTe[:,2]
        lnTe_xy = Dxy*lnTe[:,2]
        lnTe_xx = Dxx*lnTe[:,2]
        lnTe_yy = Dyy*lnTe[:,2]
        lnTe_s = partial_s(lnTe[:,2], lnTe_x, lnTe_y, psi_x, psi_y, Ds, am)
        lnTe_ss = partial_ss(lnTe[:,2], lnTe_s, lnTe_x, lnTe_y, lnTe_xy, lnTe_xx, lnTe_yy, psi_x, psi_y, psi_xy, psi_xx, psi_yy, Dss, Ds, Dx, Dy, am)

        lnTe_t = lnT_partial_thermal(lnTe_s, lnTe_ss, Te, n, ke) + STe
        lnTe[:,3] = forward_euler(lnTe[:,2], lnTe_t, bval_hom, K1, K2, dt/N_subcycle, GC_nmann)

        # shift indices
        lnTe[:,1:2] = lnTe[:,2:3]
    end
    lnTe[:,2] = 0.5*(lnTe0 + lnTe[:,2])

    # subcycle ion thermal conduction terms
    lnTi0 = lnTi[:,2]
    for _=1:N_subcycle

        # derivatives
        lnTi_x = Dx*lnTi[:,2]
        lnTi_y = Dy*lnTi[:,2]
        lnTi_xy = Dxy*lnTi[:,2]
        lnTi_xx = Dxx*lnTi[:,2]
        lnTi_yy = Dyy*lnTi[:,2]
        lnTi_s = partial_s(lnTi[:,2], lnTi_x, lnTi_y, psi_x, psi_y, Ds, am)
        lnTi_ss = partial_ss(lnTi[:,2], lnTi_s, lnTi_x, lnTi_y, lnTi_xy, lnTi_xx, lnTi_yy, psi_x, psi_y, psi_xy, psi_xx, psi_yy, Dss, Ds, Dx, Dy, am)

        lnTi_t = lnT_partial_thermal(lnTi_s, lnTi_ss, Te, n, ke) + STi
        lnTi[:,3] = forward_euler(0.5*(lnTi[:,1] + lnTi[:,2]), lnTi_t, bval_hom, K1, K2, dt/N_subcycle, GC_nmann)

        # derivatives again
        lnTi_x = Dx*lnTi[:,2]
        lnTi_y = Dy*lnTi[:,2]
        lnTi_xy = Dxy*lnTi[:,2]
        lnTi_xx = Dxx*lnTi[:,2]
        lnTi_yy = Dyy*lnTi[:,2]
        lnTi_s = partial_s(lnTi[:,2], lnTi_x, lnTi_y, psi_x, psi_y, Ds, am)
        lnTi_ss = partial_ss(lnTi[:,2], lnTi_s, lnTi_x, lnTi_y, lnTi_xy, lnTi_xx, lnTi_yy, psi_x, psi_y, psi_xy, psi_xx, psi_yy, Dss, Ds, Dx, Dy, am)

        lnTi_t = lnT_partial_thermal(lnTi_s, lnTi_ss, Te, n, ke) + STi
        lnTi[:,3] = forward_euler(lnTi[:,2], lnTi_t, bval_hom, K1, K2, dt/N_subcycle, GC_nmann)

        # shift indices
        lnTi[:,1:2] = lnTi[:,2:3]
    end
    lnTi[:,2] = 0.5*(lnTi0 + lnTi[:,2])

    # take other derivatives
    phi_x = Dx*phi
    phi_y = Dy*phi
    phi_xy = Dxy*phi
    phi_xx = Dxx*phi
    phi_yy = Dyy*phi
    phi_xxx = Dxxx*phi
    phi_yyy = Dyyy*phi
    phi_xxy = Dxxy*phi
    phi_xyy = Dxyy*phi
    phi_c = -2*phi_y[:]
    phi_s = partial_s(phi, phi_x, phi_y, psi_x, psi_y, Ds, am)

    Pe_x = Dx*Pe
    Pe_y = Dy*Pe
    Pe_xy = Dxy*Pe
    Pe_xx = Dxx*Pe
    Pe_yy = Dyy*Pe
    Pe_xxx = Dxxx*Pe
    Pe_yyy = Dyyy*Pe
    Pe_xxy = Dxxy*Pe
    Pe_xyy = Dxyy*Pe
    Pe_c = -2*Pe_y[:]
    Pe_s = partial_s(Pe, Pe_x, Pe_y, psi_x, psi_y, Ds, am)

    Pi_x = Dx*Pi
    Pi_y = Dy*Pi
    Pi_c = -2*Pi_y[:]
    Pi_s = partial_s(Pi, Pi_x, Pi_y, psi_x, psi_y, Ds, am)

    lnn_x = Dx*lnn[:,2]
    lnn_y = Dy*lnn[:,2]
    lnn_xy = Dxy*lnn[:,2]
    lnn_xx = Dxx*lnn[:,2]
    lnn_yy = Dyy*lnn[:,2]
    lnn_s = partial_s(lnn[:,2], lnn_x, lnn_y, psi_x, psi_y, Ds, am)

    lnTe_x = Dx*lnTe[:,2]
    lnTe_y = Dy*lnTe[:,2]
    lnTe_c = -2*lnTe_y[:]
    lnTe_s = partial_s(lnTe[:,2], lnTe_x, lnTe_y, psi_x, psi_y, Ds, am)

    lnTi_x = Dx*lnTi[:,2]
    lnTi_y = Dy*lnTi[:,2]
    lnTi_c = -2*lnTi_y[:]
    lnTi_s = partial_s(lnTi[:,2], lnTi_x, lnTi_y, psi_x, psi_y, Ds, am)

    u_x = Dx*u[:,2]
    u_y = Dy*u[:,2]
    u_c = -2*u_y[:]
    u_s = partial_s(u[:,2], u_x, u_y, psi_x, psi_y, Ds, am)

    j_x = Dx*j
    j_y = Dy*j
    j_s = partial_s(j, j_x, j_y, psi_x, psi_y, Ds, am)

    jn_x = Dx*jn
    jn_y = Dy*jn
    jn_s = partial_s(jn, jn_x, jn_y, psi_x, psi_y, Ds, am)

    G = G_def(Ti, phi_c, Pi_c, u_s, er, ev, ad)
    G_x = Dx*G
    G_y = Dy*G
    G_c = -2*G_y[:]
    G_s = partial_s(G, G_x, G_y, psi_x, psi_y, Ds, am)

    # evaluate total derivatives
    lnn_Dt = lnn_total_ion_derivative(phi_c, Pe_c, n, j_s, u_s, er, ev, ad)
    lnTe_Dt = lnTe_total_electron_derivative_adiabatic(lnn_Dt, j, lnn_s, n, Te, lnTe_c, er, ad)
    lnTi_Dt = lnTi_total_ion_derivative_adiabatic(lnn_Dt, Ti, lnTi_c, er, ad)
    u_Dt = u_total_ion_derivative(Pe_s, Pi_s, G_s, n, Ti, u_c, ev, eg, er, ad)

    # evaluate partial derivatives
    lnn_t = general_partial_derivative(lnn_Dt, lnn_x, lnn_y, lnn_s, phi_x, phi_y, u[:,2]) + Sn
    lnTe_t = general_partial_derivative(lnTe_Dt, lnTe_x, lnTe_y, lnTe_s, phi_x, phi_y, u[:,2] - jn / ev) + STe
    lnTi_t = general_partial_derivative(lnTi_Dt, lnTi_x, lnTi_y, lnTi_s, phi_x, phi_y, u[:,2]) + STi
    u_t = general_partial_derivative(u_Dt, u_x, u_y, u_s, phi_x, phi_y, u[:,2])
    w_t = w_partial_t(n, lnn_x, lnn_y, lnn_xy, lnn_xx, lnn_yy,
                         Pe_x, Pe_y, Pe_xy, Pe_xx, Pe_yy, Pe_xxx, Pe_yyy, Pe_xxy, Pe_xyy, Pe_c,
                         phi_x, phi_y, phi_xy, phi_xx, phi_yy, phi_xxx, phi_yyy, phi_xxy, phi_xyy,
                         G_c, Pi_c, j_s, ad, eg)
    A_t = A_partial_t(phi_x, phi_y, phi_s, j, jn, jn_x, jn_y, jn_s, Pe_s, u[:,2], n, Te, de2, ev, ad, am, eta)

    # compute inhomogeneous boundary values
    u_b = bval_u(Te, Ti, trgt_sgn, H3)
    phi_b = bval_phi(phi, Te, lcfs_avg, H1, H2, H3)

    # compute predictor step
    lnn[:,3] = forward_euler(0.5*(lnn[:,1] + lnn[:,2]), lnn_t, bval_hom, K1, K2, dt, GC_nmann)
    lnTe[:,3] = forward_euler(0.5*(lnTe[:,1] + lnTi[:,2]), lnTe_t, bval_hom, K1, K2, dt, GC_nmann)
    lnTi[:,3] = forward_euler(0.5*(lnTi[:,1] + lnTi[:,2]), lnTi_t, bval_hom, K1, K2, dt, GC_nmann)
    u[:,3] = forward_euler(0.5*(u[:,1] + u[:,2]), u_t, u_b, K1, K2, dt, GC_u)
    w[:,3] = forward_euler(0.5*(w[:,1] + w[:,2]), w_t, bval_hom, K1, K2, dt, GC_dchlt)
    A[:,3] = forward_euler(0.5*(A[:,1] + A[:,2]), A_t, bval_hom, K1, K2, dt, GC_dchlt)

    # solve diffusion problems
    lnn[:,3] = solve_diffusion(lnn[:,3], bval_hom, kdiff_lnn, Dxx, Dyy, GC_nmann)
    lnTe[:,3] = solve_diffusion(lnTe[:,3], bval_hom, kdiff_lnTe, Dxx, Dyy, GC_nmann)
    lnTi[:,3] = solve_diffusion(lnTi[:,3], bval_hom, kdiff_lnTi, Dxx, Dyy, GC_nmann)
    u[:,3] = solve_diffusion(u[:,3], u_b, kdiff_u, Dxx, Dyy, GC_u)
    w[:,3] = solve_diffusion(w[:,3], bval_hom, kdiff_w, Dxx, Dyy, GC_nmann)
    A[:,3] = solve_diffusion(A[:,3], bval_hom, kdiff_A, Dxx, Dyy, GC_nmann)

    # exponentiate log terms
    n[:] = exp.(lnn[:,3])
    Te[:] = exp.(lnTe[:,3])
    Ti[:] = exp.(lnTi[:,3])
    Pe[:] = n .* Te
    Pi[:] = n .* Ti

    # solve vorticity and helmholtz equations
    lnn_x = Dx * lnn[:,3]
    lnn_y = Dy * lnn[:,3]
    Pi_xx = Dxx * Pi
    Pi_yy = Dyy * Pi
    phi = solve_vorticity_eqn(phi, phi_b, w[:,3], n, lnn_x, lnn_y, Pi_xx, Pi_yy, ad, Dx, Dy, Dxx, Dyy, GC_dchlt)
    psi = solve_helmholtz_eqn(psi, A[:,3], de2, Dxx, Dyy, GC_dchlt)

    # compute current
    j[:] = (Dxx + Dyy) * psi
    jn[:] = j ./ n

    # take derivatives
    psi_x = Dx*psi
    psi_y = Dy*psi
    psi_xy = Dxy*psi
    psi_xx = Dxx*psi
    psi_yy = Dyy*psi

    phi_x = Dx*phi
    phi_y = Dy*phi
    phi_xy = Dxy*phi
    phi_xx = Dxx*phi
    phi_yy = Dyy*phi
    phi_xxx = Dxxx*phi
    phi_yyy = Dyyy*phi
    phi_xxy = Dxxy*phi
    phi_xyy = Dxyy*phi
    phi_c = -2*phi_y[:]
    phi_s = partial_s(phi, phi_x, phi_y, psi_x, psi_y, Ds, am)

    Pe_x = Dx*Pe
    Pe_y = Dy*Pe
    Pe_xy = Dxy*Pe
    Pe_xx = Dxx*Pe
    Pe_yy = Dyy*Pe
    Pe_xxx = Dxxx*Pe
    Pe_yyy = Dyyy*Pe
    Pe_xxy = Dxxy*Pe
    Pe_xyy = Dxyy*Pe
    Pe_c = -2*Pe_y[:]
    Pe_s = partial_s(Pe, Pe_x, Pe_y, psi_x, psi_y, Ds, am)

    Pi_x = Dx*Pi
    Pi_y = Dy*Pi
    Pi_c = -2*Pi_y[:]
    Pi_s = partial_s(Pi, Pi_x, Pi_y, psi_x, psi_y, Ds, am)

    lnn_x = Dx*lnn[:,3]
    lnn_y = Dy*lnn[:,3]
    lnn_xy = Dxy*lnn[:,3]
    lnn_xx = Dxx*lnn[:,3]
    lnn_yy = Dyy*lnn[:,3]
    lnn_s = partial_s(lnn[:,3], lnn_x, lnn_y, psi_x, psi_y, Ds, am)

    lnTe_x = Dx*lnTe[:,3]
    lnTe_y = Dy*lnTe[:,3]
    lnTe_c = -2*lnTe_y[:]
    lnTe_s = partial_s(lnTe[:,3], lnTe_x, lnTe_y, psi_x, psi_y, Ds, am)

    lnTi_x = Dx*lnTi[:,3]
    lnTi_y = Dy*lnTi[:,3]
    lnTi_c = -2*lnTi_y[:]
    lnTi_s = partial_s(lnTi[:,3], lnTi_x, lnTi_y, psi_x, psi_y, Ds, am)

    u_x = Dx*u[:,3]
    u_y = Dy*u[:,3]
    u_c = -2*u_y[:]
    u_s = partial_s(u[:,3], u_x, u_y, psi_x, psi_y, Ds, am)

    j_x = Dx*j
    j_y = Dy*j
    j_s = partial_s(j, j_x, j_y, psi_x, psi_y, Ds, am)

    jn_x = Dx*jn
    jn_y = Dy*jn
    jn_s = partial_s(jn, jn_x, jn_y, psi_x, psi_y, Ds, am)

    G = G_def(Ti, phi_c, Pi_c, u_s, er, ev, ad)
    G_x = Dx*G
    G_y = Dy*G
    G_c = -2*G_y[:]
    G_s = partial_s(G, G_x, G_y, psi_x, psi_y, Ds, am)

    # evaluate total derivatives
    lnn_Dt = lnn_total_ion_derivative(phi_c, Pe_c, n, j_s, u_s, er, ev, ad)
    lnTe_Dt = lnTe_total_electron_derivative_adiabatic(lnn_Dt, j, lnn_s, n, Te, lnTe_c, er, ad)
    lnTi_Dt = lnTi_total_ion_derivative_adiabatic(lnn_Dt, Ti, lnTi_c, er, ad)
    u_Dt = u_total_ion_derivative(Pe_s, Pi_s, G_s, n, Ti, u_c, ev, eg, er, ad)

    # evaluate partial derivatives
    lnn_t = general_partial_derivative(lnn_Dt, lnn_x, lnn_y, lnn_s, phi_x, phi_y, u[:,3])
    lnTe_t = general_partial_derivative(lnTe_Dt, lnTe_x, lnTe_y, lnTe_s, phi_x, phi_y, u[:,3] - jn / ev)
    lnTi_t = general_partial_derivative(lnTi_Dt, lnTi_x, lnTi_y, lnTi_s, phi_x, phi_y, u[:,3])
    u_t = general_partial_derivative(u_Dt, u_x, u_y, u_s, phi_x, phi_y, u[:,3])
    w_t = w_partial_t(n, lnn_x, lnn_y, lnn_xy, lnn_xx, lnn_yy,
                         Pe_x, Pe_y, Pe_xy, Pe_xx, Pe_yy, Pe_xxx, Pe_yyy, Pe_xxy, Pe_xyy, Pe_c,
                         phi_x, phi_y, phi_xy, phi_xx, phi_yy, phi_xxx, phi_yyy, phi_xxy, phi_xyy,
                         G_c, Pi_c, j_s, ad, eg)
    A_t = A_partial_t(phi_x, phi_y, phi_s, j, jn, jn_x, jn_y, jn_s, Pe_s, u[:,3], n, Te, de2, ev, ad, am, eta)

    # compute inhomogeneous boundary values
    u_b = bval_u(Te, Ti, trgt_sgn, H3)
    phi_b = bval_phi(phi, Te, lcfs_avg, H1, H2, H3)

    # compute corrector step
    lnn[:,3] = forward_euler(lnn[:,2], lnn_t, bval_hom, K1, K2, dt, GC_nmann)
    lnTe[:,3] = forward_euler(lnTi[:,2], lnTe_t, bval_hom, K1, K2, dt, GC_nmann)
    lnTi[:,3] = forward_euler(lnTi[:,2], lnTi_t, bval_hom, K1, K2, dt, GC_nmann)
    u[:,3] = forward_euler(u[:,2], u_t, u_b, K1, K2, dt, GC_u)
    w[:,3] = forward_euler(w[:,2], w_t, bval_hom, K1, K2, dt, GC_dchlt)
    A[:,3] = forward_euler(A[:,2], A_t, bval_hom, K1, K2, dt, GC_dchlt)

    # solve diffusion problems
    lnn[:,3] = solve_diffusion(lnn[:,3], bval_hom, kdiff_lnn, Dxx, Dyy, GC_nmann)
    lnTe[:,3] = solve_diffusion(lnTe[:,3], bval_hom, kdiff_lnTe, Dxx, Dyy, GC_nmann)
    lnTi[:,3] = solve_diffusion(lnTi[:,3], bval_hom, kdiff_lnTi, Dxx, Dyy, GC_nmann)
    u[:,3] = solve_diffusion(u[:,3], u_b, kdiff_u, Dxx, Dyy, GC_u)
    w[:,3] = solve_diffusion(w[:,3], bval_hom, kdiff_w, Dxx, Dyy, GC_nmann)
    A[:,3] = solve_diffusion(A[:,3], bval_hom, kdiff_A, Dxx, Dyy, GC_nmann)

    # exponentiate log terms
    n[:] = exp.(lnn[:,3])
    Te[:] = exp.(lnTe[:,3])
    Ti[:] = exp.(lnTi[:,3])
    Pe[:] = n .* Te
    Pi[:] = n .* Ti

    # solve vorticity and helmholtz equations
    lnn_x = Dx * lnn[:,3]
    lnn_y = Dy * lnn[:,3]
    Pi_xx = Dxx * Pi
    Pi_yy = Dyy * Pi
    phi = solve_vorticity_eqn(phi, phi_b, w[:,3], n, lnn_x, lnn_y, Pi_xx, Pi_yy, ad, Dx, Dy, Dxx, Dyy, GC_dchlt)
    psi = solve_helmholtz_eqn(psi, A[:,3], de2, Dxx, Dyy, GC_dchlt)

    # compute current
    j[:] = (Dxx + Dyy) * psi
    jn[:] = j ./ n

    # shift indices
    lnn[:,1:2] = lnn[:,2:3]
    lnTe[:,1:2] = lnTe[:,2:3]
    lnTi[:,1:2] = lnTi[:,2:3]
    u[:,1:2] = u[:,2:3]
    w[:,1:2] = w[:,2:3]
    A[:,1:2] = A[:,2:3]
end


text:
function pythagorean_triplets(N::Int)::Array{Tuple{Int,Int,Int},1}
    # Return all (a,b,c) such that:
    # a + b + c = N
    # a^2 + b^2 = c^2
    # a < b < c

    # Pythagorean triplets: Euclid did this a while ago
    # a = m^2 - n^2, b = 2mn, c = m^2 + n^2
    # so a + b + c = 2m^2 + 2mn = N
    # m(m + n) = N/2: not nice for generation, so let's give up and filter by it instead
    #
    # Need to take care of int multiples so also k.*(a,b,c) => km(m+n) = N/2 too
    # => N/2 / m(m+n) is an integer: divrem(N/2,m(m+n))

    # TODO / exercise for the reader: ensure bounds on N/2 are as tight as they can be
    sort!(triple.(Iterators.filter(x->conditions(x[1],x[2],N),Iterators.product(1:N/2,1:N/2)),N) |> unique; by=x->x[1])
end

function conditions(m,n,N)
    m > n > 0 &&
    rem(N/2,m*(m+n)) == 0
end

triple(t,N) = triple(t[1],t[2],N)

# We only make the triple of answers that we know will work so there's no point optimising this
triple(m,n,N) = div(N/2,m*(m+n)).*sort!((m^2 - n^2, 2m*n, m^2 + n^2) |> collect) |> Tuple{Int,Int,Int}

# This is vastly faster than mine, it's not fair
# https://exercism.io/tracks/julia/exercises/pythagorean-triplet/solutions/369f511a364d4ce0b9c054c55b49d8d8
function pythagorean_triplets(x::Int)

	triplets = Array{Tuple{Int64,Int64,Int64}}(undef, 0)
	for a in 1:(x ÷ 3)
		for b in (a+1):(x ÷ 2 + 1)

			c = x - a - b
			if (a*a + b*b == c*c)
				push!(triplets, (a, b, c))
			end
		end
	end
	return triplets
end


text:
using AbstractStorage, Documenter, FolderStorage

makedocs(sitename="FolderStorage", modules=[FolderStorage])

deploydocs(
    repo = "github.com/ChevronETC/FolderStorage.jl.git",
)


text:
# ------------------------------------------------------------------------------------------
# # Introduction to DataFrames
# **[Bogumił Kamiński](http://bogumilkaminski.pl/about/), Apr 21, 2018**
# ------------------------------------------------------------------------------------------

using DataFrames # load package

# ------------------------------------------------------------------------------------------
# ## Working with CategoricalArrays
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Constructor
# ------------------------------------------------------------------------------------------

x = categorical(["A", "B", "B", "C"]) # unordered

y = categorical(["A", "B", "B", "C"], ordered=true) # ordered, by default order is sorting order

z = categorical(["A","B","B","C", missing]) # unordered with missings

c = cut(1:10, 5) # ordered, into equal counts, possible to rename labels and give custom breaks

by(DataFrame(x=cut(randn(100000), 10)), :x, d -> DataFrame(n=nrow(d)), sort=true) # just to make sure it works right

v = categorical([1,2,2,3,3]) # contains integers not strings

Vector{Union{String, Missing}}(z) # sometimes you need to convert back to a standard vector

# ------------------------------------------------------------------------------------------
# ### Managing levels
# ------------------------------------------------------------------------------------------

arr = [x,y,z,c,v]

isordered.(arr) # chcek if categorical array is orderd

ordered!(x, true), isordered(x) # make x ordered

ordered!(x, false), isordered(x) # and unordered again

levels.(arr) # list levels

unique.(arr) # missing will be included

y[1] < y[2] # can compare as y is ordered

v[1] < v[2] # not comparable, v is unordered although it contains integers

levels!(y, ["C", "B", "A"]) # you can reorder levels, mostly useful for ordered CategoricalArrays

y[1] < y[2] # observe that the order is changed

levels!(z, ["A", "B"]) # you have to specify all levels that are present

levels!(z, ["A", "B"], allow_missing=true) # unless the underlying array allows for missings and force removal of levels

z[1] = "B"
z # now z has only "B" entries

levels(z) # but it remembers the levels it had (the reason is mostly performance)

droplevels!(z) # this way we can clean it up
levels(z)

# ------------------------------------------------------------------------------------------
# ### Data manipulation
# ------------------------------------------------------------------------------------------

x, levels(x)

x[2] = "0"
x, levels(x) # new level added at the end (works only for unordered)

v, levels(v)

v[1] + v[2] # even though underlying data is Int, we cannot operate on it

Vector{Int}(v) # you have either to retrieve the data by conversion (may be expensive)

get(v[1]) + get(v[2]) # or get a single value

get.(v) # this will work for arrays witout missings

get.(z) # but will fail on missing values

Vector{Union{String, Missing}}(z) # you have to do the conversion

z[1]*z[2], z.^2 # the only exception are CategoricalArrays based on String - you can operate on them normally

recode([1,2,3,4,5,missing], 1=>10) # recode some values in an array; has also in place recode! equivalent

recode([1,2,3,4,5,missing], "a", 1=>10, 2=>20) # here we provided a default value for not mapped recodings

recode([1,2,3,4,5,missing], 1=>10, missing=>"missing") # to recode Missing you have to do it explicitly

t = categorical([1:5; missing])
t, levels(t)

recode!(t, [1,3]=>2)
t, levels(t) # note that the levels are dropped after recode

t = categorical([1,2,3], ordered=true)
levels(recode(t, 2=>0, 1=>-1)) # and if you introduce a new levels they are added at the end in the order of appearance

t = categorical([1,2,3,4,5], ordered=true) # when using default it becomes the last level
levels(recode(t, 300, [1,2]=>100, 3=>200))

# ------------------------------------------------------------------------------------------
# ### Comparisons
# ------------------------------------------------------------------------------------------

x = categorical([1,2,3])
xs = [x, categorical(x), categorical(x, ordered=true), categorical(x, ordered=true)]
levels!(xs[2], [3,2,1])
levels!(xs[4], [2,3,1])
[a == b for a in xs, b in xs] # all are equal - comparison only by contents

signature(x::CategoricalArray) = (x, levels(x), isordered(x)) # this is actually the full signature of CategoricalArray
# all are different, notice that x[1] and x[2] are unordered but have a different order of levels
[signature(a) == signature(b) for a in xs, b in xs]

x[1] < x[2] # you cannot compare elements of unordered CategoricalArray

t[1] < t[2] # but you can do it for an ordered one

isless(x[1], x[2]) # isless works within the same CategoricalArray even if it is not ordered

y = deepcopy(x) # but not across categorical arrays
isless(x[1], y[2])

isless(get(x[1]), get(y[2])) # you can use get to make a comparison of the contents of CategoricalArray

x[1] == y[2] # equality tests works OK across CategoricalArrays

# ------------------------------------------------------------------------------------------
# ### Categorical columns in a DataFrame
# ------------------------------------------------------------------------------------------

df = DataFrame(x = 1:3, y = 'a':'c', z = ["a","b","c"])

categorical!(df) # converts all eltype(AbstractString) columns to categorical

showcols(df)

categorical!(df, :x) # manually convert to categorical column :x

showcols(df)

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using Test
using DelimitedFiles
import ClausenFunctions

include("Cl1.jl")
include("Cl2.jl")
include("Cl3.jl")
include("Cl4.jl")
include("Cl5.jl")
include("Cl6.jl")
include("Cl.jl")
include("Sl.jl")


text:
using aITaskDispatcher
using Test

@testset "Generator Tests" begin
    @test isa(aI_generate_00(), Array{Float64,2}) == true
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule GAP_pkg_ferret_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("GAP_pkg_ferret")
JLLWrappers.@generate_main_file("GAP_pkg_ferret", UUID("d6652ff9-bdf9-5f19-9e6f-c10727807222"))
end  # module GAP_pkg_ferret_jll


text:
using VideoIO

# imgstack = []
# for i in 1: 144
# 	img = load("./lotofdata/screenshot_$(i).jpeg")
#     img = convert(Array{Float64}, Gray.(img)) .> 0.5
#     first, last = firstlast(img)
#     center, radius = centerradius(first, last)
# 	img = draw!(RGB.(img), Polygon(RectanglePoints(Point(center[2]-radius-5, center[1]-radius-5), Point(center[2]+radius+5, center[1]+radius+5))), RGB(0,1,0))
#     save("./lotofdata/results/$(i).png", img)
# end

# framestack = map(x->rand(UInt8, 100, 100), 1:100) #vector of 2D arrays
framestack = map(x->load("./lotofdata/screenshot_$(x).jpeg"), 1:1432)
encoder_options = (crf=23, preset="medium")
framerate=24
open_video_out("video.mp4", framestack[1], framerate=framerate, encoder_options=encoder_options) do writer
    for frame in framestack
        write(writer, frame)
    end
end

max_stars_repo_name: ashwani-rathee/Sigil

max_issues_repo_name: ashwani-rathee/Sigil

max_forks_repo_name: ashwani-rathee/Sigil


text:
ROS.init("testTime")
nh = ROS.NodeHandle()
r = ROS.Rate(2)
d = ROS.Duration(2)
t = ROS.Time(2)
ROS.Rate(d)
ROS.reset(r)
ROS.sleep(r)
ROS.cycleTime(r)
ROS.expectedCycleTime(r)

ROS.Duration()
ROS.Duration(r)
ROS.Duration(2,3)
ROS.sleep(d)
@test ROS.toSec(d) == 2

ROS.Time()
ROS.Time(2,3)
@test ROS.toSec(t) == 2
ROS.init(t)
@test !ROS.isSimTime(t)
@test ROS.isSystemTime(t)
@test ROS.isValid(t)
ROS.now(t)
ROS.setNow(t, ROS.now())
ROS.setNow(t, ROS.Time(0))
ROS.sleepUntil(t, ROS.now())
ROS.useSystemTime(t)
ROS.shutdown(t)
t1 = ROS.now()
t.sec = t1.sec + 2
t.nsec = 0
ROS.toSec(t1 + d)
@test ROS.toSec(d + d) == 4
@test ROS.toSec(d - d) == 0
ROS.waitForValid(t)

println("All $(basename(@__FILE__)) tests passed.")


text:
using Mosek

try
    Pkg.installed("Convex")
    using Convex
    using FactCheck

    set_default_solver(MosekSolver(QUIET=true));
    include(joinpath(Pkg.dir("Convex"),"test","test_sdp.jl"))
    FactCheck.exitstatus()
catch
    # If convex not installed, skip this test
end


text:
# Julia wrapper for header: /usr/include/scip/scip_nodesel.h
# Automatically generated using Clang.jl wrap_c

function SCIPincludeNodesel(scip, name, desc, stdpriority, memsavepriority, nodeselcopy, nodeselfree, nodeselinit, nodeselexit, nodeselinitsol, nodeselexitsol, nodeselselect, nodeselcomp, nodeseldata)
    ccall((:SCIPincludeNodesel, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Cstring, Cstring, Cint, Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_NODESELDATA}), scip, name, desc, stdpriority, memsavepriority, nodeselcopy, nodeselfree, nodeselinit, nodeselexit, nodeselinitsol, nodeselexitsol, nodeselselect, nodeselcomp, nodeseldata)
end

function SCIPincludeNodeselBasic(scip, nodesel, name, desc, stdpriority, memsavepriority, nodeselselect, nodeselcomp, nodeseldata)
    ccall((:SCIPincludeNodeselBasic, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{Ptr{SCIP_NODESEL}}, Cstring, Cstring, Cint, Cint, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{SCIP_NODESELDATA}), scip, nodesel, name, desc, stdpriority, memsavepriority, nodeselselect, nodeselcomp, nodeseldata)
end

function SCIPsetNodeselCopy(scip, nodesel, nodeselcopy)
    ccall((:SCIPsetNodeselCopy, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODESEL}, Ptr{Cvoid}), scip, nodesel, nodeselcopy)
end

function SCIPsetNodeselFree(scip, nodesel, nodeselfree)
    ccall((:SCIPsetNodeselFree, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODESEL}, Ptr{Cvoid}), scip, nodesel, nodeselfree)
end

function SCIPsetNodeselInit(scip, nodesel, nodeselinit)
    ccall((:SCIPsetNodeselInit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODESEL}, Ptr{Cvoid}), scip, nodesel, nodeselinit)
end

function SCIPsetNodeselExit(scip, nodesel, nodeselexit)
    ccall((:SCIPsetNodeselExit, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODESEL}, Ptr{Cvoid}), scip, nodesel, nodeselexit)
end

function SCIPsetNodeselInitsol(scip, nodesel, nodeselinitsol)
    ccall((:SCIPsetNodeselInitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODESEL}, Ptr{Cvoid}), scip, nodesel, nodeselinitsol)
end

function SCIPsetNodeselExitsol(scip, nodesel, nodeselexitsol)
    ccall((:SCIPsetNodeselExitsol, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODESEL}, Ptr{Cvoid}), scip, nodesel, nodeselexitsol)
end

function SCIPfindNodesel(scip, name)
    ccall((:SCIPfindNodesel, libscip), Ptr{SCIP_NODESEL}, (Ptr{SCIP_}, Cstring), scip, name)
end

function SCIPgetNodesels(scip)
    ccall((:SCIPgetNodesels, libscip), Ptr{Ptr{SCIP_NODESEL}}, (Ptr{SCIP_},), scip)
end

function SCIPgetNNodesels(scip)
    ccall((:SCIPgetNNodesels, libscip), Cint, (Ptr{SCIP_},), scip)
end

function SCIPsetNodeselStdPriority(scip, nodesel, priority)
    ccall((:SCIPsetNodeselStdPriority, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODESEL}, Cint), scip, nodesel, priority)
end

function SCIPsetNodeselMemsavePriority(scip, nodesel, priority)
    ccall((:SCIPsetNodeselMemsavePriority, libscip), SCIP_RETCODE, (Ptr{SCIP_}, Ptr{SCIP_NODESEL}, Cint), scip, nodesel, priority)
end

function SCIPgetNodesel(scip)
    ccall((:SCIPgetNodesel, libscip), Ptr{SCIP_NODESEL}, (Ptr{SCIP_},), scip)
end


text:
function generate_n_samples(model, grad;
  epsilon = 0.03, # Step size
  L = 11, # No of leapfrog steps
  pr = 0.3, # Plot range
  n_samples = 4, # No of samples
  q = [-0.1, 0.2]) # Initial position

  fig = plot(ylab="muy", xlab="mux", xlim=(-pr, pr), ylim=(-pr, pr), leg=false)

  for i in 1:n_samples
    q, ptraj, qtraj, accept, dH = HMC(model, grad, 0.03, 11, q)
    if n_samples < 10
      for j in 1:L
        K0 = sum(ptraj[j, :] .^ 2) / 2
        plot!(fig, [qtraj[j, 1], qtraj[j+1, 1]], [qtraj[j, 2], qtraj[j+1, 2]], leg=false)
        scatter!(fig, [(qtraj[j, 1], qtraj[j, 2])])
      end
    end
  end
  fig

end


text:
"""
Sample a value For each parameter uniformly at random from the candidate vectors. Log-uniform sampling available by providing a log-spaced candidate vector.
"""
struct RandomSampler <: Sampler end
function (s::RandomSampler)(ho)
    [list[rand(1:length(list))] for list in ho.candidates]
end

# Blue noise ===================================================================

"""
Try to spread out parameters as blue noise (only high frequency). Should sample the space better than random sampling. Use this if you intend to run less than, say, 2000.
"""
@with_kw mutable struct BlueNoiseSampler <: Sampler
    samples = zeros(0,0)
    orders = Int[]
end
function init!(s::BlueNoiseSampler, ho)
    s.samples != zeros(0,0) && return # Already initialized
    dims = length(ho.candidates)
    all(length(c) == ho.iterations for c in ho.candidates) ||
        throw(ArgumentError("BlueNoiseSampler requires all candidate vectors to have the same length as the number of iterations, got lengths $(repr.(collect(zip(ho.params, length.(ho.candidates))))) with $(ho.iterations) iterations"))
    s.samples = bluenoise(dims = dims, nsamples = ho.iterations)
    s.orders = [sortperm(s.samples[dim,:]) for dim = 1:dims]
end

function (s::BlueNoiseSampler)(ho)
    init!(s, ho)
    iter = length(ho.history)+1
    [list[s.orders[dim][iter]] for (dim,list) in enumerate(ho.candidates)]
end

function bluenoise(;
    dims,
    nsamples,
    iters = 100,
    points = rand(dims, nsamples),
    α = 0.005
    )

    n          = nsamples = size(points,2)
    dims       = size(points,1)
    rdmax      = [1/(2n), sqrt(1/(2sqrt(3)*n)), (1/(4sqrt(2)*n))^(1/3), (1/(8n))^(1/4)]
    w          = 1 ./ rdmax
    w        ./= sum(w)
    fadelength = iters ÷ 2
    fadein(i)  = min(i,fadelength)*w[1]/fadelength
    points2    = copy(points)
    grad       = zeros(dims)
    d          = zeros(dims)
    for i = 1:iters
        i == iters ÷ 2 && (α *= 0.2)
        perms = [sortperm(points[dim,:]) for dim = 1:dims]
        for p1 = 1:nsamples
            grad .*= 0
            for dim = 1:dims
                optimal = (-0.5 + perms[dim][p1])/n
                grad[dim] += fadein(i)*sign(optimal - points[dim,p1])
            end
            for p2 = 1:nsamples
                d .= points[:,p1] .- points[:,p2]
                nd = norm(d)
                # nd > 0.2 && continue
                grad .+= w[dims]*d./n/(0.001 + nd)^3
            end
            points2[:,p1] .+= α*grad
        end
        clamp!(points2,0.,1.)
        copyto!(points, points2)
    end
    points
end

# Latin Hypercube Sampler ======================================================

"""
Sample from a latin hypercube
"""
@with_kw mutable struct LHSampler <: Sampler
    samples = zeros(0,0)
    iters = -1
end
function init!(s::LHSampler, ho)
    s.samples != zeros(0,0) && return # Already initialized
    ndims = length(ho.candidates)
    all(length(c) == ho.iterations for c in ho.candidates) ||
        throw(ArgumentError("Latin hypercube sampling requires all candidate vectors to have the same length as the number of iterations, got lengths $(repr.(collect(zip(ho.params, length.(ho.candidates))))) with $(ho.iterations) iterations"))
        s.iters == -1 && (s.iters = (1000*100*2)÷ho.iterations÷ndims)
    X, fit = LHCoptim(ho.iterations,ndims,s.iters)
    s.samples = copy(X')
end

function (s::LHSampler)(ho)
    init!(s, ho)
    iter = length(ho.history)+1
    [list[s.samples[dim,iter]] for (dim,list) in enumerate(ho.candidates)]
end

"""
    CLHSampler(dims=[Continuous(), Categorical(2), ...])
Sample from a categorical/continuous latin hypercube. All continuous variables must have the same length of the candidate vectors.
"""
@with_kw mutable struct CLHSampler <: Sampler
    samples = zeros(0,0)
    dims = []
end
function init!(s::CLHSampler, ho)
    s.samples != zeros(0,0) && return # Already initialized
    all(zip(s.dims, ho.candidates)) do (d,c)
        d isa Categorical || length(c) == ho.iterations
    end || throw(ArgumentError("Latin hypercube sampling requires all candidate vectors for Continuous variables to have the same length as the number of iterations, got lengths $(repr.(collect(zip(ho.params, length.(ho.candidates)))))"))
    ndims = length(ho.candidates)
    initialSample = randomLHC(ho.iterations,s.dims)
    X,_ = LHCoptim!(initialSample, 500, dims=s.dims)
    s.samples = copy(X')
end

function (s::CLHSampler)(ho)
    init!(s, ho)
    iter = length(ho.history)+1
    [list[s.samples[dim,iter]] for (dim,list) in enumerate(ho.candidates)]
end

# GaussianProcesses sampler ====================================================

"""
Sample using Bayesian optimization. `GPSampler(Min)/GPSampler(Max)` fits a Gaussian process to the data and tries to use this model to figure out where the best point to sample next is (using expected improvement). Underneath, the package [BayesianOptimization.jl](https://github.com/jbrea/BayesianOptimization.jl/) is used. We try to provide reasonable defaults for the underlying model and optimizer and we do not provide any customization options for this sampler. If you want advanced control, use BayesianOptimization.jl directly.
"""
@with_kw mutable struct GPSampler <: Sampler
    sense
    model = nothing
    # opt = nothing
    modeloptimizer = nothing
    logdims = nothing
    candidates = nothing
end
GPSampler() = error("The GPSampler needs to know if you want to maximize or minimize. Choose between `GPSampler(Max)/GPSampler(Min)`")
GPSampler(sense) = GPSampler(sense=sense)

function islogspace(x)
    all(x->x > 0, x) || return false
    length(x) > 4    || return false
    std(diff(log.(x))) < sqrt(eps(eltype(x)))
end

function to_logspace(x,logdims)
    map(x,logdims) do x,l
        l ? log10.(x) : x
    end
end

function from_logspace(x,logdims)
    map(x,logdims) do x,l
        l ? exp10.(x) : x
    end
end

function init!(s::GPSampler, ho)
    s.model === nothing || return # Already initialized
    ndims = length(ho.candidates)
    logdims = islogspace.(ho.candidates)
    candidates = to_logspace(ho.candidates, logdims)
    lower_bounds = [minimum.(candidates)...]
    upper_bounds = [maximum.(candidates)...]
    widths = upper_bounds - lower_bounds
    kernel_widths = widths/10
    log_function_noise = 0.
    hypertuning_interval = max(9, ho.iterations ÷ 9)
    model = ElasticGPE(ndims, mean = MeanConst(0.),
                       kernel = SEArd(log.(kernel_widths), log_function_noise), logNoise = 0.)
    modeloptimizer = MLGPOptimizer(every = hypertuning_interval, noisebounds = [-5, 0], # log bounds on the function noise?
                                    maxeval = 40) # max iters for optimization of the GP hyperparams
    # opt = BOpt(f, model, ExpectedImprovement(),
    #            modeloptimizer, lower_bounds, upper_bounds,
    #            maxiterations = ho.iterations, sense = Min, repetitions = 1,
    #            acquisitionoptions = (maxeval = 4000, restarts = 50),
    #            verbosity = Progress)
    # result = boptimize!(opt)

    # s.opt = opt
    s.model = model
    s.modeloptimizer = modeloptimizer
    s.logdims = logdims
    s.candidates = candidates
end

function (s::GPSampler)(ho)
    init!(s, ho)
    iter = length(ho.history)+1
    if iter == 1
        return [rand(list) for (dim,list) in enumerate(ho.candidates)]
    else
        input = reshape(to_logspace(float.(ho.history[end]), s.logdims), :, 1)
        try
            BayesianOptimization.update!(s.model, input, [Int(s.sense)*ho.results[end]])
            BayesianOptimization.optimizemodel!(s.modeloptimizer, s.model) # This determines whether to run or not internally?
        catch ex
            @warn("BayesianOptimization failed at iter $iter: error: ", ex)
        end
    end

    acqfunc = BayesianOptimization.acquisitionfunction(ExpectedImprovement(Int(s.sense)*maximum(s.model.y)), s.model)
    ho2 = Hyperoptimizer(iterations=min(1000, prod(length, s.candidates)), params=ho.params, candidates=s.candidates)
    for params in ho2
        params2 = [params[2:end]...]
        res = -Inf
        try
            res = acqfunc(params2)
        catch ex
            @warn("BayesianOptimization acqfunc failed at iter $iter: error: ", ex)
        end
        push!(ho2.results, res)
        push!(ho2.history, params2)
    end
    return from_logspace(ho2.minimizer, s.logdims)

end

# Hyperband ====================================================================

"""

"""
@with_kw mutable struct Hyperband <: Sampler
    R
    η::Int = 3
    minimum = (Inf,)
    inner = RandomSampler()
end
Hyperband(R) = Hyperband(R=R)

function macrobody(ex, params, candidates, sampler::Hyperband)
    quote
        iters = $(esc(candidates[1]))
        if $sampler.inner isa LHSampler
            smax = floor(Int, log($sampler.η,$sampler.R))
            B = (smax + 1)*$sampler.R
            iters = floor(Int,$sampler.R*$sampler.η^smax)
            ss = string($sampler.inner)
            @info "Starting Hyperband with inner sampler $(ss). Setting the number of iterations to R*η^log(η,R)=$(iters), make sure all candidate vectors have this length as well!"
        end
        ho = Hyperoptimizer(iterations = iters, params = $(esc(params[2:end])), candidates = $(Expr(:tuple, esc.(candidates[2:end])...)), sampler=$sampler)

        costfun = $(Expr(:tuple, esc.(params)...)) -> begin
            $(esc(:(state = nothing)))
            $(esc(ex.args[2]))
        end
        (::$typeof(costfun))($(esc(params[1])), $(esc(:state))) = $(esc(ex.args[2]))
        hyperband($sampler, ho, costfun)
        ho
    end
end

function hyperband(hb::Hyperband, ho, costfun)
    R,η = hb.R, hb.η
    hb.minimum = (Inf,)
    smax = floor(Int, log(η,R))
    B = (smax + 1)*R
    # ho.iterations >= R*η^smax || error("The number of iterations must be larger than R*η^log(η,R) = $(R*η^smax)")
    Juno.progress() do id
        for s in smax:-1:0
            n = ceil(Int, (B/R)*((η^s)/(s+1)))
            r = R / (η^s)
            minᵢ = successive_halving(hb, ho, costfun, n, r, s)
            if minᵢ[1] < hb.minimum[1]
                hb.minimum = minᵢ
            end
            Base.CoreLogging.@logmsg -1 "Hyperband" progress=(smax-s)+1/(smax+1)  _id=id
        end
    end
    return hb.minimum
end

function successive_halving(hb, ho, costfun, n, r=1, s=round(Int, log(hb.η, n)))
    η = hb.η
    minimum = Inf
    T = [ hb.inner(ho) for i=1:n ]
    # append!(ho.history, T)
    Juno.progress() do id
        for i in 0:s
            nᵢ = floor(Int,n/(η^i))
            rᵢ = floor(Int,r*(η^i))
            if i == 0
                LTend = [ costfun(rᵢ, t...) for t in T ]
            else
                LTend = [ costfun(rᵢ, t) for t in T ]
            end
            L, T = first.(LTend), last.(LTend)
            # if i == 0
            #     append!(ho.results, L)
            # else
            #     for t in eachindex(T) # Update results for those that were continued
            #         ho.results[findlast(x->x==T[t], ho.history)] = L[t]
            #     end
            # end
            append!(ho.history, T)
            append!(ho.results, L)

            perm = sortperm(L)
            besti = perm[1]
            if L[besti] < minimum[1]
                minimum = (L[besti], rᵢ, T[besti])
            end
            T = T[perm[1:floor(Int,nᵢ/η)]]
            # T, minimum
            # T, minimum = top_k(hb,T,L,nᵢ,minimum)
            Base.CoreLogging.@logmsg -1 "successive_halving" progress=i/s  _id=id
        end
    end
    return minimum
end

function top_k(hb,T,L,nᵢ,minimum)
    η = hb.η
end


text:
struct CGControl
    abstol :: Real
    reltol :: Real
    maxiter :: Int
    log :: Bool
    verbose :: Bool
    Pl
end

mutable struct EventStreamPredConjGrad{T, L} <: LinPred
    X :: AbstractEventStreamMatrix{T, L}
    beta0 :: Vector{T}
    delbeta :: Vector{T}
    scratchbeta :: Vector{T}
    control :: CGControl
    function EventStreamPredConjGrad{T, L}(E::AbstractEventStreamMatrix{T, L}, beta0::Vector{T}, cntrl::CGControl) where {T, L}
        n, p = size(E)
        length(beta0) == p || throw(DimensionMismatch("length(β0) ≠ size(X,2)"))
        new{T, L}(E, beta0, zeros(eltype(E.δ),p), zeros(eltype(E.δ),p), cntrl)
    end
    function EventStreamPredConjGrad{T, L}(E::AbstractEventStreamMatrix{T, L}, beta0::Vector{T}) where {T, L}
        n, p = size(E)
        length(beta0) == p || throw(DimensionMismatch("length(β0) ≠ size(X,2)"))
        cntrl = CGControl(zero(real(eltype(E.δ))), sqrt(eps(real(eltype(E.δ)))), size(E, 2), false, false, Identity())
        new{T, L}(E.δ, beta0, zeros(eltype(E.δ),p), zeros(eltype(E.\delta),p), cntrl)
    end
    function EventStreamPredConjGrad{T, L}(E::AbstractEventStreamMatrix{T, L}) where {T, L}
        n, p = size(E)
        cntrl = CGControl(zero(real(eltype(E))), sqrt(eps(real(eltype(E)))), size(E, 2), false, false, Identity())
        new{T, L}(E, zeros(eltype(E.δ), p), zeros(eltype(E.δ),p), zeros(eltype(E.δ),p), cntrl)
    end
end

EventStreamPredConjGrad(E :: AbstractEventStreamMatrix) = EventStreamPredConjGrad{eltype(fineness(E)), eltype(E.labels)}(E)

mutable struct DensePredConjGrad{T<:BlasReal} <: DensePred
    X :: Matrix{T}
    beta0 :: Vector{T}
    delbeta :: Vector{T}
    scratchbeta :: Vector{T}
    control :: CGControl
    intercept :: Bool
    function DensePredConjGrad{T}(X::Matrix{T}, beta0::Vector{T}, cntrl::CGControl) where T
        n, p = size(X)
        length(beta0) == p || throw(DimensionMismatch("length(β0) ≠ size(X,2)"))
        new{T}(X, beta0, zeros(eltype(X),p), zeros(eltype(X),p), cntrl, true)
    end
    function DensePredConjGrad{T}(X::Matrix{T}, beta0::Vector{T}) where T
        n, p = size(X)
        length(beta0) == p || throw(DimensionMismatch("length(β0) ≠ size(X,2)"))
        cntrl = CGControl(zero(real(eltype(X))), sqrt(eps(real(eltype(X)))), size(X, 2), false, false, Identity())
        new{T}(X, beta0, zeros(eltype(X),p), zeros(eltype(X),p), cntrl, true)
    end
    function DensePredConjGrad{T}(X::Matrix{T}) where T
        n, p = size(X)
        cntrl = CGControl(zero(real(eltype(X))), sqrt(eps(real(eltype(X)))), size(X, 2), false, false, Identity())
        new{T}(X, zeros(eltype(X), p), zeros(eltype(X),p), zeros(eltype(X),p), cntrl, true)
    end
end

DensePredConjGrad(X::Matrix) = DensePredConjGrad{eltype(X)}(X)

#TODO: Make these efficient
cholesky!(p::DensePredConjGrad) = cholesky(p.X'*p.X)
cholesky(p::DensePredConjGrad) = cholesky(p.X'*p.X)

cholesky(p :: EventStreamPredConjGrad) = Identity() #cholesky(Hermitian(XtWX(p.X, ones(p.X.nbins))))
cholesky!(p :: EventStreamPredConjGrad) = Identity() #cholesky(Hermitian(XtWX(p.X, ones(p.X.nbins))))

# Standard errors are way off
invchol(p :: EventStreamPredConjGrad) = inv(cholesky(p))

function linpred!(out, p :: EventStreamPredConjGrad, f::Real=1.0)
    out[:] = XWb(p.X, ones(size(p.X)[2]), p.beta0 + p.delbeta)
end

function delbeta!(p :: DensePredConjGrad{T}, r::Vector{T}) where T<: BlasReal
    delbeta!(p, r, ones(eltype(r), length(r)))
end

function delbeta!(p :: EventStreamPredConjGrad{T, L}, r :: Vector{T}, wt :: Vector{T}) where {T, L}
    G = WeightedGramMatrix(p.X, wt)
    p.delbeta = cg(G, XtWy(p.X, wt, r),
        abstol=p.control.abstol,
        reltol=p.control.reltol,
        maxiter=p.control.maxiter,
        log=p.control.log,
        verbose=p.control.verbose,
        Pl=p.control.Pl
    )
end

function delbeta!(p :: DensePredConjGrad{T}, r :: Vector{T}, wt::Vector{T}) where T <: BlasReal
    G = WeightedNormGramMatrix(p.X, wt)
    p.delbeta = cg(G, G.X' * Diagonal(wt) * r ;
        abstol=p.control.abstol,
        reltol=p.control.reltol,
        maxiter=p.control.maxiter,
        log=p.control.log,
        verbose=p.control.verbose,
        Pl=p.control.Pl
        )
    return p
end


text:
using SimulacoesTemplate
using Test

@testset "Distance function" begin

   # Testing rpbc
   x = [ 0, 0 ]
   y = [ 12, 12 ]
   side = 10
   @test SimulacoesTemplate.rpbc(x,y,side) ≈ sqrt(8)
   y = [ 17, 17 ]
   @test SimulacoesTemplate.rpbc(x,y,side) ≈ sqrt(18)

end

@testset "Potential Energy" begin

  data = Data()

  @test SimulacoesTemplate.upair(π,data) ≈ -0.0003250451775148634

  p = initial_point(data.N, data.side, 0.9)
  @test utotal(p,data) ≈ -7472.164849146468

end


text:
    function estimate_numberofeigenvalues(A,x,radius,origin,Z2randomize!,randomize!;L0=10,Nq=64,M=20,κ=2,δ=1e-14,MaxCGstep=3000,epsCG=1e-16,α=1)
        γ = origin
        ρ = radius

        println("Estimating number of eigenvalues...")
        println("------------------------------------------------------------")
        xtype = typeof(x)
        V = Vector{xtype}(undef,L0)
        for i=1:L0
            V[i] = zero(x)
            Z2randomize!(V[i])
        end

        Sk = Vector{xtype}(undef,L0*M)
        for i=1:L0*M
            Sk[i] = zero(x)
        end

        Y = Vector{xtype}(undef,Nq)
        for j=1:Nq
            Y[j] = zero(x)
        end

        evaluate_Eq22!(Sk,V,A,x,radius,origin,α,L0,Nq,M,MaxCGstep,epsCG)

        S0 = Vector{xtype}(undef,L0)
        for j=1:L0
            S0[j] = zero(x)
            axpy!(1,Sk[j],S0[j])
        end

        mst = evaluate_Eq26(S0,V,L0)

        println("done.")
        println("------------------------------------------------------------")
        println("Estimated number of eigenvalues: ",mst)

        return mst
    end

    function eigensystem(A,x,radius,origin,Z2randomize!,randomize!;L0=10,Nq=64,M=10,κ=2,δ=1e-14,MaxCGstep=3000,epsCG=1e-16,α=1)
        mst = estimate_numberofeigenvalues(A,x,radius,origin,Z2randomize!,randomize!,L0=L0,Nq=Nq,M=M,κ=κ,δ=δ,MaxCGstep=MaxCGstep,epsCG=epsCG,α=α)
        xtype = typeof(x)
        Nrh = Int64(round(κ*mst/M,RoundNearestTiesAway))
        if Nrh < 2
            Nrh = 2
        end
        Nint = Nq
        Nmm = 2*M
        N = length(x)

        z = zeros(ComplexF64,Nint)
        w = zero(z)

        for j = 1:Nint
            θ = (2π/Nint)*(j-1/2)
            z[j] = origin + radius*(cos(θ)+im*α*sin(θ))
            w[j] = α*cos(θ) + im*sin(θ)
        end

        V = Vector{xtype}(undef,Nrh)
        for i=1:Nrh
            V[i] = zero(x)
            #for j=1:N
            #    V[i][j] = rand()
            #end
            randomize!(V[i])
        end

        Sk = Vector{xtype}(undef,Nrh*Nmm*2)
        for i=1:Nrh*Nmm*2
            Sk[i] = zero(x)
        end

        #S = zeros(ComplexF64,N,Nrh*Nmm*2)

        for i=1:Nrh
            println("i = ",i,"/",Nrh)
            #(zj I - H) Yj = V
            #(-zj I + H) (-Yj) = V
            Y = shiftedbicg_2003(-z,A,V[i],maxsteps = MaxCGstep,eps=epsCG,verboselevel = 2)

            for k in 0:Nmm*2-1
                ik = k*Nrh+i
                for j in 1:Nint
                    coeff = radius*w[j]*z[j]^k/Nint
                    axpy!(-coeff,Y[j],Sk[ik])
                    #vec_Sk[:,ik] += ρ*vec_w[j]*vec_z[j]^k*vec_y[:,j]/Nq
                end
            end
        end

        μ = zeros(ComplexF64,Nrh,Nrh,2*Nmm)
        S = zeros(ComplexF64,N,Nrh*Nmm*2)
        for i=1:Nrh*Nmm*2
            S[:,i] = convert_to_normalvector(Sk[i])
        end

        Vmatrix = zeros(ComplexF64,N,Nrh)
        for i=1:Nrh
            Vmatrix[:,i] = convert_to_normalvector(V[i])
        end

        μ = zeros(ComplexF64,Nrh,Nrh,2*Nmm)
        for k=0:Nmm*2-1
            μ[:,:,k+1] = Vmatrix'*S[:,1+k*Nrh:Nrh+k*Nrh]
        end
        Tv = zeros(ComplexF64,Nmm*Nrh,Nmm*Nrh)
        Tt = zeros(ComplexF64,Nmm*Nrh,Nmm*Nrh)

        for i=1:Nmm
            for j=1:Nmm
                Tv[1+(i-1)*Nrh:Nrh+(i-1)*Nrh,1+(j-1)*Nrh:Nrh+(j-1)*Nrh] = μ[:,:,i+j-1+1]
                Tt[1+(i-1)*Nrh:Nrh+(i-1)*Nrh,1+(j-1)*Nrh:Nrh+(j-1)*Nrh] = μ[:,:,i+j-2+1]
            end
        end

        F = svd(Tt)
        U = F.U
        Σ = F.S
        W = F.V

        eps = 1e-8

        δ = 1e-10
        j = 0
        ratio = 1
        while ratio > δ
            j += 1
            ratio = Σ[j]/Σ[1]
            if j == length(Σ)
                ratio = 0.0
            end
        end
        ms = j
        Σinv = zeros(ms,ms)
        for i=1:ms
            Σinv[i,i] = 1/Σ[i]
        end

        H = U[:,1:ms]'*Tv*W[:,1:ms]*Σinv
        ε,v = eigen(H)
        println(ε)

        vec_xvec = zeros(ComplexF64,N,ms)
        for is=1:ms
            for i=1:Nmm
                a = Σinv*v[:,is]
                b = W[:,1:ms]*a
                c = S[:,1:Nrh*Nmm]*b
                vec_xvec[:,is] += c
            end
        end

        vec_x = Array{typeof(x),1}(undef,ms)
        for i=1:ms
            vec_x[i] = zero(x)
            for j=1:N
                vec_x[i][j] = vec_xvec[j,i]
            end
        end

        Ax = similar(x)
        r = similar(x)

        is = 0
        eigenvalues = ComplexF64[]
        residuals = ComplexF64[]
        eigenvectors = []

        for i in 1:ms
            rr = real( ε[i]-origin)^2+imag( ε[i]-origin)^2
            if rr > radius^2
                continue
            end

            mul!(Ax,A,vec_x[i])
            add!(0,r,-ε[i],vec_x[i])
            add!(r,1,Ax)

            residual = sqrt(real(dot(r,r)))/sqrt(real(dot(vec_x[i],vec_x[i])))
            #println(residual)
            println("$(ε[i]) residual: $residual")

            if residual  < eps
                println("$(ε[i]) residual: $residual")
                is += 1
                push!(eigenvalues,ε[i])
                push!(residuals,residual)
                push!(eigenvectors,vec_x[i])
            end

        end

        #println("-----------------------------------------")
        println("done.")
        println("------------------------------------------------------------")
        return eigenvalues,residuals,eigenvectors,is

    end

    function eigensystem_old2(A,x,radius,origin,Z2randomize!,randomize!;L0=10,Nq=64,M=20,κ=2,δ=1e-14,MaxCGstep=3000,epsCG=1e-16,α=1)

        γ = origin
        ρ = radius

        println("Estimating number of eigenvalues...")
        println("------------------------------------------------------------")
        xtype = typeof(x)
        V = Vector{xtype}(undef,L0)
        for i=1:L0
            V[i] = zero(x)
            Z2randomize!(V[i])
        end

        Sk = Vector{xtype}(undef,L0*M)
        for i=1:L0*M
            Sk[i] = zero(x)
        end

        Y = Vector{xtype}(undef,Nq)
        for j=1:Nq
            Y[j] = zero(x)
        end

        evaluate_Eq22!(Sk,V,A,x,radius,origin,α,L0,Nq,M,MaxCGstep,epsCG)

        S0 = Vector{xtype}(undef,L0)
        for j=1:L0
            S0[j] = zero(x)
            axpy!(1,Sk[j],S0[j])
        end

        mst = evaluate_Eq26(S0,V,L0)

        println("done.")
        println("------------------------------------------------------------")
        println("Estimated number of eigenvalues: ",mst)
        L = Int64(round(κ*mst/M,RoundNearestTiesAway))
        if L < 2
            L = 2
        end
        println("L is ",L)

        println("Calculating eigenvalues and eigenvectors...")
        println("------------------------------------------------------------")

        V = Vector{xtype}(undef,L)
        for i=1:L
            V[i] = zero(x)
            randomize!(V[i])
        end

        ε,vec_x,ms = evaluate_SS(Sk,V,A,x,radius,origin,α,L,Nq,M,MaxCGstep,epsCG,δ)
        println(ε)
        Ax = similar(x)
        r = similar(x)

        is = 0
        eigenvalues = ComplexF64[]
        residuals = ComplexF64[]
        eigenvectors = []

        for i in 1:ms
            mul!(Ax,A,vec_x[i])
            add!(0,r,-ε[i],vec_x[i])
            add!(r,1,Ax)

            resiup = sqrt(real(dot(r,r)))
            println(resiup)

            if resiup < 0.001  && γ-ρ <= real(ε[i]) <= ρ+γ
                is += 1
                push!(eigenvalues,ε[i])
                push!(residuals,resiup)
                push!(eigenvectors,vec_x[i])
            end

            #=
            residown1 = sqrt(dot(Ax,Ax))
            residown2 = sqrt(dot(vec_x[i],vec_x[i]))
            resi = real(resiup/(residown1+abs(ε[i])*residown2))
            if resi < 0.1 # && γ-ρ <= ε[i] <= ρ+γ
                is += 1
                push!(eigenvalues,ε[i])
                push!(residuals,resi)
                push!(eigenvectors,vec_x[i])
            end
            =#

        end

        #println("-----------------------------------------")
        println("done.")
        println("------------------------------------------------------------")
        return eigenvalues,residuals,eigenvectors,is

    end

    function evaluate_Eq22!(Sk,V,A,x,radius,origin,α,L,Nq,M,MaxCGstep,epsCG)
        vec_z = zeros(ComplexF64,Nq)
        vec_w = zeros(ComplexF64,Nq)
        for j in 1:Nq
            θj = (2π/Nq)*(j-1/2)
            vec_z[j] = origin + radius*(cos(θj)+im*α*sin(θj))
            vec_w[j] = α*cos(θj)+im*sin(θj)
        end

        for i in 1:L
            println("i = ",i,"/",L)
            #(zj I - H) Yj = V
            #(-zj I + H) (-Yj) = V
            Y = shiftedbicg_2003(-vec_z,A,V[i],maxsteps = MaxCGstep,eps=epsCG,verboselevel = 2)

            for k in 0:M-1
                ik = k*L+i
                for j in 1:Nq
                    coeff = radius*vec_w[j]*vec_z[j]^k/Nq
                    axpy!(-coeff,Y[j],Sk[ik])
                    #vec_Sk[:,ik] += ρ*vec_w[j]*vec_z[j]^k*vec_y[:,j]/Nq
                end
            end

        end
    end

    function evaluate_Eq26(S0,V,L0)
        m = 0.0im
        for i in 1:L0
            m += dot(V[i],S0[i])
            #m += hatV[:,i]'*vec_s0[:,i]
        end
        ms = Int64(round(real(m)/L0))
        if ms < 0
            ms = 1
        end
        return ms
    end

    function evaluate_SS(Sk,V,A,x,radius,origin,α,L,Nq,M,MaxCGstep,epsCG,δ)
        xvec = convert_to_normalvector(x)
        N = length(xvec)

        evaluate_Eq22!(Sk,V,A,x,radius,origin,α,L,Nq,M,MaxCGstep,epsCG)
        Skmatrix = zeros(ComplexF64,N,L*M)
        for i=1:L*M
            Skmatrix[:,i] = convert_to_normalvector(Sk[i])
        end
        Vmatrix = zeros(ComplexF64,N,L)
        for i=1:L
            Vmatrix[:,i] = convert_to_normalvector(V[i])
        end
        μk = zeros(ComplexF64,L,L,M)
        for k=0:M-1
            μk[:,:,k+1] = Vmatrix'*Skmatrix[:,1+k*L:L+k*L]
            #μk[:,:,i] = Vmatrix'*Skmatrix[:,1+(i-1)*L:L+(i-1)*L]
        end
        M2 = div(M,2)
        T1matrix = zeros(ComplexF64,L*M2,L*M2)
        T2matrix = zeros(ComplexF64,L*M2,L*M2)
        for i=1:M2
            for j=1:M2
                T1matrix[1+(i-1)*L:L+(i-1)*L,1+(j-1)*L:L+(j-1)*L] = μk[:,:,i+j-1+1]
                T2matrix[1+(i-1)*L:L+(i-1)*L,1+(j-1)*L:L+(j-1)*L] = μk[:,:,i+j-2+1]
            end
        end
        F =  svd(T2matrix)
        U = F.U
        Σ = F.S
        W = F.V
        #(U,Σ,W ) = svd(T2matrix)
        println( Σ)
        j = 0
        ratio = 1
        while ratio > δ
            j += 1
            ratio = Σ[j]/Σ[1]
            if j == length(Σ)
                ratio = 0.0
            end
        end
        ms = j
        println("ms = $ms")
        Σmatrixinv = zeros(ComplexF64,ms,ms)
        for i=1:ms
            Σmatrixinv[i,i] = 1/Σ[i]
        end
        H = U[:,1:ms]'*T1matrix*W[:,1:ms]*Σmatrixinv

        ε,vec_w = eigen(H)
        vec_xvec = zeros(ComplexF64,N,ms)
        for is=1:ms
            for i=1:M2
                a = Σmatrixinv*vec_w[:,is]
                b = W[:,1:ms]*a
                c = Skmatrix[:,1:L*M2]*b
                vec_xvec[:,is] += c
                #vec_xvec[:,is] += Skmatrix[:,1+(i-1)*L:L+(i-1)*L]*W[:,1:ms]*Σmatrixinv*vec_w[:,is]
            end
        end

        vec_x = Array{typeof(x),1}(undef,ms)
        for i=1:ms
            vec_x[i] = zero(x)
            for j=1:N
                vec_x[i][j] = vec_xvec[j,i]
            end
        end

        return ε,vec_x,ms

    end

function eigensystem_old(A,x,N,ρ,γ;α=0.1,L0=10,Nq=64,M=10,κ=2,δ=1e-14,numr = 0,cg=true,MaxCGstep=3000,epsCG=1e-16)
        println("Estimating number of eigenvalues...")
        println("------------------------------------------------------------")

        V = Array{typeof(x),1}(undef,L0)
        for i=1:L0
            temp = rand([-1,1],N)
            V[i] = similar(x)
            for j=1:N
                V[i][j] = temp[j]
            end
            #=
            if typeof(x) <: Abstractfermion
                V[i] = similar(x)
                Z2_distribution_fermion!(V[i])
            else
                V[i] = rand([-1,1],N)
            end
            =#
        end

        Sk = Array{typeof(x),1}(undef,L0*M)
        for i=1:L0*M
            Sk[i] = similar(x)
        end

        Y = Array{typeof(x),1}(undef,Nq)
        for j=1:Nq
            Y[j] = similar(x)
        end

        #V = rand([-1,1],N,L0)
        calculate_Eq22!(Sk,Y,V,A,x,N,ρ,γ,α,L0,Nq,M,cg,MaxCGstep,epsCG)

        S0 = Array{typeof(x),1}(undef,L0)
        for i=1:L0
            S0[i] = similar(x)
            if typeof(x) <: Abstractfermion
                add!(0,S0[i],1,Sk[i]) #S0[i] = S[i]
            else
                S0[i] .= Sk[i]
            end
        end

        mst = calc_msinEq26(S0,V,L0)

        println("done.")
        println("------------------------------------------------------------")
        println("Estimated number of eigenvalues: ",mst)
        L = Int64(round(κ*mst/M,RoundNearestTiesAway))
        if L < 2
            L = 2
        end
        println("L is ",L)

        println("Calculating eigenvalues and eigenvectors...")
        println("------------------------------------------------------------")

        V = Array{typeof(x),1}(undef,L)
        for i=1:L
            temp = rand(N).*2 .- 1
            V[i] = similar(x)
            for j=1:N
                V[i][j] = temp[j]
            end

            #=
            if typeof(x) <: Abstractfermion
                V[i] = similar(x)
                uniform_distribution_fermion!(V[i])
            else
                V[i] = rand(N).*2 .- 1
            end
            =#

        end

        Sk = Array{typeof(x),1}(undef,L*M)
        for i=1:L*M
            Sk[i] = similar(x)
        end

        ε,vec_x,ms = calc_SS(Sk,Y,V,A,x,N,ρ,γ,α,L,Nq,M,δ,cg,MaxCGstep,epsCG)
        println(ε)
        Ax = similar(x)
        r = similar(x)

        is = 0
        eigenvalues = ComplexF64[]
        residuals = ComplexF64[]
        eigenvectors = []

        for i in 1:ms
            mul!(Ax,A,vec_x[i])
            add!(0,r,-ε[i],vec_x[i])
            add!(r,1,Ax)

            resiup = sqrt(real(dot(r,r)))
            println(resiup)

            if resiup < 0.001  && γ-ρ <= real(ε[i]) <= ρ+γ
                is += 1
                push!(eigenvalues,ε[i])
                push!(residuals,resiup)
                push!(eigenvectors,vec_x[i])
            end

            #=
            residown1 = sqrt(dot(Ax,Ax))
            residown2 = sqrt(dot(vec_x[i],vec_x[i]))
            resi = real(resiup/(residown1+abs(ε[i])*residown2))
            if resi < 0.1 # && γ-ρ <= ε[i] <= ρ+γ
                is += 1
                push!(eigenvalues,ε[i])
                push!(residuals,resi)
                push!(eigenvectors,vec_x[i])
            end
            =#

        end

        #println("-----------------------------------------")
        println("done.")
        println("------------------------------------------------------------")
        return eigenvalues,residuals,eigenvectors,is

    end

    function calc_msinEq26(S0,V,L0)
        m = 0.0im
        for i in 1:L0
            m += dot(V[i],S0[i])
        end
        #m = real(m)/L0
        #println(m)
        ms = Int64(round(real(m)/L0))
        if ms < 0
            ms = 1
        end

        return ms
    end

    function calc_SS(Sk,Y,mV,A,x,N,ρ,γ,α,L,Nq,M,δ,cg,MaxCGstep,epsCG)
        calculate_Eq22!(Sk,Y,mV,A,x,N,ρ,γ,α,L,Nq,M,cg,MaxCGstep,epsCG)

        Skmatrix = zeros(ComplexF64,N,L*M)
        if typeof(x) <: Abstractfermion
            fermions2vectors!(Skmatrix,Sk)
        else
            for i=1:L*M
                for j=1:N
                    Skmatrix[j,i] = Sk[i][j]
                end
            end
        end

        (U,Σ,W ) = svd(Skmatrix)
        j = 0
        ratio = 1
        while ratio > δ
            j += 1
            ratio = Σ[j]/Σ[1]
            if j == length(Σ)
                ratio = 0.0
            end
        end

        ms = j
        Q = Array{typeof(x),1}(undef,ms)
        for i=1:ms
            Q[i] = similar(x)
            for j=1:N
                Q[i][j]=U[j,i]
            end
            if typeof(x) <: Abstractfermion
                set_wing_fermion!(Q[i])
            end
        end

        temp = Array{typeof(x),1}(undef,ms)
        for i=1:ms
            temp[i] = similar(x)
            mul!(temp[i],A,Q[i])
        end

        mat_Ht =zeros(ComplexF64,ms,ms)
        for i=1:ms
            for j=1:ms
                mat_Ht[j,i] = dot(Q[j],temp[i])
            end
        end

        #mat_Ht = (mat_Ht'+mat_Ht)/2

        ε,vec_w = eigen(mat_Ht)
        vec_x = Array{typeof(x),1}(undef,ms)
        for i=1:ms
            vec_x[i] = similar(x)
            if typeof(x) <: Abstractfermion
            else
                vec_x[i] .= 0
            end

            for j=1:N
                for k=1:ms
                    vec_x[i][j] += Q[k][j]*vec_w[k,i]
                end
            end

            if typeof(x) <: Abstractfermion
                set_wing_fermion!(vec_x[i])
            end
        end
        #vec_x = mat_Q*vec_w

        return ε,vec_x,ms
    end

    function calculate_Eq22!(Sk,Y,mV,A,x,N,ρ,γ,α,L,Nq,M,cg,MaxCGstep,epsCG) #mV = -V
        if typeof(x) <: Abstractfermion
            for i in 1:L
                for k = 0:M-1
                    ik = k*L+i
                    clear_fermion!(Sk[ik])
                end
            end
        else
            for i in 1:L
                for k = 0:M-1
                    ik = k*L+i
                    Sk[ik] .= 0
                end
            end
        end

        z = zeros(ComplexF64,Nq)
        w = zero(z)

        for j = 1:Nq
            θ = (2π/Nq)*(j-1/2)
            z[j] = γ + ρ*(cos(θ)+im*α*sin(θ))
            w[j] = α*cos(θ) + im*sin(θ)
        end

        #(zI - A) Y = V
        #(A - z I) (-Y) = V
        for i in 1:L
            #shiftedcg(Y,-z,x,A,mV[i])
            Y = shiftedbicg_2003(-z,A,mV[i],maxsteps = MaxCGstep,eps=epsCG,verboselevel = 1)
            #shiftedbicg(Y,-z,x,A,mV[i],maxsteps = MaxCGstep,eps=epsCG)
            for k = 0:M-1
                ik = k*L+i
                for j = 1:Nq
                    coef = ρ*w[j]*z[j]^k/Nq
                    add!(Sk[ik],-coef,Y[j])
                end
            end
        end

    end


text:
module ModZ

using CustomUnitRanges: filename_for_zerorange
include(filename_for_zerorange)

end

using ModZ: ZeroRange

@testset "ZeroRange" begin
    r = ZeroRange(-5)
    @test isempty(r)
    @test length(r) == 0
    @test size(r) == (0,)
    r = ZeroRange(3)
    @test !isempty(r)
    @test length(r) == 3
    @test size(r) == (3,)
    @test step(r) == 1
    @test first(r) == 0
    @test last(r) == 2
    @test minimum(r) == 0
    @test maximum(r) == 2
    @test r[1] == 0
    @test r[2] == 1
    @test r[3] == 2
    @test_throws BoundsError r[4]
    @test_throws BoundsError r[-1]
    @test r[1:3] === 0:2
    @test_throws BoundsError r[r]
    @test r+1 === 1:3
    @test 2*r === 0:2:4
    k = -1
    for i in r
        @test i == (k+=1)
    end
    @test k == length(r)-1
    @test intersect(r, ZeroRange(2)) === intersect(ZeroRange(2), r) === ZeroRange(2)
    @test intersect(r, -1:5) === intersect(-1:5, r) === 0:2
    @test intersect(r, 2:5) === intersect(2:5, r) === 2:2
    io = IOBuffer()
    show(io, r)
    str = takebuf_string(io)
    @test str == "ModZ.ZeroRange(3)"

    r = ZeroRange(5)
    @test checkindex(Bool, r, 4)
    @test !checkindex(Bool, r, 5)
    @test checkindex(Bool, r, :)
    @test checkindex(Bool, r, 1:4)
    @test !checkindex(Bool, r, 1:5)
    @test !checkindex(Bool, r, trues(4))
    @test !checkindex(Bool, r, trues(5))
    @test convert(UnitRange, r) == 0:4
    @test convert(StepRange, r) == 0:1:4
    @test !in(-1, r)
    @test in(0, r)
    @test in(4, r)
    @test !in(5, r)
    @test issorted(r)
    @test maximum(r) == 4
    @test minimum(r) == 0
    @test sortperm(r) == 1:5
    @test r == 0:4
    @test r+r == 0:2:8
    @test (5:2:13)-r == 5:9
    @test -r == 0:-1:-4
    @test reverse(r) == 4:-1:0
    @test r./2 == 0:0.5:2

    r = ZeroRange{Int16}(5)
    @test length(r) === 5
    @test start(r) === 0
    k = -1
    for i in r
        @test i == (k+=1)
    end
    @test k == length(r)-1
    x, y = promote(ZeroRange(5), ZeroRange{Int16}(8))
    @test x === ZeroRange(5)
    @test y === ZeroRange(8)
    x, y = promote(ZeroRange(5), 0:7)
    @test x === 0:4
    @test y === 0:7
    @test convert(ZeroRange{Int16}, ZeroRange(5)) === ZeroRange{Int16}(5)
    @test convert(ZeroRange{Int}, ZeroRange(5)) === ZeroRange(5)
    @test convert(UnitRange, ZeroRange(4)) === 0:3
    r = ZeroRange(Int128(10))
    @test length(r) === Int128(10)
end


text:
using BoardGames
using Test

include("TicTacToe.jl")
include("RandomStrategy.jl")

b = TicTacToeBoard()
b = play(b, 5) #x

function somemoves()
    b = TicTacToeBoard()
    b = play(b, 5) #x
    b = play(b, 1) #o
    b = play(b, 2) #x
    b = play(b, 7) #o
    b = play(b, 8) #x

    b
end

@testset "BoardGames.jl" begin
    # Write your tests here.

    @test b.v == [0,0,0,0,1,0,0,0,0]
    @test getmoves(b) == [1,2,3,4,6,7,8,9]
    @test somemoves().winner == 1

    @test typeof(playgame(TicTacToe(), RandomStrategy(), RandomStrategy())) ==
          TicTacToeBoard
    @test boardtype(TicTacToe) == TicTacToeBoard
    @test movetype(TicTacToe) == Int

end


text:
# ---
# title: 1662. Check If Two String Arrays are Equivalent
# id: problem1662
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/check-if-two-string-arrays-are-equivalent/description/>
# hidden: true
# ---
#
# Given two string arrays `word1` and `word2`, return __`true` _if the two
# arrays **represent** the same string, and _`false` _otherwise._
#
# A string is **represented** by an array if the array elements concatenated
# **in order** forms the string.
#
#
#
# **Example 1:**
#
#
#
#     Input: word1 = ["ab", "c"], word2 = ["a", "bc"]
#     Output: true
#     Explanation:
#     word1 represents string "ab" + "c" -> "abc"
#     word2 represents string "a" + "bc" -> "abc"
#     The strings are the same, so return true.
#
# **Example 2:**
#
#
#
#     Input: word1 = ["a", "cb"], word2 = ["ab", "c"]
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input: word1  = ["abc", "d", "defg"], word2 = ["abcddefg"]
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * `1 <= word1.length, word2.length <= 103`
#   * `1 <= word1[i].length, word2[i].length <= 103`
#   * `1 <= sum(word1[i].length), sum(word2[i].length) <= 103`
#   * `word1[i]` and `word2[i]` consist of lowercase letters.
#
#
## @lc code=start
using LeetCode

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


text:
"""
Binding state that must be set in order for
drawing commands to render correctly.
"""
struct BindRequirements
    pipeline::Pipeline
    push_data::Any
end

"""
Describes the current binding state.
"""
struct BindState
    pipeline::Optional{Pipeline}
    push_data::Any
end

BindState() = BindState(nothing, nothing)

function Base.bind(cbuffer::CommandBuffer, reqs::BindRequirements, state::BindState = BindState())
    (;pipeline, push_data) = reqs
    current_push_data = state.push_data

    if pipeline ≠ state.pipeline
        Vk.cmd_bind_pipeline(cbuffer, Vk.PipelineBindPoint(pipeline.type), pipeline)
        if isnothing(state.pipeline) || pipeline.layout.push_constant_ranges ≠ state.pipeline.layout.push_constant_ranges
            # push data gets invalidated
            current_push_data = nothing
        end
    end

    if !isnothing(push_data) && push_data ≠ current_push_data
        ref = Ref(push_data)
        GC.@preserve ref Vk.cmd_push_constants(cbuffer, pipeline.layout, Vk.SHADER_STAGE_ALL, 0, sizeof(push_data), Base.unsafe_convert(Ptr{Cvoid}, ref))
    end
end


text:
"""
This module provides:

- An interface for physical fields, basis of physical fields and related objects.
- Helpers functions to work with fields and arrays of fields.
- Helpers functions to create lazy operation trees from fields and arrays of fields

The exported names are:

$(EXPORTS)
"""
module Fields

using InplaceArrays.Helpers
using InplaceArrays.Inference
using InplaceArrays.TensorValues
using InplaceArrays.Arrays
using InplaceArrays.Arrays: BCasted
using InplaceArrays.Arrays: NumberOrArray
using InplaceArrays.Arrays: AppliedArray
using InplaceArrays.Arrays: Contracted

using Test
using DocStringExtensions
using FillArrays

export Point
export field_gradient
export evaluate_field!
export evaluate_field
export field_cache
export field_return_type
export evaluate
export evaluate!
export gradient
export ∇
export Field
export test_field
export apply_kernel_to_field
export apply_to_field_array
export test_array_of_fields
export compose
export lincomb
export varinner
export attachmap
export integrate
export field_caches
export field_return_types
export evaluate_fields
export evaluate_fields!
export field_gradients
export field_array_cache
export evaluate_field_array
export field_array_gradient
export gradient_type
export curl
export grad2curl
export laplacian
export divergence
export Δ
export ε

import InplaceArrays.Arrays: kernel_cache
import InplaceArrays.Arrays: apply_kernel!
import InplaceArrays.Arrays: kernel_return_type
import InplaceArrays.TensorValues: outer
import InplaceArrays.TensorValues: inner
import InplaceArrays.TensorValues: symmetic_part
import Base: +, - , *
import LinearAlgebra: cross
import LinearAlgebra: tr
import Base: transpose
import Base: adjoint

include("FieldInterface.jl")

include("MockFields.jl")

include("ConstantFields.jl")

include("FieldApply.jl")

include("FieldArrays.jl")

include("Lincomb.jl")

include("Compose.jl")

include("Varinner.jl")

include("Attachmap.jl")

include("Integrate.jl")

include("FieldOperations.jl")

include("DiffOperators.jl")

end # module


text:
using CUDA
using Random
using Test
using VUMPS
using VUMPS: vumps, vumps_env

@testset "$(Ni)x$(Nj) VUMPSRuntime with $atype{$dtype}" for Ni = [1,2,3], Nj = [1,2,3], atype = [Array, CuArray], dtype = [Float64, ComplexF64]
    @test SquareLattice <: AbstractLattice

    M = Array{atype{dtype,4},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        M[i,j] = atype(rand(dtype,2,2,2,2))
    end
    rt = SquareVUMPSRuntime(M, Val(:random), 2)
end

@testset "$(Ni)x$(Nj) vumps with $atype{$dtype}" for Ni = [2], Nj = [2], atype = [Array, CuArray], dtype = [Float64, ComplexF64]
    M = Array{atype{dtype,4},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        M[i,j] = atype(rand(dtype,2,2,2,2))
    end
    rt = SquareVUMPSRuntime(M, Val(:random), 2)
    env = vumps(rt)
    @test env !== nothing
end

@testset "$(Ni)x$(Nj) vumps_env with $atype{$dtype}" for Ni = [2], Nj = [2], atype = [Array, CuArray], dtype = [Float64, ComplexF64]
    M = Array{atype{dtype,4},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        M[i,j] = atype(rand(dtype,2,2,2,2))
    end
    env = vumps_env(M; χ = 10)
    @test env !== nothing
end

@testset "$(Ni)x$(Nj) obs_env" for Ni = [2], Nj = [2], atype = [Array, CuArray], dtype = [Float64, ComplexF64]
    M = Array{atype{dtype,4},2}(undef, Ni, Nj)
    for j = 1:Nj, i = 1:Ni
        M[i,j] = atype(rand(dtype,2,2,2,2))
    end
    env = obs_env(M; χ = 10, verbose = true, savefile = true, infolder = "./data1", outfolder = "./data2", show_every = 3)
    @test env !== nothing
end


text:
using NCDatasets
using Distributed
using SharedArrays
using Formatting

s_map_file = "b.e11.B1850C5CN.f09_g16.005.pop.h.SST.100001-109912.nc"
d_map_file = "domain.ocn.gx3v7.120323.nc"

wgt_file = "wgt_gx1v6_to_gx3v7.nc"

const NNN = 9
missing_value = 1e20

Dataset(s_map_file, "r") do ds
        global s_mask, s_Nx, s_Ny, s_N, s_lat, s_lon

        s_mask = reshape(ds["REGION_MASK"][:], :)
        s_Nx = ds.dim["nlon"]
        s_Ny = ds.dim["nlat"]
        s_N  = s_Nx * s_Ny
        s_lon = reshape(ds["TLONG"][:] .|> deg2rad, :)
        s_lat = reshape(ds["TLAT"][:] .|> deg2rad, :)
end

Dataset(d_map_file, "r") do ds
        global d_mask, d_Nx, d_Ny, d_N, d_lat, d_lon

        d_mask = reshape(ds["mask"][:], :)
        d_Nx = ds.dim["ni"]
        d_Ny = ds.dim["nj"]
        d_N = d_Nx * d_Ny
        d_lon = reshape(ds["xc"][:] .|> deg2rad, :)
        d_lat = reshape(ds["yc"][:] .|> deg2rad, :)
end

#trans = spzeros(Float64, d_N, s_N)
trans = SharedArray{Float64}(NNN, d_N)

# s_coord and d_coord are the coordinates of grid points
# in 3-dimensional cartesian coordinate

s_coord = zeros(Float64, s_N, 3)
d_coord = zeros(Float64, d_N, 3)

s_coord[:, 1] .= cos.(s_lat) .* cos.(s_lon)
s_coord[:, 2] .= cos.(s_lat) .* sin.(s_lon)
s_coord[:, 3] .= sin.(s_lat)

d_coord[:, 1] .= cos.(d_lat) .* cos.(d_lon)
d_coord[:, 2] .= cos.(d_lat) .* sin.(d_lon)
d_coord[:, 3] .= sin.(d_lat)

s_NaN_idx = (s_mask .== 0)

println("Start making transform matrix... ")
#idx_arr = zeros(Integer, s_N)
#dist2   = zeros(Float64, s_N)

@time @sync @distributed for i = 1:d_N

    # For every point find its nearest-neighbors

    print("\r", i, "/", d_N)

    if d_mask[i] == 0
        trans[:, i] .= missing_value
        continue
    end

    dist2 = (  (s_coord[:, 1] .- d_coord[i, 1]).^2
             + (s_coord[:, 2] .- d_coord[i, 2]).^2
             + (s_coord[:, 3] .- d_coord[i, 3]).^2 )

    # Decided not to apply this condition because in
    # extreme cases there might be a small area of water
    # that is surrounded by lands.

    dist2[s_NaN_idx] .= NaN

    idx_arr = 1:s_N |> collect
    sort!(idx_arr; by=(k)->dist2[k])
    trans[:, i] = idx_arr[1:NNN]

end
println("done.")

s_lat = reshape(s_lat, s_Nx, :)
s_lon = reshape(s_lon, s_Nx, :)
d_lat = reshape(d_lat, d_Nx, :)
d_lon = reshape(d_lon, d_Nx, :)

Dataset(wgt_file, "c") do ds

    defDim(ds, "s_Nx", s_Nx)
    defDim(ds, "s_Ny", s_Ny)
    defDim(ds, "s_N",  s_N)

    defDim(ds, "d_Nx", d_Nx)
    defDim(ds, "d_Ny", d_Ny)
    defDim(ds, "d_N",  d_N)

    defDim(ds, "NNN", NNN)

    for (varname, vardata, dims) in (
        ("NN_idx", trans, ("NNN", "d_N")),
        ("s_lat", s_lat, ("s_Nx", "s_Ny")),
        ("s_lon", s_lon, ("s_Nx", "s_Ny")),
        ("d_lat", d_lat, ("d_Nx", "d_Ny")),
        ("d_lon", d_lon, ("d_Nx", "d_Ny")),
    )

        print(format("Output data: {} ...", varname))
        v = defVar(ds, varname, Float64, dims)
        v.attrib["_FillValue"] = missing_value
        v[:] = vardata
        println("done.")
    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Luna
Luna.set_fftw_mode(:estimate)
import FFTW
import NumericalIntegration: integrate, SimpsonEven
import Luna.PhysData: wlfreq

gas = :Ar
pres = 4

τ = 30e-15
λ0 = 800e-9

w0 = 2e-3
energy = 1.5e-9
L = 2

R = 6e-3
N = 128

grid = Grid.RealGrid(L, 800e-9, (400e-9, 2000e-9), 0.2e-12)
xygrid = Grid.FreeGrid(R, N)

x = xygrid.x
y = xygrid.y
energyfun, energyfunω = Fields.energyfuncs(grid, xygrid)

densityfun = let dens0=PhysData.density(gas, pres)
    z -> dens0
end

ionpot = PhysData.ionisation_potential(gas)
ionrate = Ionisation.ionrate_fun!_ADK(ionpot)
responses = (Nonlinear.Kerr_field(PhysData.γ3_gas(gas)),)
#  Nonlinear.PlasmaCumtrapz(grid.to, grid.to, ionrate, ionpot))

linop = LinearOps.make_const_linop(grid, xygrid, PhysData.ref_index_fun(gas, pres))
normfun = NonlinearRHS.const_norm_free(grid, xygrid, PhysData.ref_index_fun(gas, pres))

inputs = Fields.GaussGaussField(λ0=λ0, τfwhm=τ, energy=energy, w0=w0)

Eω, transform, FT = Luna.setup(grid, xygrid, densityfun, normfun, responses, inputs)

# statsfun = Stats.collect_stats(grid, Eω, Stats.ω0(grid))
output = Output.MemoryOutput(0, grid.zmax, 21)

Luna.run(Eω, grid, linop, transform, FT, output, max_dz=Inf, init_dz=1e-1)

ω = grid.ω
t = grid.t

zout = output.data["z"]
Eout = output.data["Eω"]

println("Transforming...")
Eωyx = FFTW.ifft(Eout, (2, 3))
Etyx = FFTW.irfft(Eout, length(grid.t), (1, 2, 3))
println("...done")

Ilog = log10.(Maths.normbymax(abs2.(Eωyx)))

Iωyx = abs2.(Eωyx);

Iyx = zeros(Float64, (length(y), length(x), length(zout)));
energy = zeros(length(zout));
for ii = 1:size(Etyx, 4)
    energy[ii] = energyfun(Etyx[:, :, :, ii]);
    Iyx[:, :, ii] = (grid.ω[2]-grid.ω[1]) .* sum(Iωyx[:, :, :, ii], dims=1);
end

ω0idx = argmin(abs.(grid.ω .- 2π*PhysData.c/λ0))

E0ωyx = FFTW.ifft(Eω[ω0idx, :, :], (1, 2));

Iωyx = abs2.(Eωyx)
Iωyxlog = log10.(Maths.normbymax(Iωyx));

import PyPlot:pygui, plt
pygui(true)
plt.figure()
plt.pcolormesh(x, y, (abs2.(E0ωyx)))
plt.colorbar()
plt.xlabel("X")
plt.ylabel("Y")
plt.title("I(ω=ω0, x, y, z=0)")

plt.figure()
plt.pcolormesh(x, y, (abs2.(Eωyx[ω0idx, :, :, end])))
plt.colorbar()
plt.xlabel("X")
plt.ylabel("Y")
plt.title("I(ω=ω0, x, y, z=L)")

plt.figure()
plt.pcolormesh(zout, ω.*1e-15/2π, Iωyxlog[:, N÷2+1, N÷2+1, :])
plt.xlabel("Z (m)")
plt.ylabel("f (PHz)")
plt.title("I(ω, x=0, y=0, z)")
plt.clim(-6, 0)
plt.colorbar()

plt.figure()
plt.pcolormesh(x.*1e3, y.*1e3, Iyx[:, :, 1])
plt.xlabel("X (mm)")
plt.ylabel("Y (mm)")
plt.title("\$\\int I(\\omega, x, y, z=0) d\\omega\$")
plt.colorbar()

plt.figure()
plt.pcolormesh(x.*1e3, y.*1e3, Iyx[:, :, end])
plt.xlabel("X (mm)")
plt.ylabel("Y (mm)")
plt.title("\$\\int I(\\omega, x, y, z=L) d\\omega\$")
plt.colorbar()

plt.figure()
plt.plot(zout.*1e2, energy.*1e6)
plt.xlabel("Distance [cm]")
plt.ylabel("Energy [μJ]")


text:
module MonadFunctions

export Just, None, just, NONE, fmap, cata, or_else
export Either, left, right, either, is_left, is_right, left_value, right_value, result
export list, flatten

include("maybe.jl")
include("either.jl")
include("cross.jl")
include("list.jl")

end # module