text:
using AQFED, Test
using AQFED.Basket
using Printf

@testset "KornTable1" begin
    spot = 80.0
    strike = 100.0
    r = -0.005
    q = -0.05
    tte = 1.0
    vol = 0.2

    refPrices = [21.6921, 25.0293, 28.0074, 30.7427, 32.0412, 33.9187]
    weight = [0.25, 0.25, 0.25, 0.25]
    spot = [100.0, 100.0, 100.0, 100.0]
    strike = 100.0
    r = 0.0
    sigma = 0.4
    rhos = [0.1, 0.3, 0.5, 0.7, 0.8, 0.95]
    tte = 5.0
    p = DeelstraBasketPricer(3, 3)
    pl = DeelstraLBBasketPricer(3, 3)
    for (i, rho) in enumerate(rhos)
        correlation = [
            1.0 rho rho rho
            rho 1.0 rho rho
            rho rho 1.0 rho
            rho rho rho 1.0
        ]
        tvar = [sigma, sigma, sigma, sigma] .^ 2 .* tte
        q = 0.0
        forward = spot .* exp((r - q) * tte)
        discountFactor = exp(-r * tte)
        price = priceEuropean(pl, true, strike, discountFactor, spot, forward, tvar, weight, correlation)
        @printf("%.2f DeelstraLB %.4f %.2e\n", rho, price, price - refPrices[i])
        price = priceEuropean(p, true, strike, discountFactor, spot, forward, tvar, weight, correlation)
        @printf("%.2f Deelstra %.4f %.2e\n", rho, price, price - refPrices[i])
        @test isapprox(refPrices[i], price, atol = 1e-1)
    end
end

@testset "KornTable5" begin
    spot = 80.0
    strike = 100.0
    r = -0.005
    q = -0.05
    tte = 1.0
    vol = 0.2

    refPrices = [19.450, 20.959, 22.995, 25.370, 30.593]
    weight = [0.25, 0.25, 0.25, 0.25]
    spot = [100.0, 100.0, 100.0, 100.0]
    strike = 100.0
    r = 0.0
    sigmas = [0.05, 0.1, 0.15, 0.2, 0.3]
    rho = 0.5
    tte = 5.0
    p = DeelstraBasketPricer(3, 3)
    for (i, sigma) in enumerate(sigmas)
        correlation = [
            1.0 rho rho rho
            rho 1.0 rho rho
            rho rho 1.0 rho
            rho rho rho 1.0
        ]
        tvar = [1.0^2, sigma^2, sigma^2, sigma^2] .* tte
        q = 0.0
        forward = spot .* exp((r - q) * tte)
        discountFactor = exp(-r * tte)
        price = priceEuropean(p, true, strike, discountFactor, spot, forward, tvar, weight, correlation)
        @printf("%.2f Deelstra %.4f %.2e\n", rho, price, price - refPrices[i])
        @test isapprox(refPrices[i], price, atol = 1e-1)
    end
end

@testset "JuAsianTable4" begin
    #Table 4 of Ju "Pricing Asian and Basket Options via Taylor Expansions"
    spot = 100.0
    strikes = [95.0, 100.0, 105.0]
    r = 0.09
    q = 0.0
    tte = 3.0
    vol = 0.3
    nWeights = Int(tte * 52)
    weights = zeros(Float64, nWeights + 1)
    tvar = zeros(Float64, nWeights + 1)
    forward = zeros(Float64, nWeights + 1)
    for i = 1:length(weights)
        weights[i] = 1.0 / (nWeights + 1)
        ti = (i-1) / (nWeights) * tte
        tvar[i] = vol^2 * ti
        forward[i] = spot * exp((r - q) * ti)
    end
    discountFactor = exp(-r * tte)
    refPrices = [19.0144, 16.5766, 14.3830] #TS values
    p = DeelstraBasketPricer(1, 3)
    pl = DeelstraLBBasketPricer(1, 3)  #1,3 best results = GC of Ju paper.
    for (refPrice,strike) in zip(refPrices,strikes)
        price = priceAsianFixedStrike(p, true, strike, discountFactor, spot, forward, tvar, weights)
        @printf("%6.2f Deelstra    %.4f %.2e\n", strike, price, price - refPrice)
        @test isapprox(refPrice, price, atol = 1e-4)
        price = priceAsianFixedStrike(pl, true, strike, discountFactor, spot, forward, tvar, weights)
        @printf("%6.2f Deelstra-LB %.4f %.2e\n", strike, price, price - refPrice)
        @test isapprox(refPrice, price, atol = 1e-2)
    end

    vol = 0.5
    for i = 1:length(weights)
        ti = (i-1) / (nWeights) * tte
        tvar[i] = vol^2 * ti
    end
    refPrices = [24.5526, 22.6115, 20.8241] #TS values
    for (refPrice,strike) in zip(refPrices,strikes)
        price = priceAsianFixedStrike(p, true, strike, discountFactor, spot, forward, tvar, weights)
        @printf("%6.2f Deelstra    %.4f %.2e\n", strike, price, price - refPrice)
        @test isapprox(refPrice, price, atol = 2e-4)
        price = priceAsianFixedStrike(pl, true, strike, discountFactor, spot, forward, tvar, weights)
        @printf("%6.2f Deelstra-LB %.4f %.2e\n", strike, price, price - refPrice)
        @test isapprox(refPrice, price, atol = 3e-2)
    end
    #Deelstra is much more accurate than TE6, as expected with vols relatively high
end


text:
"""
    struct ScaledMeasure{R,M} <: AbstractMeasure
        logscale :: R
        base :: M
    end

"""
struct ScaledMeasure{R,M} <: AbstractMeasure
    logscale :: R
    base :: M
end

function logdensity(sm::ScaledMeasure{R,M}, x::X) where {X, R, M <: AbstractMeasure}
    logdensity(sm.base, x) + sm.logscale
end

function Base.:*(k::Number, m::AbstractMeasure)
    if iszero(k)
        return TrivialMeasure
    else
        return ScaledMeasure{typeof(k), typeof(m)}(log(k),m)
    end
end

Base.:*(m::AbstractMeasure, k::Real) = k * m

≪(::M, ::ScaledMeasure{R,M}) where {R,M} = true
≪(::ScaledMeasure{R,M}, ::M) where {R,M} = true

basemeasure(μ::ScaledMeasure{R,M}) where {R,M} = μ.base

sampletype(μ:: ScaledMeasure) = sampletype(μ.base)

ext: jl


text: include("svds/SVD.jl")


text:
using InvertibleNetworks, NNlib, LinearAlgebra, Test
using PyPlot, Flux, Random
import Flux.Optimise.update!

# Random seed
#Random.seed!(100)

# Data
nx = 1
ny = 1
n_in = 2
batchsize = 100
hidden_factor = 64   # use this times as many weights for the hidden layer
k = 1   # kernel size
s = 1   # stride
p = 0   # padding
depth = 40

###################################################################################################
# Hyperbolic layer

# Layer
HL = Array{HyperbolicLayer}(undef, depth)
AN = ActNorm(n_in; logdet=true)
Params = get_params(AN)
for j=1:depth
    HL[j] = HyperbolicLayer(nx, ny, Int(n_in/2), batchsize, k, s, p; action="same",  α=1f0, hidden_factor=4)
    global Params = cat(Params, get_params(HL[j]); dims=1)
end

# Forward pass
function forward(X)
    X, logdet = AN.forward(X)
    X_prev, X_curr = tensor_split(X)
    for j=1:depth
        X_prev, X_curr = HL[j].forward(X_prev, X_curr)
    end
    X = tensor_cat(X_prev, X_curr)
    return X, logdet
end

# Inverse pass
function inverse(X)
    X_curr, X_new = tensor_split(X)
    for j=depth:-1:1
        X_curr, X_new = HL[j].inverse(X_curr, X_new)
    end
    X = tensor_cat(X_curr, X_new)
    X = AN.inverse(X)
    return X
end

# Backward pass
function backward(ΔX, X)
    ΔX_curr, ΔX_new = tensor_split(ΔX)
    X_curr, X_new = tensor_split(X)
    for j=depth:-1:1
        ΔX_curr, ΔX_new, X_curr, X_new = HL[j].backward(ΔX_curr, ΔX_new, X_curr, X_new)
    end
    ΔX = tensor_cat(ΔX_curr, ΔX_new)
    X = tensor_cat(X_curr, X_new)
    ΔX, X = AN.backward(ΔX, X)
    return ΔX, X
end

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

# Loss
function loss(X)
    Y, logdet = forward(X)
    f = -log_likelihood(Y) - logdet
    ΔY = -∇log_likelihood(Y)
    backward(ΔY, Y)
    return f
end

# Training
maxiter = 4000
opt = Flux.ADAM(1f-2)
fval = zeros(Float32, maxiter)
for j=1:maxiter

    # Evaluate objective and gradients
    X = sample_swirl(batchsize)
    fval[j] = loss(X)
    mod(j, 100) == 0 && (print("Iteration: ", j, "; f = ", fval[j], "\n"))

    # Update params
    for p in Params
        update!(opt, p.data, p.grad)
    end
    clear_grad!(Params)
end

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

# Evaluate network
test_size = 2000
X = sample_swirl(test_size)
Y = forward(X)[1]
Y_ = randn(Float32, nx, ny, n_in, test_size)
X_ = inverse(Y_)

# Plot
figure(figsize=[8,8])
ax1 = subplot(2,2,1); plot(X[1, 1, 1, :], X[1, 1, 2, :], "."); title(L"Data space: $x \sim \hat{p}_X$")
ax1.set_xlim([-15,15]); ax1.set_ylim([-15,18])
ax2 = subplot(2,2,2); plot(Y[1, 1, 1, :], Y[1, 1, 2, :], "g."); title(L"Latent space: $z = f(x)$")
ax2.set_xlim([-5, 5]); ax2.set_ylim([-5, 5])
ax3 = subplot(2,2,3); plot(X_[1, 1, 1, :], X_[1, 1, 2, :], "g."); title(L"Data space: $x = f^{-1}(z)$")
ax3.set_xlim([-15,15]); ax3.set_ylim([-15,18])
ax4 = subplot(2,2,4); plot(Y_[1, 1, 1, :], Y_[1, 1, 2, :], "."); title(L"Latent space: $z \sim \hat{p}_Z$")
ax4.set_xlim([-5, 5]); ax4.set_ylim([-5, 5])

max_stars_repo_name: slimgroup/LoopUnrolledSeismicImaging

max_issues_repo_name: slimgroup/LoopUnrolledSeismicImaging

max_forks_repo_name: slimgroup/LoopUnrolledSeismicImaging


text:
using Encoref
using Documenter

makedocs(;
    modules=[Encoref],
    authors="Endre Márk Borza",
    repo="https://github.com/endremborza/Encoref.jl/blob/{commit}{path}#L{line}",
    sitename="Encoref.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://endremborza.github.io/Encoref.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/endremborza/Encoref.jl",
)


text:
export note!, clear_notes!

# ============================================================================
# Annotation

