text:
using JuMP

"""
    MathProgNLPModel(model; kwargs...)

Construct a `MathProgNLPModel` from a JuMP `Model`.
"""
function MathProgNLPModel(jmodel :: JuMP.Model; kwargs...)
  JuMP.setsolver(jmodel, ModelReader())
  jmodel.internalModelLoaded || JuMP.build(jmodel)
  return MathProgNLPModel(jmodel.internalModel; kwargs...)
end

"""
    MathProgNLSModel(cmodel, F)

Construct a `MathProgNLSModel` from a JuMP `Model` and a vector of NLexpression.
"""
function MathProgNLSModel(cmodel :: JuMP.Model,
                          F :: Vector{JuMP.NonlinearExpression}; kwargs...)
  # Hessian-Vector products don't work:
  # https://github.com/JuliaOpt/JuMP.jl/issues/1204
  @NLobjective(cmodel, Min, 0.5 * sum(Fi^2 for Fi in F))
  JuMP.setsolver(cmodel, ModelReader())
  cmodel.internalModelLoaded || JuMP.build(cmodel)
  ev = cmodel.internalModel.eval

  Fmodel = JuMP.Model()
  @NLobjective(Fmodel, Min, 0.0)
  Fmodel.nlpdata.user_operators = cmodel.nlpdata.user_operators
  @variable(Fmodel, x[1:MathProgBase.numvar(cmodel)])
  for Fi in F
    expr = ev.subexpressions_as_julia_expressions[Fi.index]
    replace!(expr, x)
    expr = :($expr == 0)
    JuMP.addNLconstraint(Fmodel, expr)
  end
  JuMP.setsolver(Fmodel, ModelReader())
  Fmodel.internalModelLoaded || JuMP.build(Fmodel)

  return MathProgNLSModel(Fmodel.internalModel,
                          cmodel.internalModel; kwargs...)
end

function replace!(ex, x)
  if isa(ex, Expr)
    for (i, arg) in enumerate(ex.args)
      if isa(arg, Expr)
        if arg.head == :ref && arg.args[1] == :x
          ex.args[i] = x[arg.args[2]]
        else
          replace!(arg, x)
        end
      end
    end
  end
end


text:
using SaddleSearch
using Base.Test
using SaddleSearch: numE, numdE, res_trans, res_rot, maxres
include("testsets.jl"); using TestSets

function heading1(str)
   dashes = "-"^length(str)
   print_with_color(:magenta, dashes); println()
   print_with_color(:magenta, str, bold=true); println()
   print_with_color(:magenta, dashes);println()
end

heading2(str) = print_with_color(:green, str*"\n", bold=true)

# change this to 2 to get full information or to 0 to get virtually no information
verbose=1

# test files are test_*** where *** is the string in this list
tests = [
   "staticvsbb",
   "lsdimer",
   "vacancy",
   "string",
]

println("≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡")
println(" Starting `SaddleSearch` tests")
println("≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡≡")
@testset "SaddleSearch" begin
   for test in tests
      include("test_$(test).jl")
   end
end

ext: jl


text:
@testset "Monomial selection for certificate" begin
    @polyvar x y z
    @testset "Multipartite error not disjoint: $parts" for parts in [([x], [x]),
                                                                     ([x], [x, y]),
                                                                     ([x], [y], [x, y])]
        err = ArgumentError("Parts are not disjoint in multipartite Newton polytope estimation: $parts")
        @test_throws err SumOfSquares.Certificate.monomials_half_newton_polytope([x*y, y^2], parts)
    end
    @testset "Unipartite" begin
        @test SumOfSquares.Certificate.monomials_half_newton_polytope([x*y, y^2], tuple()) == [y]
        @test isempty(SumOfSquares.Certificate.monomials_half_newton_polytope([x, y], tuple()))
        @test SumOfSquares.Certificate.monomials_half_newton_polytope([x^2, y^2], tuple()) == [x, y]
        @test SumOfSquares.Certificate.monomials_half_newton_polytope([x^2, y^2], ([x, y],)) == [x, y]
        @test SumOfSquares.Certificate.monomials_half_newton_polytope([x^2, y^2], ([y, x],)) == [x, y]
    end
    @testset "Multipartite" begin
        # In the part [y, z], the degree is between 0 and 2
        X = [x^4, x^2*y^2, x^2*z^2, x^2*y*z, y*z]
        @test SumOfSquares.Certificate.monomials_half_newton_polytope(X, tuple()) == [x^2, x*y, x*z, y*z, x, y, z]
        function full_test(X, Y, part1, part2)
            @test SumOfSquares.Certificate.monomials_half_newton_polytope(X, (part1,)) == Y
            @test SumOfSquares.Certificate.monomials_half_newton_polytope(X, (part2,)) == Y
            a = SumOfSquares.Certificate.monomials_half_newton_polytope(X, (part2,))
            @test SumOfSquares.Certificate.monomials_half_newton_polytope(X, (part1, part2)) == Y
            @test SumOfSquares.Certificate.monomials_half_newton_polytope(X, (part2, part1)) == Y
        end
        full_test(X, monovec([x^2, x*y, x*z,      x, y, z]), [x], [y, z])
        full_test(X, monovec([x^2, x*y, x*z, y*z, x,    z]), [y], [x, z])
        full_test(X, monovec([x^2, x*y, x*z, y*z, x, y   ]), [z], [x, y])
        # FIXME: With recursive merging, it should give [x^2, x*y, x*z, x]
        @test SumOfSquares.Certificate.monomials_half_newton_polytope([x^4, x^2*y^2, x^2*z^2, x^2*y*z, y*z], ([x], [y], [z])) == [x^2, x*y, x*z, y*z, x, y, z]
    end
end

@testset "Random SOS should be SOS with $(factory.constructor)" for factory in sdp_factories
    @polyvar x y
    x = [1, x, y, x^2, y^2, x*y]
    @test_throws ArgumentError randsos(x, monotype=:Unknown)
    for i in 1:10
        for monotype in [:Classic, :Gram]
            p = randsos(x, monotype=monotype)

            m = SOSModel(factory)

            @constraint m p >= 0

            JuMP.optimize!(m)

            @test JuMP.primal_status(m) == MOI.FEASIBLE_POINT
        end
    end
end


text:
# Parameters controlling how a plot appears

type Theme
    # If the color aesthetic is not mapped to anything, this is the color that
    # is used.
    default_color::ColorOrNothing

    # Default size when the size aesthetic is not mapped.
    default_point_size::Measure

    # Width of lines in the line geometry.
    line_width::Measure

    # Background color of the plot.
    panel_fill::ColorOrNothing

    # Border color of the plot panel.
    panel_stroke::ColorOrNothing

    # Grid line color.
    grid_color::ColorOrNothing

    # Grid lines for focused item.
    grid_color_focused::ColorOrNothing

    # Width of grid lines
    grid_line_width::Measure

    # Font name, size, and color used for tick labels, entries in keys, etc.
    minor_label_font::String
    minor_label_font_size::Measure
    minor_label_color::ColorOrNothing

    # Font name, size and color used for axis labels, key title, etc.
    major_label_font::String
    major_label_font_size::Measure
    major_label_color::ColorOrNothing

    # Font name, size and color used for labels on plot elements.
    point_label_font::String
    point_label_font_size::Measure
    point_label_color::ColorOrNothing

    # Spacing between bars for Geom.bar.
    bar_spacing::Maybe(Measure)

    # Spacing between boxplots in Geom.boxplot.
    boxplot_spacing::Maybe(Measure)

    # Points, etc, are highlighted by stroking in slightly different color. This
    # is the stroke width.
    highlight_width::Maybe(Measure)

    # A function mapping fill color to stoke color for highlights.
    highlight_color::Maybe(Function)

    # A function mapping base fill color to the color of the median marker in a
    # boxplot.
    middle_color::Maybe(Function)

    # Number of annealing iterations.
    label_placement_iterations::Maybe(Int)

    # Penalty for a label not being contained within the plot frame.
    label_out_of_bounds_penalty::Maybe(Float64)

    # Penalty for making a label hidden to avoid overlaps.
    label_hidden_penalty::Maybe(Float64)

    # Probability of proposing a visibility flip during label layout.
    label_visibility_flip_pr::Maybe(Float64)
end

# Choose highlight color by darkening the fill color
function default_highlight_color(fill_color::ColorValue)
    fill_color = convert(LCHab, fill_color)
    c = LCHab(fill_color.l, fill_color.c, fill_color.h)
    LCHab(c.l - 15, c.c, c.h)
end

# Choose a middle color by darkening the fill color
function default_middle_color(fill_color::ColorValue)
    fill_color = convert(LCHab, fill_color)
    LCHab(fill_color.l + 40, fill_color.c, fill_color.h)
end

const default_font_desc = "PT Sans,Helvetica Neue,Helvetica,sans"

const default_theme =
    Theme(color("steel blue"),      # default_color
          0.6mm,                    # default_point_size
          0.2mm,                    # line_width
          color("#fcfcfc"),         # panel_fill
          color("#f4f4f8"),         # panel_stroke
          color("#f4f4f8"),         # grid_color
          color("#f0f0f0"),         # grid_color_focused
          0.2mm,                    # grid line width
          default_font_desc,        # minor_label_font
          9pt,                      # minor_label_font_size
          color("#4c404b"),         # minor_label_color
          default_font_desc,        # major_label_font
          11pt,                     # major_label_font_size
          color("#362a35"),         # major_label_color
          default_font_desc,        # point_label_font
          8pt,                      # point_label_font_size
          color("#4c404b"),         # point_label_color
          0.0mm,                    # bar_spacing
          1mm,                      # boxplot_spacing
          0.3mm,                    # highlight_width
          default_highlight_color,  # highlight_color
          default_middle_color,     # middle_color
          1000,                     # label_placement_iterations
          10.0,                     # label_out_of_bounds_penalty
          0.5,                      # label_hidden_penalty
          0.2)                      # label_visibility_flip_pr


text:
import ..Types: move!
import ..Utils: acceleration

@inline function _move_half_step!(r::Vector3s, v::Vector3s, a::Vector3s, dt::Float64)
    N = length(r)
    v .+= a * dt / 2
    r .+= v * dt
    r, v
end

@inline function _move_full_step!(r::Vector3s, v::Vector3s, a::Vector3s, dt::Float64)
    N = length(r)
    v .+= a * dt / 2
    r, v
end

function move!(integrator::VerletIntegrator)
    r = integrator.positions
    v = integrator.velocities
    dt = integrator.timestep
    m = integrator.mass_function
    f = integrator.force_function
    fv = integrator.forces

    a = acceleration(fv, m)
    _move_half_step!(r, v, a, dt)

    fv .= f(r)

    a = acceleration(fv, m)
    _move_full_step!(r, v, a, dt)

    integrator, dt
