text:
function  NCLUSTERS = GetNegativeSpanningTrees(NCLUSTERS,Lx,Ly)
clusternames = fieldnames(NCLUSTERS);
numclusters = length(clusternames);
for zz = 1:numclusters

  positions = eval(['NCLUSTERS.cluster' num2str(zz) '.positions']);
  xi = positions(:,1);
  yi = positions(:,2);
  nc = length(xi);
  %Get the spanning tree
  xdiff= repmat(xi,1,nc) - repmat(xi',nc,1);
  ydiff= repmat(yi,1,nc) - repmat(yi',nc,1);
  rij = sqrt(xdiff.^2 + ydiff.^2);
  rij = tril(rij);
  rij = sparse(rij);

  [T,~] = graphminspantree(rij);

  %Try wrapping all positions, save the spanning tree which is
  %smallest
  T_old = T;
  [~,I] = sort(xi);
  xi_new = xi;
  yi_new = yi;

  %Wrapping in x...
  if ( any(xi<-Lx/4) && any(xi>Lx/4) )
    for ii =  1:length(xi)
      test_xi = xi;
      test_xi(I(1:ii)) = test_xi(I(1:ii)) + Lx;
      xdiff= repmat(test_xi,1,nc) - repmat(test_xi',nc,1);
      ydiff= repmat(yi,1,nc) - repmat(yi',nc,1);
      rij = sqrt(xdiff.^2 + ydiff.^2);
      rij = tril(rij);
      rij = sparse(rij);
      [T,~] = graphminspantree(rij);
      if sum(T(:)) < sum(T_old(:))
        disp('Found smaller spanning tree by wrapping in x')
        T_old = T;
        xi_new = test_xi;

      end
    end
  end

  %Wrapping in y...
  [~,I] = sort(yi);
  if (any(yi<-Ly/4) && any(yi>Ly/4)) %Changed to any() from sum()>0 (faster) ..
    for ii =  1:length(yi)
      test_yi = yi;
      test_yi(I(1:ii)) = test_yi(I(1:ii)) + Ly;
      ydiff= repmat(test_yi,1,nc) - repmat(test_yi',nc,1);
      xdiff= repmat(xi_new,1,nc) - repmat(xi_new',nc,1);
      rij = sqrt(xdiff.^2 + ydiff.^2);
      rij = tril(rij);
      rij = sparse(rij);
      [T,~] = graphminspantree(rij);
      if sum(T(:)) < sum(T_old(:))
        disp('Found smaller spanning tree by wrapping in y')
        T_old = T;
        yi_new = test_yi;
      end
    end
  end

  %Get centre-of-mass position and radius
  XC = mean(xi_new);
  YC = mean(yi_new);
  Xrad = mean(abs(xi_new-XC));
  Yrad = mean(abs(yi_new-YC));
  Radius = sqrt(Xrad.^2 + Yrad.^2); %#ok
  XC = XC -Lx*(XC > Lx/2);
  YC = YC -Ly*(YC > Ly/2);
  XC = XC +Lx*(XC < -Lx/2);
  YC = YC +Ly*(YC < -Ly/2);
  RC = [XC YC]; %#ok

  clear negindex
  [I,J] = find(T_old); %#ok
  eval(['NCLUSTERS.cluster' num2str(zz) '.spanningtreepointers = [I J];']);
  eval(['NCLUSTERS.cluster' num2str(zz) '.RC = RC;']);
  eval(['NCLUSTERS.cluster' num2str(zz) '.Radius = Radius;']);

 end


text:
println("There is more than one way to say hello:")

# This is an array consisting of three strings
helloArray = ["Hello","G'day","Shalom"]

for i in 1:3
    println("\t", helloArray[i], " World!")
end

println("\nThese squares are just perfect:")

# This construct is called a `comprehension' (or 'list comprehension')
squares = [i^2 for i in 0:10]

# You can loop on elements of arrays without having to use indexing
for s in squares
    print("  ",s)
end

# The last line of every code snippet is also evaluated as output (in addition to
# any figures and printing output generated previously).
sqrt.(squares)

ext: jl

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Test, Random, CUDA, Knet, AutoGrad

using CUDA.CUDNN:
    cudnnMultiHeadAttnForward,
    cudnnMultiHeadAttnForward!,
    cudnnMultiHeadAttnBackwardData,
    cudnnMultiHeadAttnBackwardWeights,
    cudnnGetMultiHeadAttnBuffers,
    cudnnGetMultiHeadAttnWeights,
    cudnnAttnDescriptor,
        cudnnAttnDescriptor_t,
        cudnnCreateAttnDescriptor,
        cudnnDestroyAttnDescriptor,
        cudnnSetAttnDescriptor,
        cudnnGetAttnDescriptor,
        cudnnDataType_t,
        cudnnDropoutDescriptor_t,
    cudnnAttnQueryMap_t,
        CUDNN_ATTN_QUERYMAP_ALL_TO_ONE, # 0         /* multiple Q-s map to a single (K,V) set when beam size > 1, beam sizes for (K,V) = 1 */
        CUDNN_ATTN_QUERYMAP_ONE_TO_ONE, # (1U << 0) /* multiple Q-s map to multiple (K,V) sets when beam size > 1, beam sizes for (K,V) = beam size for (Q) */
        CUDNN_ATTN_DISABLE_PROJ_BIASES, # 0         /* no biases in attention input and output projections */
        CUDNN_ATTN_ENABLE_PROJ_BIASES,  # (1U << 1) /* use biases in attention input and output projections */
    cudnnMultiHeadAttnWeightKind_t,
        CUDNN_MH_ATTN_Q_WEIGHTS, # 0, /* input projection weights for 'queries' */
        CUDNN_MH_ATTN_K_WEIGHTS, # 1, /* input projection weights for 'keys' */
        CUDNN_MH_ATTN_V_WEIGHTS, # 2, /* input projection weights for 'values' */
        CUDNN_MH_ATTN_O_WEIGHTS, # 3, /* output projection weights */
        CUDNN_MH_ATTN_Q_BIASES,  # 4, /* input projection bias tensor for 'queries' */
        CUDNN_MH_ATTN_K_BIASES,  # 5, /* input projection bias for 'keys' */
        CUDNN_MH_ATTN_V_BIASES,  # 6, /* input projection bias for 'values' */
        CUDNN_MH_ATTN_O_BIASES,  # 7, /* output projection biases */
    cudnnMathType_t,
        CUDNN_DEFAULT_MATH,                    # 0,
        CUDNN_TENSOR_OP_MATH,                  # 1,
        CUDNN_TENSOR_OP_MATH_ALLOW_CONVERSION, # 2,
        CUDNN_FMA_MATH,                        # 3,
    cudnnWgradMode_t,
        CUDNN_WGRAD_MODE_ADD,  # 0,
        CUDNN_WGRAD_MODE_SET,  # 1,
    cudnnSeqDataDescriptor,
        cudnnSeqDataDescriptor_t,
        cudnnCreateSeqDataDescriptor,
        cudnnDestroySeqDataDescriptor,
        cudnnSetSeqDataDescriptor,
        cudnnGetSeqDataDescriptor,
    cudnnSeqDataAxis_t,
        CUDNN_SEQDATA_TIME_DIM,  # 0, /* index in time */
        CUDNN_SEQDATA_BATCH_DIM, # 1, /* index in batch */
        CUDNN_SEQDATA_BEAM_DIM,  # 2, /* index in beam */
        CUDNN_SEQDATA_VECT_DIM,  # 3  /* index in vector */
        CUDNN_SEQDATA_DIM_COUNT, # 4
    cudnnDataType,
    cudnnSeqDataDefaultAxes,
    math_mode,
    sdim,
    handle

@testset "cudnn/multiheadattn" begin

    Q,K,V,B,T,F = 6,6,5,4,3,Float64
    w, q, k, v = (Param(CUDA.randn(x...)) for x in ((F,400),(F,Q,B,T),(F,K,B,T),(F,V,B,T)))

    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; )

    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; seqLengthsQO = Cint[1,2,3,1])
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; seqLengthsKV = Cint[1,2,3,1])
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; nHeads = 2)
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; smScaler = 2)
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; mathType = CUDNN_DEFAULT_MATH)
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; mathType = CUDNN_TENSOR_OP_MATH)
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; mathType = CUDNN_TENSOR_OP_MATH_ALLOW_CONVERSION)
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; mathType = CUDNN_FMA_MATH)
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; kProjSize = 7, qProjSize = 7) # k and q have to match
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; vProjSize = 7)
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; oProjSize = 7)
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; qoMaxSeqLength = 7)
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; kvMaxSeqLength = 7)
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; maxBatchSize = 7)
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; maxBeamSize = 7)
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; loWinIdx = fill(Cint(1),T))
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; hiWinIdx = fill(Cint(1),T))
    # currIdx >= 0 can only be used in inference mode
    # @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; currIdx = 0)

    # Test residuals: residuals and output (and thus v unless oProjSize>0) must match q in vector size
    v, residuals = (Param(CUDA.randn(x...)) for x in ((F,Q,B,T),(F,Q,B,T)))
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; residuals = residuals)

    # PROJ_BIASES not supposed to be supported according to docs, it works for some cases
    # (attnDesc.attnMode & CUDNN_ATTN_ENABLE_PROJ_BIASES == 0) || (@warn "The CUDNN_ATTN_ENABLE_PROJ_BIASES option is not supported in the multi-head attention gradient functions." maxlog=1)
    @test_skip @gcheck cudnnMultiHeadAttnForward(w, q, k, v; attnMode = CUDNN_ATTN_QUERYMAP_ALL_TO_ONE | CUDNN_ATTN_ENABLE_PROJ_BIASES |> Cuint, qProjSize=9, kProjSize=9, vProjSize=6, oProjSize=4)

    # Test nonstandard axes order
    w, q, k, v = (Param(CUDA.randn(x...)) for x in ((F,100),(F,Q,T,B),(F,K,T,B),(F,V,T,B)))
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; axes = [CUDNN_SEQDATA_VECT_DIM, CUDNN_SEQDATA_TIME_DIM, CUDNN_SEQDATA_BATCH_DIM, CUDNN_SEQDATA_BEAM_DIM])

    # Test beam handling
    w, q, k, v = (Param(CUDA.randn(x...)) for x in ((F,100),(F,Q,B,T,2),(F,K,B,T,1),(F,V,B,T,1)))
    @test @gcheck cudnnMultiHeadAttnForward(w, q, k, v; )
    # CUDNN_ATTN_QUERYMAP_ONE_TO_ONE does not seem to be supported
    w, q, k, v = (Param(CUDA.randn(x...)) for x in ((F,100),(F,Q,B,T,2),(F,K,B,T,2),(F,V,B,T,2)))
    @test_skip @gcheck cudnnMultiHeadAttnForward(w, q, k, v; attnMode = CUDNN_ATTN_QUERYMAP_ONE_TO_ONE | CUDNN_ATTN_DISABLE_PROJ_BIASES |> Cuint) ## Not supported

end


text:
using Test

if !@isdefined RandomNumbers
    include("../common.jl")
end

@testset "Xorshifts" begin

    using RandomNumbers.Xorshifts

    @info "Testing Xorshifts"
    stdout_ = stdout
    pwd_ = pwd()
    cd(dirname(@__FILE__))
    rm("./actual"; force=true, recursive=true)
    mkpath("./actual")

    seeds = [0x1ed73cdc948901ab, 0x6d76dcf952161e22, 0x8d0eba5421bc695b, 0xa16065e26c36cce2,
            0x944867a9aacc1bfa, 0x173d478b3ac9de91, 0x97becbc60ac9c1f2, 0x1b1cfeff0505b998,
            0x815a28dc8372d7ab, 0x5c6746d45ae2c5c4, 0x17f22f9aa839c717, 0x1fd7e8194bb9c598,
            0x442ffe4a57ed7987, 0xfd686ec417788eb6, 0x1ec96940807a8749, 0x64d93a8608a988ef]

    for (rng_name, seed_t) in (
            (:SplitMix64,           UInt64),
            (:Xoroshiro64Star,      NTuple{2, UInt32}),
            (:Xoroshiro64StarStar,  NTuple{2, UInt32}),
            (:Xoroshiro128Plus,     NTuple{2, UInt64}),
            (:Xoroshiro128StarStar, NTuple{2, UInt64}),
            (:Xoshiro128Plus,       NTuple{4, UInt32}),
            (:Xoshiro128StarStar,   NTuple{4, UInt32}),
            (:Xoshiro256Plus,       NTuple{4, UInt64}),
            (:Xoshiro256StarStar,   NTuple{4, UInt64}),

            (:Xorshift64,       UInt64),
            (:Xorshift64Star,   UInt64),
            (:Xorshift128,      NTuple{2, UInt64}),
            (:Xorshift128Star,  NTuple{2, UInt64}),
            (:Xorshift128Plus,  NTuple{2, UInt64}),
            (:Xorshift1024,     NTuple{16, UInt64}),
            (:Xorshift1024Star, NTuple{16, UInt64}),
            (:Xorshift1024Plus, NTuple{16, UInt64}),
            (:Xoroshiro128,     NTuple{2, UInt64}),
            (:Xoroshiro128Star, NTuple{2, UInt64}),
        )

        outfile = open(string(
            "./actual/check-$(lowercase("$rng_name")).out"
        ), "w")
        redirect_stdout(outfile)

        rng_type = @eval Xorshifts.$rng_name

        x = rng_type()
        x = rng_type(1234)
        if rng_type == Xorshifts.SplitMix64
            rng_type(0)
        else
            @test_throws(
                ErrorException("0 cannot be the seed"),
                seed!(x, seed_t <: Tuple ? (zeros(seed_t.types[1], length(seed_t.types))...,) : 0)
            )
        end
        seed!(x)
        @test seed_type(x) == seed_t
        @test copyto!(copy(x), x) == x

        seed = seed_t <: Tuple ? ((seeds[i] % seed_t.types[i] for i in 1:length(seed_t.types))...,) : seeds[1] % seed_t
        x = rng_type(seed)

        for i in 1:100
            @printf "%.9f\n" rand(x)
        end

        close(outfile)
    end
    redirect_stdout(stdout_)

    compare_dirs("expected", "actual")

    # Package content should not be modified.
    # rm("./actual"; force=true, recursive=true)

    cd(pwd_)

end