# Adding a string to SeisData writes a note; if the string mentions a channel
# name or ID, the note is restricted to the given channels(s), else it's
# added to all channels
"""
    note!(S::SeisData, i::Int64, s::String)

Append `s` to channel `i` of `S` and time stamp.

    note!(S::SeisData, id::String, s::String)

As above for the first channel in `S` whose id is an exact match to `id`.
"""
note!(S::T, i::Int64, s::String) where {T<:GphysData} = push!(S.notes[i], tnote(s))

"""
    note!(S::SeisData, s::String)

Append `s` to `S.notes` and time stamp. If `txt` contains a channel name or ID, only the channel mentioned is annotated; otherwise, all channels are annotated.
"""
function note!(S::GphysData, s::String)
    J = [occursin(i, s) for i in S.name]
    K = [occursin(i, s) for i in S.id]
    j = findall(max.(J,K) .== true)
  if !isempty(j)
    [push!(S.notes[i], tnote(s)) for i in j]
  else
    for i = 1:S.n
      push!(S.notes[i], tnote(s))
    end
  end
  return nothing
end

function note!(S::GphysData, id::String, s::String)
  i = findid(id, S)
  (i == 0) && error(string("id = ", id, " not found in S!"))
  push!(S.notes[i], tnote(s))
  return nothing
end

function note!(S::GphysData, chans::Array{Int64,1}, s::String)
  for c in chans
    note!(S, c, s)
  end
  return nothing
end

note!(S::GphysChannel, s::String) = push!(S.notes, tnote(s))
# note!(R::SeisSrc, s::String) = push!(R.notes, tnote(s))
# note!(H::SeisHdr, s::String) = push!(H.notes, tnote(s))

# DND, these methods prevent memory reuse
"""
    clear_notes!(U::Union{SeisData,SeisChannel,SeisHdr})

Clear all notes from `U` and leaves a note about this.

    clear_notes!(S::SeisData, i::Int64, s::String)

Clear all notes from channel `i` of `S` and leaves a note about this.

    clear_notes!(S::SeisData, id::String, s::String)

As above for the first channel in `S` whose id is an exact match to `id`.
"""
function clear_notes!(S::GphysData)
  cstr = tnote("notes cleared.")
  for i = 1:S.n
    empty!(S.notes[i])
    push!(S.notes[i], identity(cstr))
  end
  return nothing
end

function clear_notes!(S::GphysData, i::Int64)
  empty!(S.notes[i])
  push!(S.notes[i], tnote("notes cleared."))
  return nothing
end

function clear_notes!(S::GphysData, id::String)
  i = findid(id, S)
  (i == 0) && error(string("id = ", id, " not found in S!"))
  empty!(S.notes[i])
  push!(S.notes[i], tnote("notes cleared."))
  return nothing
end

clear_notes!(U::SeisChannel) = (U.notes = Array{String,1}(undef,1); U.notes[1] = tnote("notes cleared."); return nothing)


text:
module SimpleValueGraphs

using LightGraphs
using LightGraphs.SimpleGraphs: AbstractSimpleGraph, AbstractSimpleEdge,
    SimpleEdge, IsDirected

using SimpleTraits

using Base: OneTo

import Base: eltype, show, reverse, iterate, length, replace_in_print_matrix,
    getindex, size, zero, Tuple

import SparseArrays: blockdiag

import LightGraphs:
    src, dst, edgetype, nv, ne, vertices, edges, is_directed,
    add_vertex!, add_edge!, rem_vertex!, rem_edge!,
    has_vertex, has_edge, inneighbors, outneighbors,
    indegree, outdegree, degree, has_self_loops, num_self_loops,
    add_vertex!, adjacency_matrix, weights,
    SimpleGraph, SimpleDiGraph, IsDirected,

    # operators
    complement

export AbstractValGraph,
    AbstractEdgeValGraph,
    EdgeValGraph, EdgeValOutDiGraph, EdgeValDiGraph,
    AbstractValEdge, ValEdge, ValDiEdge,
    get_edgeval, get_edgeval_or, get_edgevals, get_edgevals_or,
    set_edgeval!, set_edgevals!,
    val, vals,
    outedgevals, inedgevals,
    outedgevals_for_key, inedgevals_for_key,
    default_edgeval, edgevals_type,
    transform_edgevals!,
    KeyView, keyview,
    ValMatrix, AdjacencyMatrix

# ==== Includes ===========================

include("utils.jl")

include("abstractvaluegraph.jl")
include("valueedge.jl")
include("valuegraph.jl")
#include("keyview.jl")
include("matrices.jl")

#include("operators.jl")
#include("modified_functions.jl")

#include("export.jl") # Temporary disabled

end # module


text:
using DataFrames
using DecisionTree
using JLD

# Read the iris data set.
df = readtable(ARGS[1], header = false)

# Get the features and labels.
features = convert(Array, df[:, 1:4])
labels = convert(Array, df[:, 5])

# Train decision tree classifier.
model = DecisionTreeClassifier(pruning_purity_threshold=0.9, maxdepth=6)
DecisionTree.fit!(model, features, labels)

# Save the model.
save(ARGS[2], "model", model)

max_stars_repo_name: carlsonp/pachyderm

max_issues_repo_name: carlsonp/pachyderm

max_forks_repo_name: carlsonp/pachyderm


text:
using AlgebraPDF
using Test

@makefuntype Amazing2D(x;p) = (x[1]-p.x0)^2+(x[2]-p.y0)^2-p.R0^2
a = Amazing2D((x0=1.1, y0=2.1, R0=0.0))

@testset "Two-dimensional function" begin
    @test a((1.1,2.1)) == 0
    @test a((1.1+1,2.1+1)) ≈ 2
    data = collect(zip(rand(10), rand(10)))
    @test length(a(data)) == 10
end

@makefuntype Amazing3D(x;p) = (x[1]-p.x0)^2+(x[2]-p.y0)^2+(x[3]-p.z0)^2-p.R0^2
a = Amazing3D((x0=1.1, y0=2.1, z0=3.1, R0=0.0))

@testset "Two-dimensional function" begin
    @test a((1.1,2.1,3.1)) == 0
    @test a((1.1+1,2.1+1,3.1+1)) ≈ 3
    data = collect(zip(rand(10), rand(10), rand(10)))
    @test length(a(data)) == 10
end


text:
@testset "simple workflow" begin
    env = CartPoleEnv{Float32}() |> StateOverriddenEnv(deepcopy)
    agent = Agent(;
        policy = RandomPolicy(env),
        trajectory = CircularCompactSARTSATrajectory(;
            capacity = 10_000,
            state_type = Float32,
            state_size = (4,),
        ),
    )
    N_EPISODE = 10000
    hook = TotalRewardPerEpisode()
    run(agent, env, StopAfterEpisode(N_EPISODE), hook)

    @test isapprox(sum(hook.rewards) / N_EPISODE, 21; atol = 2)
end


text:
function make_vec(format_j)

  b = PetscVec(sys_size_global, format_j, comm)
  low, high = VecGetOwnershipRange(b)
  global_indices = collect(low:PetscInt(high - 1))

  idxm = global_indices + PetscInt(1)  # 1 based indexing
  vals = Array{PetscScalar}(sys_size_local)
  for i=1:sys_size_local
    vals[i] = rhs[i]
  end
  set_values1!(b, idxm, vals)
  PETSc2.VecAssemble(b)

  return b
end

#------------------------------------------------------------------------------
# functions that run tests

# test indexing
function test_indexing(format_j)
  println("testing indexing")

  b = make_vec(format_j)
  #=
  b = PetscVec(sys_size_global, format_j, comm)
  idxm = global_indices + PetscInt(1)  # 1 based indexing
  vals = Array{PetscScalar}(sys_size_local)
  for i=1:sys_size_local
    vals[i] = i
  end
  set_values1!(b, idxm, vals)
  PETSc2.VecAssemble(b)
  =#
  vals2 = zeros(rhs)

  low, high = VecGetOwnershipRange(b)
  global_indices = collect(low:PetscInt(high - 1))
  idxm = global_indices + PetscInt(1)  # 1 based indexing

  VecGetValues(b, global_indices, vals2)

  @test isapprox( rhs, vals2) atol=1e-13
  fill!(vals2, 0.0)
  get_values1!(b,  idxm, vals2)
  @test isapprox( rhs, vals2) atol=1e-13

  b2s = zeros(3)
  idxm = collect(PetscInt, 1:3)
  vals = rand(3)
  set_values1!(b2s, idxm, vals)

  @test isapprox( b2s, vals) atol=1e-13

  set_values1!(b2s, idxm, vals, ADD_VALUES)

  @test isapprox( b2s, 2*vals) atol=1e-13
  vals2 = zeros(vals)
  get_values1!(b2s, idxm, vals2)

  @test isapprox( vals2, 2*vals) atol=1e-13

  for i=1:sys_size_local
    idxm = [global_indices[i]]   # index
    val = [ rhs[i] ]  # value
    VecSetValues(b, idxm, val, INSERT_VALUES)
  end

  PETSc2.VecAssemble(b)

  PetscDestroy(b)
  return nothing
end

# test copy, getArray
function test_copy(format_j)
  println("testing copy")
  b = make_vec(format_j)

  low, high = VecGetOwnershipRange(b)
  global_indices = collect(low:PetscInt(high - 1))

  # check that the vector was set/assembled correctly
  # check all the methods of copying/accesing a vector work
  b_copy = zeros(PetscScalar, sys_size_local)
  b2_copy = zeros(PetscScalar, sys_size_local)
#  idx = collect(0:2)
#  idx = Array{PetscInt}(3)
#  for i=1:sys_size_local
#    idx[i] = global_indices[i]
#  end

  b_arr = VecGetArray(b)
  b_arr_ro = VecGetArrayRead(b)
  b2 = VecDuplicate(b)
  VecCopy(b, b2)

  VecGetValues(b, sys_size_local, global_indices, b_copy)
  VecGetValues(b2, sys_size_local, global_indices, b2_copy)
  for i=1:sys_size_local
     @test isapprox( b_copy[i], rhs[i])
     @test isapprox( b_arr[i], rhs[i])
     @test isapprox( b_arr_ro[i], rhs[i])
     @test isapprox( b2_copy[i], rhs[i])
  end

  VecRestoreArray(b, b_arr)
  VecRestoreArrayRead(b, b_arr_ro)

  @test ( VecGetSize(b) )== comm_size*sys_size_local

  PetscDestroy(b)
  PetscDestroy(b2)
  return nothing
end

# test norms/BLAS stuff
function test_linalg(format_j)

  b = make_vec(format_j)
  b2 = copy(b)