end


text:
@testset "save_tetgen.jl" begin

end


text:
export TraceNorm_dual

function TraceNorm_dual(x::Tx,weights::Tw, params::Dict{String,Any}) where {ETx<:Number, Tx<:AbstractVector{ETx}, ETw<:Number, Tw<:AbstractArray{ETw}}

numr::Int = params["numr"]
numc::Int= params["numc"]
E = reshape(x, numr, numc)

#DEVNOTE# -Performance: Very expensive line
tmp = svds(E; nsv = 1, ritzvec = false)[1]

d = tmp.S[1]

return d

end


text:
n = 10 # rand(1:10)
@test matrixdepot("binomial", n) == matrixdepot("binomial", Float64, n)

A = matrixdepot("binomial", n)

@test A*A ≈ 2^(n-1)*Matrix(1.0I, n, n)
println("'binomial' passed test...")


text:
#=
############ Streams ############
This file contains several random stream types. Any new streams type should implement the following interface:

create_rng(::Stream) creates a (mutable) rng of an appropriate type to be emitted from the stream
fill_random_streams!(::Stream, rand_max::Int) is called at the beginning of action selection
set_rng_state!(rng::AbstractRNG, ::Stream, scenario::Int, depth::Int) is called at the beginning of each simulation step to set the rng state to a value appropriate for the given scenario and depth
###############################
=#

# RandomStreams
# This type encapsulates the streams of random numbers used in state
# transitions during simulations. It also provides random-number seeds
# for different components of the system.

mutable struct RandomStreams
    num_streams::Int64
    len_streams::Int64
    streams::Array{Float64,2}     # each particle is associated with a single stream of numbers
    seed::UInt32

  # default constructor
  function RandomStreams(num_streams::Int64,
                len_streams::Int64,
                seed::UInt32)
          this = new()

          this.num_streams = num_streams
          this.len_streams = len_streams
          this.streams = Array{Float64}(num_streams, len_streams)
          this.seed = seed

          return this
    end
end

get_stream_seed(streams::RandomStreams, streamId::UInt32) =
    streams.seed ⊻ streamId # bitwise XOR

get_world_seed(streams::RandomStreams) =
    streams.seed ⊻ streams.num_streams

get_model_seed(streams::RandomStreams) =
    streams.seed ⊻ (streams.num_streams + 2)

function fill_random_streams!(empty_streams::RandomStreams, rand_max::Int64)
    # Populate random streams
    if is_linux()
        ccall((:srand, "libc"), Void, (Cuint,), 1)
        for i in 1:empty_streams.num_streams
            seed = Cuint[get_stream_seed(empty_streams, convert(UInt32, i-1))]
            ccall( (:rand_r, "libc"), Int, (Ptr{Cuint},), seed)
            for j = 1:empty_streams.len_streams
                empty_streams.streams[i,j] = ccall((:rand_r, "libc"), Int, (Ptr{Cuint},), seed) / rand_max
            end
        end
    else #Windows, etc
        for i in 1:empty_streams.num_streams
            seed = get_stream_seed(empty_streams, convert(UInt32, i-1))
            srand(seed)
            empty_streams.streams[i,:] = Base.rand(convert(Int64, empty_streams.len_streams))
        end
    end
end

function set_rng_state!(rng::DESPOTRandomNumber, rs::RandomStreams, scenario::Int, depth::Int)
    rng.number = rs.streams[scenario, depth]
end

create_rng(::RandomStreams) = DESPOTRandomNumber(-1)

### MersenneTwister Streams ###

mutable struct MersenneStreamArray
    rng::AbstractRNG
    streams::Vector{Vector{MersenneTwister}}
end
MersenneStreamArray(rng::AbstractRNG) = MersenneStreamArray(rng, Vector{MersenneTwister}[])

function set_rng_state!(rng::MersenneTwister, sa::MersenneStreamArray, scenario::Int, depth::Int)
    while length(sa.streams) < scenario
        push!(sa.streams, MersenneTwister[])
    end
    stream = sa.streams[scenario]
    while length(stream) < depth
        push!(stream, MersenneTwister(rand(sa.rng, UInt32))) # maybe just advance a large number of steps and copy, rather than seed?
    end
    copy!(rng, stream[depth])
end

fill_random_streams!(streams::MersenneStreamArray, rand_max::Int64) = nothing

create_rng(::MersenneStreamArray) = copy(Base.GLOBAL_RNG)


text:
module CapacityExpansionData
    include(joinpath(@__DIR__,"utils","data_path.jl"))
end


text:
bigfib(n) = ((BigInt[1 1; 1 0])^n)[2,1]
bigfib(20)


text:
# Common Utils
# ============
#
# Common utilities shared among several algorithms (internal use only).
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# k: gap length
function affinegap_score(k, gap_open_penalty, gap_extend_penalty)
    return -(gap_open_penalty + gap_extend_penalty * k)
end

# Trace type
# ----------

# trace type for pairwise alignment
typealias Trace UInt8

# trace bitmap
const TRACE_NONE   = 0b00000
const TRACE_MATCH  = 0b00001
const TRACE_DELETE = 0b00010
const TRACE_INSERT = 0b00100
const TRACE_EXTDEL = 0b01000
const TRACE_EXTINS = 0b10000

# Traceback
# ---------

macro start_traceback()
    esc(quote
        anchor_point = (i, j)
        op = OP_INVALID
    end)
end

macro finish_traceback()
    quote
        push!(anchors, AlignmentAnchor(anchor_point, op))
        push!(anchors, AlignmentAnchor((i, j), OP_START))
        reverse!(anchors)
        pop!(anchors)  # remove OP_INVALID
    end
end

macro anchor(ex)
    esc(quote
        if op != $ex
            push!(anchors, AlignmentAnchor(anchor_point, op))
            op = $ex
            anchor_point = (i, j)
        end
        if ismatchop(op)
            i -= 1
            j -= 1
        elseif isinsertop(op)
            i -= 1
        elseif isdeleteop(op)
            j -= 1
        else
            @assert false
        end
    end)
end


text:
using GalacticOptimJL, GalacticOptimJL.Optim, GalacticOptim, ForwardDiff, Zygote, Random, ModelingToolkit
using Test

@testset "GalacticOptimJL.jl" begin
    rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
    x0 = zeros(2)
    _p = [1.0, 100.0]
    l1 = rosenbrock(x0, _p)
    f = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff())
    prob = OptimizationProblem(f, x0, _p)
    Random.seed!(1234)
    sol = solve(prob, SimulatedAnnealing())
    @test 10 * sol.minimum < l1

    Random.seed!(1234)
    prob = OptimizationProblem(f, x0, _p, lb=[-1.0, -1.0], ub=[0.8, 0.8])
    sol = solve(prob, SAMIN())
    @test 10 * sol.minimum < l1

    prob = OptimizationProblem(rosenbrock, x0, _p)
    sol = solve(prob, Optim.NelderMead(;initial_simplex=Optim.AffineSimplexer(;a = 0.025, b = 0.5)))
    @test 10 * sol.minimum < l1

    cons = (x, p) -> [x[1]^2 + x[2]^2]
    optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff(); cons=cons)
    optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoModelingToolkit(); cons=cons)

    prob = OptimizationProblem(optprob, x0, _p)

    sol = solve(prob, Optim.BFGS())
    @test 10 * sol.minimum < l1

    sol = solve(prob, Optim.Newton())
    @test 10 * sol.minimum < l1

    sol = solve(prob, Optim.KrylovTrustRegion())
    @test 10 * sol.minimum < l1

    prob = OptimizationProblem(optprob, x0, _p, lcons=[-Inf], ucons=[Inf])
    sol = solve(prob, IPNewton())
    @test 10 * sol.minimum < l1

    prob = OptimizationProblem(optprob, x0, _p, lcons=[-5.0], ucons=[10.0])
    sol = solve(prob, IPNewton())
    @test 10 * sol.minimum < l1

    prob = OptimizationProblem(optprob, x0, _p, lcons=[-Inf], ucons=[Inf], lb=[-500.0, -500.0], ub=[50.0, 50.0])
    sol = solve(prob, IPNewton())
    @test sol.minimum < l1

    function con2_c(x, p)
        [x[1]^2 + x[2]^2, x[2] * sin(x[1]) - x[1]]
    end

    optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff(); cons=con2_c)
    prob = OptimizationProblem(optprob, x0, _p, lcons=[-Inf, -Inf], ucons=[Inf, Inf])
    sol = solve(prob, IPNewton())
    @test 10 * sol.minimum < l1

    cons_circ = (x, p) -> [x[1]^2 + x[2]^2]
    optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff(); cons=cons_circ)
    prob = OptimizationProblem(optprob, x0, _p, lcons=[-Inf], ucons=[0.25^2])
    sol = solve(prob, IPNewton())
    @test sqrt(cons(sol.u, nothing)[1]) ≈ 0.25 rtol = 1e-6

    optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoZygote())

    prob = OptimizationProblem(optprob, x0, _p, lb=[-1.0, -1.0], ub=[0.8, 0.8])
    sol = solve(prob, Optim.Fminbox())
    @test 10 * sol.minimum < l1

    Random.seed!(1234)
    prob = OptimizationProblem(optprob, x0, _p, lb=[-1.0, -1.0], ub=[0.8, 0.8])
    sol = solve(prob, Optim.SAMIN())
    @test 10 * sol.minimum < l1

    optprob = OptimizationFunction((x, p) -> -rosenbrock(x, p), GalacticOptim.AutoZygote())
    prob = OptimizationProblem(optprob, x0, _p; sense=GalacticOptim.MaxSense)

    sol = solve(prob, NelderMead())
    @test 10 * sol.minimum < l1

    sol = solve(prob, BFGS())
    @test 10 * sol.minimum < l1

    function g!(G, x)
        G[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
        G[2] = 200.0 * (x[2] - x[1]^2)
    end
    optprob = OptimizationFunction((x, p) -> -rosenbrock(x, p), GalacticOptim.AutoZygote(), grad=g!)
    prob = OptimizationProblem(optprob, x0, _p; sense=GalacticOptim.MaxSense)
    sol = solve(prob, BFGS())
    @test 10 * sol.minimum < l1

    optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoModelingToolkit())
    prob = OptimizationProblem(optprob, x0, _p)
    sol = solve(prob, Optim.BFGS())
    @test 10 * sol.minimum < l1
end


text:
# undefined fixture
let
  @test_throws UndefVarError @pytest function(g) end
end