text:
# Test RepresentativeStateTable and DiagonalizationSector with reflection
# symmetry at k=0 and k=π
let L = 10
    hs = SpinHalfHilbertSpace(ChainLattice([L]))
    seed_state!(hs, N_up=0)
    apply_hamiltonian = spin_half_hamiltonian(J1_z=1, h_x=0.45, h_z=0.87)
    rst = RepresentativeStateTable(hs, apply_hamiltonian, [reflection_symmetry])

    full_ham = operator_matrix(hs, apply_hamiltonian)
    @test length(hs.indexer) == 2^L

    k_values = [0]
    rem(L, 2) == 0 && push!(k_values, div(L, 2)) # momentum π, if available
    for k_idx in k_values .+ 1
        for reflection_idx in 0:1
            diagsect = DiagonalizationSector(rst, 1, k_idx, [reflection_idx])
            reduced_ham = Matrix(construct_reduced_hamiltonian(diagsect))
            fact = eigen(Hermitian((reduced_ham + reduced_ham') / 2))
            evals, evecs = fact.values, fact.vectors
            for i in 1:length(evals)
                eval = evals[i]
                evec = evecs[:,i]
                ψ = get_full_psi(diagsect, evec)
                check_eigenstate(full_ham, eval, ψ)
            end
        end
    end
end


text:
reset_default_graph(); sess = Session()
function test_flow(sess, flows, dim)
    prior = ADCME.MultivariateNormalDiag(loc=zeros(dim))
    model = NormalizingFlowModel(prior, flows)
    x = rand(3, dim)
    z, _, J = model(x)
    xs, iJ = ADCME.backward(model, z[end])
    init(sess)
    x_, J_, iJ_ = run(sess, [xs[end], J, iJ])
    x, x_, J_, iJ_
end

@testset "LinearFlow" begin
    dim = 3
    flows = [LinearFlow(dim)]
    x1, x2, j1, j2 = test_flow(sess, flows, dim)
    @test x1≈x2
    @test j1≈-j2
end

@testset "AffineConstantFlow" begin
    dim = 3
    flows = [AffineConstantFlow(dim)]
    x1, x2, j1, j2 = test_flow(sess, flows, dim)
    @test x1≈x2
    @test j1≈-j2
end

@testset "ActNorm" begin
    dim = 3
    flows = [ActNorm(dim)]
    x1, x2, j1, j2 = test_flow(sess, flows, dim)
    @test x1≈x2
    @test j1≈-j2
end

@testset "SlowMAF" begin
    dim = 3
    flows = [SlowMAF(dim, false, [x->fc(x, [20,20,2], "fc1"); x->fc(x,[20,20,2],"fc2")])]
    x1, x2, j1, j2 = test_flow(sess, flows, dim)
    @test x1≈x2
    @test j1≈-j2
end

@testset "MAF" begin
    dim = 3
    flows = [MAF(dim, false, [20,20,20])]
    x1, x2, j1, j2 = test_flow(sess, flows, dim)
    @test x1≈x2
    @test j1≈-j2
end

@testset "IAF" begin
    dim = 3
    flows = [IAF(dim, false, [20,20,20])]
    x1, x2, j1, j2 = test_flow(sess, flows, dim)
    @test x1≈x2
    @test j1≈-j2
end

@testset "Invertible1x1Conv" begin
    dim = 3
    flows = [Invertible1x1Conv(dim)]
    x1, x2, j1, j2 = test_flow(sess, flows, dim)
    @test x1≈x2
    @test j1≈-j2
end

@testset "AffineHalfFlow" begin
    dim = 3
    flows = [AffineHalfFlow(dim, true)]
    x1, x2, j1, j2 = test_flow(sess, flows, dim)
    @test x1≈x2
    @test j1≈-j2
end

@testset "NeuralCouplingFlow" begin
    @test_skip begin
        dim = 6
        K = 8
        n = dim÷2
        flows = [NeuralCouplingFlow(dim, x->fc(x, [20,20,n*(3K-1)], "ncf1"), x->fc(x,[20,20,n*(3K-1)],"ncf2"), K)]
        x1, x2, j1, j2 = test_flow(sess, flows, dim)
        @test x1≈x2
        @test j1≈-j2
    end
end

@testset "Permute" begin
    dim = 6
    flows = [Permute(dim, randperm(6) .- 1)]
    x1, x2, j1, j2 = test_flow(sess, flows, dim)
    @test x1≈x2
    @test j1≈-j2
end

@testset "composite" begin
    function create_linear_transform(dim)
        permutation = randperm(dim) .- 1
        flow = [Permute(dim, permutation);LinearFlow(dim)]
    end
    dim = 5
    flows = create_linear_transform(dim)
    x1, x2, j1, j2 = test_flow(sess, flows, dim)
    @test x1≈x2
    @test j1≈-j2
end

# @testset "composite 2" begin
#     function create_base_transform(dim, K=8)
#         n1 = dim÷2
#         n2 = dim - n1
#         r1 = x->Resnet1D(n1 * (3K-1), 256, dropout_probability=0.0, use_batch_norm=false)(x)
#         r2 = x->Resnet1D(n2 * (3K-1), 256, dropout_probability=0.0, use_batch_norm=false)(x)
#         NeuralCouplingFlow(dim, r1, r2)
#     end
#     dim = 6
#     flows = [create_base_transform(dim)]
#     x1, x2, j1, j2 = test_flow(sess, flows, dim)
#     @test x1≈x2
#     @test j1≈-j2
# end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    KBallSearcher(object, k, ball)

A method that searches `k` nearest neighbors and then filters
these neighbors using a norm `ball`.
"""
struct KBallSearcher{O,B,K} <: AbstractBoundedNeighborSearcher
  # input fields
  object::O
  k::Int
  ball::B

  # state fields
  kdtree::K
end

function KBallSearcher(object::O, k::Int, ball::B) where {O,B}
  kdtree = KDTree(coordinates(object), metric(ball))
  KBallSearcher{O,B,typeof(kdtree)}(object, k, ball, kdtree)
end

maxneighbors(searcher::KBallSearcher) = searcher.k

function search!(neighbors, xₒ::AbstractVector,
                 searcher::KBallSearcher; mask=nothing)
  k = searcher.k
  r = radius(searcher.ball)

  inds, dists = knn(searcher.kdtree, xₒ, k, true)

  # keep neighbors inside ball
  keep = dists .≤ r

  # possibly mask some of the neighbors
  isnothing(mask) || (keep .*= mask[inds])

  nneigh = 0
  @inbounds for i in 1:k
    if keep[i]
      nneigh += 1
      neighbors[nneigh] = inds[i]
    end
  end

  nneigh
end


text:
abstract type AbstractModels end

struct Model1 <: AbstractModels end
struct Model2 <: AbstractModels end

m1 = Model1()
m2 = Model2()

data_union = Union{Missing, AbstractString, Dict, Array{T, 1} where T, NamedTuple}
init_union = Union{Missing, AbstractString, Dict, Array{T, 1} where T, NamedTuple}

function method2(model, data::Union{NamedTuple, Dict})

  println("Method 2 called.\n")
end

function method1(model::AbstractModels; kwargs...)
    println("Method 1 called.")
    :init in keys(kwargs) && println("$(kwargs[:init])")
    :init in keys(kwargs) && println("$(kwargs[:init].v)")
    :data in keys(kwargs) && println("$(kwargs[:data])")
    :data in keys(kwargs) && println("$(kwargs.data)")
    println()
end

method1(m1)
method1(m2, data=(z=12,))
method1(m1; data = (x=2, z=13))
method1(m1; init = (v=2,))
method1(m1; init = (v=2,), data = (y=3,))

method2(m2, (x=2,))
method2(m1, Dict(:x=>2))


text:
# A collection of simple programs that more fully test aspects of the processor.

# This function will read numbers from input fifo 0. It will store a single
# number at dmem 0
#
# If a value it reads is less than the value stored in DMEM, it will forward
# the newly read value to output 0
#
# Otherwise, it will swap the two numbers and write the old number stored in
# DMEM to output 0
@asap4asm function dumb_sort()
    @label start
        # Read from FIFO 0 and store the result into DMEM 1
        MOVE(dmem[1], ibuf[0])
        # Subtract DMEM 1 from DMEM 0. This will set the "negative" alu flag if
        # dmem[1] > dmem[0]. Branch based on the outcome of this comparison.
        SUBU(null, dmem[0], bypass[1])
        BRL(swap, N, j)
        # Write out DMEM[1]
        NOP()
        NOP()
        NOP()
        MOVE(output[0], dmem[1])
        BRL(start)

    @label swap
        # Take advantage of latency in the pipeline to swap the two entries in
        # DMEM
        MOVE(dmem[1], dmem[0])
        MOVE(dmem[0], dmem[1], nop3)
        # Add NOPs to ensure writeback happens.
        NOP()
        BRL(back)
end


text:
using DataFrames
using CSV
using DataStructures: SortedDict, Reverse
using StatsBase
using CategoricalArrays
using HypothesisTests
using Plots
using StatsPlots
using FreqTables

DIR_OUT = "/Users/olivierlabayle/Documents/JoubertSyndrome"

DATAPATH = joinpath(DIR_OUT, "AA_base_joubert_statistiques_OL_11621.csv")

PHEN_OF_INTEREST = ("Neonatal breathing dysregulations",
                    "Epilepsy",
                    "Obesity or overweight",
                    "Macrocephaly",
                    "Retinal dystrophy",
                    "Coloboma",
                    "Kidney defect",
                    "Liver defect",
                    "Orofacial anomalies",
                    "Multiple pulmonary infections",
                    "Endocrinal and/or genital features ",
                    "polydactyly",
                    "pure JS")

GENES_OF_INTEREST = ("CC2D2A",
                    "CPLANE1",
                    "AHI1",
                    "CEP290",
                    "TMEM67",
                    "KIAA0586",
                    "INPP5E",
                    "CSPP1",
                    "RPGRIP1L",
                    "OFD1",
                    "MKS1",
                    "NPHP1")

data = DataFrame(CSV.File(DATAPATH; missingstrings=["na"]))

# ###########################################################################
#         PVALUES                                                           #
# ###########################################################################

function getfigures(x, y)
    ct = freqtable(x, y)
    odds_phen_if_gen = (ct[2,2]/ct[2,1]) / (ct[1,2]/ct[1,1])
    se = sqrt(sum(1 ./ ct))
    lower_bound = exp(log(odds_phen_if_gen) - 1.96*se)
    upper_bound = exp(log(odds_phen_if_gen) + 1.96*se)
    min_freq = minimum(ct)
    reschi = ChisqTest(ct)
    resfish = FisherExactTest(ct[1,1], ct[1,2], ct[2,1], ct[2,2])
    return min_freq, reschi, resfish, odds_phen_if_gen, lower_bound, upper_bound
end

function gene_phenotype_it(data::DataFrame, gene, phenotype)
    data_restrict = dropmissing(data[:, ["GENE", phenotype]])
    x = categorical(data_restrict.GENE .== gene)
    y = categorical(data_restrict[:, phenotype] .== "oui" )
    return getfigures(x, y)
end

function phenotype_phenotype_it(data::DataFrame, phenotype1, phenotype2)
    data_restrict = dropmissing(data[:, [phenotype1, phenotype2]])
    x = categorical(data_restrict[:, phenotype1] .== "oui")
    y = categorical(data_restrict[:, phenotype2] .== "oui")
    return getfigures(x, y)
end

function gene_phenotype_it_results(data::DataFrame, genes, phenotypes)
    results = DataFrame(Gene = String[], Phenotype = String[],
                PVAL_CHI2 = Float64[], PVAL_FISH = Float64[], MIN_FREQ = Int[],
                OR = Float64[], OR_LB = Float64[], OR_UB = Float64[],
                RES_CHI2 = Any[], RES_FISH = Any[])
    for gene in genes
        for phenotype in phenotypes
            min_freq, reschi, resfish, odds_phen_if_gen, lower_bound, upper_bound = gene_phenotype_it(data, gene, phenotype)
            push!(results,
            (gene, phenotype, pvalue(reschi), pvalue(resfish), min_freq, odds_phen_if_gen, lower_bound, upper_bound, reschi, resfish))
        end
    end
    return results
end

function phenotype_phenotype_it_results(data::DataFrame, phenotypes)
    results = DataFrame(Phenotype1 = String[], Phenotype2 = String[],
                        PVAL_CHI2 = Float64[], PVAL_FISH = Float64[],
                        MIN_FREQ = Int[], OR = Float64[], OR_LB = Float64[], OR_UB = Float64[],
                        RES_CHI2 = Any[], RES_FISH = Any[])
    for phenotype1 in phenotypes
        for phenotype2 in phenotypes
            if phenotype1 != phenotype2
                min_freq, reschi, resfish, odds_phen_if_gen, lower_bound, upper_bound = phenotype_phenotype_it(data, phenotype1, phenotype2)
                push!(results,
                (phenotype1, phenotype2, pvalue(reschi), pvalue(resfish), min_freq, odds_phen_if_gen, lower_bound, upper_bound, reschi, resfish))
            end
        end
    end
    return results
end

gene_phen_results = gene_phenotype_it_results(data, GENES_OF_INTEREST, PHEN_OF_INTEREST)

phen_phen_results = phenotype_phenotype_it_results(data, PHEN_OF_INTEREST)

CSV.write(joinpath(DIR_OUT, "gene_phen.csv"), gene_phen_results[:, Not(["RES_CHI2", "RES_FISH"])])
CSV.write(joinpath(DIR_OUT, "phen_phen.csv"), phen_phen_results[:, Not(["RES_CHI2", "RES_FISH"])])

# ###########################################################################
#         PLOTS                                                             #
# ###########################################################################

genes_to_plot = ["CC2D2A", "CPLANE1", "AHI1", "CEP290", "TMEM67", "OFD1", "NPHP1", "CSPP1"]

phens_to_plot = ["Neonatal breathing dysregulations", "Epilepsy", "Retinal dystrophy", "Coloboma",
"Kidney defect", "Liver defect", "Orofacial anomalies", "polydactyly", "Multiple pulmonary infections",
"Endocrinal and/or genital features ", "pure JS"]

cols_to_xticks_fr = Dict(
    "Multiple pulmonary infections" => "Infections respiratoires \nrépétées",
    "polydactyly" => "Polydactylie",
    "Endocrinal and/or genital features " => "Anomalies \nendocrino-génitales",
    "Orofacial anomalies" => "Anomalies orofaciales",
    "Liver defect" => "Atteinte hépatique",
    "Kidney defect" =>  "Atteinte rénale",
    "Retinal dystrophy" => "Atteinte rétinienne",
    "Coloboma" => "Colobome",
    "Epilepsy" => "Epilepsie",
    "Neonatal breathing dysregulations" => "Troubles respiratoires \nnéonataux",
    "pure JS" => """JS \"pur\""""
 )

 cols_to_xticks_en = Dict(
    "Multiple pulmonary infections" => "Multiple pulmonary \ninfections",
    "polydactyly" => "Polydactyly",
    "Endocrinal and/or genital features " => "Endocrinal and/or \ngenital features",
    "Orofacial anomalies" => "Orofacial anomalies",
    "Liver defect" => "Liver disease",
    "Kidney defect" =>  "Kidney disease",
    "Retinal dystrophy" => "Retinal disease",
    "Coloboma" => "Coloboma",
    "Epilepsy" => "Epilepsy",
    "Neonatal breathing dysregulations" => "Neonatal breathing \ndysregulations",
    "pure JS" => """\"Pure\" JS"""
    )

function myfreq(x)
    nomiss_x = skipmissing(x)
    nbyes = sum(nomiss_x .== "oui")
    nb = length(collect(nomiss_x))
    return round(100*nbyes/nb, digits=2)
end

function getstars(gene_phen_results::DataFrame, gene, phen)
    pvals = filter(row->(row.Gene == gene && row.Phenotype == phen),
            gene_phen_results)[1, [:PVAL_FISH, :PVAL_CHI2]]
    pvals = Vector(pvals)
    if all(pvals .< 0.001)
        return " (***)"
    elseif all(pvals .< 0.01)
        return " (**)"
    elseif all(pvals .< 0.05)
        return " (*)"
    else
        return ""
    end
end

function Base.unique(ctg::CategoricalArray)
    l = levels(ctg)
    newctg = CategoricalArray(l)
    levels!(newctg, l)
end

function save_plots(data, genes_to_plot, phens_to_plot, gene_phen_results, cols_to_xticks;
                     dirout=DIR_OUT,
                     lang="en",
                     titlebase="Phenotype affection for gene ",
                     xlabel="Phenotypes",
                     ylabel="Percentage of affected individuals",
                     groups=["Mutated", "Not mutated"])
    export_df = DataFrame()
    for gene in genes_to_plot
        data[!, gene] = data[!, "GENE"] .== gene
        grouped = groupby(data, gene)
        freqs = combine(grouped, phens_to_plot .=> myfreq)
        mutation_true_freqs = Vector(filter(row->row[gene] === true, freqs)[1, Not(gene)])
        mutation_false_freqs = Vector(filter(row->row[gene] === false, freqs)[1, Not(gene)])
        xticks = [cols_to_xticks[phen] * getstars(gene_phen_results, gene, phen) for phen in phens_to_plot]
        xticks = categorical(repeat(xticks, outer=2), levels=xticks)
        p = groupedbar(xticks, [mutation_true_freqs mutation_false_freqs],
                    bar_position=:dodge,
                    bar_width=0.7,
                    # xlabel=xlabel,
                    xtickfontsize=7,
                    title=titlebase * gene,
                    rotation=45,
                    fontfamily="Times",
                    ylims=(0, 100),
                    legend=:best,
                    ylabel=ylabel,
                    group=repeat(groups, inner=length(phens_to_plot)))
        savefig(p, joinpath(dirout, "$(gene)_barplot_$lang.png"))
        export_df[!, "$(gene)_mutation_true_freqs"] = mutation_true_freqs
        export_df[!, "$(gene)_mutation_false_freqs"] = mutation_false_freqs
    end
    return export_df
end

data_to_plot = data[:, vcat(phens_to_plot, ["GENE"])]

save_plots(data_to_plot, genes_to_plot, phens_to_plot, gene_phen_results, cols_to_xticks_en;
                     dirout=DIR_OUT,
                     lang="en",
                     titlebase="Phenotype affection for gene ",
                     xlabel="Phenotypes",
                     ylabel="Percentage of affected individuals",
                     groups=["Mutated", "Not mutated"])

export_df = save_plots(data_to_plot, genes_to_plot, phens_to_plot, gene_phen_results, cols_to_xticks_fr;
                     dirout=DIR_OUT,
                     lang="fr",
                     titlebase="Phénotype du gène ",
                     xlabel="Phénotypes",
                     ylabel="Pourcentage d'individus affectés",
                     groups=["Muté", "Non muté"])


text:
"""
subcommands for managing large Julia package
that contains sub-packages in its `lib` directory.
"""
@cast module Package

using Pkg
using TOML
using Comonicon
using ..Internal: collect_lib_deps, foreach_subpackage, develop_local_deps

"""
    initall(;root_path::String=root_dir(), no_docs::Bool=false, no_examples::Bool=false)

Init all pacakge dependencies by developing
local dependencies in `lib`.

# Intro

This command will init the local dependencies
of folders such as `lib`, `examples` and `docs`.

# Options

- `--root-path=<path/to/package>`: path to the root of the package,
    default will try to find the closest git repository.

# Flags

- `--no-docs`: exclude `docs` folder.
- `--no-examples`: exclude `examples` folder.
"""
@cast function initall(;root_path::String=root_dir(), no_docs::Bool=false, no_examples::Bool=false)
    proj_path = Pkg.project().path

    foreach_subpackage(root_path) do package, path
        develop_local_deps(root_path, path)
    end

    # setup docs env
    docs_dir = joinpath(root_path, "docs")
    isdir(docs_dir) && develop_local_deps(root_path, docs_dir)

    # setup examples env
    examples_dir = joinpath(root_path, "examples")
    isdir(examples_dir) && develop_local_deps(root_path, examples_dir)

    # dev libs for root package
    develop_local_deps(root_path, root_path)

    # activate current project
    Pkg.activate(proj_path; io=devnull)
    return
end

"""
    init(path::String; root_path::String=root_dir())

Init a project environment at `path` by developing local dependencies from `root_path`.

# Args

- `path`: path to the project to init.

# Options

- `--root-path=<path/to/package>`: path to the root of the package,
    default will try to find the closest git repository.

"""
@cast function init(path::String; root_path::String=root_dir())
    proj_path = Pkg.project().path
    ispath(path) || cmd_error("$path does not exists")
    develop_local_deps(root_path, path)
    # activate current project
    Pkg.activate(proj_path; io=devnull)
    return
end

end


text:
module Pheno

using Cropbox
using TimeZones
import Dates

@system Estimator begin
    year ~ preserve::int(parameter)

    # 270(+1)th days of the first year (around end of September)
    Ds: start_date_offset => 270 ~ preserve::int(u"d")

    # 150 days after the second year (around end of May)
    De: end_date_offset => 150 ~ preserve::int(u"d")

    tz: timezone => tz"UTC" ~ preserve::TimeZone(parameter)

    t0(year, tz, Ds): start_date => begin
        ZonedDateTime(year-1, 1, 1, tz) + Dates.Day(Ds)
    end ~ preserve::datetime

    t1(year, tz, De): end_date => begin
        ZonedDateTime(year, 1, 1, tz) + Dates.Day(De)
    end ~ preserve::datetime

    calendar(context, init=t0') ~ ::Calendar
    t(calendar.time): current_date ~ track::datetime

    s: store ~ provide(init=t, parameter)

    match => false ~ flag
    stop(m=match, t, t1) => (m || t >= t1) ~ flag

    T: temperature ~ drive(from=s, by=:tavg, u"°C")
end

estimate(S::Type{<:Estimator}, years; config, index=[:year, "calendar.time"], target=[:match], stop=:stop, kwargs...) = begin
    configs = @config config + !(S => :year => years)
    simulate(S; index, target, configs, stop, snap=:match, kwargs...)
end

@system BetaFuncEstimator(BetaFunction, Estimator, Controller) <: Estimator begin
    Rg: growth_requirement ~ preserve(parameter)
    Cg(ΔT): growth_cumulated ~ accumulate
    match(Cg, Rg) => Cg >= Rg ~ flag
end

end

using DataFrames
using TimeZones
import Dates
@testset "pheno" begin
    t0 = ZonedDateTime(2016, 9, 1, tz"UTC")
    t1 = ZonedDateTime(2018, 9, 30, tz"UTC")
    dt = Dates.Hour(1)
    T = collect(t0:dt:t1);
    df = DataFrame(index=T, tavg=25.0);

    Rg = 1000
    config = (
        :Estimator => (
            :tz => tz"UTC",
            :store => df,
        ),
        :BetaFuncEstimator => (
            :To => 20,
            :Tx => 35,
            :Rg => Rg,
        )
    )

    @testset "single" begin
        r = Pheno.estimate(Pheno.BetaFuncEstimator, 2017; config, target=[:ΔT, :Cg, :match, :stop])
        @test nrow(r) == 1
        r1 = r[1, :]
        @test r1.match == true
        @test r1.stop == true
        @test r1.Cg >= Rg
    end

    @testset "double" begin
        r = Pheno.estimate(Pheno.BetaFuncEstimator, [2017, 2018]; config, target=[:ΔT, :Cg, :match, :stop])
        @test nrow(r) == 2
        @test all(r.match)
        @test all(r.stop)
        @test all(r.Cg .>= Rg)
    end
end


text:
"""A package for graph computations related to graph Laplacians

Graphs are represented by sparse adjacency matrices, etc.
"""
module Laplacians

  function __init__()

    #=
    if !isdefined(Main, :LAPLACIANS_NOPLOT) & !haskey(ENV,"LAPLACIANS_NOPLOT")
        eval(Expr(:using, :PyPlot))
    end
    =#

    if isdefined(Main, :LAPLACIANS_AMG)
        eval(Expr(:using, :PyAMG))

    end
  end

  using RandomV06
  V06 = RandomV06.V06
  V07 = RandomV06.V07
  Vcur = RandomV06.Vcur

  using Arpack
  using SuiteSparse
  using DataStructures
  using SparseArrays
  using Random
  using LinearAlgebra
  using Statistics
  using Printf

  using DelimitedFiles

  using Plots

  include("IJV.jl")

  include("fastCSC.jl")
  export symPermuteCSC
  export symTransposeCSC
  export submatrixCSC

  include("graphUtils.jl")
  export deg
  export nbri
  export weighti
  export nbrs
  export wdeg
  export setValue
  export backIndices
  export flipIndex
  export findEntries
  export compConductance
  export getVolume
  export getObound

  include("graphGenerators.jl")
  export ring_graph
  export generalized_ring
  export rand_matching
  export rand_regular
  export grown_graph
  export grown_graph_d
  export pref_attach
  export hypercube
  export complete_binary_tree
  export complete_graph
  export empty_graph
  export path_graph

  #export wgrid2

  export grid2
  export grid2coords
  export grid3

  export rand_gen_ring
  export randperm
  export ErdosRenyi
  export ErdosRenyiCluster
  export ErdosRenyiClusterFix
  export pure_random_graph

  export chimera
  export rand_weight
  export wted_chimera, semiwted_chimera

  include("latinSquares.jl")
  export latin_square, latin_square_graph

  include("IO.jl")
  export writeIJV, read_graph

  include("graphOps.jl")

  export unweight, unweight!
  export mapweight
  export uniformWeight, uniformWeight!

  export edgeVertexMat, wtedEdgeVertexMat

  export power, thicken_once, thicken

  export product_graph
  export generalized_necklace
  export subsampleEdges

  export two_lift
  export join_graphs, join_graphs!, disjoin

  export plot_graph

  export shortIntGraph, floatGraph

  export lap
  export adj
  export spectral_coords
  export spectral_drawing

  export diagmat

  include("graphAlgs.jl")

  export components
  export biggestComp
  export vecToComps
  export isConnected

  export shortestPaths, shortestPathTree, pathFromParents
  export kruskal, prim

  include("treeAlgs.jl")

  export RootedTree
  export matToTree
  export matToTreeDepth
  export tarjanStretch
  export compDepth
  export comp_stretches
  export dfsOrder

  include("pcg.jl")

  export cg, cgSolver, cgLapSolver
  export pcg, pcgSolver, pcgLapSolver

  include("flow.jl")

  export maxflow

  include("akpw.jl")

  export akpw, akpwU

  include("localClustering.jl")

  export prn
  export apr
  export localImprove

  include("cutHeuristics.jl")

  export refineCut
  export dumbRefineCut

  include("randTrees.jl")
  export randishKruskal, randishPrim

  include("sampler.jl")
  include("fastSampler.jl")
  export FastSampler, sample, sampleMany
  export blockSample

  include("solverInterface.jl")
  export chol_sddm, chol_lap, lapWrapSDDM

  include("augTreeSolver.jl")

  export augmentTree, augTreePrecon, augTreeSddm
  export augTreeLapPrecon, augTreeLap, AugTreeParams, AugTreeParamsOld

  include("KMPSolver.jl")
  export KMPSDDMSolver
  export KMPLapSolver
  export KMPParams

  include("approxCholTypes.jl")
  include("approxChol.jl")
  export approxchol_lap, ApproxCholParams, approxchol_sddm

  include("fiedler.jl")
  export fiedler

  include("complexSolvers.jl")
  export SDDMSolvers
  export LapSolvers

  include("compare_solvers.jl")
  export SolverTest, speedTestLapSolvers

  include("conditionNumber.jl")
  export support, approxQual, conditionNumber

  include("sparsify.jl")
  export sparsify

  include("johnlind.jl")
  export johnlind

  include("toposort.jl")
  export toposort, dirEdgeVertexMat

  # include("isotonicIPM.jl")
  # export isotonicIPM, isotonicIPMrelEps

  include("harmonic.jl")
  export harmonic_interp

  include("from_cholmod.jl")
  export cholmod_perm, ask_cholmod

  include("deprecated.jl")

end # module Laplacians.jl


text:
import Term: inspect

println("\nTesting logging, stdout temporarily disabled")
@suppress_out begin
    @testset "\e[34mINSPECT test" begin

        # define some types
        abstract type T1 end

        abstract type T2 <: T1 end

        """§§
            MyType

        It's just a useless type we've defined to provide an example of
        Term's awesome `inspect` functionality!
        """
        struct MyType <: T2
            x::Int
            y::String
            fn::Any
        end

        #  Now, what is MyType like?
        @test_nowarn inspect(MyType)

        # Let's define some constructors and methods using MyType

        """
        constructors!
        """
        MyType(x::Int) = MyType(x, "no string", nothing)

        MyType(x::Int, y) = MyType(x, y, nothing)

        # methods
        useless_method(m::MyType) = m

        """
        these methods don't do much, just an example
        """
        another_method(m1::MyType, m2::MyType) = print(m1, m2)

        # let's inspect MyType again!
        @test_nowarn inspect(MyType)

        """
        Not bad huh !?
        It also works with instances of types
        """

        @test_nowarn inspect(1)

        """
        And for abstract types of course1
        """

        @test_nowarn inspect(T2)

        """
        finally, you can also use inspect to look at functions.
        Let's get a bit meta
        """

        @test_nowarn inspect(inspect)

        """
        as you can see inspect prints the docstrings of the methods it can find and lists some of
        the methods available.
        """
    end
end


text:
"""
    get_funcdef_expr(ex::Expr[, throws=true]) -> def::Expr, isinside::Bool, container::Expr

Given an expression `ex` that contains a single function definition,
return the "inner" expression `def` defining the function.  `def` may
be `ex`, or it may be `ex` stripped of "trivial" (single-expression)
`Expr(:block, ...)` wrappers or any macro calls such as `@inline` . If
container expressions are found, `isinside` will be `true` and
`container` will be the immediate "parent expression" of `def`---if
you rewrite `def`, then you should stash the result in the last of
`container`'s `args`.

If `throws = true` and the expression does not correspond to a
single function definition, an ArgumentError will be thrown.  When
`throws = false`, you should check `def` manually using
`is_funcdef_expr`.

# Examples:

```julia
julia> ex = :(foo(x) = 1)

julia> def, isinside, container = get_funcdef_expr(ex);

julia> def
:(foo(x) = begin
            1
        end)

julia> isinside
false

julia> container
quote
    nothing
end

julia> ex = quote
Base.@pure @inline function foo{T}(::Type{T})
    false
end
end

julia> def, isinside, container = get_funcdef_expr(ex);

julia> def
:(function foo{T}(::Type{T})
        false
    end)

julia> isinside
true

julia> container
:(@inline function foo{T}(::Type{T})
            false
        end)
"""
function get_funcdef_expr(ex::Expr, throws::Bool=true)
    def = ex
    isinside = false
    container = Expr(:block, nothing)
    while (def.head == :macrocall && isa(def.args[end], Expr)) || is_trivial_block_wrapper(def)
        container = def
        def = def.args[end]::Expr
        isinside = true
    end
    if throws && !is_funcdef_expr(def)
        throw(ArgumentError(string("expected function expression, got ", def)))
    end
    def, isinside, container
end

"""
    is_funcdef_expr(ex)

Return `true` if `ex` is an expression defining a function, either the
"long form" (`function`) or "short form" (`f(x) = 1`).

`ex` must be "stripped" of any macro calls or `Expr(:block, ...)`
wrappers; see `get_funcdef_expr`.
"""
function is_funcdef_expr(ex::Expr)
    ex.head == :function ||
        (ex.head == :(=) && isa(ex.args[1], Expr) && (ex.args[1]::Expr).head == :call)
end
is_funcdef_expr(arg::ANY) = false

"""
    funcdef_longform(ex)

Cannonicalizes a function-definition expression to "long form," like
```julia
function f(x)
    x^2
end
```
rather than
```julia
f(x) = x^2
```
Such expressions have head `:function` and two `args`, corresponding
to the signature and body, respectively.  In the returned expression,
the body is guaranteed to be a `:block` expression.
"""
function funcdef_longform(ex::Expr)
    if ex.head == :function
        return ex
    end
    if ex.head == :(=) && isa(ex.args[1], Expr) && (ex.args[1]::Expr).head == :call
        sig::Expr = ex.args[1]
        if !isa(ex.args[2], Expr)
            return Expr(:function, sig, Expr(:block, ex.args[2]))
        end
        body::Expr = ex.args[2]
        if body.head != :block
            body = Expr(:block, body)
        end
        return Expr(:function, sig, body)
    end
    throw(ArgumentError(string("expected a function definition, got ", ex)))
end

"""
    funcdef_name(sig)

Return the function's Symbol given a signature expression `sig`. `sig`
is the first argument of the function definition expression.
"""
function funcdef_name(sig::Expr)
    sig.head == :call || throw(ArgumentError(string("expected call expression, got ", sig)))
    decl = sig.args[1]
    isa(decl, Symbol) && return decl::Symbol
    if isa(decl, Expr) && (decl::Expr).head == :curly
        return ((decl::Expr).args[1])::Symbol
    end
    throw(ArgumentError(string("unexpected call expression ", sig)))
end

"""
    funcdef_params(sig)

Return a Vector{Any} of the function's parameters (the part inside the
curlies) given a signature expression `sig`. `sig` is the first
argument of the function definition expression.

The elements of the return may be Symbols or Exprs, where the latter
is used for expressions like `T<:Integer`.
"""
function funcdef_params(sig::Expr)
    sig.head == :call || throw(ArgumentError(string("expected call expression, got ", sig)))
    ret = []
    decl = sig.args[1]
    isa(decl, Symbol) && return ret
    if isa(decl, Expr) && (decl::Expr).head == :curly
        for i = 2:length((decl::Expr).args)
            push!(ret, (decl::Expr).args[i])
        end
        return ret
    end
    throw(ArgumentError(string("unexpected call expression ", sig)))
end

"""
    funcdef_args(sig)

Return a Vector{Any} of the function's argument expressions, which
includes both the variable name and any type declaration, if any.

The elements of the return may be Symbols or Exprs, where the latter
is used for expressions like `x::Integer`.

See also `funcdef_argnames` and `funcdef_argtypes`.
"""
function funcdef_args(sig::Expr)
    sig.head == :call || throw(ArgumentError(string("expected call expression, got ", sig)))
    return sig.args[2:end]
end

"""
    funcdef_argnames(sig)

Return a vector of the function's argument names.  These are typically
Symbols, except for non-canonical syntax like that of Traitor.jl.

See also `funcdef_args` and `funcdef_argtypeexprs`.
"""
function funcdef_argnames(sig::Expr)
    sig.head == :call || throw(ArgumentError(string("expected call expression, got ", sig)))
    return [argname(a) for a in sig.args[2:end]]
end

argname(s::Symbol) = s
function argname(ex::Expr)
    ex.head == :(::) || throw(ArgumentError("expected :(::) expression, got ", ex))
    ex.args[1]
end

"""
    funcdef_argtypeexprs(sig)

Return a vector of the function's argument type expressions.  These
may be Symbols (`:Int`) or expressions (`AbstractArray{T}`).

See also `funcdef_args` and `funcdef_argnames`.
"""
function funcdef_argtypeexprs(sig::Expr)
    sig.head == :call || throw(ArgumentError(string("expected call expression, got ", sig)))
    return [argtypeexpr(a) for a in sig.args[2:end]]
end

argtypeexpr(s::Symbol) = :Any
function argtypeexpr(ex::Expr)
    ex.head == :(::) || throw(ArgumentError("expected :(::) expression, got ", ex))
    ex.args[2]
end


text:
using SimilaritySearch
export fftraversal, fftclustering

function _ignore3(a, b, c)
end

"""
Selects a number of farthest points in `X`, using a farthest first traversal

- The callback function is called on each selected far point as `callback(centerID, dmax)` where `dmax` is the distance to the nearest previous reported center (the first is reported with typemax)
- The selected objects are far under the `dist` distance function with signature (T, T) -> Float64
- The number of points is determined by the stop criterion function with signature (Float64[], T[]) -> Bool
    - The first argument corresponds to the list of known distances (far objects)
    - The second argument corresponds to the database
- Check `criterions.jl` for basic implementations of stop criterions
- The callbackdist function is called on each distance evaluation between pivots and items in the dataset
    `callbackdist(index-pivot, index-item, distance)`
"""
function fftraversal(callback::Function, dist::Function, X::AbstractVector{T}, stop, callbackdist=_ignore3) where {T}
    N = length(X)
    D = Vector{Float64}(undef, N)
    dmaxlist = Float64[]
    dset = [typemax(Float64) for i in 1:N]
    imax::Int = rand(1:N)
    dmax::Float64 = typemax(Float64)
    if N == 0
        return
    end

    k::Int = 0

    @inbounds while k <= N
        k += 1
        pivot = X[imax]
        push!(dmaxlist, dmax)
        callback(imax, dmax)
        println(stderr, "computing fartest point $k, dmax: $dmax, imax: $imax, stop: $stop")
        dmax = 0.0
        ipivot = imax
        imax = 0

        D .= 0.0
        Threads.@threads for i in 1:N
            D[i] = dist(X[i], pivot)
        end

        for i in 1:N
            # d = dist(X[i], pivot)
            d = D[i]
            callbackdist(i, ipivot, d)

            if d < dset[i]
                dset[i] = d
            end

            if dset[i] > dmax
                dmax = dset[i]
                imax = i
            end
        end

        if dmax == 0 || stop(dmaxlist, X)
            break
        end
    end
end

"""
    fftclustering(dist::Function, X::Vector{T}, numcenters::Int, k::Int) where T

Clustering algorithm based on Farthest First Traversal (an efficient solution for the K-centers problem).
- `dist` distance function
- `X` contains the objects to be clustered
- `numcenters` number of centers to be computed
- `k` number of nearest references per object (k=1 makes a partition)

Returns a named tuple ``(NN, irefs, dmax)``.

- `NN` contains the ``k`` nearest references for each object in ``X``.
- `irefs` contains the list of centers (indexes to ``X``)
- `dmax` smallest distance among centers

"""
function fftclustering(dist::Function, X::Vector{T}, numcenters::Int; k::Int=1) where T
    # refs = Vector{Float64}[]
    irefs = Int[]
    NN = [KnnResult(k) for i in 1:length(X)]
    dmax = 0.0

    function callback(c, _dmax)
        push!(irefs, c)
        dmax = _dmax
    end

    function capturenn(i, refID, d)
        push!(NN[i], refID, d)
    end

    fftraversal(callback, dist, X, size_criterion(numcenters), capturenn)
    return (NN=NN, irefs=irefs, dmax=dmax)
end


text:
const out_dir = "../../docs/src/integrators/projection_methods/lotka_volterra_2d/"
const run_id  = "lotka_volterra_2d_vprk_radau"

# create list of figures to put in list
figures_custom = ("_momentum_error",)

include("lotka_volterra_2d_settings.jl")

# solver settings
set_config(:nls_atol, 8eps())
set_config(:nls_rtol, 2eps())
set_config(:nls_stol, 2eps())
set_config(:nls_nmax, 20)


text:
#=

Example demonstrating sparsity control by means of the `lazy_tolerance`-factor passed to the lazy AFW variant

A larger lazy_tolerance >= 1 favors sparsity by favoring optimization over the current active set rather than
adding a new FW vertex.

The default for AFW is lazy_tolerance = 2.0

=#

include("activate.jl")

using LinearAlgebra
using Random

n = Int(1e4)
k = 1000

s = rand(1:100)
@info "Seed $s"
Random.seed!(s)

xpi = rand(n);
total = sum(xpi);

# here the optimal solution lies in the interior if you want an optimal solution on a face and not the interior use:
# const xp = xpi;

const xp = xpi ./ total;

f(x) = norm(x - xp)^2
function grad!(storage, x)
    @. storage = 2 * (x - xp)
end

const lmo = FrankWolfe.KSparseLMO(100, 1.0)

## other LMOs to try
# lmo_big = FrankWolfe.KSparseLMO(100, big"1.0")
# lmo = FrankWolfe.LpNormLMO{Float64,5}(1.0)
# lmo = FrankWolfe.ProbabilitySimplexOracle(1.0);
# lmo = FrankWolfe.UnitSimplexOracle(1.0);

const x00 = FrankWolfe.compute_extreme_point(lmo, rand(n))

## example with BirkhoffPolytopeLMO - uses square matrix.
# const lmo = FrankWolfe.BirkhoffPolytopeLMO()
# cost = rand(n, n)
# const x00 = FrankWolfe.compute_extreme_point(lmo, cost)

function build_callback(trajectory_arr)
    return function callback(state)
        return push!(trajectory_arr, (Tuple(state)[1:5]..., state.active_set_length))
    end
end

FrankWolfe.benchmark_oracles(f, grad!, () -> randn(n), lmo; k=100)

x0 = deepcopy(x00)
@time x, v, primal, dual_gap, trajectorySs = FrankWolfe.frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Shortstep(),
    L=2,
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
);

trajectoryAda = []
callback = build_callback(trajectoryAda)

x0 = deepcopy(x00)
@time x, v, primal, dual_gap, _ = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    trajectory=true,
    callback=callback,
);

println("\n==> Lazy AFW.\n")

trajectoryAdaLoc15 = []
callback = build_callback(trajectoryAdaLoc15)

x0 = deepcopy(x00)
@time x, v, primal, dual_gap, _ = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    lazy=true,
    lazy_tolerance=1.5,
    trajectory=true,
    callback=callback,
);