#    b2 = VecDuplicate(b)
#    VecCopy(b, b2)
  # create 3rd vector to store results in
  b3 = PetscVec(sys_size_global, format_j, comm)

  low, high = VecGetOwnershipRange(b)
  global_indices = collect(low:PetscInt(high - 1))

  b_copy = zeros(PetscScalar, sys_size_local)
  b2_copy = zeros(PetscScalar, sys_size_local)
  rhs_tmp2 = zeros(PetscScalar, sys_size_local)
  for i=1:length(vec_norms)
    norm_petsc = VecNorm(b, vec_norms[i])
    norm_julia = norm(rhs, julia_vec_norms[i])
    @test isapprox( VecNorm(b, vec_norms[i]), norm(rhs_global, julia_vec_norms[i]))
    print("\n")
  end

  # check for non zero exit status is all we can really do here
  @test ( PetscView(b, 0) )== 0

  # test math functions
  VecSqrtAbs(b)
  rhs_tmp = deepcopy(rhs)  # don't modifiy original rhs
  VecGetValues(b, sys_size_local, global_indices, b_copy)
  for i=1:sys_size_local
    rhs_tmp[i] = sqrt(abs(rhs_tmp[i]))
    @test isapprox( b_copy[i], rhs_tmp[i])
  end

  VecLog(b)
  VecGetValues(b, sys_size_local, global_indices, b_copy)
  for i=1:sys_size_local
    rhs_tmp[i] = log(rhs_tmp[i])
    @test isapprox( b_copy[i], rhs_tmp[i])
  end

  VecExp(b)
  VecGetValues(b, sys_size_local, global_indices, b_copy)
  for i=1:sys_size_local
    rhs_tmp[i] = exp(rhs_tmp[i])
    @test isapprox( b_copy[i], rhs_tmp[i])
  end

  VecAbs(b)
  VecGetValues(b, sys_size_local, global_indices, b_copy)
  for i=1:sys_size_local
    rhs_tmp[i] = abs(rhs_tmp[i])
    @test isapprox( b_copy[i], rhs_tmp[i])
  end

  petsc_r, petsc_idx = VecMax(b)
  julia_r = maximum(real(b_copy))
  julia_idx = indmax(real(b_copy))

  @test isapprox( petsc_r, julia_r)
  @test ( petsc_idx )== (julia_idx - 1)

  petsc_r, petsc_idx = VecMin(b)
  julia_r = minimum(real(b_copy))
  julia_idx = indmin(real(b_copy))

  @test isapprox( petsc_r, julia_r)
  @test ( petsc_idx )== (julia_idx - 1)

  VecReciprocal(b)
  VecGetValues(b, sys_size_local, global_indices, b_copy)
  for i=1:sys_size_local
    rhs_tmp[i] = 1/rhs_tmp[i]
    @test isapprox( b_copy[i], rhs_tmp[i])
  end

  VecShift(b, PetscScalar(2.0))
  VecGetValues(b, sys_size_local, global_indices, b_copy)
  for i=1:sys_size_local
    rhs_tmp[i] = rhs_tmp[i] + 2.0
    @test isapprox( b_copy[i], rhs_tmp[i])
  end

  VecPointwiseMult(b3, b, b2)
  VecGetValues(b3, sys_size_local, global_indices, b_copy)
  for i=1:sys_size_local
    rhs_tmp2[i] = rhs_tmp[i]*rhs[i]
    @test isapprox( b_copy[i], rhs_tmp2[i])
  end

  VecPointwiseDivide(b3, b, b2)
  VecGetValues(b3, sys_size_local, global_indices, b_copy)
  for i=1:sys_size_local
    rhs_tmp2[i] = rhs_tmp[i]/rhs[i]
    @test isapprox( b_copy[i], rhs_tmp2[i])
  end

  # test vector multiplication, addition functions
  # use b, b2 as the Petsc vectors
  # use rhs_tmp, rhs as the julia vectors

  alpha = PetscScalar(2.3)
  beta = PetscScalar(3.5)
  gamma = PetscScalar(4.8)

  VecAXPY(b, alpha, b2)
  rhs_tmp += alpha*rhs
  VecGetValues(b, sys_size_local, global_indices, b_copy)
  for i=1:sys_size_local
    @test isapprox( b_copy[i], rhs_tmp[i])
  end

  #=
  VecAXPBY(b, alpha, beta, b2)
  rhs_tmp = alpha*rhs + beta*rhs_tmp
  VecGetValues(b, sys_size_local, global_indices, b_copy)
  for i=1:sys_size_local
    @test isapprox( b_copy[i], rhs_tmp[i])
  end
=#

  VecAYPX(b, alpha, b2)
  rhs_tmp = alpha*rhs_tmp + rhs
  VecGetValues(b, sys_size_local, global_indices, b_copy)
  for i=1:sys_size_local
    @test isapprox( b_copy[i], rhs_tmp[i])
  end

  VecWAXPY(b3, alpha, b, b2)
  rhs_tmp2 = alpha*rhs_tmp + rhs
  VecGetValues(b3, sys_size_local, global_indices, b_copy)
  for i=1:sys_size_local
    @test isapprox( b_copy[i], rhs_tmp2[i])
  end

  # MAXPY
  scalar_arr = [alpha, beta]
  vec_arr = [b2.pobj, b3.pobj]
  VecMAXPY(b, PetscInt(2), scalar_arr, vec_arr)
  rhs_tmp += alpha*rhs + beta*rhs_tmp2
  VecGetValues(b, sys_size_local, global_indices, b_copy)
  for i=1:sys_size_local
    @test isapprox( b_copy[i], rhs_tmp[i])
  end

  VecAXPBYPCZ(b, alpha, beta, gamma, b2, b3)
  rhs_tmp = alpha*rhs + beta*rhs_tmp2 + gamma*rhs_tmp
  VecGetValues(b, sys_size_local, global_indices, b_copy)
  for i=1:sys_size_local
    @test isapprox( b_copy[i], rhs_tmp[i])
  end

  VecScale(b, alpha)
  rhs_tmp *= alpha
   VecGetValues(b, sys_size_local, global_indices, b_copy)
  for i=1:sys_size_local
    @test isapprox( b_copy[i], rhs_tmp[i])
  end

  petsc_r = VecDot(b, b2)
  julia_r = comm_size*rhs'*rhs_tmp  # note reversed b, b2

  @test isapprox( petsc_r, julia_r[1])

#=
  petsc_r = VecTDot(b, b2)
  julia_r = rhs_tmp.'*rhs

  @test isapprox( petsc_r, julia_r[1])

  println("finished testing VecTDot")

  petsc_r =  VecSum(b2)
  julia_r = sum(rhs_global)

  @test isapprox( petsc_r, julia_r)

  println("finished testing VecSum")

  VecSwap(b, b2)
  b2_copy = zeros(PetscScalar, sys_size_local)
  VecGetValues(b, sys_size_local, global_indices, b_copy)
  VecGetValues(b2, sys_size_local, global_indices, b2_copy)
  for i=1:sys_size_local
    @test isapprox( b_copy[i], rhs[i])
    @test isapprox( b2_copy[i], rhs_tmp[i])
  end

  println("finished testing VecSwap")

  VecSwap(b, b2)  # swap back

=#

  PetscDestroy(b)
  PetscDestroy(b2)
  PetscDestroy(b3)
end

@testset "\n   ---Testing vector functions---" begin

for j=1:length(vec_formats)
#  format_j = "standard"
  format_j = vec_formats[j]
  println("testing vector format ", format_j)

  test_indexing(format_j)
  test_copy(format_j)
  test_linalg(format_j)

  end  # end loop over formats
end # end fact check block


text:
include("utils.jl")

""" Cutting plane facility location model.
    Note this is just the nominal problem, but with linear policies y(z) = u + Vz when z = 0. """
function CP_facility_model(c::Matrix, f::Vector)
    n, m = size(c)
    @assert length(f) == n
    model = Model(Gurobi.Optimizer)

    # VARIABLES
    @variable(model, x[1:n], Bin)     # Facility locations
    @variable(model, V[1:n, 1:m, 1:m])  # Linear decision rule w.r.t. z
    @variable(model, u[1:n, 1:m])  # Fixed policy

    # CONSTRAINTS
    for j=1:m # Demand constraints
        dt = zeros(n, m);
        [dt[i,j] = 1 for i=1:n];
        @constraint(model, sum(dt .* u) >= d[j])
    end
    for i=1:n # Capacity constraints
        dt = zeros(n, m);
        [dt[i,j] = 1 for j=1:m];
        @constraint(model, sum(dt .* u) <= s[i] * x[i])
    end
    # OBJECTIVE
    @variable(model, F)
    @constraint(model, F >= sum(f[j] * x[j] for j = 1:n) + sum(c[i, j] * u[i, j] for i=1:n, j=1:m))
    @objective(model, Min, F)
    return model, x, u, V
end

function apply_heuristic(model, x, u, V)
    @constraint(model, u .>= 0)
    @constraint(model, V .>= 0)
    return
end

""" Finds and adds worst case cuts for the facility location problem. """
function find_wc_cuts(model, x, u, V, xvals, uvals, Vvals, rho, Gamm, atol = 1e-5)
    obj_value = objective_value(model)
    wc_model = Model(Gurobi.Optimizer)
    set_optimizer_attribute(wc_model, "OutputFlag", 0) # suppressing printouts.
    @variable(wc_model, -rho <= z[1:m] <= rho)
    @variable(wc_model, normdummy[1:m] >= 0)
    @constraint(wc_model, [i=1:m], normdummy[i] >= z[i])
    @constraint(wc_model, [i=1:m], normdummy[i] >= -z[i])
    @constraint(wc_model, sum(normdummy) <= Gamm)

    count = 0
    for i = 1:n # nonnegativity constraints
        for j=1:m
            @objective(wc_model, Min, uvals[i,j] + sum(Vvals[i,j,:] .* z))
            optimize!(wc_model)
            if objective_value(wc_model) < -atol # if constraint is violated
                new_z = value.(z)
                @constraint(model, u[i,j] + sum(V[i,j,:] .* new_z) >= 0)
                count += 1
            end
        end
    end
    @info("$(count) nonnegativity cuts added. ")
    ct = count
    for j = 1:m # demand constraints
        dt = zeros(n, m);
        [dt[i,j] = 1 for i=1:n];
        @objective(wc_model, Min, sum(dt .* uvals) + sum(Vvals[:,j,:] * z) - d[j] - (P*z)[j])
        optimize!(wc_model)
        if objective_value(wc_model) < -atol # if constraint is violated
            new_z = value.(z)
            @constraint(model, sum(dt .* u) + sum(V[:,j,:] * new_z) >= d[j] + (P*new_z)[j])
            count += 1
        end
    end
    @info("$(count - ct) demand constraint cuts added.")
    ct = count
    for i = 1:n # capacity constraints
        dt = zeros(n, m);
        [dt[i,j] = 1 for j=1:m];
        @objective(wc_model, Min, - sum(dt .* uvals) - sum(Vvals[i,:,:] * z) + s[i] * xvals[i])
        optimize!(wc_model)
        if objective_value(wc_model) < -atol # if constraint is violated
            new_z = value.(z)
            @constraint(model, sum(dt .* u) + sum(V[i,:,:] * new_z) <= s[i] * x[i])
            count += 1
        end
    end
    @info("$(count - ct) capacity constraint cuts added.")
    ct = count
    # objective
    @objective(wc_model, Max, sum(f[j] * xvals[j] for j = 1:n) + sum(c[i, j] * (uvals[i, j] + sum(Vvals[i,j,:] .* z)) for i=1:n, j=1:m))
    optimize!(wc_model)
    if objective_value(wc_model) > obj_value + atol # if constraint is violated
        new_z = value.(z)
        @constraint(model, objective_function(model) >= sum(f[j] * x[j] for j = 1:n) + sum(c[i, j] * (u[i, j] + sum(V[i,j,:] .* new_z)) for i=1:n, j=1:m))
        count += 1
        @info("Objective cut added.")
    end
    return count