# unnamed fixture
let
  @test_throws ArgumentError eval(macroexpand(:(@fixture function() end)))
end

# incorrect middle arguments to fixture
let
  @test_throws ArgumentError eval(macroexpand(:(@fixture wrooooong function f() end)))
end


text:
@testset "TaylorExpansion" begin
    s0 = [0.2, 0.4, 1.1]
    x0 = [1.2, 0.9]

    N = 1000
    points = rand(N, 3)

    X_s = rand(2, 3)
    X_ss = rand(2, 3, 3)
    X_sss = rand(2, 3, 3, 3)
    dr1 = TaylorExpansion(s0, x0, X_s)
    dr2 = TaylorExpansion(s0, x0, X_s, X_ss)
    dr3 = TaylorExpansion(s0, x0, X_s, X_ss, X_sss)

    out1 = dr1(points)
    out2 = dr2(points)
    out3 = dr3(points)

    out1_1d = dr1(vec(points[1, :]))
    out2_1d = dr2(vec(points[1, :]))
    out3_1d = dr3(vec(points[1, :]))

    @test maximum(abs, out1_1d - vec(out1[1, :])) == 0.0
    @test maximum(abs, out2_1d - vec(out2[1, :])) == 0.0
    @test maximum(abs, out3_1d - vec(out3[1, :])) == 0.0

    ds = points .- s0'
    verif1 = x0 .+ X_s*ds'
    @test maximum(abs, out1 - verif1') < 1e-12

    # TODO: write verification tests for 2d and 3d
end


text:
using KaHyPar
using SparseArrays
using DataStructures

export KaHyPar_contraction_plan, KaHyPar_contraction!

# *************************************************************************** #
#                             KaHyPar functions
# *************************************************************************** #

"""
    function KaHyPar_contraction_plan(network::TensorNetworkCircuit,
                                      k::Int, ϵ::Real, V::Int)

This function uses KaHyPar to partition a TensorNetworkCircuit into k parts of
roughly equal size (how roughly is determined by ϵ) and then finds a contraction
plan for the partitions. This function is then recursively called on the
partitions to compute a complete contraction plan for the whole Tensor Network.
If the number of nodes in the network is lower than the threshold V then netcon
is used to find the optimal contraction plan.

network - The tensor network to build a contraction plan for.

k - The number of partitions to divide the network into at each recursive step.

ϵ - Imbalance parameter for KaHyPar

V - If the number of nodes in the given network is below this threshold then
    netcon is used to find the optimal contraction plan for it.
"""
function KaHyPar_contraction_plan(network::TensorNetworkCircuit,
                                  k::Int=2, ϵ::Real=0.1, V::Int=20,
                                  index_dims::Union{Dict{Symbol, Int}, Nothing}=nothing)

    # TODO: This won't be needed when the bgreedy method is added as an option
    # for finding contraction plans
    if k > 36
        error("The number of partitions should be less than 36 for netcon")
    end

    # Base case: If the number of nodes is less than the threshold V
    # use netcon to get the optimal contraction plan.
    if length(network.nodes) < V
        return netcon(network, :tree)
    end

    # If index_dims dictionary not given then create one. This is used to create
    # the network of partitions with the correct dimensions for edges connecting
    # different partitions (needed by netcon to find the contraction plan).
    if index_dims === nothing
        index_dims = Dict{Symbol, Int}()
        for node in values(network.nodes)
            for (index, dim) in zip(node.indices, node.dims)
                index_dims[index] = dim
            end
        end
    end
    # TODO: something to think about/consider:
    # The above could be made more efficient or avoided if dims were matched
    # with edge indices in the data structure somehow. eg if dims were stored
    # in a dictionary, mapping edges to dims, inside the node struct instead of
    # an array or if dims were stored in the edge structs.

    # Use KaHyPar to find k partitions of the tensor network.
    A = incidence_matrix(network)
    h = KaHyPar.HyperGraph(A)
    partitions = KaHyPar.partition(h, k, configuration = :edge_cut)

    # Find a contraction tree for the partitions.
    partitioned_net = network_of_partitions(network, partitions, index_dims)
    contraction_plan = netcon(partitioned_net, :tree)

    # Create a dictionary mapping partitions to subnetworks of the Tensor
    # Network. These subnetworks are to be partitioned in the next recursive
    # step.
    subnets = Dict{Symbol, TensorNetworkCircuit}()
    nodes = [keys(network.nodes)...]
    for i = 0:maximum(partitions)
        partition_i = findall(x -> x==i, partitions)
        label = Symbol("partition_$i")
        subnets[label] = TensorNetworkCircuit(network, nodes[partition_i])
    end

    # Recursively call KaHyPar_contraction_plan on the next layer of the
    # contraction tree.
    contraction_tree = Array{Any, 1}(undef, 2)
    contraction_tree[1] = KaHyPar_contraction_plan(contraction_plan[1],
                                                   subnets, k, ϵ, V, index_dims)
    contraction_tree[2] = KaHyPar_contraction_plan(contraction_plan[2],
                                                   subnets, k, ϵ, V, index_dims)

    contraction_tree
end

"""
    function KaHyPar_contraction_plan(contraction_step::Array{Any, 1},
                                      subnets::Dict{Symbol, TensorNetworkCircuit},
                                      k::Int, ϵ::Real, V::Int)

While using KaHyPar_contraction_plan to find a contraction plan for a network,
if an element of the next layer of a contraction tree happens to be a
contraction of two other objects, this function recursively calls
KaHyPar_contraction_plan on both of those elements.
"""
function KaHyPar_contraction_plan(contraction_step::Array{Any, 1},
                                  subnets::Dict{Symbol, TensorNetworkCircuit},
                                  k::Int, ϵ::Real, V::Int,
                                  index_dims::Dict{Symbol, Int})
    contraction_plan[1] = KaHyPar_contraction_plan(contraction_plan[1],
                                                   subnets, k, ϵ, V, index_dims)
    contraction_plan[2] = KaHyPar_contraction_plan(contraction_plan[2],
                                                   subnets, k, ϵ, V, index_dims)

    contraction_plan
end

"""
    function KaHyPar_contraction_plan(partition::Symbol,
                                      subnets::Dict{Symbol, TensorNetworkCircuit},
                                      k::Int, ϵ::Real, V::Int)

While using KaHyPar_contraction_plan to find a contraction plan for a network,
if an element of the next layer of a contraction tree is a symbol representing
a partition of the original tneor network, this function calls
KaHyPar_contraction_plan on the corresponding subnetwork.
"""
function KaHyPar_contraction_plan(partition::Symbol,
                                  subnets::Dict{Symbol, TensorNetworkCircuit},
                                  k::Int, ϵ::Real, V::Int,
                                  index_dims::Dict{Symbol, Int})
    KaHyPar_contraction_plan(subnets[partition], k, ϵ, V, index_dims)
end

"""
    function incidence_matrix(tn::TensorNetworkCircuit)

Create an incidence matrix for the given tensor network.
"""
function incidence_matrix(tn::TensorNetworkCircuit)
    vertices, edges = [], []

    for (node_label, node) in tn.nodes
        node_number = parse(Int, String(node_label)[6:end])
        for edge_label in node.indices
            edge_number = parse(Int, String(edge_label)[7:end])

            append!(vertices, node_number)
            append!(edges, edge_number)
        end
    end

    weigths = ones(Int, length(vertices))
    sparse(vertices, edges, weigths)
end

# TODO: The following function should probably be in layer3. However, a general
# 'subnetwork function' could return a network independent of the parent network
# (i.e. so changes to the subnetwork are not reflected in the parent network and
# should probably have its own backend etc.) which would be inefficient for this
# particular use case. Shoud look at other possible use cases before deciding
# on the best implementation.
"""
    function TensorNetworkCircuit(tn::TensorNetworkCircuit,
                                  nodes::Array{Symbol, 1})

Outer constructor to create an instance of TensorNetworkCircuit representing
a sub network of a given TensorNetworkCircuit. The given array of nodes are
copied from the original network to create the subnetwork.
"""
function TensorNetworkCircuit(tn::TensorNetworkCircuit, nodes::Array{Symbol, 1})

    # Copy the given nodes from the given tensor network.
    subnet_nodes = OrderedDict{Symbol, Node}(node_label => tn.nodes[node_label]
                   for node_label in nodes)

    # Copy the edges connected to the copied nodes.
    subnet_edges = OrderedDict{Symbol, Edge}()
    for node in values(subnet_nodes)
        for index in node.indices
            index in keys(subnet_edges) && continue
            edge = deepcopy(tn.edges[index])
            if !(edge.dst in nodes); edge.dst = nothing end
            if !(edge.src in nodes); edge.src = nothing end
            subnet_edges[index] = edge
        end
    end

    # Return a TensorNetworkCircuit representing the subnetwork.
    TensorNetworkCircuit(tn.backend, tn.number_qubits,
                         tn.input_qubits, tn.output_qubits,
                         subnet_nodes, subnet_edges,
                         tn.qubit_ordering, tn.counters, tn.node_layers)
end

"""
    function network_of_partitions(tn::TensorNetworkCircuit,
                                   partitioning::Array{Int, 1})

Given a Tensor Network and a partitioning of it, this function returns
a tensor network whose nodes represent the partitions and whose edges
correspond to the edges connecting partitions.
"""
function network_of_partitions(network::TensorNetworkCircuit,
                               partitions::Array{Int, 1},
                               index_dims::Dict{Symbol, Int})

    # Get the number of partitions and create a dictionary mapping nodes
    # to the partition to which they belong. (this construction depends
    # on tn.nodes being an ordered dictionary.)
    num_nodes = maximum(partitions)
    node_partition = Dict{Symbol, Int}(keys(network.nodes) .=> partitions)

    # Create a node for each partition.
    nodes = OrderedDict{Symbol, Node}()
    for i = 0:num_nodes
        label = Symbol("partition_$i")
        nodes[label] = Node(label)
    end

    # Create an edge for any edge of the original tensor network that connects
    # two partitions, or is an open edge connected to a particular partition.
    edges = OrderedDict{Symbol, Edge}()
    for (index, edge) in network.edges
        src, dst = edge.src, edge.dst
        src === nothing || (src = node_partition[src])
        dst === nothing || (dst = node_partition[dst])
        src == dst && continue # Skip this edge if it is a loop.

        if !(src === nothing)
            src = Symbol("partition_$src")
            append!(nodes[src].indices, [index])
            append!(nodes[src].dims, [index_dims[index]])
        end

        if !(dst === nothing)
            dst = Symbol("partition_$dst")
            append!(nodes[dst].indices, [index])
            append!(nodes[dst].dims, [index_dims[index]])
        end

        edges[index] = Edge(src, dst, edge.qubit, edge.virtual)
    end

    TensorNetworkCircuit(network.backend, 0, [], [], nodes, edges,
                         [], network.counters, network.node_layers)
end

"""
    function convert_tree_to_plan(tree::Union{Array{<:Any, 1}, Int},
                                  last_node_label::Int,
                                  tensor_labels::Array{Symbol, 1},
                                  plan::Array{Array{Symbol, 1}, 1})

This function converts a contraction tree, found by the netcon method, to a
contraction plan (array of node-label pairs) for PicoQuant contract functions.

tree - The contraction tree found by netcon.

last_node_label - Initially, this should be the number used to create the last
                  node added to the network being contracted, it is subsequantly
                  incremented to account for nodes created during contraction.

converted_plan - Initially, this argument should be an empty array. It will hold
                 the sequence of pairs when the function completes.
"""
function convert_tree_to_plan(tree::Union{Array{<:Any, 1}, Symbol},
                              last_node_label::Int,
                              plan::Array{Array{Symbol, 1}, 1}
                                    =Array{Array{Symbol, 1}, 1}())

    # Recursively convert the plan into a sequence of pairs.
    if typeof(tree) == Symbol
        # Base case: If the plan is a Symbol corresponding to a single tensor,
        # return that Symbol as the label for this tensor and don't change the
        # number of nodes or add anything to the converted plan.
        return tree, last_node_label, plan

    else
        # First, get the label for the tensor created by contracting the tree
        # in plan[1] and append its contraction plan to converted plan. Then
        # do the same for plan[2]. Finally, append the contraction of these
        # trees to the converted plan.
        a, last_node_label = convert_tree_to_plan(tree[1], last_node_label,
                                                  plan)
        b, last_node_label = convert_tree_to_plan(tree[2], last_node_label,
                                                  plan)
        append!(plan, [[a, b]])
        last_node_label = last_node_label + 1
        return Symbol("node_$last_node_label"), last_node_label, plan
    end
end

"""
    function KaHyPar_contraction!(network::TensorNetworkCircuit,
                                  k::Int, ϵ::Real, V::Int,)

Function to contract a tensor network according to a contraction plan found
using the KaHyPar method.
"""
function KaHyPar_contraction!(network::TensorNetworkCircuit,
                              k::Int=2, ϵ::Real=1, V::Int=7,
                              output_shape::Union{String, Array{Int, 1}}="")
    contraction_tree = KaHyPar_contraction_plan(network, k, ϵ, V)
    contraction_plan = convert_tree_to_plan(contraction_tree,
                                            network.counters["node"])

    contract_network!(network, contraction_plan[3], output_shape)
end


text:
import JSON
import Glob
import GZip

# include("textmodel.jl")
include("io.jl")

separator='\t'
emptydict = Dict()

function normspaces(text)
    replace(replace(text, "\n", " "), "\t", " ")
end

function formattext(tweet)
    key = tweet["key"]
    country = "UNK"
    place = get(tweet, "place", emptydict)
    country = get(place, "country_code", "UNK")
    text = normspaces(tweet["text"])
    date = replace(tweet["created_at"], "+", "")
    date = string(Dates.DateTime(date, "e u d H:M:S s y"))
    coordinates = get(tweet, "coordinates", emptydict)
    string(key, separator, date, separator, country, separator, JSON.json(coordinates), separator, text, '\n')
end

if length(ARGS) > 0
    for filename in ARGS
        itertweets(filename) do tweet
            write(STDOUT, tweet |> formattext)
        end
    end
else
    while !eof(STDIN)
        line = readline(STDIN)
        write(STDOUT, line |> parsetweet |> formattext)
    end
end

max_stars_repo_name: INGEOTEC/PatternAnalysisTools

max_issues_repo_name: INGEOTEC/PatternAnalysisTools

max_forks_repo_name: INGEOTEC/PatternAnalysisTools


text:
"A funtion to create an `OMPVariables` object from a `BEDFile` object `x`, a `SharedVector` `y`, and a desired sparsity level `k`."
function OMPVariables{T <: Float}(
    x :: BEDFile{T},
    y :: SharedVector{T},
    k :: Int;
    pids :: Vector{Int} = procs(x)
)
    n,p  = size(x)
    r    = SharedArray(T, n, pids=pids) :: SharedVector{T}
    copy!(r, y)
    idxs = ones(Int, k)
    a    = zeros(T, n)
    b    = SharedArray(T, p, pids=pids) :: SharedVector{T}
    dots = SharedArray(T, p, pids=pids) :: SharedVector{T}

    OMPVariables{T}(r, idxs, a, b, dots)
end

"""
    omp!(w::OMPVariables, x::BEDFile, y, k)

Greedily computes the vector `b` of sparsity level `k` that minimizes `norm(y - x*b), where `x` is a `BEDFile` object`.
This is the kernel of the function `omp`.
"""
function omp!{T <: Float}(
    w    :: OMPVariables{T},
    x    :: BEDFile{T},
    y    :: SharedVector{T},
    k    :: Int;
    m    :: Vector{Int} = ones(Int, length(y)), # bitmask of 0s and 1s
    pids :: Vector{Int} = procs(x)
)
    # dot_p = x' * r
    # need this for basis pursuit, since OMP selects dot product of largest magnitude
    At_mul_B!(w.dots, x, w.r, m, pids=pids)

    # compute index of next index to add to support
    λ = indmax_abs(w.dots)

    # expand active set indices Λ = [Λ i]
    w.idxs[k] = λ
    idx       = w.idxs[1:k]

    # get subset of x corresponding to active set
    temp = zeros(T, (n,k))
    decompress_genotypes!(temp, x, idx)

    # z = argmin_{b} ( norm(y - temp*b) )
    z, = lsqr(temp, y)

    # r = y - temp*z
    #BLAS.gemv!('N', -one(T), temp, z, zero(T), w.r)
    #BLAS.axpy!(one(T), y, w.r)
    A_mul_B!(w.r, temp, z)
    BLAS.axpy!(-one(T), y, w.r)

    # save current model to b
    # no need to erase w.b beforehand since OMP builds models stepwise,
    # and previous model (sparsity level k-1) is subset of current model of size k
    w.b[idx] = z

    return nothing
end

"""
    omp(x::BEDFile, y, k) -> SparseMatrixCSC

Perform *o*rthogonal *m*atching *p*ursuit using a BEDFile object `x`, a response vector `y`, and a sparsity level `k`.
This function will compute in greedy fashion all sparsity levels `1, 2, ..., k`.
The result is a matrix of type `SparseMatrixCSC` containing the sparse vector that minimizes `sumsq(y - x*b)` for each model size.

Arguments:
- `x` is the design matrix (dictionary)
- `y` is the response vector (signal)
- `k` is the desired sparsity level

Output:
- `B` is the sparse matrix whose columns provide estimates of `y` with sparisty level `1, 2, ..., k`.
"""
function omp{T <: Float}(
    x :: BEDFile{T},
    y :: DenseVector{T},
    k :: Int;
    m :: Vector{Int} = ones(Int, length(y)),
    quiet :: Bool = true,
)
    # size of problem?
    n,p = size(x)

    # initialize all temporary arrays
    w = OMPVariables(x, y, k)

    # initialize sparse matrix of models
    # will fill this iteratively as we grow our beta
    # result is a "path" of betas
    B = spzeros(T, p, k)

    # compute models 1, 2, ..., k
    for i = 1:k

        # this computes model size i
        omp!(w, x, y, i)

        # output progress if desired
        quiet || @printf("Sparsity level: %d, norm(residuals): %f\n", i, norm(w.r));

        # save model for sparsity level i
        B[:,i] = sparsevec(w.b)
    end

    # return path
    return B
end

### NOT DONE: CROSSVALIDATION ROUTINES


text: include("../server/server.jl").serve(6666)


text:
using Revise # reduce recompilation time
using Plots
# using Documenter
using LinearAlgebra
using SparseArrays
using BenchmarkTools
using UnPack
using ToeplitzMatrices

push!(LOAD_PATH, "./src")
using CommonUtils
using Basis1D
using Basis2DTri
using UniformTriMesh

push!(LOAD_PATH, "./examples/EntropyStableEuler")
using EntropyStableEuler

using SetupDG

S_N(x) = @. sin(pi*x/h)/(2*pi/h)/tan(x/2)
"""
Vandermonde matrix of sinc basis functions determined by h,
evaluated at r
"""
function vandermonde_Sinc(h,r)
    N = convert(Int, 2*pi/h)
    V = zeros(length(r),N)
    for n = 1:N
        V[:,n] = S_N(r.-n*h)
    end
    V[1,1] = 1
    V[end,end] = 1
    return V
end

# For testing TODO: redundant
γ = 1.4
function Unorm(U)
    unorm = zeros(size(U[1]))
    for u in U
        @. unorm += u^2
    end
    return unorm
end
"pressure as a function of ρ,u,v,E"
function pfun(rho,rhoU,E,rhounorm)
    return @. (γ-1)*(E-.5*rhounorm)
end
function pfun(rho,rhoU,E)
    rhounorm = Unorm(rhoU)./rho
    return pfun(rho,rhoU,E,rhounorm)
end

"Constants"
const sp_tol = 1e-12

"Program parameters"
compute_L2_err = false

"Approximation Parameters"
N_P   = 2;    # The order of approximation in polynomial dimension
Np_P  = Int((N_P+1)*(N_P+2)/2)
Np_F  = 8;    # The order of approximation in Fourier dimension
K1D   = 10;   # Number of elements in polynomial (x,y) dimension
CFL   = 1.0;
T     = 5.0;  # End time

"Time integration Parameters"
rk4a,rk4b,rk4c = rk45_coeffs()
CN = (N_P+1)*(N_P+2)*3/2  # estimated trace constant for CFL
dt = CFL * 2 / CN / K1D
Nsteps = convert(Int,ceil(T/dt))
dt = T/Nsteps

"Initialize Reference Element in Fourier dimension"
h = 2*pi/Np_F
column = [0; .5*(-1).^(1:Np_F-1).*cot.((1:Np_F-1)*h/2)]
column2 = [-pi^2/3/h^2-1/6; -((-1).^(1:Np_F-1)./(2*(sin.((1:Np_F-1)*h/2)).^2))]
Dt = Toeplitz(column,column[[1;Np_F:-1:2]])
D2t = Toeplitz(column2,column2[[1;Np_F:-1:2]])
t = LinRange(h,2*pi,Np_F)

"Initialize Reference Element in polynomial dimension"
rd = init_reference_tri(N_P);
@unpack fv,Nfaces,r,s,VDM,V1,Dr,Ds,rf,sf,wf,nrJ,nsJ,rq,sq,wq,Vq,M,Pq,Vf,LIFT = rd
Nq_P = length(rq)
Nfp_P = length(rf)
Nh_P = Nq_P+Nfp_P # Number of hybridized points
Lq = LIFT

"Mesh related variables"
# First initialize 2D triangular mesh
VX,VY,EToV = uniform_tri_mesh(K1D,K1D)
@. VX = 1+VX
@. VY = 1+VY
md = init_mesh((VX,VY),EToV,rd)
VX = repeat(VX,2)
VY = repeat(VY,2)
VZ = [2/Np_F*ones((K1D+1)*(K1D+1),1); 2*ones((K1D+1)*(K1D+1),1)]
EToV = [EToV EToV.+(K1D+1)*(K1D+1)]

# Make domain periodic
@unpack Nfaces,Vf = rd
@unpack xf,yf,K,mapM,mapP,mapB = md
LX,LY = (x->maximum(x)-minimum(x)).((VX,VY)) # find lengths of domain
mapPB = build_periodic_boundary_maps(xf,yf,LX,LY,Nfaces*K,mapM,mapP,mapB)
mapP[mapB] = mapPB
@pack! md = mapP

# Initialize 3D mesh
@unpack x,y,xf,yf,xq,yq,rxJ,sxJ,ryJ,syJ,J,nxJ,nyJ,sJ,mapM,mapP,mapB = md
x,y,xf,yf,xq,yq,rxJ,sxJ,ryJ,syJ,J,nxJ,nyJ,sJ = (x->reshape(repeat(x,inner=(1,Np_F)),size(x,1),Np_F*K)).((x,y,xf,yf,xq,yq,rxJ,sxJ,ryJ,syJ,J,nxJ,nyJ,sJ))
z,zq,zf = (x->reshape(repeat(collect(2/Np_F:(2/Np_F):2),inner=(1,x),outer=(K,1))',x,Np_F*K)).((Np_P,Nq_P,Nfp_P))
mapM = reshape(1:Nfp_P*Np_P*K,Nfp_P,Np_P*K)
mapP_2D = (x->mod1(x,Nfp_P)+div(x-1,Nfp_P)*Nfp_P*Np_F).(mapP)
mapP = reshape(repeat(mapP_2D,inner=(1,Np_F)),Nfp_P,Np_F,K)
for j = 1:Np_F
    mapP[:,j,:] = mapP[:,j,:].+(j-1)*Nfp_P
end
mapP = reshape(mapP,Nfp_P,Np_F*K)

# scale by Fourier dimension
M = h*M
wq = h*wq
wf = h*wf

# Hybridized operators
Vh = [Vq;Vf]
rxJ,sxJ,ryJ,syJ = (x->mapslices((y->Vh*y),x,dims=(1,2))).((rxJ,sxJ,ryJ,syJ))
Ef = Vf*Pq
Br,Bs = (x->diagm(wf.*x)).((nrJ,nsJ))
Qr,Qs = (x->Pq'*M*x*Pq).((Dr,Ds))
Qrh,Qsh = (x->1/2*[x[1]-x[1]' Ef'*x[2];
                   -x[2]*Ef   x[2]]).(((Qr,Br),(Qs,Bs)))
Qrh_skew,Qsh_skew = (x->1/2*(x-x')).((Qrh,Qsh))
Qt = Dt
Qth = Qt # Not the SBP operator, weighted when flux differencing

# TODO: assume mesh uniform affine, so Jacobian are constants
# TODO: fix other Jacobian parts
JP = 1/K1D^2
JF = 1/pi
J = JF*JP
wq = J*wq
wf = JF*wf
Lq = 1/JP*Lq
Qrh = JF*Qrh
Qsh = JF*Qsh
Qth = JP*Qth
Qrh_skew = 1/2*(Qrh-Qrh')
Qsh_skew = 1/2*(Qsh-Qsh')

# TODO: refactor
ops = (Vq,Vf,wq,wf,Pq,Lq,Qrh_skew,Qsh_skew,Qth)
mesh = (rxJ,sxJ,ryJ,syJ,sJ,nrJ,nsJ,nxJ,nyJ,JP,JF,J,h,mapM,mapP)
param = (K,Np_P,Nfp_P,Np_F,Nq_P,Nh_P)
function rhs(Q,ops,mesh,param,compute_rhstest)

    Vq,Vf,wq,wf,Pq,Lq,Qrh_skew,Qsh_skew,Qth = ops
    rxJ,sxJ,ryJ,syJ,sJ,nrJ,nsJ,nxJ,nyJ,JP,JF,J,h,mapM,mapP = mesh
    K,Np_P,Nfp_P,Np_F,Nq_P,Nh_P = param
    w = [wq;wf]

    # Entropy projection
    VU = v_ufun(Q...)
    vh = (x->[Vq;Vf]*Pq*x).(VU)
    (ρ,ρu,ρv,ρw,E) = u_vfun(vh...)
    Uf = (x->x[Nq_P+1:end,:]).((ρ,ρu,ρv,ρw,E))

    # Convert to rho,u,v,beta vars
    β = betafun(ρ,ρu,ρv,ρw,E)
    Qh = (ρ,ρu./ρ,ρv./ρ,ρw./ρ,β)

    # Compute face values
    QM = (x->x[Nq_P+1:end,:]).(Qh)
    QP = (x->x[mapP]).(QM)

    # Lax Friedrichs Dissipation flux
    # TODO: Fix Jacobian. incorrect sJ?
    (ρM,ρuM,ρvM,ρwM,EM) = Uf
    ρuM_n = @. (ρuM*nxJ+ρvM*nyJ)/sJ # TODO: 3D lax-friedrichs?
    lam = abs.(wavespeed(ρM,ρuM_n,EM))
    LFc = .5*max.(lam,lam[mapP]).*sJ
    # LFc = 0 #TODO: TEST

    fSx,fSy,_ = euler_fluxes(QM,QP)
    normal_flux(fx,fy,u) = fx.*nxJ + fy.*nyJ - LFc.*(u[mapP]-u)
    flux = normal_flux.(fSx,fSy,Uf)
    rhsQ = (x->Vq*Lq*x).(flux)

    # Flux differencing
    ∇fh = [zeros(size(Qh[1])) for i in eachindex(Q)]

    # TODO: fix Jacobian
    for k = 1:K
        for nf = 1:Np_F
            j = nf+Np_F*(k-1)
            Qhj = [(x->x[i,j]).(Qh) for i = 1:Nh_P]
            Fxj = [euler_fluxes(UL,UR)[1] for UL in Qhj, UR in Qhj]
            Fyj = [euler_fluxes(UL,UR)[2] for UL in Qhj, UR in Qhj]
            for i = 1:5
                ∇fh[i][:,j] += 2*((rxJ[1,j]*Qrh_skew+sxJ[1,j]*Qsh_skew).*(x->x[i]).(Fxj)
                                 +(ryJ[1,j]*Qrh_skew+syJ[1,j]*Qsh_skew).*(x->x[i]).(Fyj))*ones(Nh_P,1)
            end
        end
    end

    for k = 1:K
        for nh = 1:Nh_P
            j_idx = (k-1)*Np_F+1:k*Np_F
            Qhi = [(x->x[nh,j]).(Qh) for j = j_idx]
            Fzi = [euler_fluxes(UL,UR)[3] for UL in Qhi, UR in Qhi]
            for i = 1:5
                # TODO: Why multiply by 2?
                ∇fh[i][nh,(k-1)*Np_F+1:k*Np_F] += 2*[wq/J;zeros(Nfp_P,1)][nh]*(Qth.*(x->x[i]).(Fzi))*ones(Np_F,1)
            end
        end
    end

    ∇f = (x->Vq*[Pq Lq]*diagm(1 ./ w)*x).(∇fh)
    rhsQ = @. -(∇f+rhsQ)

    rhstest = 0
    if compute_rhstest
        for fld in eachindex(rhsQ)
            rhstest += sum(wq.*VU[fld].*rhsQ[fld])
        end
    end

    return rhsQ,rhstest
end

xq,yq,zq = (x->reshape(x,Nq_P,Np_F*K)).((xq,yq,zq))
# All directions
println(" ")
println("======= All directions =======")
ρ_exact(x,y,z,t) = @. 1+0.2*sin(pi*(x+y+z-3/2*t))
ρ = @. 1+0.2*sin(pi*(xq+yq+zq))
u = ones(size(xq))
v = -1/2*ones(size(xq))
w = ones(size(xq))
p = ones(size(xq))
Q_exact(x,y,z,t) = (ρ_exact(x,y,z,t),u,v,w,p)

# # x direction test case
# println(" ")
# println("======= x direction =======")
# ρ_exact(x,y,z,t) = @. 2+ .5*sin(pi*(x-t))
# ρ = ρ_exact(xq,yq,zq,0)
# u = ones(size(xq))
# v = zeros(size(xq))
# w = zeros(size(xq))
# p = ones(size(xq))

# # y direction test case
# println(" ")
# println("======= y direction =======")
# ρ_exact(x,y,z,t) = @. 2+ .5*sin(pi*(y-t))
# ρ = ρ_exact(xq,yq,zq,0)
# u = zeros(size(xq))
# v = ones(size(xq))
# w = zeros(size(xq))
# p = ones(size(xq))

# # z direction test case
# println(" ")
# println("======= z direction =======")
# ρ_exact(x,y,z,t) = @. 2+ .5*sin(pi*(z-t))
# ρ = ρ_exact(xq,yq,zq,0)
# u = zeros(size(xq))
# v = zeros(size(xq))
# w = ones(size(xq))
# p = ones(size(xq))

Q = primitive_to_conservative(ρ,u,v,w,p)
Q = collect(Q)
resQ = [zeros(size(Q[1])) for _ in eachindex(Q)]
# rhs(Q,ops,mesh,param,false)
# @btime rhs(Q,ops,mesh,param,false)
@btime begin

for i = 1:Nsteps
    rhstest = 0
    for INTRK = 1:5
        compute_rhstest = INTRK==5
        rhsQ,rhstest = rhs(Q,ops,mesh,param,compute_rhstest)
        @. resQ = rk4a[INTRK]*resQ + dt*rhsQ
        @. Q += rk4b[INTRK]*resQ
    end

    if i%10==0 || i==Nsteps
        println("Time step: $i out of $Nsteps with rhstest = $rhstest")
    end
end

end # end time

rq2,sq2,wq2 = quad_nodes_2D(N_P+2)
Vq2 = vandermonde_2D(N_P,rq2,sq2)/VDM
xq2,yq2,zq2 = (x->Vq2*x).((x,y,z))
ρ = Vq2*Pq*Q[1]
ρ_ex = ρ_exact(xq2,yq2,zq2,T)
Q = (x->Vq2*Pq*x).(Q)
p = pfun(Q[1],(Q[2],Q[3],Q[4]),Q[5])
Q = (Q[1],Q[2]./Q[1],Q[3]./Q[1],Q[4]./Q[1],p)
Q_ex = Q_exact(xq2,yq2,zq2,T)

# # Check velocity and pressure
# println("Velocity (Avg by entries)")
# @show sum(Q[2]./Q[1])/size(Q[1],1)/size(Q[1],2)
# @show sum(Q[3]./Q[1])/size(Q[1],1)/size(Q[1],2)
# @show sum(Q[4]./Q[1])/size(Q[1],1)/size(Q[1],2)
# println("Pressure")
# p = pfun(Q[1],(Q[2],Q[3],Q[4]),Q[5])
# @show sum(p/size(Q[1],1))/size(Q[1],2)
# println("L2 error density")
# @show L2_err_ρ = sum(h*J*wq2.*(ρ-ρ_ex).^2)
#

L2_err = 0.0
for fld in eachindex(Q)
    global L2_err
    L2_err += sum(h*J*wq2.*(Q[fld]-Q_ex[fld]).^2)
end
println("L2err at final time T = $T is $L2_err\n")

max_stars_repo_name: yiminllin/ESDG-Fourier

max_issues_repo_name: yiminllin/ESDG-Fourier

max_forks_repo_name: yiminllin/ESDG-Fourier


text:
const LARGE_GFILE = "https://drive.google.com/file/d/1OiX6gEWRm57kb1H8L0K_HWN_pzc-sk8y/view?usp=sharing"
const SMALL_GDRIVE = "https://drive.google.com/file/d/1BENwlCOlGEMF8zsM4rC-SYVxyL7f8xw0/view?usp=sharing"
const GFILE_NAME = "FetchTest"

@testset "large file" begin
    f = gdownload(LARGE_GFILE, pwd())
    DataDeps.unpack(f)

    open(joinpath(pwd(), GFILE_NAME, "$GFILE_NAME.txt"), "r") do file
        @test readline(file) == "Test"
    end

    rm(GFILE_NAME, recursive=true, force=true)
end

@testset "small file" begin
    gdownload(SMALL_GDRIVE, pwd())
    open(joinpath(pwd(), "$GFILE_NAME.txt"), "r") do file
        @test readline(file) == "Test"
    end

    rm("$GFILE_NAME.txt", recursive=true, force=true)
end

const GSHEET = "https://docs.google.com/spreadsheets/d/1rwoDt5HrcP6TTgBe3BWxp7kPdVMWuS1cGVZp7BXhPwc/edit?usp=sharing"
const GSHEET_NAME = "FetchTest"

@testset "sheet" begin
    gdownload(GSHEET, pwd())
    open(joinpath(pwd(), "$GSHEET_NAME-1.csv"), "r") do file
        @test readline(file) == "Test,1"
    end

    rm("$GSHEET_NAME-1.csv", recursive=true, force=true)
end

@testset "unshortlink" begin
    @test Fetch.unshortlink("https://bit.ly/3ir0gYu") == LARGE_GFILE
    @test Fetch.unshortlink("https://bit.ly/3yqqdwK") == SMALL_GDRIVE
    @test Fetch.unshortlink("https://bit.ly/3yDDi69") == GSHEET
end


text:
using IntervalRootFinding
using Base.Test

include("bisect.jl")
include("findroots.jl")
include("roots.jl")
include("newton1d.jl")


text:
const CMDSTAN_HOME_VAR = if haskey(ENV, "CMDSTAN")
    "CMDSTAN"
else
    "JULIA_CMDSTAN_HOME"
end

"""

Get the path to the `CMDSTAN_HOME` environment variable.

$(SIGNATURES)

# Extended help

Example: `get_cmdstan_home()`
"""
function get_cmdstan_home()
    get(ENV, CMDSTAN_HOME_VAR) do
        throw(ErrorException("The environment variable $CMDSTAN_HOME_VAR needs to be set."))
    end
end

"""

Set the path for `CMDSTAN_HOME`.

$(SIGNATURES)

# Extended help

Example: `set_cmdstan_home!(homedir() * "/Projects/Stan/cmdstan/")`
"""
set_cmdstan_home!(path) = global CMDSTAN_HOME = path


text:
# Autogenerated wrapper script for Pinchot_jll for x86_64-w64-mingw32-cxx11
export libpinchot

JLLWrappers.@generate_wrapper_header("Pinchot")
JLLWrappers.@declare_library_product(libpinchot, "libpinchot.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libpinchot,
        "bin\\libpinchot.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Resources

const resources = Dict{String, String}()

resource(f, name = basename(f)) = (@assert isfile(f); resources[name] = f)

const resroute =
  branch(req -> length(req[:path]) == 1 && haskey(resources, req[:path][1]),
         req -> d(:body => open(read, resources[req[:path][1]]),
                  :headers => Mux.fileheaders(req[:path][1])))

# Server setup

const maintp = Mustache.template_from_file(joinpath(dirname(@__FILE__), "main.html"))

app(f) = req -> render(maintp, d("id"=>Page(f).id))

function page_handler(req)
  id = try parse(Int, req[:params][:id]) catch e @goto fail end
  haskey(pool, id) || @goto fail
  active(pool[id].value) && @goto fail

  callback_id = try split(req[:query], "=")[2] catch e nothing end
  callback_script = callback_id != nothing ? """var callback_id = $callback_id""" : ""
  return render(maintp, d("id"=>id,
                          "optional_create_window_callback"=>callback_script))

  @label fail
  return d(:body => "Not found",
           :status => 404)
end

function ws_handler(req)
  id = try parse(Int, req[:path][end]) catch e @goto fail end
  client = req[:socket]
  haskey(pool, id) || @goto fail
  p = pool[id].value
  active(p) && @goto fail

  p.sock = client
  @async @errs get(handlers(p), "init", identity)(p)
  put!(p.cb, true)
  while active(p)
    local data
    try
      data = read(client)
    catch e
      if (isa(e, ArgumentError) && contains(e.msg, "closed")) || isa(e, WebSockets.WebSocketClosedError)
        handle_message(p, d("type"=>"close", "data"=>nothing))
        yield() # Prevents an HttpServer task error (!?)
        return
      else
        rethrow()
      end
    end
    @errs handle_message(p, JSON.parse(String(data)))
  end
  return

  @label fail
  close(client)
end

http_default =
  mux(Mux.defaults,
      resroute,
      page(":id", page_handler),
      Mux.notfound())

ws_default =
  mux(Mux.wdefaults,
      ws_handler)

const serving = Ref(false)

function serve()
  serving[] && return
  serving[] = true
  http = Mux.http_handler(Mux.App(http_default))
  ws = Mux.ws_handler(Mux.App(ws_default))
  @async begin
    # Suppress log output from HTTP.jl unless the user has already
    # configured a custom logger. See:
    # https://github.com/JuliaLang/julia/pull/28229#issuecomment-410229612
    if current_task().logstate === nothing
      with_logger(NullLogger()) do
        WebSockets.serve(WebSockets.ServerWS(http, ws), ip"127.0.0.1", port[], false)
      end
    else
      WebSockets.serve(WebSockets.ServerWS(http, ws), ip"127.0.0.1", port[], false)
    end
  end
end


text:
function generate(file)
    dir = joinpath(@__DIR__, "src", "examples")
    name = first(splitext(last(splitdir(file))))
    md = julia_example_to_markdown(file)
    # replace some things
    md = replace(md, "@__NAME__" => name)
    md = replace(md, "@__DIR__" => dir)
    code = """
        ```@setup $(name)
        include("$(@__FILE__)")
        write("$(name).jl", extract_code("$(file)"))
        ```
        ```julia
        $(extract_code(file))
        ```
        """
    md = replace(md, "@__CODE__" => code)
    open(joinpath(dir, name*".md"), "w") do f
        write(f, md)
        write(f,
            """

            *This page was rendered automatically based on
            [`$(name).jl`](https://github.com/KristofferC/JuAFEM.jl/blob/docs/examples/$(name).jl)*
            """)
    end
    return "examples/$(name).md"
end

# strip all comments
function extract_code(file)
    io = IOBuffer() # output
    for line in eachline(file, chomp = false)
        (!startswith(strip(line), "#") || isempty(strip(line))) && write(io, line)
    end
    strip(String(take!(io)))
end

function julia_example_to_markdown(file)
    io = IOBuffer()
    for line in eachline(file)
        if isempty(line)
            write(io, '\n')
        elseif startswith(strip(line), "#")
            write(io, replace(line, r"\s*#\s" => "", count = 1), '\n')
        else # code
            write(io, "```julia\n", line, "\n```\n")
        end
    end
    # collapse all adjacent code blocks
    replace(String(take!(io)), r"\n```\n+```julia" => "")
end


text:
## Drop in LaTeX mode.
## Default TeX engine is XeLaTeX
## Load font package

## Typeset tick labels in math font by escaping them in math mode $...$
## Either use a formatter, e.g
## or convert ticks using

module MathTicks

export ticksAsMath!, labelsAsMath!, toggleTeX, preamble

import Plots
import PyPlot
import PyCall

try
    global rcParams = PyCall.PyDict(PyPlot.matplotlib["rcParams"]);
catch
    error("Error importing rcParams from matplotlib. Check your PyCall and PyPlot!")
end

global preamble = MathTicks.rcParams["text.latex.preamble"]

function toggleTeX(;rcParams=MathTicks.rcParams)
    rcParams["text.usetex"] ⊻= true
end

function ticksAsMath!(A::Plots.Axis)
    values,labels = Plots.get_ticks(A)
    A[:ticks] = (values, map(x->"\$$x\$", labels))
end

function ticksAsMath!(P::Plots.Plot{Plots.PyPlotBackend})
    for p in P.subplots
        ticksAsMath!(Plots.get_axis(p,:x))
        ticksAsMath!(Plots.get_axis(p,:y))
        try
            ticksAsMath!(Plots.get_axis(p,:z))
        catch
        end
    end
end

function labelsAsMath!(P::Plots.Plot{Plots.PyPlotBackend})
    for p in P.subplots
        for s in p.series_list
            s[:label] = "\$$(s[:label])\$"
        end
    end
end

end


text:
module Day25

end

max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
include("7.1.1 (a+b arcsinh(c x))^n.jl")
include("7.1.2 (d x)^m (a+b arcsinh(c x))^n.jl")
include("7.1.3 (d+e x^2)^p (a+b arcsinh(c x))^n.jl")
include("7.1.4 (f x)^m (d+e x^2)^p (a+b arcsinh(c x))^n.jl")
include("7.1.5 u (a+b arcsinh(c x))^n.jl")
include("7.1.6 Miscellaneous inverse hyperbolic sine.jl")

max_stars_repo_path: src/rules/7 Inverse hyperbolic functions/7.1 Inverse hyperbolic sine/.jl


text:
using ArgParse
using DrWatson
@quickactivate
using BSON
using Random
using FileIO
using Statistics
using DataFrames

# pkgs which come from deserialized BSONs
# have to be present in the Main module
using ValueHistories
using LinearAlgebra

s = ArgParseSettings()
@add_arg_table! s begin
    "modelname"
		arg_type = String
		default = "sptn"
		help = "Name of the model from which to compute ensemble."
	"dataset"
		arg_type = String
		default = "iris"
		help = "Dataset name."
    "dataset_type"
		arg_type = String
		default = "tabular"
        help = "images | tabular"
    "max_seed"
        arg_type = Int
        default = 5
        help = "Number of seeds to go through with each dataset."
    "anomaly_classes"
        arg_type = Int
        default = 10
        help = "Number of anomaly_classes to go through with each dataset."
end

function ensemble_experiment(eval_directory, exp_directory, out_directory)
    eval_files = isdir(eval_directory) ? readdir(eval_directory, join=true) : String[]
    if length(eval_files) == 0
        @warn "There are no valid files in $eval_directory"
        return
    end

    # load the dataframe
    df = reduce(vcat, map(f -> load(f)[:df], eval_files))

    # select best based on criterion
    for criterion in [:val_auc, :val_tpr_5, :val_pat_10]
        best = sortperm(df, order(criterion, rev=true))
        for select_top in [5, 10] # 0 ... automatic ensemble size
            if (select_top > 0) && (length(best) > select_top)
                top = best[1:select_top]
            elseif (length(best) <= select_top)
                top = best
            else
                @warn "Automatic ensemble size is not yet implemented."
                continue
            end

            # get rid of the "eval_" prefix
            exp_files = replace.(basename.(eval_files[top]), "eval_" => "")

            # this assumes that evaluation and experiments files are in sync
            # more precisely for every evaluation file there exist it's exp.
            files_to_load = joinpath.(exp_directory, exp_files)
            results = if all(isfile.(files_to_load))
                load.(files_to_load)
            else
                @warn "Discrepancy between $exp_directory and $eval_directory"
                return
            end

            scores, ensemble = _init_ensemble(results)
            ensemble[:ensemble_files] = exp_files
            for method in [:max, :mean]
                for ignore in [true, false]
                    eagg, fit_t, _, _, _ = @timed aggregate_score!(deepcopy(ensemble), scores;
                                        method=method, ignore_nan=ignore)
                    parameters = (
                        modelname=ensemble[:modelname],
                        criterion=criterion,
                        size=select_top,
                        method=method,
                        ignore_nan=ignore)
                    eagg[:parameters] = parameters
                    eagg[:fit_t] = fit_t
                    eagg[:tr_eval_t] = 0.0
                    eagg[:tst_eval_t] = 0.0
                    eagg[:val_eval_t] = 0.0

                    savef = joinpath(out_directory, savename("ensemble", parameters, "bson"))
                    @info "Saving ensemble experiment to $savef"
                    tagsave(savef, eagg, safe = true)
                end
            end
        end
    end
end

const SPLITS = ["tr", "tst", "val"]

_prefix_symbol(prefix, s) = Symbol("$(prefix)_$(s)")

function _init_ensemble(results)
    ensemble = Dict{Symbol, Any}()      # new ensemble experiment dictionary
    r = first(results)                  # reference dictionary

    for key in vcat([:dataset, :modelname, :seed, :model], _prefix_symbol.(SPLITS, :labels))
        ensemble[key] = r[key]
    end
    # add anomaly class if present
    if Symbol("anomaly_class") in keys(r)
        ensemble[:anomaly_class] = r[:anomaly_class]
	elseif Symbol("ac") in keys(r)
        ensemble[:anomaly_class] = r[:ac]
	end

    ensemble[:ensemble_phash] = [hash(r[:parameters]) for r in results]

    # sum training times
    ensemble[:ensemble_fit_t] = sum([r[:fit_t] for r in results])
    ensemble[:ensemble_eval_t] = sum([r[:tr_eval_t] + r[:tst_eval_t] + r[:val_eval_t] for r in results])

    scores = Dict()
    for key in _prefix_symbol.(SPLITS, :scores)
        kscores = similar(r[key], (length(r[key]), length(results)))
        for (i, rr) in enumerate(values(results))
            kscores[:,i] .= rr[key][:]
        end
        scores[key] = kscores
    end
    scores, ensemble
end

function aggregate_score!(ensemble, scores, weights=nothing;
                            method=:max, ignore_nan=true)
    for (k, s) in scores
        if method == :max
            mask_nan_max = (x) -> (isnan(x) ? -Inf : x)
            ensemble[k] = ignore_nan ?
                        maximum(mask_nan_max, s, dims=2)[:] :
                        maximum(s, dims=2)[:]
        elseif method == :mean
            ensemble[k] = ignore_nan ?
                        [mean(filter(x -> !isnan(x), s[i,:])) for i in 1:size(s, 1)] :
                        mean(s, dims=2)[:]
        elseif method == :wsum
            ensemble[k] = sum(weights' .* s, dims=2)[:]
        else
            error("Unsupported ensemble aggregation.")
        end
    end
    ensemble
end

function main(args)
    @unpack modelname, dataset, dataset_type, max_seed, anomaly_classes = args
    if dataset_type == "tabular"
        for s in 1:max_seed
            eval_directory = datadir("evaluation/$(dataset_type)/$(modelname)/$(dataset)/seed=$(s)/")
            exp_directory = datadir("experiments/$(dataset_type)/$(modelname)/$(dataset)/seed=$(s)/")
            # for now the ensemble outputs are kept in separate directory
            out_directory = datadir("experiments_ensembles/$(dataset_type)/$(modelname)/$(dataset)/seed=$(s)/")

            # run experiment
            @info "Generating ensemble scores of $(modelname) on $(dataset):$(s)"
            ensemble_experiment(eval_directory, exp_directory, out_directory)
        end
    elseif dataset_type == "images"
        for s in 1:max_seed
            for ac in 1:anomaly_classes
                eval_directory = datadir("evaluation/$(dataset_type)/$(modelname)/$(dataset)/ac=$(ac)/seed=$(s)/")
                exp_directory = datadir("experiments/$(dataset_type)/$(modelname)/$(dataset)/ac=$(ac)/seed=$(s)/")
                # for now the ensemble outputs are kept in separate directory
                out_directory = datadir("experiments_ensembles/$(dataset_type)/$(modelname)/$(dataset)/ac=$(ac)/seed=$(s)/")

                # run experiment
                @info "Generating ensemble scores of $(modelname) on $(dataset):$(ac):$(s)"
                ensemble_experiment(eval_directory, exp_directory, out_directory)
            end
        end
    else
        @error "Unsupported dataset type."
    end
end

main(parse_args(ARGS, s))


text:
module MockBlockSolver
  import Base.convert, Base.show, Base.copy, Base.pointer

  using JuMP, BlockDecomposition

  importall MathProgBase.SolverInterface
  importall BlockDecomposition.BlockSolverInterface

  include("MockBlockSolverInterface.jl")

end # module


text:
function init_sketch(n, R; rseed=0)
    # Random.seed!(rseed)
    return randn(n, R), zeros(n, R)
end

# cache should have length R
function rank_one_update!(S, Ω, v, η; cache::AbstractVector=nothing)
    @assert length(cache) == size(S, 2)
    S .= (1 - η) .* S
    if isnothing(cache)
        tmp = Ω'*v
        BLAS.ger!(η, v, tmp, S)
    else
        mul!(cache, Ω', v)
        BLAS.ger!(η, v, cache, S)
    end
end

# TODO: maybe should make this unallocating for experiments?
# Numerical stable version of eq (5.3)
#  X̂ = S*pinv(Ω'*S)*S' = U*Λ*U'
function reconstruct(Ω, S; correction=false)
    n = size(S, 1)
    σ = sqrt(n)*eps()*norm(S)
    Sσ = S + σ .* Ω
    M = Ω'*S
    M .= 0.5 .* (M .+ M')
    F = nothing

    # Some error handling -- PosDefException on cholesky for first few iters
    # TODO: figure out why this occurs
    try
        F = cholesky(M)
    catch e
        !isa(e, PosDefException) && error("Sketch reconstruction error")
        M[diagind(M)] .+= minimum(eigvals(M)) .+ σ
        F = cholesky(M)
    end
    U, Σ, _ = svd(Sσ / F.U)
    Λ = Diagonal(max.(0, Σ.^2 .- σ))

    # Implements trace correction (Remark 6.1)
    # Note: correction doubles the error at worse; usually works better in practice
    if correction
        R = size(S, 2)
        correction_factor = (1.0 - tr(Λ))/R
        Λ.diag .= Λ.diag .+ correction_factor
    end

    return U, Λ
end

function construct_Xhat(soln::SCGALResults)
    return soln.UT*soln.ΛT*soln.UT'
end

# Avoids overhead of constructing Xhat fully
# cache is a n x R matrix
function compute_objective(C, U, Λ; cache=nothing)
    n = size(C, 1)
    if isnothing(cache)
        cache = C*U
    else
        mul!(cache, C, U)
    end

    obj_val = 0
    @views for i in 1:size(Λ, 1)
        obj_val += Λ.diag[i] * dot(U[:,i], cache[:,i])
    end

    return obj_val
end

# Avoids overhead of constructing Xhat fully
# cache is a n x R matrix
function compute_primal_infeas_mc(U, Λ; cache=nothing)
    n, R = size(U)
    @assert size(cache) == size(U)
    if isnothing(cache)
        cache = U*Λ
    else
        mul!(cache, U, Λ)
    end

    infeas_val = 0
    @views for i in 1:n
        infeas_val += (dot(cache[i,:], U[i,:]) - 1/n)^2
    end
    infeas_val = sqrt(infeas_val)
    return infeas_val / (1 + sqrt(1/n))
end

# MAXCUT objective with GW randomized rounding
# C = -0.25*(Diagonal(G*ones(n)) - G)
function max_cut_val(C, U, Λ)
    n, R = size(U)
    hp = randn(R)
    x = sign.(U*sqrt.(Λ)*hp)
    return -dot(x, C*x)
end


text:
thisdir = dirname(@__FILE__())
any(path -> path==thisdir, LOAD_PATH) || push!(LOAD_PATH, thisdir)
using Variogram
#using Plots
#gr()

const INDIR = ARGS[1]
const OUTDIR = ARGS[2]

@everywhere function makevariogram(indir, outdir, file)

    local decorrelationlength::Int64 = -1

    try

        # Load the dataset
        local datasetname::String = file[1:end - 4]
        print("$datasetname\n")

        local variogramfile::String = "$outdir/$datasetname.csv"
        local plotfile::String = "$outdir/$datasetname.png"
        local dataset::Array{Float64, 2} = readdlm("$indir/$file", '\t', Float64, '\n', header=true)[1]

        # Calculate the variogram
        local variogram = Variogram.variogram(dataset, 25)
        writedlm("$outdir/$datasetname.csv", variogram, ',')

        #Plots.scatter(variogram[:,1], variogram[:,2])

        # Get the decorrelation length of the variogram
        decorrelationlength = Variogram.getdecorrelationlength(variogram, true)

        #plotmarker = [[decorrelationlength, decorrelationlength] [0, maximum(variogram[:,2])]]
        #Plots.plot!(plotmarker[:,1], plotmarker[:,2], linewidth=1, color="red", linestyle=:dot)

        #Plots.png(plotfile)
    catch ex
        println(typeof(ex))
        showerror(STDOUT, ex)
        throw(ex)
    end

    return decorrelationlength
end

files = readdir(INDIR)

decorrelationlengths = pmap(f -> makevariogram(INDIR, OUTDIR, f), files, retry_delays=zeros(1000))

open("decorrelation_lengths.csv", "w") do f
    for i in 1:length(files)
        write(f, "$(files[i]),$(decorrelationlengths[i])\n")
    end
end


text:
# ##############################################################
# ### optimisation modelling sets
# ##############################################################
# ECH, TS, S, NAMES = Workflow.get_ech_ts_s_name(launcher);
# BUSES =  Workflow.get_bus(launcher, NAMES);
# units_by_kind = Workflow.get_units_by_kind(launcher);
# units_by_bus =  Workflow.get_units_by_bus(launcher, BUSES);

# println("NAMES : ", NAMES);
# println("BUSES : ", BUSES);
# println("units_by_bus : ", units_by_bus);

# println("Number of scenario ", length(S));
# println("Number of time step is ", length(TS));
# println("ECH : ", ECH);
# K_IMPOSABLE = Workflow.K_IMPOSABLE;
# K_LIMITABLE = Workflow.K_LIMITABLE;

# # netloads = Dict([(bus, ts, s) => launcher.uncertainties[bus, s, ts, ech] for bus in BUSES, ts in TS, s in S]);
# eod_slack = Dict([(ts,s) => 0.0 for s in S, ts in TS]);
# factor = Dict([name => 1 for name in NAMES]);
# for bus in BUSES
#     factor[bus] = -1;
# end
# for ((name, s, ts, ech), v) in launcher.uncertainties
#     eod_slack[ts,  s] += factor[name] * v;
# end
# println("eod_slack is $eod_slack");

function worse_case(launcher::Launcher, ech::DateTime)
    #######################################################################
    # for a given ts, at a given ech, find the worst case within
    # uncertainties for why the economical balance maximise the constraints
    #######################################################################
    # max rho
    # u in [lb, ub]
    # sum p = u
    # p empilement économique
    # rho est la surcharge minimale

    is_on = Dict([kvp[1] => true for kvp in launcher.units])
    units_cost = Dict();
    # name->P-maxP-startCost-propCost
    for unit in launcher.units
        name = unit[1];
        minP = unit[2][1];
        maxP = unit[2][2];
        startCost = unit[2][3];
        propCost = unit[2][4];

    end

    units_by_kind = Workflow.get_units_by_kind(launcher);
    println(is_on);
end

max_stars_repo_name: azized/pscopf

max_issues_repo_name: azized/pscopf

max_forks_repo_name: azized/pscopf


text:
"""
    -------------------------------------------------------------------
    Build the fundamental data structures to define vertices, edges and
    the graph itself.
    -------------------------------------------------------------------
"""

"""
    Defines the structure for an edge of a graph. It is defined by its
    'index'{Int}, 'source'{Int} and 'target'{Int}.

    'source' and 'target' refers to the node at the tail of the edge and
    the node of the head of the edge, respectively.
"""
mutable struct Edge
    index::Int
    source::Int
    target::Int
    function Edge(source::Int, target::Int)
        index = -1
        new(index, source, target)
    end
end

"""
    Defines the structure for a vertex/node of a graph. It is defined by its
    'index'{Int}, 'edges_source'{Array{Edge}} and 'edges_target'{Array{Edge}}.
"""
mutable struct Vertex
    index::Int
    edges_source::Array{Edge}
    edges_target::Array{Edge}
    function Vertex(index::Int)
        edges_source = Edge[]
        edges_target = Edge[]
        new(index, edges_source, edges_target)
    end
end

"""
    Graph struct to define graph structures.

    is_directed:
        {true, false} -> If true the graph is defined as directed,
        otherwise it is defined as directed.
"""
mutable struct Graph
    N::Int
    M::Int
    is_directed::Bool
    vertices::Array{Vertex}
    edges::Array{Edge}
    # ---------- Vertex properties ----------- #
    int_vproperties::Dict{String, Array{Int}}
    float_vproperties::Dict{String, Array{Float64}}
    string_vproperties::Dict{String, Array{String}}
    # ----------- Edge properties ------------ #
    int_eproperties::Dict{String, Array{Int}}
    float_eproperties::Dict{String, Array{Float64}}
    string_eproperties::Dict{String, Array{String}}
    # ---------------------------------------- #
    function Graph(is_directed::Bool)
        N = 0
        M = 0
        vertices = Vertex[]
        edges = Edge[]
        int_vproperties = Dict{String, Array{Int}}()
        float_vproperties = Dict{String, Array{Float64}}()
        string_vproperties = Dict{String, Array{String}}()
        int_eproperties = Dict{String, Array{Int}}()
        float_eproperties = Dict{String, Array{Float64}}()
        string_eproperties = Dict{String, Array{String}}()
        new(N, M, is_directed, vertices, edges, int_vproperties,
            float_vproperties, string_vproperties, int_eproperties,
            float_eproperties, string_eproperties)
    end
end


text:
loopSym = [:while, :for]

function addLoopInvariant!(loopBody::Expr, agreement::Agreement)
    check = createCheckExpressions(agreement)
    newExpr = Expr(:block)
    newExpr.args = [loopBody.args[2].args; check]

    innerExpr = copy(loopBody)
    innerExpr.args = [loopBody.args[1], newExpr]
    loopBody.head = :block
    loopBody.args = [check..., innerExpr]
end

"""
@loopinvariant

Takes in a loop expression and fills it with invariant checks.

```julia-repl
julia> a = 0
0

julia> @loopinvariant (a >= 0) for i in 1:10
           a -= 10
       end
ERROR: ContractBreachException: Breach on Loop Invariant Expression 'a >= 0'
...

julia> a = 100
100

julia> @loopinvariant (a >= 0) for i in 1:10
           a -= 10
       end

julia> a = 100
100

julia> i = 1
1

julia> @loopinvariant (a >= 0) (i % 2 == 0) while(i <= 11)
           a -= 10
           i += 2
       end
ERROR: ContractBreachException: Breach on Loop Invariant Expression 'i % 2 == 0'
...
```
"""
macro loopinvariant(expr...)
    loopAgreement = Agreement(loopInvariant, addLoopInvariant!, nothing)
    loopBody = nothing
    for item in expr
        if typeof(item) == Expr
            item.head in loopSym ? (loopBody = item; break) :
            push!(loopAgreement.expressions, item)
        else
            newExpr = Expr(:block)
            newExpr.args = [item]
            push!(loopAgreement.expressions, newExpr)
        end
    end
    loopAgreement(loopBody)
    return loopBody |> esc
end


text:
using Test,
      DataFrames,
      Jib

include("client.jl")
include("decoder.jl")
include("requests.jl")
include("roundtrip.jl")

@testset "TickType" begin

  @test Jib.Reader.Decoder.tickname( 0) == "BID_SIZE"

  @test Jib.Reader.Decoder.tickname(90) == "DELAYED_HALTED"

  @test Jib.Reader.Decoder.tickname(99) == "ETF_NAV_LOW"

end


text:
# # Macros and expression interpolation

using JuliaSyntax: SourceFile, SyntaxNode, parseall, child, setchild!

# The following shows that SyntaxNode works nicely for simple macros which
# just interpolate expressions into one another. In particular it shows how
# precise source information from multiple files can coexist within the same
# syntax tree.

# First, here's the functionality that we're going to implement as a normal
# Julia macro. It's similar to the standard @show macro.
macro show2(ex)
    name = sprint(Base.show_unquoted, ex)
    quote
        value = $(esc(ex))
        println($name, " = ", value)
        value
    end
end

# Now, let's implement the same expression interpolation but using SyntaxNode
# (and with a normal Julia function which we need to use, absent any deeper
# integration with the Julia runtime)
function at_show2(ex::SyntaxNode)
    name = sprint(show, MIME"text/x.sexpression"(), ex)
    quote
        value = $(esc(ex))
        println($name, " = ", value)
        value
    end
    # The following emulates the expression interpolation lowering which is
    # usually done by the compiler.
    # 1. Extract the expression literal as `block`
    tree = parseall(SyntaxNode, String(read(@__FILE__)), filename=@__FILE__)
    block = child(tree, 3, 2, 2, 1)
    # 2. Interpolate local variables into the block at positions of $'s
    # Interpolating a SyntaxNode `ex` is simple:
    setchild!(block, (1, 2), ex)
    # The interpolation of a Julia *value* should inherit the source location
    # of the $ interpolation expression. This is different to when substituting
    # in a SyntaxNode which should just be inserted as-is.
    setchild!(block, (2, 2),
              JuliaSyntax.interpolate_literal(child(block, 2, 2), name))
    block
end

# Usage of at_show2()

# Let's have some simple expression to pass to at_show2. This will be
# attributed to a different file foo.jl
s2 = parseall(SyntaxNode, "foo +\n42", filename="foo.jl", rule=:statement)

# Calling at_show2, we see that the precise source information is preserved for
# both the surrounding expression and the interpolated fragments.
println("\nInterpolation example")
s3 = at_show2(s2)


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct _Any <: ProtoType
    type_url::AbstractString
    value::Array{UInt8,1}
    _Any(; kwargs...) = (o=new(); fillunset(o); isempty(kwargs) || ProtoBuf._protobuild(o, kwargs); o)
end #mutable struct _Any

export _Any