trajectoryAdaLoc2 = []
callback = build_callback(trajectoryAdaLoc2)

x0 = deepcopy(x00)
@time x, v, primal, dual_gap, trajectoryAdaLoc2 = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    lazy=true,
    lazy_tolerance=2.0,
    trajectory=true,
    callback=callback,
);

trajectoryAdaLoc4 = []
callback = build_callback(trajectoryAdaLoc4)

x0 = deepcopy(x00)
@time x, v, primal, dual_gap, _ = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    lazy_tolerance=4.0,
    lazy=true,
    trajectory=true,
    callback=callback,
);

trajectoryAdaLoc10 = []
callback = build_callback(trajectoryAdaLoc10)

x0 = deepcopy(x00)
@time x, v, primal, dual_gap, _ = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=k,
    line_search=FrankWolfe.Adaptive(),
    print_iter=k / 10,
    emphasis=FrankWolfe.memory,
    lazy=true,
    lazy_tolerance=10.0,
    verbose=true,
    trajectory=true,
    callback=callback,
);

# Reduction primal/dual error vs. sparsity of solution

dataSparsity =
    [trajectoryAda, trajectoryAdaLoc15, trajectoryAdaLoc2, trajectoryAdaLoc4, trajectoryAdaLoc10]
labelSparsity = ["AFW", "LAFW-K-1.5", "LAFW-K-2.0", "LAFW-K-4.0", "LAFW-K-10.0"]