end

rho = 1
Gamm = 5
model, x, u, V = CP_facility_model(c, f)
unset_binary.(x) # First solve a relaxation of the problem.
apply_heuristic(model, x, u, V)
optimize!(model)

# Do 25 iterations of the relaxed problem, and then tighten binary variables.
for i=1:30
    @info("Iteration $(i).")
    xvals, uvals, Vvals = value.(x), value.(u), value.(V)
    count = find_wc_cuts(model, x, u, V, xvals, uvals, Vvals , rho, Gamm)
    if i == 25
        set_binary.(x)
    end
    if count == 0
        @info("Optimum reached.")
        @info("Optimal cost:$(objective_value(model)).")
        break
    end
    optimize!(model)
end

# Problem should converge in ~100 adversarial iterations.

plot_solution(model, x, u)


text:
#= -------------------------------------------------------------------------
# @file solution.jl
#
# @date 12/16/15 22:43:38
# @author Martin Noblia
# @email martin.noblia@openmailbox.org
#
# @brief
#
# @detail
#
  Licence:
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or (at
# your option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# General Public License for more details.
#
# You should have received a copy of the GNU General Public License

---------------------------------------------------------------------------=#

module Credit
"""
Type for encapsulate a test cases
"""
type Cases
    c::Int
    I::Int
    L::Vector{Int}
end

"""
Function for returning the lines of text file
"""
function read_file(str::AbstractString)
    file = open(str)
    lines = eachline(file)
    close(file)
    return lines
end

function lines(io::IO)
    lines_array = Array{UTF8String, 1}()
    for line in eachline(io)
        push!(lines_array, chomp(line))
    end
    return lines_array
end

lines(path::AbstractString) = open(lines, path)

"""
Parse the text file and create the test cases

Input: file_name(ASCIIString)
------
Ouput: cases_vec(Vector{Cases})
------

Example:
--------
`cases_vec = create_cases(file_name)`
"""
function create_cases(file_name::String)

    data = Credit.lines(file_name)
    data_size = length(data)

    re = r"\+?\d+" # regular expresion for find numbers
    number_of_cases = parse(Int, data[1])
    cases_vec = Vector{Credit.Cases}()
    for i = 2:3:data_size - 1
        c = parse(Int, data[i])
        I = parse(Int, data[i+1])
        L = [parse(Int, s)  for s in matchall(re, data[i+2])]
        push!(cases_vec, Credit.Cases(c, I, L))

    end
    return cases_vec
end

"""
Function for find the cases in that the credit is a sum of two items

Input: cases(Array{Credit.Cases, 1})
"""

function filter_cases(cases::Array{Credit.Cases, 1})
   for case in 1:length(cases)
      size = cases[case].I
      s = cases[case].c
      for i in 1:size
         for j in i:size
            if((cases[case].L[i] + cases[case].L[j] == s) && i != j)
               #println("index $i + index $j = credit $s")
               println("Case #$case: $i  $j")
            end
         end
      end
   end
end

end

max_stars_repo_name: elsuizo/CodeJam

max_issues_repo_name: elsuizo/CodeJam

max_forks_repo_name: elsuizo/CodeJam


text:
mutable struct Curriculum
    name::AbstractString        # Name of the curriculum (can be used as an identifier)
    terms::Array{Term}          # Array of terms in order
    courses::Array{Course}      # Array of courses in the curriculum
    numCourses::Int             # Number of courses in curriculum
    creditHours::Int            # Number of credit hours

    complexity::Int             # Sum of all course crucialities
    delay::Int                  # Sum of course delay factors
    blocking::Int               # Sum of course blocking factors
    centrality::Int             # Total Centrality
    reachability::Int           # Total Reachability

    passrate::Float64           # Average passrate of courses
    stopoutModel::Dict          # Dictionary that stores the the model for predicint stopouts

    graph::Any                  # LightGraph representation

    function Curriculum(curriculumName::AbstractString, terms::Array{Term})
        this = new()

        this.name = curriculumName
        this.terms = terms
        this.courses = Course[]

        for (i, term) in enumerate(terms)
            for course in term.courses
                # Add courses to array
                push!(this.courses, course)

                # Add postreq links
                for prereq in course.prereqs
                    push!(prereq.postreqs, course)
                end

                # Add co reqs
                for coreq in course.coreqs
                    push!(coreq.postreqs, course)
                end

                # Mark term the course belongs to
                course.term = i
            end
        end

        for (id, course) in enumerate(this.courses)
            # Compute complexity
            course.id = id
            cruciality(course)
        end

        # Create graph
        this.graph = curriculumGraph(this)

        # Compute centrality
        this.centrality = curriculumCentrality!(this)

        # Compute reachability
        this.reachability = curriculumReachability!(this)

        this.complexity = sum(map(x->x.cruciality, this.courses))
        this.delay = sum(map(x->x.delay, this.courses))
        this.blocking = sum(map(x->x.blocking, this.courses))
        this.numCourses = length(this.courses)
        this.passrate = sum(map(x->x.passrate, this.courses)) / this.numCourses
        this.stopoutModel = Dict()

        this.creditHours = sum(map(x->x.credits, this.courses))

        return this
    end

    function Curriculum(curriculumName::AbstractString, data::Dict)
        # Sort courses by term
        sort!(data["courses"], by = x -> x["term"])

        # Get number of terms
        terms = data["terms"]

        # Array of course arrays, one for each term
        courses = Array{Array{Course}}(undef, terms)
        for i=1:terms
            courses[i] = Course[]
        end

        allCourses = Course[]

        # Create all course objects
        for course in data["courses"]
            name = course["name"]
            credits = course["credits"]
            passrate = course["passrate"]
            haskey(data, "termReq") ? termReq = data["termReq"] : termReq = 0

            c = Course(name, credits, termReq, passrate, Course[], Course[])
            push!(courses[course["term"]], c)
            push!(allCourses, c)
        end

        # Assign pre and coreqs
        for (i, course) in enumerate(data["courses"])
            c = allCourses[i]

            # Prereqs
            prereqs = Course[]
            for prereq in course["prerequisites"]
                ind = findfirst(x -> x.name == prereq, allCourses)
                if ind != 0
                    push!(prereqs, allCourses[ind])
                end
            end

            # Coreqs
            coreqs = Course[]
            for coreq in course["corequisites"]
                ind = findfirst(x -> x.name == coreq, allCourses)
                if ind != 0
                    push!(coreqs, allCourses[ind])
                end
            end

            c.prereqs = prereqs
            c.coreqs = coreqs
        end

        terms = Term[]
        for courseArray in courses
            push!(terms, Term(courseArray))
        end

        this = Curriculum(curriculumName, terms)
        return this
    end

    function Curriculum(curriculumName::AbstractString, path::AbstractString)
        f = open(path)
        data = JSON.parse(read(f, String))
        close(f)

        this = Curriculum(curriculumName, data)
        return this
    end
end

function curriculumGraph(curriculum)
    numCourses = length(curriculum.courses)

    graph = DiGraph(numCourses)

    for course in curriculum.courses
        for prereq in course.prereqs
            add_edge!(graph, prereq.id, course.id)
        end
    end

    return graph
end

function curriculumCentrality!(curriculum)
    measures = zeros(length(curriculum.courses))

    adj = adjacency_matrix(curriculum.graph)

    for i=1:length(curriculum.courses)
        for j=1:length(curriculum.courses)
            paths = path_enumeration(i, j, adj)
            all_nodes = []
            for p in paths
                all_nodes = vcat(all_nodes, p.path)
            end

            for v in unique(all_nodes)
                measures[v] += 1
            end
        end
    end

    for (i, measure) in enumerate(measures)
        curriculum.courses[i].centrality = measure
    end

    return sum(measures)
end

function curriculumReachability!(curriculum)
    newC = deepcopy(curriculum)

    for course in newC.courses
        course.prereqs = Course[]
        course.postreqs = Course[]
    end

    for (ci, course) in enumerate(curriculum.courses)
        for prereq in course.prereqs
            push!(newC.courses[ci].postreqs, newC.courses[prereq.id])
        end
    end

    for (i, course) in enumerate(newC.courses)
        cruciality(course)
        curriculum.courses[i].reachability = course.blocking
    end

    return sum(map(x->x.blocking, newC.courses))
end

max_stars_repo_name: M-S-H/CASL

max_issues_repo_name: M-S-H/CASL

max_forks_repo_name: M-S-H/CASL


text:
# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# Parse a scene

"""
    parse_scene(stream::InputStream, scene::Scene = Scene())

Return the [`Scene`](@ref) instance resulting parsing the SceneLang script associated with the given [`InputStream`](@ref).
"""
function parse_scene(stream::InputStream, scene::Scene = Scene())
    while !eof(stream)
        command_token = expect_command(stream, (USING, SET, UNSET, SPAWN, DUMP))
        unread_token(stream, command_token)
        command = command_token.value
        if command == USING
            parse_using_command(stream, scene)
        elseif command == SET
            parse_set_command(stream, scene)
        elseif command == UNSET
            parse_unset_command(stream, scene)
        elseif command == SPAWN
            parse_spawn_command(stream, scene)
        elseif command == DUMP
            parse_dump_command(stream, scene)
        else
            @assert false "Got unparsable command '$command' from expect_command."
        end
    end
    scene
end

function print_subsequent_lexer_exceptions(stream::InputStream, except::InterpreterException)
    is_first_exception = true
    while true
        try
            iterate(stream) |> isnothing && break
        catch e
            if is_first_exception
                is_first_exception = false
                printstyled(stderr, "### OTHER LEXER ERRORS ###\n\n", color = :magenta)
            end
            isa(e, InterpreterException) || rethrow(e)
            showerror(stderr, e)
            println(stderr)
        end
    end
    if !is_first_exception
        println(stderr)
        printstyled(stderr, "### MAIN PARSER ERROR ###\n\n", color = :magenta)
    end
    except
end


text:
function get_segment_id(distances, cumulative_length)
    segment_id = []
    for i in distances
        index = findfirst(x -> x >= i, cumulative_length)
        push!(segment_id, index)
    end
    return segment_id
end

function rotate_by(x, theta)
    rx = cos(theta) * x[1] - sin(theta) * x[2]
    ry = sin(theta) * x[1] + cos(theta) * x[2]
    return [rx, ry]
end

function compute_position(distance, fr, to, sensor_distance)
    dx = to[1] - fr[1]
    dy = to[2] - fr[2]
    theta = atan(dy, dx)
    D = sqrt(dx * dx + dy * dy)
    l_s = [distance, sensor_distance]
    r_s = [distance, -sensor_distance]
    l_r = rotate_by(l_s, theta)
    r_r = rotate_by(r_s, theta)
    return (l_r + collect(fr), r_r + collect(fr))
end

function get_straight_placements(td::TunnelData, drop_distance)
    num_sensors = 2
    distances = Vector()
    l_positions = Vector()
    r_positions = Vector()
    segment_id = Dict()
    num_sensors_in_segment = []

    push!(distances, 0.0)

    while true
        last_placement_distance = distances[end]
        if (last_placement_distance + drop_distance > td.c_distances[end])
            break
        end
        push!(distances, last_placement_distance + drop_distance)
        num_sensors += 2
    end

    for i in 1:length(distances)
        d = distances[i]
        index = findfirst(x -> d < x, td.c_distances)
        segment_id[i] = index
    end

    for i in 1:td.num_segments
        count = length([k for (k, v) in segment_id if v == i])
        push!(num_sensors_in_segment, count)
    end

    for i in 1:length(distances)
        d = distances[i]
        seg_id = segment_id[i]
        (seg_id > 1) && (d -= td.c_distances[seg_id - 1])
        fr = td.wp[seg_id]
        to = td.wp[seg_id + 1]
        l_pos, r_pos = compute_position(d, fr, to, 10.0)
        push!(l_positions, l_pos)
        push!(r_positions, r_pos)
    end

    return Placement(
        num_sensors, distances,
        l_positions, r_positions,
        segment_id, num_sensors_in_segment,
        Vector()
    )
end

function get_over_placements(td::TunnelData,
    sp::Placement, drop_distance)::Placement
    num_sensors = 2
    distances = Vector()
    l_positions = Vector()
    r_positions = Vector()
    segment_id = Dict()
    num_sensors_in_segment = []

    num_sensors = deepcopy(sp.num_sensors)
    distances = deepcopy(sp.distances)

    for i in 1:(length(td.c_distances)-1)
        turn_at = td.c_distances[i]
        if !(turn_at in distances)
            push!(distances, turn_at)
            num_sensors += 2
        end
    end

    distances = sort(distances)

    for i in 1:length(distances)
        d = distances[i]
        index = findfirst(x -> d < x, td.c_distances)
        segment_id[i] = index
    end

    for i in 1:td.num_segments
        count = length([k for (k, v) in segment_id if v == i])
        push!(num_sensors_in_segment, count)
    end

    for i in 1:length(distances)
        d = distances[i]
        seg_id = segment_id[i]
        (seg_id > 1) && (d -= td.c_distances[seg_id - 1])
        fr = td.wp[seg_id]
        to = td.wp[seg_id + 1]
        l_pos, r_pos = compute_position(d, fr, to, 10.0)
        push!(l_positions, l_pos)
        push!(r_positions, r_pos)
    end

    return Placement(
        num_sensors, distances,
        l_positions, r_positions,
        segment_id, num_sensors_in_segment,
        Vector()
    )
end

function update_placement!(p, td)
    num_sensors = length(p.distances) * 2
    distances = p.distances
    p.l_positions = Vector()
    p.r_positions = Vector()
    p.segment_id = Dict()
    p.num_sensors_in_segment = []

    for i in 1:length(distances)
        d = distances[i]
        index = findfirst(x -> d < x, td.c_distances)
        p.segment_id[i] = index
    end

    for i in 1:td.num_segments
        count = length([k for (k, v) in p.segment_id if v == i])
        push!(p.num_sensors_in_segment, count)
    end

    for i in 1:length(distances)
        d = distances[i]
        seg_id = p.segment_id[i]
        (seg_id > 1) && (d -= td.c_distances[seg_id - 1])
        fr = td.wp[seg_id]
        to = td.wp[seg_id + 1]
        l_pos, r_pos = compute_position(d, fr, to, 10.0)
        push!(p.l_positions, l_pos)
        push!(p.r_positions, r_pos)
    end
end

function compute_last_sensor(tunnel_data, placement, segment_id)
    distances = []
    sensor_ids = []
    for i in 1:Int(placement.num_sensors/2)
        d = placement.distances[i]
        id = findfirst(x -> x > d, tunnel_data.c_distances)
        if (id == segment_id)
            push!(distances, d)
            push!(sensor_ids, i)
        end
    end
    max_element = findmax(distances)
    return max_element[1], sensor_ids[max_element[2]]
end

function get_segment_ids(td::TunnelData, segment_id,
    fr::Float64, to_turn::Float64, drop_distance)

    segment_ids = []
    additional_distance = drop_distance - to_turn
    for i in (segment_id+1):td.num_segments
        if (td.c_distances[i] > drop_distance)
            push!(segment_ids, i)
            break
        else
            push!(segment_ids, i)
        end
    end
    return segment_ids
end

function get_los_point(s_l, s_r, w_l, w_r, seg_fr, seg_to)

    r_ray = (w_r - s_r)./dist(w_l - s_l, zeros(2))
    l_ray = (w_l - s_l)./dist(w_l - s_l, zeros(2))

    m_r = JuMP.Model()
    m_l = JuMP.Model()

    @variable(m_r, mu_r, lower_bound = 0.0)
    @variable(m_r, 0.0 <= lambda_r <= 1.0)

    @variable(m_l, mu_l, lower_bound = 0.0)
    @variable(m_l, 0.0 <= lambda_l <= 1.0)

    @constraint(m_r, s_r + mu_r * r_ray .==
        lambda_r * seg_fr + (1-lambda_r) * seg_to)
    @constraint(m_l, s_l + mu_l * l_ray .==
        lambda_l * seg_fr + (1-lambda_l) * seg_to)

    set_optimizer(m_r, CPLEX.Optimizer)
    set_silent(m_r)
    set_optimizer(m_l, CPLEX.Optimizer)
    set_silent(m_l)

    optimize!(m_l)
    optimize!(m_r)

    p_l = [NaN, NaN]
    d_l = NaN
    p_r = [NaN, NaN]
    d_r = NaN

    if (termination_status(m_l) != MOI.INFEASIBLE)
        lambda = JuMP.value(lambda_l)
        p_l = lambda * seg_fr + (1-lambda) * seg_to
        d_l = dist(p_l, seg_fr)
    end
    if (termination_status(m_r) != MOI.INFEASIBLE)
        lambda = JuMP.value(lambda_r)
        p_r = lambda * seg_fr + (1-lambda) * seg_to
        d_r = dist(p_r, seg_fr)
    end

    (isnan(d_r) && isnan(d_l)) && (return nothing)
    (isnan(d_r)) && (return p_l)
    (isnan(d_l)) && (return p_r)
    (d_r < d_l) && (return p_r)
    (d_r > d_l) && (return p_l)

end

function get_optimal_placements(td::TunnelData,
    sp::Placement, drop_distance)::Placement

    op = deepcopy(sp)

    for i in 1:(td.num_segments-1)
        lw_wp = td.lw_wp[i+1]
        rw_wp = td.rw_wp[i+1]
        (op.num_sensors_in_segment[i] == 0) && (continue)
        prev_sensor_distance, prev_sensor_id = compute_last_sensor(td, op, i)
        distance_to_turn = td.c_distances[i] - prev_sensor_distance
        (prev_sensor_id * 2 == op.num_sensors) && (break)
        next_sensor_distance = op.distances[prev_sensor_id + 1]
        next_sensor_id = prev_sensor_id + 1
        prev_sensor_l_position = op.l_positions[prev_sensor_id]
        prev_sensor_r_position = op.r_positions[prev_sensor_id]
        inter_sensor_distance = drop_distance
        segments_to_consider = get_segment_ids(
            td, i, prev_sensor_distance, distance_to_turn, drop_distance)
        for s in segments_to_consider
            fr = td.wp[s]
            to = td.wp[s+1]
            point = get_los_point(prev_sensor_l_position,
                prev_sensor_r_position,
                lw_wp, rw_wp, fr, to)
            (isa(point, Nothing)) && (continue)
            d = dist(point, fr)
            point_distance = td.c_distances[s-1] + d
            push!(op.los_points, point)
            (point_distance > next_sensor_distance) && (continue)
            offset = next_sensor_distance - point_distance
            for id in next_sensor_id:Int(op.num_sensors/2)
                op.distances[id] -= offset
            end
            update_placement!(op, td)
        end
    end

    while true
        last_placement_distance = op.distances[end]
        if (last_placement_distance + drop_distance > td.c_distances[end])
            break
        end
        push!(op.distances, last_placement_distance + drop_distance)
        op.num_sensors += 2
    end

    update_placement!(op, td)

    return op

end


text:
module NeuralQuantum

# Packages required to construct the module
using Reexport: @reexport
using Requires: @require
using MacroTools: @forward

# Packages required to work-around performance issues
using SpecializeVarargs
using UnsafeArrays

# Packages to support (mostly) GPUs
using GPUArrays
using Adapt: adapt

# Packages for Distributed/Parallelization
using MPI
include("External/TPI/TPI.jl")
using .TPI

# Math stuff
using LinearMaps: LinearMap
using Strided # for gradient of batches of networks
using NNlib
using Zygote: Zygote, Params, @adjoint, pullback

# Standard library
using Random: Random, AbstractRNG, MersenneTwister, GLOBAL_RNG, rand!
using Random: shuffle!, randn!
using LinearAlgebra
using SparseArrays
using Statistics
using Printf

@reexport using QuantumOpticsBase

# Support for Colored Graphs and Lattices
@reexport using LightGraphs
include("External/ColoredGraphs/ColoredGraphs.jl")
using .ColoredGraphs
export HyperCube, translational_symm_table

# Iterative Solvers and custom minres solver
using IterativeSolvers: minres, lsqr, cg
include("External/IterativeSolvers/minresqlp.jl")
using .MinresQlp

# Optimisers, that will be split in a separate package at some point
include("Optimisers/Optimisers.jl")
using .Optimisers
import .Optimisers: update, update!
export Optimisers

# Exports

# Standard Precision used
const STD_REAL_PREC = Float32

# Abstract Types
abstract type NeuralNetwork end
abstract type Sampler end
abstract type AbstractAccumulator end

# Type describing the parallel backend used by a solver.
abstract type ParallelType end
struct NotParallel <: ParallelType end
struct ParallelThreaded <: ParallelType end
struct ParallelMPI <: ParallelType end
export NotParallel, ParallelThreaded

# Prallelization
include("Parallel/base_parallel.jl")
include("Parallel/not_parallel.jl")

# Various utility functions
include("utils/math.jl")
include("utils/views.jl")
include("utils/stats.jl")
include("utils/rng.jl")

# Basic states
include("States/States.jl")
include("States/Iterators.jl")
include("States/StateChanges.jl")
include("States/ApplyStateChanges.jl")

# Hilbert spaces
abstract type AbstractHilbert end
include("Hilbert/base_basis.jl")
include("Hilbert/iterator.jl")
include("Hilbert/DiscreteHilbert.jl")
include("Hilbert/HomogeneousFock.jl")
include("Hilbert/HomogeneousSpin.jl")
include("Hilbert/DoubledHilbert.jl")
include("Hilbert/SuperHilbert.jl")
include("Hilbert/basis_convert.jl")

# Nets
include("base_networks.jl")
include("base_cached_networks.jl")
include("functor.jl") # from flux
include("structuring.jl")
include("tuple_logic.jl")

#AD
include("AD/base_derivatives.jl")
include("AD/RealDerivatives.jl")
include("base_batched_networks.jl")

# Linear Operators
include("Operators/BaseOperators.jl")
# Connections
include("Operators/OpConnections/OpConnection.jl")
include("Operators/OpConnections/OpConnectionIdentity.jl")
include("Operators/OpConnections/OpConnectionTensor.jl")
include("Operators/OpConnections/SuperOpConnection.jl")

include("Operators/OpConnectionIndex.jl")

include("Operators/Operators/KLocalOperator.jl")
include("Operators/Operators/KLocalIdentity.jl")
include("Operators/Operators/KLocalOperatorSum.jl")
include("Operators/Operators/KLocalOperatorTensor.jl")
include("Operators/Operators/KLocalZero.jl")
include("Operators/Operators/KLocalLiouvillian.jl")

include("Operators/SimpleOperators.jl")
include("Operators/OpConversion.jl")

export OpConnection
export KLocalOperator, KLocalOperatorTensor, KLocalOperatorRow, operators
export row_valdiff, col_valdiff, sites, conn_type

# Neural Networks
include("Networks/utils.jl")
include("Networks/activation.jl")

# Mixed Density Matrices
include("Networks/MixedDensityMatrix/NDM.jl")
include("Networks/MixedDensityMatrix/NDMBatched.jl")
include("Networks/MixedDensityMatrix/NDMComplex.jl")
include("Networks/MixedDensityMatrix/NDMSymm.jl")
include("Networks/MixedDensityMatrix/NDMSymmBatched.jl")
include("Networks/MixedDensityMatrix/RBMSplit.jl")
include("Networks/MixedDensityMatrix/RBMSplitBatched.jl")

# Closed Systems
include("Networks/ClosedSystems/RBM.jl")
include("Networks/ClosedSystems/RBMBatched.jl")

# FFNN
include("Networks/Chains/Chain.jl")
include("Networks/Chains/conv.jl")
include("Networks/Chains/NQConv.jl")

include("Networks/Chains/ChainBatched.jl")

# Layers
include("Networks/Chains/Layers/Dense.jl")
include("Networks/Chains/Layers/WeightedSum.jl")
include("Networks/Chains/Layers/sum.jl")
include("Networks/Chains/Layers/DenseSplit.jl")

include("Networks/Chains/Layers/DenseBatched.jl")
include("Networks/Chains/Layers/WeightedSumBatched.jl")
include("Networks/Chains/Layers/sumBatched.jl")
include("Networks/Chains/Layers/DenseSplitBatched.jl")

include("Networks/Chains/Layers/PositiveDefR.jl")
include("Networks/Chains/Layers/PositiveDefRBatched.jl")

# Wrappers
include("Networks/ClosedSystems/PureChainWrapper.jl")
include("Networks/MixedDensityMatrix/MixedChainWrapper.jl")

# gen state
export state
include("generate_state.jl")

# Algorithms
abstract type Algorithm end
abstract type AlgorithmCache end
export precondition!

# SR
include("Algorithms/SR/SR.jl")
# Gradient
include("Algorithms/Gradient/Gradient.jl")

include("Algorithms/batched_algorithms.jl")
include("Algorithms/SR/SR_notfull.jl")
include("Algorithms/SR/SRDirect.jl")
include("Algorithms/SR/SRIterative.jl")
include("Algorithms/Gradient/Gradient_batched.jl")

# Sampling
include("Samplers/base_samplers.jl")
export cache, init_sampler!, done, samplenext!
export get_sampler, sampler_list, multithread

include("Samplers/Exact.jl")
include("Samplers/Metropolis.jl")
include("Samplers/MCMCRules/LocalRule.jl")
include("Samplers/MCMCRules/ExchangeRule.jl")
include("Samplers/MCMCRules/OperatorRule.jl")
include("Samplers/MCMCRules/Nagy.jl")

# other
include("utils/densitymatrix.jl")
#include("utils/expectation_values.jl")
include("utils/translational_symm.jl")
#include("utils/logging.jl")
include("utils/loading.jl")
# interface
include("IterativeInterface/Samplers/BaseIterativeSampler.jl")

include("IterativeInterface/Accumulators/base_accumulators.jl")
include("IterativeInterface/Accumulators/AccumulatorLogPsi.jl")
include("IterativeInterface/Accumulators/AccumulatorLogGradPsi.jl")
include("IterativeInterface/Accumulators/AccumulatorObsScalar.jl")
include("IterativeInterface/Accumulators/AccumulatorObsGrad.jl")

include("IterativeInterface/Samplers/SimpleIterativeSampler.jl")
include("IterativeInterface/Samplers/CostFun/BatchedGradSampler.jl")
include("IterativeInterface/Samplers/CostFun/BatchedValSampler.jl")
include("IterativeInterface/Samplers/Obs/BatchedObsDMSampler.jl")
include("IterativeInterface/Samplers/Obs/BatchedObsKetSampler.jl")
include("IterativeInterface/build_Batched.jl")

export sample!, add_observable!, compute_observables

include("utils/num_grad.jl")

# Parallelization types
include("Parallel/MPI/mpi.jl")
include("Parallel/Threads/threads_base.jl")
include("Parallel/Threads/threads_wrappers.jl")

# gpu stuff
include("GPU/Base/GPUstates.jl")
include("GPU/Base/JLstates.jl")
include("GPU/IterativeInterface/AccumulatorLogPsi.jl")
include("GPU/IterativeInterface/AccumulatorObsScalar.jl")
include("GPU/Samplers/LocalRuleGPU.jl")

include("GPU/gpu.jl")

# Extra
include("Extra/RecipesBase.jl")
include("Extra/ValueHistories.jl")

function __init__()
  @require TensorBoardLogger="899adc3e-224a-11e9-021f-63837185c80f" begin
    include("Extra/Optionals/TensorBoardLogger.jl")
  end

  @require Measurements="eff96d63-e80a-5855-80a2-b1b0885c5ab7" begin
    include("Extra/Optionals/Measurements.jl")
  end
end

end # module


text:
struct Transform{F <: Function}
    f::F
end

(t::Transform)(x::AV{<:Real}) = t.f.(x)
(t::Transform)(x::ColVecs) = ColVecs(hcat(t.f.(x)...)) # In slow method error :<
# (t::Transform)(x::ColVecs) = hcat(t.f.(x)...) # In slow method
# (t::Transform)(x::ColVecs) = ColVecs(mapslices(t.f, x.X, dims=1))

#=
mapslices implementation make no sense since it finally behaves like "slow" method if
    it follows "AD style" to get gradient.

So we just disable ColVecs's error
=#

Zygote.@adjoint function ColVecs(X::AbstractMatrix)
    back(Δ::NamedTuple) = (Δ.X,)
    back(Δ::AbstractMatrix) = (Δ,)
    function back(Δ::AbstractVector{<:AbstractVector{<:Real}})
        (hcat(Δ...),)
    end
    return ColVecs(X), back
end

"""
    LinearWithTransform{T<:Real} <: BaseKernel
The standardised linear kernel / dot-product kernel with an extra transform.
"""
struct LinearWithTransform{T <: Transform} <: Stheno.BaseKernel
    transform::T
end

# Binary methods
function ew(k::LinearWithTransform, x::AV{<:Real}, x′::AV{<:Real})
    x = k.transform(x)
    x′ = k.transform(x′)
    x .* x′
end
function pw(k::LinearWithTransform, x::AV{<:Real}, x′::AV{<:Real})
    x = k.transform(x)
    x′ = k.transform(x′)
    x .* x′'
end
function ew(k::LinearWithTransform, x::ColVecs, x′::ColVecs)
    x = k.transform(x)
    x′ = k.transform(x′)
    reshape(sum(x.X .* x′.X; dims=1), :)
end
function pw(k::LinearWithTransform, x::ColVecs, x′::ColVecs)
    x = k.transform(x)
    x′ = k.transform(x′)
    x.X' * x′.X
end

# Unary methods
function ew(k::LinearWithTransform, x::AV{<:Real})
    x = k.transform(x)
    x.^2
end
function pw(k::LinearWithTransform, x::AV{<:Real})
    x = k.transform(x)
    x .* x'
end
function ew(k::LinearWithTransform, x::ColVecs)
    x = k.transform(x)
    reshape(sum(abs2.(x.X); dims=1), :)
end
function pw(k::LinearWithTransform, x::ColVecs)
    x = k.transform(x)
    x.X' * x.X
end


text:
function test_item_knn()
    println("-- Testing ItemKNN rating prediction")

    m = [NaN 3 NaN 1 2 1 NaN 4
         1 2 NaN NaN 3 2 NaN 3
         NaN 2 3 3 NaN 5 NaN 1]
    data = DataAccessor(m)

    recommender = ItemKNN(data, 1)
    build!(recommender)
    @test isapprox(recommender.sim[1, 2], 0.485, atol=1e-3)
    @test isapprox(recommender.sim[5, 6], 0.405, atol=1e-3)

    rec = recommend(recommender, 1, 4, [i for i in 1:8])
    @test rec[1] == (8 => 4.0)
    @test rec[2] == (2 => 3.0)
    @test rec[3] == (5 => 2.0)
    @test rec[4] == (4 => 1.0)

    recommender = ItemKNN(data, 1)
    build!(recommender, adjusted_cosine=true)
    @test isapprox(recommender.sim[1, 2], 0.174, atol=1e-3)
    @test isapprox(recommender.sim[5, 6], 0.039, atol=1e-3)

    rec = recommend(recommender, 1, 4, [i for i in 1:8])
    @test rec[1] == (8 => 4.0)
    @test rec[2] == (2 => 3.0)
    @test rec[3] == (5 => 2.0)
    @test rec[4] == (4 => 1.0)
end

test_item_knn()


text:
using MyPkgDemoEn
using Test

@testset "MyPkgDemoEn.jl" begin
    # Write your tests here.
    @test MyPkgDemoEn.add_nums(1,3) == 4
    @test MyPkgDemoEn.add_nums(2,3) == 5
    @test MyPkgDemoEn.add_nums(3,4) == 7
    @test MyPkgDemoEn.add_nums(3,4,1) == 8
    @test MyPkgDemoEn.add_nums(3,4,2) == 9
    @test MyPkgDemoEn.add_nums(2,4,1) == 7
end


text:
function web_render(cur_observer::Observable, cur_grid::AbstractGrid)
  cur_table = make_table(cur_grid.rows, cur_grid.cols, cur_grid.hidden_rows)

  cur_events = Dict()

  cur_events["dblclick"] = @js function (cur_event)
    $cur_observer[] = "play";
  end

  cur_events["blur"] = @js function (cur_event)
    $cur_observer[] = "suspend";
  end

  cur_events["keydown"] = @js function (cur_event)

    cur_key_code = cur_event.keyCode;

    is_action_key = false;

    if cur_event.target.tagName.toLowerCase() !== "input"

      if ( cur_key_code == 77 ) ; is_action_key = true ; $cur_observer[] = "music" ; end

      if ( cur_key_code == 80 ) ; is_action_key = true ; $cur_observer[] = "pause" ; end

      if ( cur_key_code == 37 ) ; is_action_key = true ; $cur_observer[] = "slide_left" ; end
      if ( cur_key_code == 38 ) ; is_action_key = true ; $cur_observer[] = "clockwise" ; end
      if ( cur_key_code == 39 ) ; is_action_key = true ; $cur_observer[] = "slide_right" ; end
      if ( cur_key_code == 40 ) ; is_action_key = true ; $cur_observer[] = "soft_drop" ; end

      if ( cur_key_code == 32 ) ; is_action_key = true ; $cur_observer[] = "hard_drop" ; end

      if ( cur_key_code == 67 ) ; is_action_key = true ; $cur_observer[] = "hold" ; end
      if ( cur_key_code == 90 ) ; is_action_key = true ; $cur_observer[] = "counter" ; end

    end

    if is_action_key ; cur_event.preventDefault() ; end

  end

  cur_events["keyup"] = @js function (cur_event)

    cur_key_code = cur_event.keyCode;

    is_action_key = false;

    if cur_event.target.tagName.toLowerCase() !== "input"

      if ( cur_key_code == 77 ) ; is_action_key = true ; $cur_observer[] = "" ; end

      if ( cur_key_code == 80 ) ; is_action_key = true ; $cur_observer[] = "" ; end

      if ( cur_key_code == 37 ) ; is_action_key = true ; $cur_observer[] = "" ; end
      if ( cur_key_code == 38 ) ; is_action_key = true ; $cur_observer[] = "" ; end
      if ( cur_key_code == 39 ) ; is_action_key = true ; $cur_observer[] = "" ; end
      if ( cur_key_code == 40 ) ; is_action_key = true ; $cur_observer[] = "" ; end

      if ( cur_key_code == 32 ) ; is_action_key = true ; $cur_observer[] = "" ; end

      if ( cur_key_code == 67 ) ; is_action_key = true ; $cur_observer[] = "" ; end
      if ( cur_key_code == 90 ) ; is_action_key = true ; $cur_observer[] = "" ; end

    end

    if is_action_key ; cur_event.preventDefault() ; end

  end

  cur_events["keypress"] = @js function (cur_event)

    if cur_event.target.tagName.toLowerCase() !== "input"

      cur_event.preventDefault();

    end

  end

  cur_preview_tables = []

  for cur_preview in 1:cur_grid.player.bag.previews
    push!(
      cur_preview_tables,
      make_table(2, 4, class_names="cs-preview-piece js-preview-piece js-preview-piece__$(cur_preview)")
    )
  end

  cur_side_bar = Node(
    :div,
    dom"h1"(
      "Level 1",
      attributes=Dict(:class => "cs-level-text js-level-text")
    ),
    make_table(2, 4, class_names="cs-hold-piece js-hold-piece"),
    cur_preview_tables...,
    dom"h2"(
      "00000000",
      attributes=Dict(:class => "cs-score-text js-score-text")
    ),
    attributes=Dict(:class => "col-sm-5 cs-col hidden-xs cs-side-bar")
  )

  cur_container = dom"div[tabindex=1]"(
    dom"div"(
      dom"div"(
        cur_table,
        attributes=Dict(:class => "col-sm-7 cs-col cs-main-area")
      ),
      cur_side_bar,
      attributes=Dict(:class => "row")
    ),
    dom"div"(
      dom"div"(
        dom"h1"("tetris."),
        dom"div"("double-click to play."),
        attributes=Dict(:class => "cs-splash-text js-splash-text")
      ),
      dom"div"(
        dom"div"("⬆⬆⬇⬇⬅➡⬅➡"),
        dom"div"("c. z. space bar"),
        attributes=Dict(:class => "cs-how-to")
      ),
      attributes=Dict(:class => "cs-overlay-screen cs-splash-screen js-splash-screen cs-blue-grey")
    ),
    dom"div"(
      make_buttons(cur_observer),
      attributes=Dict(:class => "cs-overlay-screen")
    ),
    dom"div"(
      dom"div"(
        make_scores(cur_grid.player.game.board, true),
        make_scores(cur_grid.player.game.board, false),
        attributes=Dict(:class => "cs-jumbotron jumbotron")
      ),
      attributes=Dict(:class => "cs-overlay-screen cs-scores-screen $(isempty(cur_grid.player.game.board.entries) ? "hidden" : "")")
    ),
    dom"div"(
      dom"div"(
        dom"div"(
          dom"div"(
            dom"h1"(
              "New High Score!",
            ),
            dom"h2"(
              lpad(0, 8, "0")
            ),
            dom"form"(
              dom"input"(
                attributes=Dict(
                  :class => "form-control",
                  :type => "text",
                  :placeholder => "________",
                  :maxlength => 8,
                  :method => "post",
                  :required => "true"
                )
              ),
              events=Dict(
                "submit" => @js function (cur_event)
                  curName = document.querySelector(".cs-new-score input").value;

                  $cur_observer[] = "submit-" + curName;
                  cur_event.preventDefault();
                end
              ),
              attributes=Dict(:class => "js-form")
            ),
            attributes=Dict(:class => "cs-new-score")
          ),
          attributes=Dict(:class => "row")
        ),
        attributes=Dict(:class => "cs-jumbotron jumbotron")
      ),
      attributes=Dict(:class => "cs-overlay-screen cs-scores-screen js-new-score hidden")
    ),
    dom"div"(
      attributes=Dict(:class => "cs-offline js-offline hidden", :title => "You are offline.")
    ),
    attributes=Dict(:class => "js-tetris-container cs-tetris-container container-fluid cs-has-splash"),
    events=cur_events
  )

  cur_container
end


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

using ParallelAccelerator
using Base.Test

include("lib.jl")
include("parfor.jl")
include("mapreduce.jl")
include("abs.jl")
include("const_promote.jl")
include("rand.jl")
include("BitArray.jl")
include("range.jl")
include("seq.jl")
include("cat.jl")
include("hcat.jl")
include("vcat.jl")
include("ranges.jl")
include("misc.jl")
include("aug_assign.jl")
include("complex.jl")
include("print.jl")
include("strings.jl")
include("test_lr.jl")
include("test_kmeans.jl")
include("gemv_test.jl")
include("transpose_test.jl")
include("vecnorm_test.jl")
include("broadcast.jl")

# Examples.  We're not including them all here, because it would take
# too long, but just including a few seems like a good compromise that
# exercises much of ParallelAccelerator.

module TestBlackScholes
using Base.Test
include("../examples/black-scholes/black-scholes.jl")
# black-scholes should have only 1 allocation and 1 parfor after optimization
if !(ParallelAccelerator.getPseMode() == ParallelAccelerator.THREADS_MODE)
@test ParallelAccelerator.get_num_acc_allocs()==1 && ParallelAccelerator.get_num_acc_parfors()==1
end
end

module TestPi
using Base.Test
include("../examples/pi/pi.jl")
# pi should have no allocation and 1 parfor after optimization
if !(ParallelAccelerator.getPseMode() == ParallelAccelerator.THREADS_MODE)
@test ParallelAccelerator.get_num_acc_allocs()==0 && ParallelAccelerator.get_num_acc_parfors()==1
end
end

module TestOptFlow
using Base.Test
include("../examples/opt-flow/opt-flow.jl")
end

module TestKMeans
using Base.Test
include("../examples/k-means/k-means.jl")
# k-means should have 3 top-level parfors after optimization
if !(ParallelAccelerator.getPseMode() == ParallelAccelerator.THREADS_MODE)
@test ParallelAccelerator.get_num_acc_parfors()==3
end
end

module TestLR
using Base.Test
include("../examples/logistic_regression/logistic_regression.jl")
# logistic_regression should have 5 top-level parfors after optimization
if !(ParallelAccelerator.getPseMode() == ParallelAccelerator.THREADS_MODE)
@test ParallelAccelerator.get_num_acc_parfors()==5
end
end

# Delete file left behind by opt-flow.
dir = pwd()
img_file = joinpath(dir, "out.flo")
rm(img_file)


text:
#=

  From http://rosettacode.org/wiki/Permutations/Derangements
  """
  Permutations/Derangements
  A derangement is a permutation of the order of distinct items in which
  no item appears in its original place.

  For example, the only two derangements of the three items
    (0, 1, 2) are (1, 2, 0), and (2, 0, 1).

  The number of derangements of n distinct items is known as the subfactorial of n,
  sometimes written as !n. There are various ways to calculate !n.
  ...
  """

  The theoretical value is 1/exp(1): 0.36787944117144233

  n: 2
  Distributions of variable total0 (num:0)
  0.00000 =>    5116  (0.511600)
  1.00000 =>    4884  (0.488400)

  n: 3
  Distributions of variable total0 (num:0)
  0.00000 =>    6821  (0.682100)
  1.00000 =>    3179  (0.317900)

  n: 4
  Distributions of variable total0 (num:0)
  0.00000 =>    6478  (0.647800)
  1.00000 =>    3522  (0.352200)

  n: 5
  Distributions of variable total0 (num:0)
  0.00000 =>    6134  (0.613400)
  1.00000 =>    3866  (0.386600)

  Cf ~/webppl/derangements.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function derangements(n=3)

    p ~ filldist(DiscreteUniform(1,n),n)

    #  Ensure unicity.
    for i in 1:n
        for j in i+1:n
            true ~ Dirac(p[i] != p[j])
        end
    end

    #  How many are in the i'th position
    total ~ Dirac(sum([i == p[i] for i in 1:n]))

    total0 ~ Dirac(total == 0)
    return p
end

function run_derangements(n)

    println("n: ", n, " theoretical: ", 1/exp(1))
    model = derangements(n)

    num_chains = 4

    # chns = sample(model, Prior(), 10_000)
    # chns = sample(model, MH(), 10_000)
    # chns = sample(model, PG(5), 10_000)
    chns = sample(model, SMC(10), 10_000)
    # chns = sample(model, IS(), 100_000)

    # display(chns)
    # display(plot(chns))

    show_var_dist_pct(chns, :total0)
    # for i in 1:n
    #     show_var_dist_pct(chns, Symbol("p[$i]"))
    # end

    # chains_params = Turing.MCMCChains.get_sections(chns, :parameters)
    # genq = generated_quantities(model, chains_params)
    # show_var_dist_pct(genq)

end

for n in 2:5
    run_derangements(n)
    println()
end

# run_derangements(7)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using PDELib
using Triangulate

function fv_laplace_circle(;Plotter=nothing,nref=0)

    function circle!(builder, center, radius; n=20*2^nref)
        points=[point!(builder, center[1]+radius*sin(t),center[2]+radius*cos(t)) for t in range(0,2π,length=n)]
        for i=1:n-1
            facet!(builder,points[i],points[i+1])
        end
        facet!(builder,points[end],points[1])
    end

    builder=SimplexGridBuilder(Generator=Triangulate)
    cellregion!(builder,1)
    maxvolume!(builder,0.1)
    regionpoint!(builder,0,0)

    facetregion!(builder,1)
    circle!(builder,(0,0),1)

    facetregion!(builder,2)
    circle!(builder,(-0.2,-0.2),0.2)
    holepoint!(builder,-0.2,-0.2)

    grid=simplexgrid(builder,maxvolume=0.1*4.0^(-nref))

    function flux(f,u0,edge)
        u=unknowns(edge,u0)
        f[1]=u[1,1]-u[1,2]
    end
    physics=VoronoiFVM.Physics(num_species=1,flux=flux)

    sys=VoronoiFVM.System(grid,physics)
    ispec=1
    enable_species!(sys,ispec,[1])
    boundary_dirichlet!(sys,ispec,1,0.0)
    boundary_dirichlet!(sys,ispec,2,1.0)
    inival=unknowns(sys,inival=0)
    solution=unknowns(sys)
    VoronoiFVM.solve!(solution,inival,sys)

    visualizer=GridVisualizer(Plotter=Plotter,layout=(1,2),resolution=(800,400))
    gridplot!(visualizer[1,1],grid)
    scalarplot!(visualizer[1,2],grid,solution[1,:])
    reveal(visualizer)
    true
end


text:
using Bridge, StaticArrays
using Test, Statistics, LinearAlgebra, Random
SV2 = SVector{2,Float64}

@test let # functions callable
    tt = range(0, stop=1, length=100)

    sample(tt, Wiener{SV2}())
    sample!(SamplePath{SV2}(collect(tt), zeros(SV2, 100)), Wiener{SV2}())
    sample(tt, WienerBridge{SV2}(1.0, SVector(2.0, -2.0)))
    sample!(SamplePath{SV2}(collect(tt), zeros(SV2, 100)), WienerBridge{SV2}(1.0, SVector(2.0, -2.0)))
    solve(Euler(), SVector(1.1, 1.0), sample(tt, Wiener{SV2}()), Wiener{SV2}())

    sample(tt, Wiener())
    sample(tt, WienerBridge(1.0, 1.0))
    true
end

# tests with alpha 0.01
# test fail 1% of the time

r = 2.576
n = 1000

# varmu(x, mu) = dot(x-mu, x-mu)/length(x)
# var0(x) = dot(x, x)/length(x)

#testing mean and variance of Brownian motion sampled at few points
X = [sample(range(0.0, stop=2.0, length=5), Wiener()).yy[end] for i in 1:n]

@test abs(mean(X)) < r * sqrt(2/n)
# see above, fails 1% of the time
# if you want to check: 99% ≈ mean([abs(mean([sample(linspace(0.0, 2.0, 5), Wiener()).yy[end] for i in 1:n]))< r*sqrt(2/n) for k in 1:10000])

chiupper = 1118.95 #upper 0.005 percentile for n = 1000
chilower = 888.56 #lower 0.005 percentile for n = 1000
@test chiupper > 1000.0 * var(X)/2 > chilower

# mean and covariance of vector brownian motion

@test norm(SVector(5.0,2.0) - mean([sample(range(0.0, stop=2.0, length=5),Wiener{SV2}(),SVector(5.0,2.0)).yy[end] for i = 1:div(n,2)])) < r *sqrt(2*2.82 / n)
@test norm(2I - cov(Bridge.mat([sample(range(0.0, stop=2.0, length=5),Wiener{SV2}(),SVector(5.0,2.0)).yy[end] for i = 1:1000]), dims=2, corrected=true))  < 0.3

@test Bridge.a(0.0, 0.0, Wiener()) == inv(Bridge.Γ(0.0, 0.0, Wiener())) == Bridge.σ(0.0, 0.0, Wiener())
@test Bridge.a(0.0, 0.0, WienerBridge(1.0, 0.0)) == inv(Bridge.Γ(0.0, 0.0, WienerBridge(1.0, 0.0))) ==
    Bridge.σ(0.0, 0.0, WienerBridge(1.0, 0.0))

@test mean(transitionprob(0.0, 0.0, 2.0, Wiener())) == 0
@test var(transitionprob(0.0, 0.0, 2.0, Wiener())) ≈ 2.0


text:
module LMMP

    using CPLEX
    using JuMP
    export LMMPSolver
    include("Types.jl")
    include("Functions.jl")

end


text:
#
## Conversion methods
#

function convert(::Type{Bool}, ::Sidef_Types_Nil_Nil)
    false
end

function convert(::Type{Bool}, v::Sidef_Types_Bool_Bool)
    v.value
end

function convert(::Type{Bool}, v::Sidef_Types_Number_Number)
    v.value != 0
end

function convert(::Type{Bool}, v::Bool)
    v
end

function convert(::Type{Bool}, v::Any)
    println("Unimplemented conversion of Any $v to Bool")
    false
end

max_stars_repo_name: trizen/Sidef-Deparse-Julia

max_issues_repo_name: trizen/Sidef-Deparse-Julia

max_forks_repo_name: trizen/Sidef-Deparse-Julia


text:
@testset "▶ set_prop/properties         " begin
  f = Figure()
  set(f, alpha=true)
  @test f.transparency

  # checkers
  @test G.id(5, :a) == 5
  @test G.id("hello", :a) == "hello"
  @test G.fl(5, :a) == 5.0
  @test G.fl((5,2),:a) == (5.0,2.0)
  @test G.not(true, :a) == false
  @test G.lc("blAh", :a) == "blah"
  @test G.lc(["bLASD", "BLad"], :a) == ["blasd", "blad"]
  @test G.col(colorant"blue", :a) == colorant"blue"

  @test G.opcol("blue", :a) == colorant"blue"
  @test G.opcol("none", :a) === nothing
  set(f, alpha=false)
  a = 0
  @test_logs (:warn, "Transparent background is only supported when the figure " *
                     "has its transparency property set to 'true'.") (a=G.opcol("none", :a))
  @test a == colorant"white"

  set(f, alpha=true)
  @test G.alpha(0.2, :a) == 0.2
  set(f, alpha=false)
  a = 0
  @test_logs (:warn, "Transparent colors are only supported when the figure " *
                     "has its transparency property set to 'true'. Ignoring α.") (a=G.alpha(.2, :a))
  @test a == 1.0

  text("blah", (0.5,0.5), font="tt")
  @test gca().objects[1].textstyle.font == "tt"
end


text:
using Documenter, BITEModel
# This does not work:
# include("../src/plotting.jl")

makedocs(
    format = :html,
    sitename = "BITE-Model",
    # pages also make the side-bar
    pages = Any[
        "index.md",
        "data-loading.md",
        "forward.md",
        "inverse.md",
        "postproc.md",
        "api.md",
        "index-index.md"]
)


text:
using Statistics,Serialization,StatsPlots, DataFrames, LaTeXStrings

filename_prefix = "..//..//data//working_paper//bailout//xi06//"
plotname_prefix = "..//..//figures//working_paper//bailout_xi06_comparison_" # store the plots here
mkpath(plotname_prefix)

show_capacity = false # for infected plot

# set intial values as in coivd par ini
datat = 7 # interval of data collection
T = 582 #582  # time periods

# define period of plotting for economic plots
start_week = 2 # cut of first week
end_week = 84
# set periods for covid plots
inf_start_day = 1
inf_end_day = 400 # cut of after period 400
inf_start_week = 1
inf_end_week = 58 # 58*7 = 406 and 42*7 = 294

datapoint = fld(T,datat)
include("$(filename_prefix)beta5//alpha1//covid_par_ini.jl")

### blue LINES
# parameter values for xi = 0.6
paras = ["beta5//alpha25","beta5//alpha50","beta5//alpha75","beta5//alpha1"]

# scatterplot
gdploss_means = []
totcas_means = []

for (i,p) in enumerate(paras)
    data_prefix = string("$(filename_prefix)", p , "//")
    # store data scatterplot
    worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    push!(gdploss_means,gdploss)
    push!(totcas_means,totcas)
    if p == "beta5//alpha1"
        global point_A_x = gdploss
        global point_A_y = totcas
    end
    if p == "beta5//alpha25"
        global point_B_x = gdploss
        global point_B_y = totcas
    end
end

ll = length(paras)
gdploss_means = reverse(gdploss_means)
totcas_means = reverse(totcas_means)

suppfig2a = quiver(gdploss_means[1:ll-1],totcas_means[1:ll-1],quiver=(gdploss_means[2:ll]-gdploss_means[1:ll-1],totcas_means[2:ll]-totcas_means[1:ll-1]), color = [:red], linewidth = 2, xlabel = "Average GDP loss [%]", ylabel = "Mortality [%]")
scatter!(suppfig2a, gdploss_means, totcas_means, markersize = 5, markercolor = [:red], legend=false)

ylims!(suppfig2a,0.008,0.065)
xlims!(suppfig2a,0,10)

### blue LINES
paras = ["beta5//alpha1","beta10//alpha1","beta30//alpha1"]

# scatterplot
gdploss_means4 = []
totcas_means4 = []

for (i,p) in enumerate(paras)
    data_prefix = string("$(filename_prefix)", p , "//")
    # store data scatterplot
    worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    push!(gdploss_means4,gdploss)
    push!(totcas_means4,totcas)
    if p == "beta30//alpha1"
        global point_C_x = gdploss
        global point_C_y = totcas
    end
end

ll = length(paras)

quiver!(suppfig2a, gdploss_means4[1:ll-1],totcas_means4[1:ll-1],quiver=(gdploss_means4[2:ll]-gdploss_means4[1:ll-1],totcas_means4[2:ll]-totcas_means4[1:ll-1]), color = [:blue], linewidth = 2)
scatter!(suppfig2a, gdploss_means4, totcas_means4, markersize = 5, markercolor = [:blue], legend=false)

### add specific points A B C
scatter!(suppfig2a, [point_A_x], [point_A_y], markersize = 6, markercolor = [:black])
annotate!(suppfig2a,[point_A_x+0.4], [point_A_y+0.0005], text.([L"A_1"]) , fontsize=60 )

scatter!(suppfig2a, [point_B_x], [point_B_y], markersize = 6, markercolor = [:black])
annotate!(suppfig2a,[point_B_x+0.4], [point_B_y+0.0015], text.([L"B_1"]) , fontsize=160 )

scatter!(suppfig2a, [point_C_x], [point_C_y], markersize = 6, markercolor = [:black])
annotate!(suppfig2a,[point_C_x+0.4], [point_C_y], text.([L"C_1"]) , fontsize=160 )

plot!(suppfig2a, size=(600,400))
savefig(suppfig2a,"$(plotname_prefix)supp_fig_2a.pdf")