FrankWolfe.plot_sparsity(dataSparsity, labelSparsity, legend_position=:topright)


text:
function seat_id(line)
    arr::Array{Int} = [x in ('B', 'R') ? 1 : 0 for x in line]
    return sum([x * 2^(n - 1) for (n, x) in enumerate(arr[end:-1:1])])
end

function run(s)
    # Your code here
    seat_ids = sort([seat_id(line) for line in split(s, '\n')])
    for idx in 1:length(seat_ids)
        if seat_ids[idx + 1] != seat_ids[idx] + 1
            return seat_ids[idx] + 1
        end
    end
    return 0
end

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

function main()
    res, time, memory = @timed run(ARGS[1])
    println("_duration:$(time * 1000)")
    println(res)
end

main()


text:
module TensorToolbox

#Tensors in Tucker format + functions
using LinearAlgebra

import LinearAlgebra: norm, normalize, normalize!
import Base: +, -, *, ==, cat, display, dropdims, isequal, kron, ndims, parent, permutedims, show, size

include("helper.jl")
include("tensor.jl")
include("ttensor.jl")
include("ktensor.jl")
include("dimtree.jl")
include("htensor.jl")
include("TTtensor.jl")

end #module


text:
include("lsssolve.jl")
"""
Implements the NILSS algorithm due to Ni and Wang JCP 2017
Below, u -> F(u) is the time-1 primal map and u_trj refers
to a primal trajectory of length n.
Inputs:
	du_trj: dxdxn n-length timeseries of Jacobian matrices
			along u_trj
	X: dxn n-length timeseries of perturbation vectors.
			X[:,i] = dF/ds(u_trj[i,:]) which is in the
			tangent space of u_trj[i+1,:].
	f: dxn n-length timeseries of the neutral CLV. By
			assumption, the only neutral CLV is the RHS of
			the primal ode: du/dt = f(u,s).
			In case of maps, u_n = F(u_{n+1}), f is zeros.
	J:mxn n-length timeseries of m scalar functions J(u)
			evaluated along u_trj.
	dJ: mxdxn n-length timeseries of DJ(u) evaluated along
			u_trj
	s: array of parameters
	d_u: integer dimension of the unstable subspace. You can 		set a guess value to obtain correct sensitivities
		as long as the guess value is greater than or equal
		to the true d_u.
Outputs:
	vsh: dxn n-length timeseries of shadowing tangent
		vectors that are bounded solutions of the
		conventional tangent equation:
			v[:,i+1] = du_trj[:,i]*v[:,i] + X[:,i]

	dJds: scalar-valued sensitivity of interest:
		d/ds sum(J) ≈ d/ds <J, mu>, where mu is
		the SRB measure.
"""
function lss(du_trj, X, f, s, d_u)
	d, n = size(X)
	lyap_exps = zeros(d_u)
    R = zeros(d_u,d_u,n)
    Q = zeros(d,d_u,n)
	v = zeros(d,1,n) #assume one parameter
	v[:,1,1] = zeros(d)
	A = qr!(rand(d,d_u))
    Q[:,:,1] = Array(A.Q)
    R[:,:,1] = A.R

	b = zeros(d_u,1,n) #1 is the # parameters
	pf_Q = zeros(d_u,1,n) #1 is the # parameters
	pf_v = zeros(n)

	ff = zeros(n)
	if all(isapprox.(f,0.))
		ff = ones(n)
	else
		[ff[i] = sum(x -> x^2, f[:,i]) for i = 1:n]
	end
	pf_Q[:,1,1] = Q[:,:,1]'*f[:,1]
	[Q[:,j,1] = Q[:,j,1] - (pf_Q[j,1,1]*
	 		f[:,1]/ff[1]) for j=1:d_u]

	for i=2:n
	    v[:,:,i] = du_trj[:,:,i-1]*v[:,:,i-1] +
					X[:,i-1]
		Q[:,:,i] = du_trj[:,:,i-1]*Q[:,:,i-1]
		pf_Q[:,:,i] = Q[:,:,i]'*f[:,i]
		[Q[:,j,i] = Q[:,j,i] - pf_Q[j,1,i]*f[:,i]/
						ff[i] for j=1:d_u]
		A = qr!(Q[:,:,i])
        Q[:,:,i] = Array(A.Q)
        R[:,:,i] = A.R
		pf_v[i] = dot(v[:,1,i],f[:,i])
		v[:,:,i] = v[:,:,i] - (pf_v[i]/
							   ff[i]*f[:,i])
        b[:,:,i] = (Q[:,:,i]')*v[:,:,i]
		v[:,:,i] = v[:,:,i] - Q[:,:,i]*b[:,:,i]
		lyap_exps .+= log.(abs.(diag(R[:,:,i])))./n
    end
	#println(lyap_exps/0.01)

	b = reshape(collect(b),d_u, n)
	#println("Solving the least squares problem... ")
	a = lsssolve(R,b)
	#println("Computing shadowing direction...")
	v = reshape(collect(v),d,n)
	vsh = zeros(d, n)
	xi = zeros(n)
	vsh[:,1] = Q[:,:,1]*reshape(a[:,1], d_u, 1) +
				  v[:,1]
	for i = 2:n
		vsh[:,i] = v[:,i] + Q[:,:,i]*reshape(a[:,i],d_u,
					1)
		xi[i] = xi[i-1] + dot(-vsh[:,i] + du_trj[:,:,i-1]*
					vsh[:,i-1] + X[:,i-1], f[:,i])/ff[i]
	end
	return vsh, xi
end
function compute_sens(vsh, xi, dJ, f)
	#println("Computing sensitivity...")
	m, d, n = size(dJ)
	dJds = zeros(m)
	for i = 1:n
		dJds .= dJds .+ dJ[:,:,i]*(vsh[:,i] .+
					 xi[i]*f[:,i])/n
	end
	return dJds
end

max_stars_repo_name: nishaChandramoorthy/energies

max_issues_repo_name: nishaChandramoorthy/energies

max_forks_repo_name: nishaChandramoorthy/energies


text:
for i in 1:length(sol.t)
    ca1 = Tnt * sol.u[i]
    cb1 = chderive(Float64(0.0), ρmax, ca1, gridnum)
    cc1 = chderive(Float64(0.0), ρmax, cb1, gridnum)
    m1 = Tn * cb1
    m2 = Tn * cc1
    m1fun=Fun(S, ApproxFun.transform(S, m1))
    m2fun=Fun(S, ApproxFun.transform(S, m2))
    minpos=optimize(Vout, 0.0, ρmax).minimizer
    mass_σ[i]=sqrt(2*minpos*m2fun(minpos)*Λ^2)
end

mass_σ=similar(sol.t)

mass_σ

89*Λ^-2

plot([sol.u[2]])
plot(ρ,sol.u[1])

myf=Fun(S, ApproxFun.transform(S, sol.u[end-10]))

plot(myf'',0.0001,0.015)

plot(m1fun,0,0.03)

plot(Λ*exp.(sol.t), mass_σ)

max_stars_repo_name: Yangyang-Tan/FRG-Code-Collection

max_issues_repo_name: Yangyang-Tan/FRG-Code-Collection

max_forks_repo_name: Yangyang-Tan/FRG-Code-Collection


text:
function parse_line_auv(line::AbstractString, Δ::Float64)
    line = strip(line)
    elts = split(line," ")
    length(elts) != 3 && throw("Format auv expects lines with three elements.")
    t,u,v = elts
    t = parse(Float64,t)
    t0 = t - Δ / 2
    t1 = t + Δ /2
    t0,t1,u,v
end

parse_line_auv(line::AbstractString)=parse_line_auv(line,1.0)

function parse_line_abuv(line::AbstractString)
    line = strip(line)
    elts = split(line," ")
    length(elts) != 4 && throw("Format abuv expects lines with four elements.")
    t0,t1,u,v = elts
    t0 = parse(Float64,t0)
    t1 = parse(Float64,t1)
    t0,t1,u,v
end

parse_line_abuv(line::AbstractString,Δ::Float64)=parse_line_abuv(line)

parsers=Dict("auv"=>parse_line_auv,
             "abuv"=>parse_line_abuv
            )

function parse_line(line::AbstractString,format::AbstractString)
    haskey(parsers,format) ? parsers[format](line) : throw("Unknown format $format")
end

function parse_line(line::AbstractString,format::AbstractString,Δ::Float64)
    haskey(parsers,format) ? parsers[format](line,Δ) : throw("Unknown format $format")
end

function load!(s::AbstractStream, f::AbstractString, format::AbstractString)
    open(f) do file
        for line in eachline(file)
            t0,t1,u,v = parse_line(line,format)
            record!(s,t0,t1,u,v)
        end
    end
end

function load!(s::AbstractStream, f::AbstractString, format::AbstractString, Δ::Float64)
    open(f) do file
        for line in eachline(file)
            t0,t1,u,v = parse_line(line,format,Δ)
            record!(s,t0,t1,u,v)
        end
    end
end

function string(e::Event)
    if typeof(e)==LinkEvent
        e.arrive ? "$(e.t) + $(e.object[1]) $(e.object[2])" : "$(e.t) - $(e.object[1]) $(e.object[2])"
    elseif typeof(e)==NodeEvent
        e.arrive ? "$(e.t) + $(e.object)" : "$(e.t) - $(e.object)"
    else
        throw("Unknown type for event object.")
    end
end

function record_event!(events::Vector{Event},
                       D::Dict{Tuple{AbstractString,AbstractString},Array{Float64,1}},
                       t0::Float64,
                       t1::Float64,
                       u::AbstractString,
                       v::AbstractString)
    if haskey(D,(u,v))
        if t0 > D[(u,v)][2]
            push!(events,LinkEvent(D[(u,v)][1],true,(u,v)))
            push!(events,LinkEvent(D[(u,v)][2],false,(u,v)))
            D[(u,v)]=[t0,t1]
        else
            D[(u,v)][2]=t1
        end
    else
        D[(u,v)]=[t0,t1]
    end
end

function parse_to_events(f::AbstractString,format::AbstractString,Δ::Float64)
    format ∉ ["auv","abuv"] && throw("Only auv and abuv formats can be parsed to events.")
    D=Dict{Tuple{AbstractString,AbstractString},Array{Float64,1}}()
    events=Event[]
    open(f) do file
        for line in eachline(file)
            t0,t1,u,v = parse_line(line,format,Δ)
            record_event!(events,D,t0,t1,u,v)
        end
    end
    for (obj,interv) in D
        push!(events,LinkEvent(interv[1],true,obj))
        push!(events,LinkEvent(interv[2],false,obj))
    end
    sort!(events, by = v -> v.t)
end

function load!(tc::TimeCursor, f::AbstractString, format::AbstractString, Δ::Float64)
    events=parse_to_events(f,format,Δ)
    load!(tc,events)
end

function dump(f::AbstractString,events::Vector{Event})
    open(f,"w") do file
        for event in events
            write(file,"$(string(event))\n")
        end
    end
end


text:
export dual_set

# Additional dual_set
function dual_set(::MOI.GreaterThan{T}) where T
    return MOI.GreaterThan(zero(T))
end

function dual_set(::MOI.LessThan{T}) where T
    return MOI.LessThan(zero(T))
end

function dual_set(::MOI.EqualTo{T}) where T
    return # Maybe return Reals in the future
end


text:
"""
    Disjunction <: Property

Type that represents a disjunction of properties.

### Fields

- `disjuncts` -- vector of properties (elements are reordered by this type)
- `reorder`   -- flag to indicate whether shuffling is allowed

### Notes

The following formula characterizes whether a set ``X`` satisfies a disjunction
``𝑃 = 𝑃_1 ∨ 𝑃_2 ∨ … ∨ 𝑃_m``:

```math
    X \\models 𝑃 \\iff X \\models 𝑃_j \\text{ for some } 1 ≤ j ≤ m
```

If the `reorder` flag is set, the disjuncts may be reordered after each call to
[`check`](@ref check(𝑃::Disjunction, X::LazySet{N}) where {N<:Real}) as a
heuristics to make subsequent checks faster.
"""
struct Disjunction <: Property
    disjuncts::Vector{Property}
    reorder::Bool
end

# default constructor with activated reordering
Disjunction(disjuncts::Vector{<:Property}) = Disjunction(disjuncts, true)

"""
    dim(𝑃::Disjunction)::Int

Return the dimension of a disjunction of properties.

### Input

- `𝑃` -- disjunction of properties

### Output

The dimension of the disjunction of properties.
"""
function dim(𝑃::Disjunction)::Int
    if isempty(𝑃.disjuncts)
        throw(ArgumentError("the dimension of an empty disjunction of " *
                            "properties is undefined"))
    end
    @inbounds return dim(𝑃.disjuncts[1])
end

"""
    check(𝑃::Disjunction, X::LazySet{N}; witness::Bool=false) where {N<:Real}

Check whether a convex set satisfies a disjunction of properties.

### Input

- `𝑃`       -- disjunction of properties
- `X`       -- convex set
- `witness` -- (optional, default: `false`) flag for returning a counterexample
               if the property is violated

### Output

* If `witness` option is deactivated: `true` iff `X` satisfies the property `𝑃`
* If `witness` option is activated:
  * `(true, [])` iff `X` satisfies the property `𝑃`
  * `(false, v)` iff `X` does not satisfy the property `𝑃` with witness `v`;
    note that `v == N[]` if 𝑃 is the empty disjunction

### Notes

By convention, the empty disjunction is equivalent to `false` and hence is
satisfied by no set.

If the `𝑃.reorder` flag is set, the disjuncts may be reordered as a heuristics
to make subsequent checks faster.
Since we check satisfaction from left to right, we move the disjunct for which
satisfaction was established to the front.

To be consistent with other propertes, the `witness` option only returns *one*
counterexample, namely for the left-most disjunct in the `disjuncts` vector.
We deactivate witness production for checking the remaining disjuncts.
"""
function check(𝑃::Disjunction, X::LazySet{N};
               witness::Bool=false) where {N<:Real}
    v = N[]
    create_witness = witness
    for (i, conjunct) in enumerate(𝑃.disjuncts)
        result = check(conjunct, X; witness=create_witness)
        if (create_witness && result[1]) || result
            _reorder!(𝑃, i)
            return witness ? (true, N[]) : true
        elseif create_witness
            v = result[2]
            # deactivate witness production for remaining checks
            create_witness = false
        end
    end
    return witness ? (false, v) : false
end

function _reorder!(𝑃::Disjunction, i::Int)
    if !𝑃.reorder || i == 1
        return nothing
    end
    first = 𝑃.disjuncts[i]
    while i > 1
        𝑃.disjuncts[i] = 𝑃.disjuncts[i-1]
        i -= 1
    end
    𝑃.disjuncts[1] = first
    return nothing
end


text:
# GLVisualize throws a lot of information in your face. I don't want that, so
# here's a Container with a show() to supress that.
mutable struct Plot3DContainer
    window::GLWindow.Screen
    cam::Camera
    hotkeys::Hotkeys
end
function show(io::IO, c::Plot3DContainer)
    print(io, "Plot3DContainer with\n")
    print(io, "\twindow\n")
    print(io, "\tcam\n")
    print(io, "\thotkeys\n")
end

function glsphere(;pos=Point3f0(0f0), r=0.1f0, quality=24)
    GLNormalMesh(HyperSphere(pos, r), quality)
end

function plot3D(
        lattice::LatticeType,
        args...;
        # Window setup
        window::GLWindow.Screen = init_window(),
        # Camera setup
        rot_v::Signal{Float32} = Signal(0.05f0),
        trans_v::Signal{Float32} = Signal(4f0),
        cam_position::Vec3f0 = Vec3f0(2),
        lookat::Vec3f0 = Vec3f0(0),
        up::Vec3f0 = Vec3f0(0., 0., 1.),
        cam::Camera = camera(
            window,
            rot_v = rot_v,
            trans_v = trans_v,
            position = cam_position,
            lookat = lookat,
            up = up
        ),
        # Visualization setup
        primitive::AbstractMesh = glsphere(),
        site_color::RGBA{Float32} = RGBA{Float32}(0.2, 0.4, 0.8, 1.0),
        bond_color::RGBA{Float32} = RGBA{Float32}(0.2, 0.2, 0.2, 1.0),
        bond_thickness::Float32 = 5f0,
        kwargs...
    ) where {LatticeType <: AbstractGraph}

    @assert dims(lattice.unitcell) == 3 "plot3D only works for 3D lattices."

    # Controls
    hotkeys = Hotkeys(window)
    init_screenshots!(hotkeys)
    angles = init_angles(hotkeys)
    init_QWEASD_controls!(hotkeys, cam, angles)

    # Recover node positions
    positions = [
        GeometryTypes.Point3f0(get_pos(
            lattice.unitcell,
            node.uvw[1],
            node.uvw[2:end]
        )) for node in lattice.nodes
    ]

    # Gather edges/bonds
    edge_positions = GeometryTypes.Point3f0[]
    for level in lattice.edges
        for edge in level
            push!(
                edge_positions,
                GeometryTypes.Point3f0(positions[edge.node1.uvw...])
            )
            push!(
                edge_positions,
                GeometryTypes.Point3f0(positions[edge.node2.uvw...])
            )
        end
    end

    # Render sites
    _view(
        visualize(
            (primitive, positions[:]),
            color = site_color
        ),
        window,
        camera = cam
    )

    # Render bonds
    _view(
        visualize(
            edge_positions,
            :linesegment,
            color = bond_color,
            thickness = bond_thickness
        ),
        window,
        camera = cam
    )

    center!(window, :custom)

    Plot3DContainer(window, cam, hotkeys)
end


text:
using VectorizationBase
using Test

const W64 = VectorizationBase.REGISTER_SIZE ÷ sizeof(Float64)
const W32 = VectorizationBase.REGISTER_SIZE ÷ sizeof(Float32)

A = randn(13, 17); L = length(A); M, N = size(A);

@testset "VectorizationBase.jl" begin
    # Write your own tests here.
@test isempty(detect_unbound_args(VectorizationBase))

@test first(A) === A[1]
@testset "Struct-Wrapped Vec" begin
@test extract_data(zero(SVec{4,Float64})) === (VE(0.0),VE(0.0),VE(0.0),VE(0.0)) === extract_data(SVec{4,Float64}(0.0))
@test extract_data(one(SVec{4,Float64})) === (VE(1.0),VE(1.0),VE(1.0),VE(1.0)) === extract_data(SVec{4,Float64}(1.0)) === extract_data(extract_data(SVec{4,Float64}(1.0)))
v = SVec((VE(1.0),VE(2.0),VE(3.0),VE(4.0)))
@test v === SVec{4,Float64}(1, 2, 3, 4) === conj(v) === v'
@test length(v) == 4 == first(size(v))
@test eltype(v) == Float64
for i in 1:4
    @test i == v[i]
    # @test i === SVec{4,Int}(v)[i] # should use fptosi (ie, vconvert defined in SIMDPirates).
end
@test zero(v) === zero(typeof(v))
@test one(v) === one(typeof(v))
# @test SVec{W32,Float32}(one(SVec{W32,Float64})) === SVec(one(SVec{W32,Float32})) === one(SVec{W32,Float32}) # conversions should be tested in SIMDPirates
    @test firstval(v) === firstval(extract_data(v)) === 1.0
    @test SVec{1,Int}(1) === SVec{1,Int}((Core.VecElement(1),))
end

@testset "alignment.jl" begin

@test all(i -> VectorizationBase.align(i) == VectorizationBase.REGISTER_SIZE, 1:VectorizationBase.REGISTER_SIZE)
@test all(i -> VectorizationBase.align(i) == 2VectorizationBase.REGISTER_SIZE, 1+VectorizationBase.REGISTER_SIZE:2VectorizationBase.REGISTER_SIZE)
@test all(i -> VectorizationBase.align(i) == 10VectorizationBase.REGISTER_SIZE, (1:VectorizationBase.REGISTER_SIZE) .+ 9VectorizationBase.REGISTER_SIZE)

@test all(i -> VectorizationBase.align(reinterpret(Ptr{Cvoid}, i)) == reinterpret(Ptr{Cvoid},   VectorizationBase.REGISTER_SIZE), 1:VectorizationBase.REGISTER_SIZE)
@test all(i -> VectorizationBase.align(reinterpret(Ptr{Cvoid}, i)) == reinterpret(Ptr{Cvoid},  2VectorizationBase.REGISTER_SIZE), 1+VectorizationBase.REGISTER_SIZE:2VectorizationBase.REGISTER_SIZE)
@test all(i -> VectorizationBase.align(reinterpret(Ptr{Cvoid}, i)) == reinterpret(Ptr{Cvoid}, 20VectorizationBase.REGISTER_SIZE), (1:VectorizationBase.REGISTER_SIZE) .+ 19VectorizationBase.REGISTER_SIZE)

@test all(i -> VectorizationBase.align(i,W32) == VectorizationBase.align(i,Float32) == VectorizationBase.align(i,Int32) == W32*cld(i,W32), 1:VectorizationBase.REGISTER_SIZE)
@test all(i -> VectorizationBase.align(i,W32) == VectorizationBase.align(i,Float32) == VectorizationBase.align(i,Int32) == W32*cld(i,W32), 1+VectorizationBase.REGISTER_SIZE:2VectorizationBase.REGISTER_SIZE)
@test all(i -> VectorizationBase.align(i,W32) == VectorizationBase.align(i,Float32) == VectorizationBase.align(i,Int32) == W32*cld(i,W32), (1:VectorizationBase.REGISTER_SIZE) .+ 29VectorizationBase.REGISTER_SIZE)

@test all(i -> VectorizationBase.align(i,W64) == VectorizationBase.align(i,Float64) == VectorizationBase.align(i,Int64) == W64*cld(i,W64), 1:VectorizationBase.REGISTER_SIZE)
@test all(i -> VectorizationBase.align(i,W64) == VectorizationBase.align(i,Float64) == VectorizationBase.align(i,Int64) == W64*cld(i,W64), 1+VectorizationBase.REGISTER_SIZE:2VectorizationBase.REGISTER_SIZE)
@test all(i -> VectorizationBase.align(i,W64) == VectorizationBase.align(i,Float64) == VectorizationBase.align(i,Int64) == W64*cld(i,W64), (1:VectorizationBase.REGISTER_SIZE) .+ 29VectorizationBase.REGISTER_SIZE)

@test reinterpret(Int, VectorizationBase.align(pointer(A))) % VectorizationBase.REGISTER_SIZE === 0

@test all(i -> VectorizationBase.aligntrunc(i) == 0, 0:VectorizationBase.REGISTER_SIZE-1)
@test all(i -> VectorizationBase.aligntrunc(i) == VectorizationBase.REGISTER_SIZE, VectorizationBase.REGISTER_SIZE:2VectorizationBase.REGISTER_SIZE-1)
@test all(i -> VectorizationBase.aligntrunc(i) == 9VectorizationBase.REGISTER_SIZE, (0:VectorizationBase.REGISTER_SIZE-1) .+ 9VectorizationBase.REGISTER_SIZE)

@test all(i -> VectorizationBase.aligntrunc(i,W32) == VectorizationBase.aligntrunc(i,Float32) == VectorizationBase.aligntrunc(i,Int32) == W32*div(i,W32), 1:VectorizationBase.REGISTER_SIZE)
@test all(i -> VectorizationBase.aligntrunc(i,W32) == VectorizationBase.aligntrunc(i,Float32) == VectorizationBase.aligntrunc(i,Int32) == W32*div(i,W32), 1+VectorizationBase.REGISTER_SIZE:2VectorizationBase.REGISTER_SIZE)
@test all(i -> VectorizationBase.aligntrunc(i,W32) == VectorizationBase.aligntrunc(i,Float32) == VectorizationBase.aligntrunc(i,Int32) == W32*div(i,W32), (1:VectorizationBase.REGISTER_SIZE) .+ 29VectorizationBase.REGISTER_SIZE)

@test all(i -> VectorizationBase.aligntrunc(i,W64) == VectorizationBase.aligntrunc(i,Float64) == VectorizationBase.aligntrunc(i,Int64) == W64*div(i,W64), 1:VectorizationBase.REGISTER_SIZE)
@test all(i -> VectorizationBase.aligntrunc(i,W64) == VectorizationBase.aligntrunc(i,Float64) == VectorizationBase.aligntrunc(i,Int64) == W64*div(i,W64), 1+VectorizationBase.REGISTER_SIZE:2VectorizationBase.REGISTER_SIZE)
@test all(i -> VectorizationBase.aligntrunc(i,W64) == VectorizationBase.aligntrunc(i,Float64) == VectorizationBase.aligntrunc(i,Int64) == W64*div(i,W64), (1:VectorizationBase.REGISTER_SIZE) .+ 29VectorizationBase.REGISTER_SIZE)

a = Vector{Float64}(undef, 0)
ptr = pointer(a)
@test UInt(VectorizationBase.align(ptr, 1 << 12)) % (1 << 12) == 0
end

    @testset "masks.jl" begin
    @test Mask{8,UInt8}(0x0f) === @inferred Mask(0x0f)
    @test Mask{16,UInt16}(0x0f0f) === @inferred Mask(0x0f0f)
    @test Mask{8,UInt8}(0xff) == mask(Val(8), 0)
    @test Mask{8,UInt8}(0xff) == mask(Val(8), 8)
    @test Mask{8,UInt8}(0xff) == mask(Val(8), 16)
    @test Mask{8,UInt8}(0xff) == mask(Val(8), VectorizationBase.Static(0))
    @test Mask{16,UInt16}(0xffff) == mask(Val(16), 0)
    @test Mask{16,UInt16}(0xffff) == mask(Val(16), 16)
    @test Mask{16,UInt16}(0xffff) == mask(Val(16), 32)
@test all(w -> VectorizationBase.mask_type(w) == UInt8, 1:8)
@test all(w -> VectorizationBase.mask_type(w) == UInt16, 9:16)
@test all(w -> VectorizationBase.mask_type(w) == UInt32, 17:32)
@test all(w -> VectorizationBase.mask_type(w) == UInt64, 33:64)
@test all(w -> VectorizationBase.mask_type(w) == UInt128, 65:128)
if VectorizationBase.REGISTER_SIZE == 64 # avx512
    @test VectorizationBase.mask_type(Float16) == UInt32
    @test VectorizationBase.mask_type(Float32) == UInt16
    @test VectorizationBase.mask_type(Float64) == UInt8
    @test VectorizationBase.max_mask(Float16) == 0xffffffff # 32
    @test VectorizationBase.max_mask(Float32) == 0xffff     # 16
    @test VectorizationBase.max_mask(Float64) == 0xff       # 8
elseif VectorizationBase.REGISTER_SIZE == 32 # avx or avx2
    @test VectorizationBase.mask_type(Float16) == UInt16
    @test VectorizationBase.mask_type(Float32) == UInt8
    @test VectorizationBase.mask_type(Float64) == UInt8
    @test VectorizationBase.max_mask(Float16) == 0xffff     # 16
    @test VectorizationBase.max_mask(Float32) == 0xff       # 8
    @test VectorizationBase.max_mask(Float64) == 0x0f       # 4
elseif VectorizationBase.REGISTER_SIZE == 16 # sse
    @test VectorizationBase.mask_type(Float16) == UInt8
    @test VectorizationBase.mask_type(Float32) == UInt8
    @test VectorizationBase.mask_type(Float64) == UInt8
    @test VectorizationBase.max_mask(Float16) == 0xff       # 8
    @test VectorizationBase.max_mask(Float32) == 0x0f       # 4
    @test VectorizationBase.max_mask(Float64) == 0x03       # 2
end
@test all(w -> bitstring(VectorizationBase.mask(Val( 8), w)) == reduce(*, ( 8 - i < w ? "1" : "0" for i in 1:8 )), 1:8 )
@test all(w -> bitstring(VectorizationBase.mask(Val(16), w)) == reduce(*, (16 - i < w ? "1" : "0" for i in 1:16)), 1:16)
@test all(w -> VectorizationBase.mask(Float64, w) === VectorizationBase.mask(VectorizationBase.pick_vector_width_val(Float64), w), 1:W64)
end

@testset "number_vectors.jl" begin
# eval(VectorizationBase.num_vector_load_expr(@__MODULE__, :(size(A)), 8)) # doesn't work?
@test VectorizationBase.length_loads(A, Val(8)) == eval(VectorizationBase.num_vector_load_expr(@__MODULE__, :((() -> 13*17)()), 8)) == eval(VectorizationBase.num_vector_load_expr(@__MODULE__, 13*17, 8)) == divrem(length(A), 8)
@test VectorizationBase.size_loads(A,1, Val(8)) == eval(VectorizationBase.num_vector_load_expr(@__MODULE__, :((() -> 13   )()), 8)) == eval(VectorizationBase.num_vector_load_expr(@__MODULE__, 13   , 8)) == divrem(size(A,1), 8)
@test VectorizationBase.size_loads(A,2, Val(8)) == eval(VectorizationBase.num_vector_load_expr(@__MODULE__, :((() ->    17)()), 8)) == eval(VectorizationBase.num_vector_load_expr(@__MODULE__,    17, 8)) == divrem(size(A,2), 8)
end

@testset "vector_width.jl" begin
@test all(VectorizationBase.ispow2, 0:1)
@test all(i -> !any(VectorizationBase.ispow2, 1+(1 << (i-1)):(1 << i)-1 ) && VectorizationBase.ispow2(1 << i), 2:9)
@test all(i ->  VectorizationBase.intlog2(1 << i) == i, 0:(Int == Int64 ? 53 : 30))
FTypes = (Float16, Float32, Float64)
Wv = ntuple(i -> VectorizationBase.REGISTER_SIZE >> i, Val(3))
for (T, N) in zip(FTypes, Wv)
    W, Wshift = VectorizationBase.pick_vector_width_shift(:IGNORE_ME, T)
    @test W == 1 << Wshift == VectorizationBase.pick_vector_width(T) == N == VectorizationBase.pick_vector_width(:IGNORE_ME, T)
    @test Vec{W,T} == VectorizationBase.pick_vector(Val(W), T) == VectorizationBase.pick_vector(T)
    @test W == VectorizationBase.pick_vector_width(Val(W), T)
    @test Val(W) === VectorizationBase.pick_vector_width_val(Val(W), T) == VectorizationBase.pick_vector_width_val(T)
    while true
        W >>= 1
        W == 0 && break
        W2, Wshift2 = VectorizationBase.pick_vector_width_shift(W, T)
        @test W2 == 1 << Wshift2 == VectorizationBase.pick_vector_width(W, T) == VectorizationBase.pick_vector_width(Val(W),T)  == W
        @test Val(W) === VectorizationBase.pick_vector_width_val(Val(W), T)
        for n in W+1:2W
            W3, Wshift3 = VectorizationBase.pick_vector_width_shift(n, T)
            @test W2 << 1 == W3 == 1 << (Wshift2+1) == 1 << Wshift3 == VectorizationBase.pick_vector_width(n, T) == VectorizationBase.pick_vector_width(Val(n),T) == W << 1
            @test VectorizationBase.pick_vector(Val(W), T) == VectorizationBase.pick_vector(W, T) == Vec{W,T}
        end
    end
end

@test all(i -> VectorizationBase.nextpow2(i) == i, 0:2)
for j in 1:10
    l, u = (1<<j)+1, 1<<(j+1)
    @test all(i -> VectorizationBase.nextpow2(i) == u, l:u)
end

end

@testset "StridedPointer" begin
A = reshape(collect(Float64(0):Float64(63)), (16, 4))
ptr_A = pointer(A)
vA = VectorizationBase.stridedpointer(A)
Att = copy(A')'
vAtt = VectorizationBase.stridedpointer(Att)
@test eltype(vA) == Float64
@test Base.unsafe_convert(Ptr{Float64}, vA) === ptr_A === pointer(vA)
@test vA == VectorizationBase.stridedpointer(vA)
@test all(i -> A[i+1] === VectorizationBase.vload(ptr_A + 8i) === VectorizationBase.vload(vA, (i,)) === Float64(i), 0:15)
VectorizationBase.vstore!(vA, 99.9, (3,))
@test 99.9 === VectorizationBase.vload(ptr_A + 8*3) === VectorizationBase.vload(vA, (VectorizationBase.Static(3),)) === VectorizationBase.vload(vA, (3,0)) === A[4,1]
VectorizationBase.vstore!(vAtt, 99.9, (3,1))
@test 99.9 === VectorizationBase.vload(vAtt, (3,1)) === VectorizationBase.vload(vAtt, (VectorizationBase.Static(3),1)) === Att[4,2]
VectorizationBase.vnoaliasstore!(ptr_A+8*4, 999.9)
@test 999.9 === VectorizationBase.vload(ptr_A + 8*4) === VectorizationBase.vload(pointer(vA), 4*sizeof(eltype(A))) === VectorizationBase.vload(vA, (4,))
@test vload(vA, (7,2)) == vload(vAtt, (7,2)) == A[8,3]
@test vload(VectorizationBase.subsetview(vA, Val(1), 7), (2,)) == vload(VectorizationBase.subsetview(vAtt, Val(1), 7), (2,)) == A[8,3]
@test vload(VectorizationBase.subsetview(vA, Val(2), 2), (7,)) == vload(VectorizationBase.subsetview(vAtt, Val(2), 2), (7,)) == A[8,3]
    @test vload(VectorizationBase.double_index(vA, Val(0), Val(1)), (2,)) == vload(VectorizationBase.double_index(vA, Val(0), Val(1)), (VectorizationBase.Static(2),)) == A[3,3]
    @test vload(VectorizationBase.double_index(vAtt, Val(0), Val(1)), (1,)) == vload(VectorizationBase.double_index(vAtt, Val(0), Val(1)), (VectorizationBase.Static(1),)) == A[2,2]
    B = rand(5, 5)
vB = VectorizationBase.stridedpointer(B)
@test vB[1, 2] == B[2, 3] == vload(VectorizationBase.stridedpointer(B, 2, 3))
@test vB[3] == B[4] == vload(VectorizationBase.stridedpointer(B, 4))
@test vload(SVec{4,Float64}, vB) == SVec{4,Float64}(ntuple(i->B[i], Val(4)))
end

end


text:
using MixedModels
using Test

using MixedModels: dataset, likelihoodratiotest
using GLM: ProbitLink
using StatsModels: lrtest, isnested

include("modelcache.jl")

@testset "isnested" begin
    slp = dataset(:sleepstudy)

    # these tests don't actually depend on the models being fit,
    # so we just construct them

    # mismatched RE terms
    m1 = LinearMixedModel(@formula(reaction ~ 1 + days + (1+days|subj)), slp)
    m2 = LinearMixedModel(@formula(reaction ~ 1 + days + (0+days|subj)), slp)
    @test !isnested(m1, m2)

    # mismatched FE
    m1 = LinearMixedModel(@formula(reaction ~ 1 + days + (1|subj)), slp)
    m2 = LinearMixedModel(@formula(reaction ~ 0 + days + (1|subj)), slp)
    @test !isnested(m1, m2)

    # mismatched grouping vars
    kb07  = dataset(:kb07)
    m1 = LinearMixedModel(@formula(rt_trunc ~ 1 + (1|subj)), kb07)
    m2 = LinearMixedModel(@formula(rt_trunc ~ 1 + (1|item)), kb07)
    @test !isnested(m1, m2)

    # fixed-effects specification in REML and
    # conversion of internal ArgumentError into @error for StatsModels.isnested
    kb07  = dataset(:kb07)
    m1 = fit(MixedModel, @formula(rt_trunc ~ 1 + prec + (1|subj)), kb07, REML=true)
    m2 = fit(MixedModel, @formula(rt_trunc ~ 1 + prec + (1+prec|subj)), kb07, REML=true)
    @test isnested(m1, m2)
    m2 = fit(MixedModel, @formula(rt_trunc ~ 1 + (1+prec|subj)), kb07, REML=true)
    @test !isnested(m1, m2)

end

@testset "likelihoodratio test" begin
    slp = dataset(:sleepstudy);

    fm0 = fit(MixedModel,@formula(reaction ~ 1 + (1+days|subj)),slp);
    fm1 = fit(MixedModel,@formula(reaction ~ 1 + days + (1+days|subj)),slp);

    lrt = likelihoodratiotest(fm0,fm1);

    @test [deviance(fm0), deviance(fm1)] == lrt.deviance
    @test deviance(fm0) - deviance(fm1) == only(lrt.tests.deviancediff)
    @test only(lrt.tests.dofdiff) == 1
    @test sum(map(length,lrt.tests)) == 3
    @test sum(map(length,lrt.pvalues)) == 1
    @test sum(map(length,lrt.models)) == 4
    @test length(lrt.formulae) == 2
    show(IOBuffer(),lrt);
    @test :pvalues in propertynames(lrt)

    # mix of REML and ML
    fm0 = fit(MixedModel,@formula(reaction ~ 1 + (1+days|subj)),slp, REML=true);
    @test_throws ArgumentError likelihoodratiotest(fm0,fm1)

    # differing FE with REML
    fm1 = fit(MixedModel,@formula(reaction ~ 1 + days + (1+days|subj)),slp, REML=true);

    @test_throws ArgumentError likelihoodratiotest(fm0,fm1)

    contra = MixedModels.dataset(:contra);
    gm0 = fit(MixedModel, @formula(use ~ 1+age+urban+livch+(1|urban&dist)), contra, Bernoulli(), fast=true);
    gm1 = fit(MixedModel, @formula(use ~ 1+age+abs2(age)+urban+livch+(1|urban&dist)), contra, Bernoulli(), fast=true);
    lrt = likelihoodratiotest(gm0,gm1);
    @test [deviance(gm0), deviance(gm1)] == lrt.deviance
    @test deviance(gm0) - deviance(gm1) == only(lrt.tests.deviancediff)
    @test first(lrt.tests.dofdiff) == 1
    @test sum(length, lrt.tests) == 3
    @test sum(length, lrt.pvalues) == 1
    @test sum(length, lrt.models) == 4
    @test length(lrt.formulae) == 2

    # mismatched links
    gm_probit = fit(MixedModel, @formula(use ~ 1+age+urban+livch+(1|urban&dist)), contra, Bernoulli(), ProbitLink(), fast=true);
    @test_throws ArgumentError likelihoodratiotest(gm0,gm_probit)

    # mismatched families
    gm_poisson = fit(MixedModel, @formula(use ~ 1+age+urban+livch+(1|urban&dist)), contra, Poisson(), fast=true);
    @test_throws ArgumentError MixedModels.likelihoodratiotest(gm0,gm_poisson)
end


text:
################################################################################
## Copyright 2016. Juri Hinz and Jeremy Yee.
################################################################################

module ConvexSwitch

## Load required packages
using Clustering
using NearestNeighbors

## Export immutable types
export CSS, Bellman, Path, PermMat, SGrid, PathPolicy, Mart, OptimalBounds

## Export functions
export disturb_mat, path_mat, subsim_mat # generate disturb matrices
export row_rearrange, expected, expected_host, expected_fast # Bellman recursion
export test_policy # test run of prescribed policy

## Include source files
include("functions_mat.jl")
include("CSS.jl")
include("Path.jl")
include("PermMat.jl")
include("functions_bellman.jl")
include("SGrid.jl")
include("Bellman.jl")
include("PathPolicy.jl")
include("test_policy.jl")
include("Mart.jl")
include("OptimalBounds.jl")

end # module


text:
using WingBiology

afl350 = Airfoil(
    0.5751,
    4.6252,
    -0.1444,
    0.0,
    [0.012971116342921363, -0.019178413993406807, 0.013755755841287076],
    -Inf,
    1.3697
)
afl350_A = 0.08181

afl1 = Airfoil(
    0.4006,
    4.2844,
    -0.05814,
    0.0,
    [0.029996785640005186, -0.021401181166492478, 0.009940000385963543],
    -Inf,
    1.1178
)
afl1_A = 0.05734

afl1_inv = Airfoil(
    -0.3732,
    4.7922,
    0.0543,
    0.0,
    [0.029996785640005186, 0.021401181166492478, 0.009940000385963543],
    -1.1178,
    0.78
)

afl216 = Airfoil(
    0.3331,
    5.9603,
    -0.0077,
    0.0,
    [0.000854, 0.0028762, 0.0064887],
    -1.1359,
    1.3482
)
afl216_A = 0.07457

GAW1 = Airfoil(
    0.5830,
    4.712,
    -0.1236,
    4.7117,
    [0.03629806248179257, -0.04975965925122239, 0.025467486439728848],
    -Inf,
    1.6642
)
GAW1_A = 0.1157

max_stars_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/WingBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/WingBiology


text:
using SEC
using Test

@testset "Hello World " begin
    @test isnothing(SEC.Calibrate.greet())
    @test isnothing(SEC.Emulate.greet())
    @test isnothing(SEC.Sample.greet())
end

max_stars_repo_name: bischtob/SEC

max_issues_repo_name: bischtob/SEC

max_forks_repo_name: bischtob/SEC


text:
using Images, ImageView, ImageMorphology

# examples on using tophat and bothat
geom_img = load("sample-images/geometrical-figures-and-noise.jpg");
geom_img_gray = Gray.(geom_img);
geom_img_th = tophat(geom_img_gray)
geom_img_bh = bothat(geom_img_gray)

geom_img_new = zeros(ColorTypes.Gray{FixedPointNumbers.Normed{UInt8,8}}, size(geom_img_gray) .* (1, 2));
geom_img_new_center = Int(size(geom_img_gray, 2))

geom_img_new[:, 1:geom_img_new_center] = geom_img_th
geom_img_new[:, geom_img_new_center:end - 1] = geom_img_bh
geom_img_new[:, geom_img_new_center] = 1

imshow(Gray.(geom_img_new))

# adjusting image contrast
minmax = scaleminmax(Gray, 0, 1)

img = load("sample-images/fabio_gray_512.png");
img_gray = Gray.(img);
img_new = minmax.(img_gray + tophat(img_gray) - bothat(img_gray))

img_center = Int(size(img_gray, 2) / 2)
img_gray[:, img_center:end] = img_new[:, img_center:end]
img_gray[:, img_center] = 0
imshow(img_gray)

img = load("sample-images/busan-night-scene.jpg");
img_gray = Gray.(img);
img_new = minmax.(img_gray + tophat(img_gray) * 0.5 - bothat(img_gray) * 0.5)

img_center = Int(size(img_gray, 2) / 4)
img_gray[:, img_center:end] = img_new[:, img_center:end]
img_gray[:, img_center] = 0
imshow(img_gray)


text:
module DoStats

include("var0.jl")
export sum0, mean0, sd0, var0

end


text:
using LinearAlgebra
using CompScienceMeshes
using SauterSchwab3D

const pI   = point(6,5,3)
const pII  = point(5,2,3)
const pIII = point(7,1,0)
const pIV  = point(3,4,0)

const qII  = point(3,8,4)
const qIII = point(7,5,6)
const qIV  = point(9,8,5)

const P = simplex(pI,pII,pIII,pIV)
const Q = simplex(pI,qII,qIII,qIV)

const sing = SauterSchwab3D.singularity_detection(P,Q)

Accuracy2 = 16
cv_ref = CommonVertex6D(sing,SauterSchwab3D._legendre(Accuracy2,0.0,1.0))

function integrand(x,y)
   ((x-pI)'*(y-pIV))*return(exp(-im*1*norm(x-y))/(4pi*norm(x-y)))
end

function INTEGRAND(u,v)
   n1 = neighborhood(P,u)
   n2 = neighborhood(Q,v)
   x = cartesian(n1)
   y = cartesian(n2)
   output = integrand(x,y)*jacobian(n1)*jacobian(n2)
   return(output)
end

print("Ref: ")
ref = sauterschwab_parameterized(INTEGRAND, cv_ref)
println(ref)
println()

res_tp =[]
res_sp =[]
res_gm = []
n1 = []
n2 = []
n3 = []

#Grauss tensor product
for i in 2:1:14
   Accuracy = i
   cv = CommonVertex6D(sing,SauterSchwab3D._legendre(Accuracy,0.0,1.0))

   int_tp = sauterschwab_parameterized(INTEGRAND, cv)

   num_pts = 2*length(cv.qps)^6
   push!(n1,num_pts)

   push!(res_tp,int_tp)
end

#ShunnHam simplex product
for i in 2:1:7
   Accuracy = i
   cv_s = CommonVertex6D_S(sing,SauterSchwab3D._shunnham3D(Accuracy))

   int_sp = sauterschwab_parameterized(INTEGRAND, cv_s)

   println(length(cv_s.qps))
   num_pts = 2*length(cv_s.qps)^2
   push!(n2,num_pts)
   push!(res_sp,int_sp)
end

#GrundmannMoeller simplex product
for i in 2:1:12
   Accuracy = i
   cv_gm = CommonVertex6D_S(sing,SauterSchwab3D._grundmannMoeller3D(2*Accuracy-1))

   int_gm = sauterschwab_parameterized(INTEGRAND, cv_gm)

   println(length(cv_gm.qps))
   num_pts = 2*length(cv_gm.qps)^2
   push!(n3,num_pts)
   push!(res_gm,int_gm)
end

err_tp = norm.(res_tp.-ref)/norm(ref)
err_sp = norm.(res_sp.-ref)/norm(ref)
err_gm = norm.(res_gm.-ref)/norm(ref)

println(err_tp)
println(err_sp)
println(err_gm)

using Plots
plot( yaxis=:log, xaxis=:log, fontfamily="Times")
plot!(n1,err_tp, label="Gauss Tensor-Product",markershape=:circle)
plot!(n2,err_sp, label="Simplex Tensor-Product",markershape=:rect)
#plot!(n3,err_gm, label="Simplex-Product GM",markershape=:x)
plot!(xlims=(1e1,10^6),ylims=(1e-7,1))
plot!(xlabel="#Quad. points/Func. evals.", ylabel="Rel. Error.", title="Common Vertex 6D",legend=:bottomleft)

using BenchmarkTools

ref = 0.003477333267758975 - 0.5997621570461863im

cv = CommonVertex6D(sing,SauterSchwab3D._legendre(6,0.0,1.0))

int_tp = sauterschwab_parameterized(INTEGRAND, cv)

num_pts = 2*length(cv.qps)^6
println(num_pts)
err_tp = norm.(int_tp.-ref)/norm(ref)
println(err_tp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, cv) end

cv_s = CommonVertex6D_S(sing,SauterSchwab3D._shunnham3D(6))

int_sp = sauterschwab_parameterized(INTEGRAND, cv_s)

num_pts = num_pts = 2*length(cv_s.qps)^2
println(num_pts)
err_sp = norm.(int_sp.-ref)/norm(ref)
println(err_sp)
@time for i in 1:100 sauterschwab_parameterized(INTEGRAND, cv_s) end

#=
cv_gm = CommonVertex6D_S(SauterSchwab3D._grundmannMoeller3D(2*6-1))

int_gm = sauterschwab_parameterized(INTEGRAND, cv_gm)

num_pts = 2*length(cv_gm.qps)^2
println(num_pts)
err_gm = norm.(int_gm.-ref)/norm(ref)
println(err_gm)
=#


text:
mutable struct SoysterMaximinObjective <: AbstractObjective
    name::String
    sense::String
    R::Int64
    cores::Vector{MaximinObjectiveCore}

    SoysterMaximinObjective(R::Int64, cores::Vector{MaximinObjectiveCore}) =
        new("soyster_maximin", "max", R, cores)
    SoysterMaximinObjective(cores::Vector{MaximinObjectiveCore}) =
        new("soyster_maximin", "max", 1, cores)
    SoysterMaximinObjective() = new("soyster_maximin", "max", 1, MaximinObjectiveCore[])
end


text:
#=

By listing the first six prime numbers: 2, 3, 5, 7, 11, and 13, we can see that the 6th prime is 13.
What is the 10 001st prime number?

=#

is_prime(n) = !any(n%i == 0 for i in 2:BigInt(floor(sqrt(n))))

function nth_prime(n)
    counter = 1
    i = 2
    while true
        if is_prime(i)
            #println(counter, ", ", i)
            if counter == n
                #println()
                return i
            end
            counter += 1
        end
        i += 1
    end
end

nth_prime(6)

nth_prime(10001)

max_stars_repo_name: pedvide/project-euler

max_issues_repo_name: pedvide/project-euler

max_forks_repo_name: pedvide/project-euler


text:
"""
    A BoundaryCondition is a Kernel... it just computes differently

    A BoundaryCondition _must_ supply a bid::Array{Int64} variable
    that specifies which boundary IDs this BC is applied to
"""
abstract BoundaryCondition <: Kernel
  # bids::Array{Int64}

" A NodalBC gets applied _at_ nodes and overwrites residual / jacobian entries "
abstract NodalBC <: BoundaryCondition

" An IntegratedBC is integrated across sidesets to compute residual / jacobians "
abstract IntegratedBC <: BoundaryCondition

" Get the boundary IDs this BoundaryCondition should be applied to "
function boundaryIDs(bc::BoundaryCondition)
    return bc.bids
end

" Generic `computeJacobian()` for Float64.  Just returns 0 "
function computeJacobian(bc::BoundaryCondition, var::Variable{Float64})
    return 0.
end

"""
    Generic `computeJacobian()` for Dual

    Uses Automatic Differentiation to compute the Jacobian.

    NOTE!  This is sub-optimal if you want both the residual and the Jacobian!
    In that case use `computeResidualAndJacobian!()` instead!
"""
function computeJacobian{N,T}(bc::BoundaryCondition, var::Variable{Dual{N,T}})
    residual = computeResidual(bc)

    return partials(residual)[var.id]
end

" Compute the Residual and Jacobian using calls to `computeJacobian()` "
function computeResidualAndJacobian!(residual::Array{Float64},
                                     var_jacobians::Array{Float64},
                                     vars::Array{Variable{Float64}},
                                     bc::BoundaryCondition)
    residual[1] = computeResidual(bc)

    for v in vars
        var_jacobians[v.id] = computeJacobian(bc, v)
    end
end

" Compute the Residual and Jacobian using Automatic Differentiation "
function computeResidualAndJacobian!{N,T}(residual::Array{Dual{N,T}},
                                          var_jacobians::Array{Float64},
                                          vars::Array{Variable{Dual{N,T}}},
                                          bc::BoundaryCondition)
    residual[1] = computeResidual(bc)

    for v in vars
        var_jacobians[v.id] = partials(residual[1])[v.id]
    end
end


text:
using Revise
using VolterraGP

D = 3 # number of outputs
C = 1 # order of volterra series
P = 2 # number of paramters in base kernel

reps = 5
nmse_l = fill(0., reps)
nlpd_l = fill(0., reps)
for i = 1:reps
    train, test = generate_toy_data(50) # genrates toy data from paper
    gp = GaussianProcess(scaledEQs, D, C, P, train)
    fit!(gp, 100, ls_lr=5e-2, σ_lr=1e-2, show_like=true)
    nmse_l[i] = NMSE(test, gp)
end

print("NMSE ")
VolterraGP.summary_stats(nmse_l)
plotgp(test.X, gp)


text:
"""
    show

"""
toString(e::True)::String = "T"
toString(e::False)::String = "F"
toString(e::Var)::String = "Var($(e.name))"
toString(e::Not)::String = "Not($(toString(e)))"

function toString(e::And)::String
    innerString = join(map(e -> toString(e), e.exprs), ", ")
    "And($innerString)"
end

function toString(e::Or)::String
    innerString = join(map(e -> toString(e), e.exprs), ", ")
    "Or($innerString)"
end

show(io::IO, e::BoolExpr) = print(io, toString(e))


text:
function get_ub_lb_lanczos(pw::PWGrid, Potentials::PotentialsT, nlancz::Int64)
    #
    Ngwx = pw.gvecw.Ngwx
    V = zeros(ComplexF64,Ngwx,nlancz)
    HV = zeros(ComplexF64,Ngwx,nlancz)
    T = zeros(Float64,nlancz,nlancz)
    f = zeros(ComplexF64,Ngwx)
    s = zeros(ComplexF64,nlancz)
    h = zeros(ComplexF64,nlancz)
    #
    V[:,1] = randn(Ngwx) + im*randn(Ngwx)
    beta = norm(V[:,1])
    V[:,1] = V[:,1] ./ beta
    #
    HV[:,1] = op_H( pw, Potentials, V[:,1] )
    h[1] = real( V[:,1]' * HV[:,1] )
    #
    T[1,1] = h[1]
    # One-step of reorthogonalization
    f[:] = HV[:,1] - V[:,1]*h[1]
    s[1] = V[:,1]' * f[:]
    h[1] = h[1] + s[1]
    f[:] = f[:] - V[:,1]*s[1]
    # MAIN LOOP
    for j = 2:nlancz
        #@printf("iter lanczos = %d\n", j)
        beta = norm(f)
        T[j,j-1] = beta
        V[:,j] = f[:]/beta
        HV[:,j] = op_H( pw, Potentials, V[:,j] )
        #
        for jj = 1:j
            h[jj] = V[:,jj]' * HV[:,j]
        end
        f[:] = HV[:,j] - V[:,1:j]*h[1:j]

        # One-step of reorthogonalization
        #s = V' * f
        #h = h + s
        for jj=1:j
            s[jj] = V[:,jj]' * f[:]
            h[jj] = h[jj] + s[jj]
        end
        f[:] = f[:] - V[:,1:j]*s[1:j]
        #
        T[1:j,j] = real(h[1:j])
    end
    #
    evalsT = eigvals(T)
    #lb = evalsT[Nstates+2]
    #ub = evalsT[2*Nstates]
    lb = evalsT[Int(nlancz/2)]
    ub = norm_matrix_induced(T) + norm(f)
    #
    return lb, ub
end


text:
export Visualization, getVisus, getVisuPath, addVisu

mutable struct Visualization
  path::String
  num::Int64
  params::Dict
end

function getVisu(d::MDFDatasetStore, study::Study, exp::Experiment, reco::Reconstruction, numVisu)

  filename = joinpath(d.path, "reconstructions", getMDFStudyFolderName(study),
			       string(exp.num), string(reco.num)*".visu")

  if isfile(filename)

  #  file = h5open(filename, "r")
 #   g = file[string(numVisu)]
    params = loadParams(filename, string(numVisu))

#    close(file)
  else
    params = Dict{Symbol,Any}()
  end
  return Visualization(filename, numVisu, params)
end

function getVisus(d::MDFDatasetStore, study::Study, exp::Experiment, reco::Reconstruction)

  visus = Visualization[]

  filename = joinpath(d.path, "reconstructions", getMDFStudyFolderName(study), string(exp.num),
			      string(reco.num)*".visu")

  if isfile(filename)

    file = h5open(filename, "r")
    g = file["/"]
    for obj in g
      key = HDF5.name(obj)
      numVisu = parse(Int64, last(splitdir(key)) )
      params = loadParams(file, key) #this is the dictionarry

      push!(visus, Visualization(filename, numVisu, params))
    end
    close(file)
  end
  return visus
end

function remove(visu::Visualization)
  h5open(visu.path, "r+") do file
    if haskey(file, string(visu.num))
      delete_object(file, string(visu.num))
    end
  end
end

function save(visu::Visualization)
  if isfile(visu.path)
    file = h5open(visu.path, "r+")
  else
    file = h5open(visu.path, "w")
  end
  if haskey(file, string(visu.num))
    delete_object(file, string(visu.num))
  end
  saveParams(file, string(visu.num), visu.params)
  close(file)
end

function addVisu(d::MDFDatasetStore, study::Study, exp::Experiment, reco::Reconstruction, visuParams)

  filename = joinpath(d.path, "reconstructions", getMDFStudyFolderName(study), string(exp.num),
			      string(reco.num)*".visu" )

  visus = getVisus(d, study, exp, reco)

  if isempty(visus)
    num = 1
  else
    num = last(visus).num + 1
  end

  visu =  Visualization(filename, num, visuParams)

  save(visu)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "libcsv"
version = v"0.0.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/rgamble/libcsv.git", "6e750805f54f81a26470486d6ed79efaf7d13805")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd libcsv/
autoreconf -fiv
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
FLAGS=()
if [[ "${target}" == *-mingw* ]]; then
    FLAGS+=(LDFLAGS="-no-undefined")
fi
make -j${nproc} "${FLAGS[@]}"
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libcsv", :libcsv)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil