text:
laptop_joltage(adapters) = maximum(adapters) + 3

parse_adapters(s) = map(x->parse(Int, x), split(s, keepempty=false))

#=
adapter_str = """
16
10
15
5
1
11
7
19
6
12
4
"""
adapters = parse_adapters(adapter_str)
start_jolt = 0
end_jolt = laptop_joltage(adapters)
=#

function diff_dist(adapters, start_jolt, end_jolt)
    v = sort([adapters..., start_jolt, end_jolt])
    diffs = diff(v)
    m = maximum(diffs)
    dist = zeros(m)
    for d in diffs
        dist[d] += 1
    end
    dist
end

function num_paths(adapters, start_jolt, end_jolt)
    v = sort([adapters..., start_jolt, end_jolt])
    n = length(v)
    acc_paths = zeros(n)
    acc_paths[1] = 1
    for i in 1:n, j in i+1:i+3
        if j > n
            continue
        end
        if v[j] - v[i] <= 3
            acc_paths[j] += acc_paths[i]
        end
    end
    last(acc_paths)
end

#=
diff_dist(adapters, start_jolt, end_jolt)
=#

function day10part1(adapters)
    start_jolt = 0
    end_jolt = laptop_joltage(adapters)
    dist = diff_dist(adapters, start_jolt, end_jolt)
    dist[1] * dist[3]
end

function day10part2(adapters)
    start_jolt = 0
    end_jolt = laptop_joltage(adapters)
    num_paths(adapters, start_jolt, end_jolt)
end

function day10(file)
    s = read(file, String)
    adapters = parse_adapters(s)
    (day10part1(adapters), day10part2(adapters))
end

#=
day10("input.txt")
=#


text:
macro test_expression(expr)
    esc(quote
            @test JuMP.isequal_canonical(@expression(m, $expr), $expr)
    end)
end

macro test_expression_with_string(expr, str)
    esc(quote
            @test string(@inferred $expr) == $str
            @test_expression $expr
    end)
end

# Test that the macro call `m` throws an error exception during pre-compilation
macro test_macro_throws(errortype, m)
    # See https://discourse.julialang.org/t/test-throws-with-macros-after-pr-23533/5878
    :(@test_throws $errortype try @eval $m catch err; throw(err.error) end)
end

# Define test data structures
struct BadSet <: AbstractInfiniteSet end
struct TestBridge{C} <: MOI.Bridges.AbstractBridge where {C} end
struct BadData <: AbstractMeasureData end
struct Bad end
struct NotASetType end

# Define test functions
function new_fn end

ext: jl


text:
vec_unwrapped = range(0, stop=5π, length=20)
mat_unwrapped = vec_unwrapped .+ vec_unwrapped'
mat_wrapped = mat_unwrapped .% (2π)
uw_test = unwrap(mat_wrapped)
difference = mat_unwrapped[1,1] - uw_test[1,1]
@test (uw_test .+ difference) ≈ mat_unwrapped
# test in-place version
unwrap!(mat_wrapped)
difference = mat_unwrapped[1,1] - mat_wrapped[1,1]
@test (mat_wrapped .+ difference) ≈ mat_unwrapped

# Test BigFloat
vec_bf = range(0, stop=5*BigFloat(π), length=20)
mat_bf_uw = vec_bf .+ vec_bf'
mat_bf_w = mat_bf_uw .% (2*BigFloat(π))
uw_bf_test = unwrap(mat_bf_w)
difference = mat_bf_uw[1,1] - uw_bf_test[1,1]
@test eltype(uw_bf_test) == BigFloat
@test (uw_bf_test .+ difference) ≈ mat_bf_uw

# Test Float32
vec_f32 = range(0, stop=5*Float32(π), length=20)
mat_f32_uw = vec_f32 .+ vec_f32'
mat_f32_w = mat_f32_uw .% (2*Float32(π))
uw_f32_test = unwrap(mat_f32_w)
difference = mat_f32_uw[1,1] - uw_f32_test[1,1]
@test eltype(uw_f32_test) == Float32
@test (uw_f32_test .+ difference) ≈ mat_f32_uw

# Test wrap_around
# after unwrapping, pixels at borders should be equal to corresponding pixels
# on other side
wrap_around = (true, true)
wa_vec = range(0, stop=5π, length=20)
wa_uw = wa_vec .+ zeros(20)'
# make periodic
wa_uw[end, :] = wa_uw[1, :]
wa_w = wa_uw .% (2π)
wa_test = unwrap(wa_w, wrap_around, 0)
difference = wa_uw[1,1] - wa_test[1,1]
# with wrap-around, the borders should be equal, but for this problem the
# image may not be recovered exactly
@test wa_test[:, 1] ≈ wa_test[:, end]
@test wa_test[end, :] ≈ wa_test[1, :]
# In this case, calling unwrap w/o wrap_around does not recover the borders
wa_test_nowa = unwrap(wa_w)
@test !(wa_test_nowa[end, :] ≈ wa_test_nowa[1, :])

max_forks_repo_name: platawiec/Unwrap


text:
using Test, PolyChaos

d, nunc = 5, 10
op = OrthoPoly("genlaguerre",d,Dict(:shape=>1.23))
opq = OrthoPolyQ(op)
mop = MultiOrthoPoly([op for i in 1:nunc],d)
mopq = MultiOrthoPoly([opq for i in 1:nunc],d)
@testset "dimenions" begin
    @test isequal(PolyChaos.dim(op), d + 1)
    @test isequal(PolyChaos.dim(opq), d + 1)
    @test isequal(PolyChaos.dim(mop), factorial(d + nunc) / (factorial(d) * factorial(nunc)))
end

@testset "degrees" begin
    @test isequal(PolyChaos.deg(op),d)
    @test isequal(PolyChaos.deg(opq),d)
    @test isequal(PolyChaos.deg(mop),d)
end

NW = [ opq.quad.nodes opq.quad.weights ]
@testset "nodes and weights" begin
    @test PolyChaos.nw(opq) == NW
    @test PolyChaos.nw([opq for i in 1:nunc])[1] == [NW[:,1] for i in 1:nunc]
    @test PolyChaos.nw([opq for i in 1:nunc])[2] == [NW[:,2] for i in 1:nunc]
    @test PolyChaos.nw(mopq)[1] == [NW[:,1] for i in 1:nunc]
    @test PolyChaos.nw(mopq)[2] == [NW[:,2] for i in 1:nunc]
end

c = [ op.α op.β ]
@testset "coefficients" begin
    @test PolyChaos.coeffs(op) == c
    @test PolyChaos.coeffs(opq) == c
    @test PolyChaos.coeffs([op for i in 1:nunc])[1] == [ c[:,1] for i in 1:nunc ]
    @test PolyChaos.coeffs([op for i in 1:nunc])[2] == [ c[:,2] for i in 1:nunc ]
    @test PolyChaos.coeffs([opq for i in 1:nunc])[1] == [ c[:,1] for i in 1:nunc ]
    @test PolyChaos.coeffs([opq for i in 1:nunc])[2] == [ c[:,2] for i in 1:nunc ]
    @test PolyChaos.coeffs(mop)[1] == [ c[:,1] for i in 1:nunc ]
    @test PolyChaos.coeffs(mop)[2] == [ c[:,2] for i in 1:nunc ]
end

@test issymmetric(op) == issymmetric(opq) == issymmetric(opq.quad)


text:
include("read_base.jl")
include("read_bedfile.jl")
include("read_alphas.jl")

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
type GrenanderEstimator <: LocalFdrEstimator
  pi0estimator:: Pi0Estimator
end

# saves everything in sorted way for now, change that later...
type GrenanderLocalFdrFit <: LocalFdrFit
  pvalues   :: Array{Float64,1}
  f         :: Array{Float64,1} #density
  F         :: Array{Float64,1} #distribution
  pi0       :: Float64
end

pvalues(fdrfit::GrenanderLocalFdrFit) = fdrfit.pvalues
pi0(fdrfit::GrenanderLocalFdrFit)  = fdrfit.pi0
distribution(fdrfit::GrenanderLocalFdrFit)    = fdrfit.F
density(fdrfit::GrenanderLocalFdrFit)    =    fdrfit.f
localfdr(fdrfit::GrenanderLocalFdrFit)   =    fdrfit.pi0./fdrfit.f #TODO make sure between [0,1]
tailfdr(fdrfit::GrenanderLocalFdrFit)    =    fdrfit.pi0.*fdrfit.pvalues./fdrfit.F #TODO make sure between [0,1]

function fit(pv:: Array{Float64,1}, fdrestimator::GrenanderEstimator)
  m = length(pv)
  pi0 = estimatepi0(pv, fdrestimator.pi0estimator)

  pv_ecdf_fun = ecdf(pv)          #going about it in quite inefficient way but this should not be the bottleneck anyway
  pv_sorted = sort(pv)
  ecdf_value = Array(Float64,m)

  # apply modification to make grenander estimate consistent with two-groups model
  for i=1:m
    ecdf_value[i] = pv_ecdf_fun(pv_sorted[i])
    if ecdf_value[i] < pi0*pv_sorted[i]
      ecdf_value[i] = pi0*pv_sorted[i]
    elseif ecdf_value[i] > 1-pi0*(1-pv_sorted[i])
      ecdf_value[i] = 1-pi0*(1-pv_sorted[i])
    end
  end

  #denote by (x,y) the coords of the modified ecdf
  Δx = diff(pv_sorted)
  Δy = float64(diff(ecdf_value))

  f = Δy./Δx #rough density estimate

  f = -isotonicregression(-f, Δx) # make previous density estimate obey monotonicity
  F  = ecdf_value[1] + vcat(0,cumsum(f .* Δx))
  f = push!(f, f[m-1])

  GrenanderLocalFdrFit(pv_sorted,f,F,pi0)
end


text:
struct VulkanError <: Exception
    msg::AbstractString
    errorcode
end
Base.showerror(io::Core.IO, e::VulkanError) = print(io, "$(e.errorcode): ", e.msg)

"""
    @vkcheck vkFunctionSomething()

Checks whether the expression returned VK_SUCCESS. Else, throw an error printing the corresponding code."""
macro vkcheck(expr)
    quote
        local expr_return_code = $(esc(expr))
        if typeof(expr_return_code) != VkResult
            throw(ErrorException("the return value is not a valid code"))
        end
        if expr_return_code != VK_SUCCESS
            local str_error = $(string(expr))
            throw(VulkanError("failed to execute $str_error", expr_return_code))
        end
    end
end


text:
@everywhere module DatasetPreloader
	using HDF5
	using Distributed
	export preload, load, getPreloadFileName

	preloadFileName = ""
	preloadFuture = 0

	function readFromHDF(filename)
		GC.gc()
		try
			return h5read(filename,"/FullSpectra/TofData")
		catch
			return 0
		end
	end

	getPreloadFileName() = preloadFileName

	function preload(filename)
		global preloadFileName = filename
		println("preloading file $filename")
		global preloadFuture = remotecall(readFromHDF,2,filename)
	end

	function load(filename)
		global preloadFileName
		global preloadFuture
		if (preloadFileName != "" && preloadFileName == filename && preloadFuture != 0)
			println("fetching $filename from future")
			ds = fetch(preloadFuture)
			preloadFuture = 0
			preloadFilename = ""
			return ds
		else
			println("reading $filename directly")
			preloadFuture = 0
			preloadFilename = ""
			return readFromHDF(filename)
		end
	end
end


text:
using YaoBlocks, YaoArrayRegister
using Random, Test

@testset "identity gate" begin
    block = igate(4)
    reg = rand_state(4)
    @test reg |> block == reg
    @test mat(block) != nothing
    @test ishermitian(block)
    @test isreflexive(block)
    @test isunitary(block)
    @test occupied_locs(block) == ()
    @test getiparams(block) == ()
    @test nqudits(igate(3; nlevel=3)) == 3
    @test mat(igate(3; nlevel=3)) == YaoBlocks.IMatrix{27}()
end


text:
include("./NAMD_example1.jl")
include("./NAMD_example2.jl")
include("./NAMD_example3.jl")
include("./NAMD_example4.jl")
include("./NAMD_example5.jl")
include("./NAMD_example6.jl")
include("./NAMD_example7.jl")

include("./ivan.jl")


text:
function get_time()
    return string(Dates.now())
end


text:
function hierNet_boot(df; nfolds::Integer = 5, nlam::Integer = 10, strong::Bool = false, nboot::Integer = 5)
# --------------------------------------------------------------------------
# Options:
# nfolds - number of folds for cross-validation
# nlam - number of values of regularization parameter to try
# strong -  whether strong hierarchy is imposed or not
# nboot - number of bootstrapped samples on which lasso procedure is run
# --------------------------------------------------------------------------

# Unpack DataFrame df
s_all = convert(Array{Float64,2}, df[:, r"shares"]);
p_all = convert(Array{Float64,2}, df[:, r"prices"]);
x_all = convert(Array{Float64,2}, df[:, r"x"]);
iv_all = convert(Array{Float64,2}, df[:, r"demand_instruments"]);

J = size(s_all,2);

@rput s_all
@rput p_all
@rput iv_all
@rput strong
@rput nlam
@rput nfolds
@rput nboot

R"
library(hierNet)
J <-  ncol(s_all)

rownames(s_all) <- seq(length=nrow(s_all))
rownames(p_all) <- seq(length=nrow(p_all))
rownames(iv_all) <- seq(length=nrow(iv_all))

for (b in 1:nboot) {
  boot_sample <-sample(1:nrow(s_all), nrow(s_all), replace=TRUE, prob=NULL)
  s <- s_all[boot_sample,]
  p <- p_all[boot_sample,]
  iv <- iv_all[boot_sample,]

  # First stage IV for market shares
  s_hat <- matrix(nrow = nrow(s),ncol = ncol(s))
  one_mat <- matrix(1,nrow(s), 1)
  iv <- as.matrix(iv)
  for (j in 1:J) {
    iv_temp <- as.matrix(cbind(matrix(1,nrow(s),1),iv))
    beta<- solve(t(iv_temp)%*%iv_temp) %*% t(iv_temp) %*% s[,j]
    s_hat[,j] <-iv_temp%*% beta
  }

  # Estimate Hierarchical Lasso
  for (j in 1:J) {
     # Estimate path of regularization parameters
    structured_fit <- hierNet.path(as.matrix(s_hat), as.vector(p[,j]),
                                  lamlist = NULL, delta=1e-8, minlam = NULL, maxlam = NULL, nlam=nlam, flmin=.05,
                                  diagonal = TRUE, strong = strong, aa = NULL, zz = NULL,
                                  stand.main = TRUE, stand.int = TRUE,
                                  rho = nrow(s), niter = 100, sym.eps = 0.001,
                                  step = 1, maxiter = 2000, backtrack = 0.2, tol = 1e-05, trace = 0)
    if (j==1) {
      matsize <- nrow(structured_fit$bp)
      coefs <- matrix(nrow = 0 ,ncol = matsize)
    }

    # Select parameter by CV, append to matrix
    a<- hierNet.cv(structured_fit, as.matrix(s_hat), as.vector(p[,j]), nfolds=nfolds,folds=NULL,trace=0)
    lambdahat <- a['lamhat']
    # lamhat lamhat.1se

    # Final Estimates
    structured_fit<- hierNet(as.matrix(s_hat), as.vector(p[,j]), as.numeric(lambdahat), delta=1e-8, strong=strong, diagonal=TRUE, aa=NULL,
                               zz=NULL,center=TRUE, stand.main=TRUE, stand.int=TRUE, rho=nrow(s), niter=100,
                               sym.eps=1e-3,step=1, maxiter=5000, backtrack=0.2, tol=1e-5, trace=0)
    coefs <- rbind(coefs,t(structured_fit$bp - structured_fit$bn))

  }
  if (b==1){
    included <- (abs(coefs)>0)
  }

  if (b>1) {
    temp <- (abs(coefs)>0)
    included <- 1*as.matrix((temp==1 & included==1))
  }

  for (j in 1:J) {
    included[j,j] = 1
  }
}
included_symmetric <- included
for (j1 in 1:J) {
  for (j2 in 1:J) {
    if (included_symmetric[j1,j2] != included_symmetric[j2,j1]){
      included_symmetric[j1,j2] = 1
      included_symmetric[j2,j1] = 1
    }
  }
}
"
@rget included
@rget included_symmetric
included = convert.(Integer, included)
included_symmetric = convert.(Integer, included_symmetric)
return included, included_symmetric
end


text:
function polyfit(x, y, n)
  A = [ float(x[i])^p for i = 1:length(x), p = 0:n ]
  A \ y
end

function polyval(p,x::Number)
	y = zero(promote_type(eltype(p), eltype(x)))
	for i=1:length(p)
		y+=p[i]*x^(i-1)
	end
	y
end
function polyval(p,x::AbstractArray)
	y = zeros(promote_type(eltype(p), eltype(x)), length(x))
	for i = 1:length(x)
		y[i] = polyval(p,x[i])
	end
	y
end
function polyresiduals(p,x,y)
	yp = polyval(p,x)
	residuals = y-yp
end

### CALIBRATION
function peakassignlog(locations_arb, locations_true)
# locations_arb is an array of possible peak locations in arbitrary units, assumed to have an unknown,
# nonlinear (but monotonic and smooth) relationship to locations_true
# locations_arb may contain more peaks than locations_true, including some spurious peaks
# the goal is to find the best possible set of assignments for entries in locations_arb to
# correspond to locations_true by examining all possible assignments and calcualting some figure of merit
@assert(issorted(locations_arb))
@assert(issorted(locations_true))
@assert(length(locations_arb)>=length(locations_true))
combos = combinations(1:length(locations_arb), length(locations_true))
length(combos)==1 && return (locations_arb[first(combos)], 0.0)

assert(!(0 in locations_true)) # don't pass zero and force zero
assert(!(0 in locations_arb))

log_locations_arb = log(locations_arb)
log_locations_true = log(locations_true)

best_combo = first(combos)
best_fom = realmax(Float64)
for c in combos
	cu = log_locations_arb[c]
	fom = peakassign_fig_of_merit2(cu, log_locations_true)
	if fom < best_fom
		best_fom = fom
		best_combo = c
	end
end
locations_arb[best_combo], best_fom
end

function peakassign2(locations_arb, locations_true)
# locations_arb is an array of possible peak locations in arbitrary units, assumed to have an unknown,
# nonlinear (but monotonic and smooth) relationship to locations_true
# locations_arb may contain more peaks than locations_true, including some spurious peaks
# the goal is to find the best possible set of assignments for entries in locations_arb to
# correspond to locations_true by examining all possible assignments and calcualting some figure of merit
@assert(issorted(locations_arb))
@assert(issorted(locations_true))
@assert(length(locations_arb)>=length(locations_true))
combos = combinations(1:length(locations_arb), length(locations_true))
length(combos)==1 && return (locations_arb[first(combos)], 0.0)

assert(!(0 in locations_true)) # don't pass zero and force zero
assert(!(0 in locations_arb))

best_combo = first(combos)
best_fom = realmax(Float64)
for c in combos
	cu = locations_arb[c]
	fom = peakassign_fig_of_merit2(cu, locations_true)
	if fom < best_fom
		best_fom = fom
		best_combo = c
	end
end
locations_arb[best_combo], best_fom
end

function peakassign2zeros(locations_arb, locations_true)
# locations_arb is an array of possible peak locations in arbitrary units, assumed to have an unknown,
# nonlinear (but monotonic and smooth) relationship to locations_true
# locations_arb may contain more peaks than locations_true, including some spurious peaks
# the goal is to find the best possible set of assignments for entries in locations_arb to
# correspond to locations_true by examining all possible assignments and calcualting some figure of merit
@assert(issorted(locations_arb))
@assert(issorted(locations_true))
@assert(length(locations_arb)>=length(locations_true))
combos = combinations(1:length(locations_arb), length(locations_true))
length(combos)==1 && return (locations_arb[first(combos)], 0.0)

assert(!(0 in locations_true)) # don't pass zero and force zero
assert(!(0 in locations_arb))
locations_true = unshift!(copy(locations_true), 0)

best_combo = first(combos)
best_fom = realmax(Float64)
for c in combos
	cu = locations_arb[c]
	unshift!(cu,0)
	fom = peakassign_fig_of_merit2(cu, locations_true)
	if fom < best_fom
		best_fom = fom
		best_combo = c
	end
end
locations_arb[best_combo], best_fom
end

function peakassign2(locations_arb, locations_true)
# locations_arb is an array of possible peak locations in arbitrary units, assumed to have an unknown,
# nonlinear (but monotonic and smooth) relationship to locations_true
# locations_arb may contain more peaks than locations_true, including some spurious peaks
# the goal is to find the best possible set of assignments for entries in locations_arb to
# correspond to locations_true by examining all possible assignments and calcualting some figure of merit
@assert(issorted(locations_arb))
@assert(issorted(locations_true))
@assert(length(locations_arb)>=length(locations_true))
combos = combinations(1:length(locations_arb), length(locations_true))
length(combos)==1 && return (locations_arb[first(combos)], 0.0)

assert(!(0 in locations_true)) # don't pass zero and force zero
assert(!(0 in locations_arb))

best_combo = first(combos)
best_fom = realmax(Float64)
for c in combos
	cu = locations_arb[c]
	fom = peakassign_fig_of_merit(cu, locations_true)
	if fom < best_fom
		best_fom = fom
		best_combo = c
	end
end
locations_arb[best_combo], best_fom
end

function peakassign_fig_of_merit(x, y)
	# average fractional error trying to predict the y value of the middle of 3 points using the slope
	# from the outer two points and the x value of the middle point
	err = 0.0
	for i = 1:length(y)-2
		yip1_prediction = y[i]+(y[i+2]-y[i])*(x[i+1]-x[i])/(x[i+2]-x[i])
		err += abs(1-y[i+1]/yip1_prediction)
	end
	err/length(y-2)
end

function peakassign_fig_of_merit2(x, y)
	n = length(x) > 2 ? 2 : 1
	p = polyfit(x,y,n)
	std(polyresiduals(p,x,y)./y)
end

# test peakassign

known_energies = [3691, 4090, 4510, 4952, 5414, 5898, 6403, 6929, 7477, 8040, 8637] # ka energies Ca-Zn

filt_values = known_energies.^0.9 # caluclate filt values with some amount of non-linearity

tested = 0
peakassignlogwrong=0
peakassign2wrong=0
peakassign2zerowrong=0
peakassignwrong=0

nc = 7
	tested+=length(combinations(1:length(known_energies),nc))
	println("testing peakassignlog $nc")
	for c in combinations(1:length(known_energies),nc)
		locations_arb = filt_values
		locations_true = known_energies[c]
		best_combo, best_fom = peakassignlog(locations_arb, locations_true)
		if !all(best_combo .== filt_values[c])
			peakassignlogwrong+=1
			# println(locations_true)
			# println(best_combo)
			# println(filt_values[c])
			# @test best_fom < 1
			# @test all(best_combo .== filt_values[c])
		end
	end

	println("testing peakassign2 $nc")
	for c in combinations(1:length(known_energies),nc)
		locations_arb = filt_values
		locations_true = known_energies[c]
		best_combo, best_fom = peakassign2(locations_arb, locations_true)
		if !all(best_combo .== filt_values[c])
			peakassign2wrong+=1

		end
	end

	println("testing peakassign2zeros $nc")
	for c in combinations(1:length(known_energies),nc)
		locations_arb = filt_values
		locations_true = known_energies[c]
		best_combo, best_fom = peakassign2zeros(locations_arb, locations_true)
		if !all(best_combo .== filt_values[c])
			peakassign2zerowrong+=1

		end
	end
	println("testing peakassign $nc")
	for c in combinations(1:length(known_energies),nc)
		locations_arb = filt_values
		locations_true = known_energies[c]
		best_combo, best_fom = peakassign2(locations_arb, locations_true)
		if !all(best_combo .== filt_values[c])
			peakassignwrong+=1

		end
	end

println("peakassignlog $(peakassignlogwrong/tested)")
println("peakassign2 $(peakassign2wrong/tested)")
println("peakassign2zeros $(peakassign2zerowrong/tested)")
println("peakassign $(peakassignwrong/tested)")


text:
function μ(gun::Gun, target::AbstractTarget)
    μx = paralax(gun,target)[1]*target.position
    μy = paralax(gun,target)[2]*target.position
    return μx,μy
end

"""
    aiming_error(aiming)

Returns the alignment error (horizontal and vertical) for the error sources defined in 'aiming'
"""
function aiming_error(aiming::AimingError)
    σ_x =sqrt(aiming.point_def[1]^2+aiming.human[1]^2+aiming.res[1]^2+aiming.opt_path[1]^2)
    σ_y =sqrt(aiming.point_def[2]^2+aiming.human[2]^2+aiming.res[2]^2+aiming.opt_path[2]^2)
    return σ_x, σ_y
end
"""
    alignment_error(align)

Returns the alignment error (horizontal and vertical) for the error sources defined in 'align'
"""
function alignment_error(align::AlignmentError)
    σ_x=sqrt(align.boresight_acc[1]^2+align.boresight_ret[1]^2+align.sync[1]^2+align.stab[1]^2+align.gun_pos[1]^2)
    σ_y=sqrt(align.boresight_acc[2]^2+align.boresight_ret[2]^2+align.sync[2]^2+align.stab[2]^2+align.gun_pos[2]^2)
    return σ_x, σ_y
end

"""
    lof_error(lof)

Returns the line of fire error (horizontal and vertical) for the error sources defined in 'lof'
"""
function lof_error(lof::LofError)
    σ_x=sqrt(lof.fcs[1]^2+lof.lead[1]^2+lof.cant[1]^2+lof.earth_rate[1]^2+lof.muzzle_vel[1]^2+lof.temp[1]^2+lof.press[1]^2+lof.cross_wind[1]^2+lof.range_wind[1]^2+lof.range[1]^2+lof.jump[1]^2+lof.zeroing[1]^2)
    σ_y=sqrt(lof.fcs[2]^2+lof.lead[2]^2+lof.cant[2]^2+lof.earth_rate[2]^2+lof.muzzle_vel[2]^2+lof.temp[2]^2+lof.press[2]^2+lof.cross_wind[2]^2+lof.range_wind[2]^2+lof.range[2]^2+lof.jump[2]^2+lof.zeroing[2]^2)
    return σ_x, σ_y
end

"""
    flight_error(flight)

Returns the projectile dispersion (horizontal and vertical)
"""
function flight_error(flight::FlightError)
    σ_x=sqrt(flight.aero[1]^2)
    σ_y=sqrt(flight.aero[2]^2)
    return σ_x, σ_y
end

"""
    σ(aiming,align,lof,flight)

Returns the total error (horizontal and vertical) with error budget defined by error sources:
* aiming: aiming error
* align: alignment error
* lof: line of fire error
* flight: projectile dispersion
"""
function σ(aiming::AimingError,align::AlignmentError,lof::LofError,flight::FlightError)
    σ_x =sqrt(aiming_error(aiming)[1]^2+alignment_error(align)[1]^2+lof_error(lof)[1]^2+flight_error(flight)[1]^2)
    σ_y =sqrt(aiming_error(aiming)[2]^2+alignment_error(align)[2]^2+lof_error(lof)[2]^2+flight_error(flight)[2]^2)
    return σ_x, σ_y
end

"""
    σ()

Returns the total error of type 'Error'. The error budget is defined by type of errors:
* σfixed : fixed bias errors
* σvariable: varible bias errors
* σrandom: random errors
"""

function σ(;σfixed=Error(0.0,0.0), σvariable=Error(0.0,0.0), σrandom=Error(0.0,0.0))
    σv = σfixed.vertical + sqrt(σvariable.vertical^2 + σrandom.vertical^2)
    σh = σfixed.horizontal + sqrt(σvariable.horizontal^2 + σrandom.horizontal^2)
    σtot = Error(σh,σv)

    return σtot
end


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

JLLWrappers.@generate_main_file_header("Attr")
JLLWrappers.@generate_main_file("Attr", UUID("1fd713ca-387f-5abc-8002-d8b8b1623b73"))
end  # module Attr_jll


text:
module NaiveNASlib

using LightGraphs, MetaGraphs
using Statistics
using Logging

import JuMP
import JuMP: @variable, @constraint, @objective, @expression, MOI, MOI.INFEASIBLE, MOI.FEASIBLE_POINT
using Cbc

#Interface
export AbstractVertex, AbstractMutationVertex, MutationOp, MutationState

# Vertex
export InputVertex, CompVertex, inputs, outputs

# Information strings
export infostr, name, RawInfoStr, NameInfoStr, InputsInfoStr, OutputsInfoStr, SizeInfoStr, MutationTraitInfoStr, ComposedInfoStr, NameAndInputsInfoStr, NinInfoStr, NoutInfoStr, NameAndIOInfoStr, FullInfoStr,MutationSizeTraitInfoStr

# Computation graph
export CompGraph, SizeDiGraph, output!,flatten, nv, vertices

# Mutation operations
#State
export InvSize, IoSize, InvIndices, IoIndices, NoOp, IoChange, nin, nout, Δnin, Δnout, clone, op, in_inds, out_inds, nin_org, nout_org

# Mutation vertex
export base, InputSizeVertex, OutputsVertex, MutationVertex

# Mutation traits
export trait, MutationTrait, DecoratingTrait, NamedTrait, Immutable, MutationSizeTrait, SizeAbsorb, SizeStack, SizeInvariant, SizeChangeLogger, SizeChangeValidation

# Size util
export minΔnoutfactor, minΔninfactor, minΔnoutfactor_only_for, minΔninfactor_only_for, findterminating, ΔSizeInfo, ΔninSizeInfo, ΔnoutSizeInfo, ΔSizeGraph, ΔninSizeGraph, ΔnoutSizeGraph, Direction, Input, Output, Both, all_in_graph, all_in_Δsize_graph, Δsize, newsizes, neighbours, fullgraph

export AbstractΔSizeStrategy, AbstractJuMPΔSizeStrategy, ΔSizeFailError, ΔSizeFailNoOp, LogΔSizeExec, DefaultJuMPΔSizeStrategy, ΔNout, ΔNoutExact, ΔNoutRelaxed, ΔNin, ΔNinExact, ΔNinRelaxed, AlignNinToNout, Exact, Relaxed

#Selection util
export AbstractSelectionStrategy, LogSelection, LogSelectionFallback, SelectionFail, NoutRevert, SelectDirection, ApplyAfter, AbstractJuMPSelectionStrategy, DefaultJuMPSelectionStrategy, OutSelect, OutSelectExact, OutSelectRelaxed, TruncateInIndsToValid, Δoutputs, solve_outputs_selection

# Connectivity mutation
export remove!, RemoveStrategy, insert!, create_edge!, remove_edge!

# Align size strategies, e.g what to do with sizes of vertices connected to a removed vertex
export AbstractAlignSizeStrategy, IncreaseSmaller, DecreaseBigger, AlignSizeBoth, ChangeNinOfOutputs, AdjustToCurrentSize, FailAlignSizeNoOp, FailAlignSizeError, FailAlignSizeWarn, FailAlignSizeRevert, NoSizeChange, CheckAligned, CheckNoSizeCycle, CheckCreateEdgeNoSizeCycle, PostAlignJuMP, SelectOutputs, ApplyMutation, PostSelectOutputs, PostApplyMutation

# Connect strategies
export AbstractConnectStrategy, ConnectAll, ConnectNone

# apply mutation
export mutate_inputs, mutate_outputs, apply_mutation

#sugar
export inputvertex, vertex, immutablevertex, absorbvertex, invariantvertex, conc, VertexConf, traitconf, mutationconf, outwrapconf

include("vertex.jl")
include("compgraph.jl")
include("prettyprint.jl")

include("mutation/op.jl")
include("mutation/vertex.jl")
include("mutation/graph.jl")
include("mutation/jumpnorm.jl")
include("mutation/size.jl")
include("mutation/apply.jl")
include("mutation/select.jl")

include("mutation/structure.jl")

include("mutation/sugar.jl")

end # module


text:
const use_stdcall = startswith(basename(libmpi), "msmpi")

# This is required in addition to __init__() so that we can call library at precompilation time.
if Sys.isunix()
    Libdl.dlopen(libmpi, Libdl.RTLD_LAZY | Libdl.RTLD_GLOBAL)
end

macro mpicall(expr)
    @assert expr isa Expr && expr.head == :call && expr.args[1] == :ccall

    # On unix systems we call the global symbols to allow for LD_PRELOAD interception
    # It can be emulated in Windows (via Libdl.dllist), but this is not fast.
    if Sys.isunix() && expr.args[2].head == :tuple &&
            (VERSION ≥ v"1.5-" || expr.args[2].args[1] ≠ :(:MPI_Get_library_version))
        expr.args[2] = expr.args[2].args[1]
    end

    # Microsoft MPI uses stdcall calling convention
    # this only affects 32-bit Windows
    # unfortunately we need to use ccall to call Get_library_version
    # so check using library name instead
    if use_stdcall
        insert!(expr.args, 3, :stdcall)
    end
    return esc(expr)
end

function Get_library_version()
    # There is no way to query at runtime what the length of the buffer should be.
    # https://github.com/mpi-forum/mpi-issues/issues/159
    # 8192 is the maximum value of MPI_MAX_LIBRARY_VERSION_STRING across known
    # implementations.
    buf = Array{UInt8}(undef, 8192)
    buflen = Ref{Cint}()

    @mpicall ccall((:MPI_Get_library_version, libmpi), Cint, (Ptr{UInt8}, Ref{Cint}), buf, buflen)
    @assert buflen[] < 8192
    resize!(buf, buflen[])
    return String(buf)
end

"""
    MPI_LIBRARY_VERSION_STRING :: String

The full version string provided by the library

# External links
$(_doc_external("MPI_Get_library_version"))
"""
const MPI_LIBRARY_VERSION_STRING = Get_library_version()

"""
    impl, version = identify_implementation()

Attempt to identify the MPI implementation based on
[`MPI_LIBRARY_VERSION_STRING`](@ref). Returns a triple of values:

- `impl`: a `String` with the name of the MPI implementation, or `"unknown"` if it cannot be determined,
- `version`: a `VersionNumber` of the library, or `nothing` if it cannot be determined.

This function is only intended for internal use. Users should use [`MPI_LIBRARY`](@ref),
[`MPI_LIBRARY_VERSION`](@ref).
"""
function identify_implementation()
    impl, version, abi = MPIPreferences.identify_implementation_version_abi(MPI_LIBRARY_VERSION_STRING)
    return impl, version
end

const MPI_LIBRARY, MPI_LIBRARY_VERSION = identify_implementation()

"""
    MPI_LIBRARY :: String

The current MPI implementation: this is determined by
"""
MPI_LIBRARY

"""
    MPI_LIBRARY_VERSION :: VersionNumber

The version of the MPI library
"""
MPI_LIBRARY_VERSION

function Get_version()
    major = Ref{Cint}()
    minor = Ref{Cint}()
    @mpicall ccall((:MPI_Get_version, libmpi), Cint,
                   (Ptr{Cint}, Ptr{Cint}), major, minor)
    VersionNumber(major[], minor[])
end

"""
    MPI_VERSION :: VersionNumber

The supported version of the MPI standard.

# External links
$(_doc_external("MPI_Get_version"))
"""
const MPI_VERSION = Get_version()


text:
# Changes from original file:
#   pyimport tinyarray
#   range(A) -> range(0,length=A)
#   A.append(...) -> append!(A,...)

tinyarray = pyimport("tinyarray")
    sigma_0 = tinyarray.array([[1, 0], [0, 1]])
    sigma_x = tinyarray.array([[0, 1], [1, 0]])
    sigma_y = tinyarray.array([[0, -1im], [1im, 0]])
    sigma_z = tinyarray.array([[1, 0], [0, -1]])

function make_system(t=1.0, alpha=0.5, e_z=0.08, W=10, L=30)
    # Start with an empty tight-binding system and a single square lattice.
    # `a` is the lattice constant (by default set to 1 for simplicity).
    lat = kwant.lattice.square()

    syst = kwant.Builder()

    #### Define the scattering region. ####
    syst[[lat(x, y) for x in range(0,length=L) for y in range(0,length=W)]] =
        4 * t * sigma_0 + e_z * sigma_z # note [...] for comprehension and 1:L/W for range(L/W) also no \ for line break
    # hoppings in x-direction
    syst[kwant.builder.HoppingKind((1, 0), lat, lat)] =
        -t * sigma_0 + 1im * alpha * sigma_y / 2 ## note 1im not 1j
    # hoppings in y-directions
    syst[kwant.builder.HoppingKind((0, 1), lat, lat)] =
        -t * sigma_0 - 1im * alpha * sigma_x / 2

    #### Define the left lead. ####
    lead = kwant.Builder(kwant.TranslationalSymmetry((-1, 0)))

    lead[[lat(0, j) for j in range(0,length=W)]] = 4 * t * sigma_0 + e_z * sigma_z
    # hoppings in x-direction
    lead[kwant.builder.HoppingKind((1, 0), lat, lat)] =
        -t * sigma_0 + 1im * alpha * sigma_y / 2
    # hoppings in y-directions
    lead[kwant.builder.HoppingKind((0, 1), lat, lat)] =
        -t * sigma_0 - 1im * alpha * sigma_x / 2

    #### Attach the leads and return the finalized system. ####
    syst.attach_lead(lead)
    syst.attach_lead(lead.reversed())

    return syst
end

function plot_conductance(syst, energies)
    # Compute conductance
    data = []
    for energy in energies
        smatrix = kwant.smatrix(syst, energy)
        append!(data,smatrix.transmission(1, 0))
    end

    pyplot.figure()
    pyplot.plot(energies, data)
    pyplot.xlabel("energy [t]")
    pyplot.ylabel("conductance [e^2/h]")
    pyplot.gcf() ## note gcf not show
end


text:
const disconnect_request_type = JSONRPC.RequestType("disconnect", DisconnectArguments, DisconnectResponseArguments)
const set_break_points_request_type = JSONRPC.RequestType("setBreakpoints", SetBreakpointsArguments, SetBreakpointsResponseArguments)
const set_exception_break_points_request_type = JSONRPC.RequestType("setExceptionBreakpoints", SetExceptionBreakpointsArguments, SetExceptionBreakpointsResponseArguments)
const set_function_exception_break_points_request_type = JSONRPC.RequestType("setFunctionBreakpoints", SetFunctionBreakpointsArguments, SetFunctionBreakpointsResponseArguments)
const stack_trace_request_type = JSONRPC.RequestType("stackTrace", StackTraceArguments, StackTraceResponseArguments)
const scopes_request_type = JSONRPC.RequestType("scopes", ScopesArguments, ScopesResponseArguments)
const source_request_type = JSONRPC.RequestType("source", SourceArguments, SourceResponseArguments)
const variables_request_type = JSONRPC.RequestType("variables", VariablesArguments, VariablesResponseArguments)
const continue_request_type = JSONRPC.RequestType("continue", ContinueArguments, ContinueResponseArguments)
const next_request_type = JSONRPC.RequestType("next", NextArguments, NextResponseArguments)
const step_in_request_type = JSONRPC.RequestType("stepIn", StepInArguments, StepInResponseArguments)
const step_in_targets_request_type = JSONRPC.RequestType("stepInTargets", StepInTargetsArguments, StepInTargetsResponseArguments)
const step_out_request_type = JSONRPC.RequestType("stepOut", StepOutArguments, StepOutResponseArguments)
const evaluate_request_type = JSONRPC.RequestType("evaluate", EvaluateArguments, EvaluateResponseArguments)
const terminate_request_type = JSONRPC.RequestType("terminate", TerminateArguments, TerminateResponseArguments)
const exception_info_request_type = JSONRPC.RequestType("exceptionInfo", ExceptionInfoArguments, ExceptionInfoResponseArguments)
const restart_frame_request_type = JSONRPC.RequestType("restartFrame", RestartFrameArguments, RestartFrameResponseResponseArguments)
const set_variable_request_type = JSONRPC.RequestType("setVariable", SetVariableArguments, SetVariableResponseArguments)
const stopped_notification_type = JSONRPC.NotificationType("stopped", StoppedEventArguments)
const threads_request_type = JSONRPC.RequestType("threads", Nothing, ThreadsResponseArguments)
const breakpointslocation_request_type = JSONRPC.RequestType("breakpointLocations", BreakpointLocationsArguments, BreakpointLocationsResponseArguments)

@dict_readable struct DebugArguments <: Outbound
    stopOnEntry::Bool
    program::String
    compiledModulesOrFunctions::Union{Missing,Vector{String}}
    compiledMode::Union{Missing,Bool}
end

@dict_readable struct ExecArguments <: Outbound
    stopOnEntry::Bool
    code::String
    file::String
    compiledModulesOrFunctions::Union{Missing,Vector{String}}
    compiledMode::Union{Missing,Bool}
end

# Our own requests
const run_notification_type = JSONRPC.NotificationType("run", NamedTuple{(:program,),Tuple{String}})
const debug_notification_type = JSONRPC.NotificationType("debug", DebugArguments)
const exec_notification_type = JSONRPC.NotificationType("exec", ExecArguments)
const finished_notification_type = JSONRPC.NotificationType("finished", Nothing)
const set_compiled_items_notification_type = JSONRPC.NotificationType("setCompiledItems", NamedTuple{(:compiledModulesOrFunctions,),Tuple{Vector{String}}})
const set_compiled_mode_notification_type = JSONRPC.NotificationType("setCompiledMode", NamedTuple{(:compiledMode,),Tuple{Bool}})


text:
@test iscyclic(permcycs((1,2,3)))
@test iscyclic(permcycs((42,100,3)))
@test !iscyclic(permcycs((1,2), (3,4)))


text:
workspace()
include("../src/LargeMatrixEquations.jl")
using LargeMatrixEquations
using Base.Test

nh = 10
n = nh^2#nh^2
T = diagm(ones(nh)*2)-diagm(ones(nh-1),-1)-diagm(ones(nh-1),1)
I = eye(nh,nh)
A = -(diagm(ones(n)*2)-diagm(ones(n-1),-1)-diagm(ones(n-1),1))#-(kron(T,I)+kron(I,T))
srand(1234)
B = ones(n,1)

E = diagm(ones(n),0)
E[1,1]=2
E[n,n]=3

Z,er2=try
        kpik(A,B,E,infoV=false)
      catch
        error("kpik failed with full matrices")
      end
@test norm(A*Z*Z'*E + E*Z*Z'*A' + B*B') < 0.001

Z,er2=try
        kpik(sparse(A),B,sparse(E))
      catch
        error("kpik failed with sparse matrices")
      end

@test norm(A*Z*Z'*E + E*Z*Z'*A' + B*B') < 0.001

Z,resnorm=try
              rksm(A,B,E,infoV=false)
          catch
              error("rksm failed with full matrices")
          end
@test norm(A*Z*Z'*E + E*Z*Z'*A' + B*B') < 0.001

Z2,resnorm2=try
            rksm(sparse(A),B,sparse(E))
          catch
              error("rksm failed with sparse matrices")
          end
@test norm(A*Z2*Z2'*E + E*Z2*Z2'*A' + B*B') < 0.001

Z,flag,res =try
              lp_lradi(A,B)
            catch
              error("lp_lradi failed with full matricies")
            end
@test norm(A*Z*Z' + Z*Z'*A' + B*B') < 0.001

Z,flag,res =try
              lp_lradi(sparse(A),B,infoV=false)
            catch
              error("lp_lradi failed with sparse matricies")
            end
@test norm(A*Z*Z' + Z*Z'*A' + B*B') < 0.001

max_stars_repo_name: garrettthomaskth/LargeMatrixEquations

max_issues_repo_name: garrettthomaskth/LargeMatrixEquations

max_forks_repo_name: garrettthomaskth/LargeMatrixEquations


text:
module CPUSummary

using Static
using Static: Zero, One, gt, lt
using IfElse: ifelse
export cache_size,
  cache_linesize,
  cache_associativity,
  cache_type,
  cache_inclusive,
  num_cache,
  num_cores,
  num_threads

# const USE_HWLOC = @load_preference("hwloc", Sys.ARCH !== :aarch64 || !Sys.isapple())
# use_hwloc(b) = @set_preferences!("hwloc" => b)

# @static if USE_HWLOC
#   try
#     script = """
#     $(Base.load_path_setup_code())
#     Hwloc = Base.require(Base.PkgId(Base.UUID("0e44f5e4-bd66-52a0-8798-143a42290a1d"), "Hwloc"))
#     Hwloc.gettopology()
#     """
#     p = run(`$(Base.julia_cmd()) -e $(script)`, wait=false)
#     wait(p)
#     if p.exitcode == 0 && p.termsignal == 0
#       include("topology.jl")
#     else
#       use_hwloc(false)
#       include("generic_topology.jl")
#     end
#   catch
#     use_hwloc(false)
#     include("generic_topology.jl")
#   end
# else
if (Sys.ARCH === :x86_64)
  include("x86.jl")
else
  include("generic_topology.jl")
end
function __init__()
  ccall(:jl_generating_output, Cint, ()) == 1 && return
  nc = _get_num_cores()
  syst = Sys.CPU_THREADS::Int
  nt = Threads.nthreads()
  if nc != num_l1cache()
    @eval num_l1cache() = static($nc)
  end
  if nc != num_cores()
    @eval num_cores() = static($nc)
  end
  if syst != sys_threads()
    @eval sys_threads() = static($syst)
  end
  if nt != num_threads()
    @eval num_threads() = static($nt)
  end
  _extra_init()
end

# end
num_cache(::Union{Val{1},StaticInt{1}}) = num_l1cache()
num_cache(::Union{Val{2},StaticInt{2}}) = num_l2cache()
num_cache(::Union{Val{3},StaticInt{3}}) = num_l3cache()
num_cache(::Union{Val{4},StaticInt{4}}) = num_l4cache()
const BASELINE_CORES = Int(num_cores()) * ((Sys.ARCH === :aarch64) && Sys.isapple() ? 2 : 1)
cache_linesize() = cache_linesize(Val(1))
function num_cache_levels()
  numl4 = num_l4cache()
  numl4 === nothing && return nothing
  ifelse(
    eq(numl4, Zero()),
    ifelse(
      eq(num_l3cache(), Zero()),
      ifelse(
        eq(num_l2cache(), Zero()),
        ifelse(eq(num_l1cache(), Zero()), Zero(), One()),
        StaticInt{2}(),
      ),
      StaticInt{3}(),
    ),
    StaticInt{4}(),
  )
end

end


text:
import TensorFlowBuilder.TFParser: greparg, grepreturn, TFFunction, findarg
import TensorFlowBuilder.TFParser
import TensorFlowBuilder.PyInspector
import TensorFlowBuilder.CoreTypes: Tensor

using PyCall
@pyimport tensorflow as tf

q = TFParser.PyInspector.pydoc(tf.random_normal)
s = "Outputs random values from a normal distribution.

  Args:
    shape: A 1-D integer Tensor or Python array. The shape of the output tensor.
    mean: A 0-D Tensor or Python value of type `dtype`. The mean of the normal
      distribution.
    stddev: A 0-D Tensor or Python value of type `dtype`. The standard deviation
      of the normal distribution.
    dtype: The type of the output.
    seed: A Python integer. Used to create a random seed for the distribution.
      See
      [`set_random_seed`](../../api_docs/python/constant_op.md#set_random_seed)
      for behavior.
    name: A name for the operation (optional).

  Returns:
    A tensor of the specified shape filled with random normal values.
  "

@test greparg(s, "shape") == "A 1-D integer Tensor or Python array. The shape of the output tensor."
@test greparg(s, "mean") == "A 0-D Tensor or Python value of type `dtype`. The mean of the normal distribution."

@test greparg(s, "seed") == "A Python integer. Used to create a random seed for the distribution. See [`set_random_seed`](../../api_docs/python/constant_op.md#set_random_seed) for behavior."

@test grepreturn(s) == "A tensor of the specified shape filled with random normal values."

tff = TFFunction(tf.random_normal)

@test tff.name == :random_normal
@test tff.returndoc == grepreturn(s)
@test length(tff.args) == 6
namearg = tff.args[6]
@test !namearg.mandatory
@test namearg.doc == greparg(s, "name")
@test namearg.doc1 == "A name for the operation (optional)."
shapearg = findarg(tff, :shape)
@test shapearg.doc1 == "A 1-D integer Tensor or Python array"

tff2 = TFFunction(tf.slice)
@test TFParser.jlname(findarg(tff2, :begin).name) == "begin_"

tff3 = TFFunction(tf.constant)
@test TFParser.guesstype(:value, findarg(tff3, :value).doc1) == Tensor


text:
module source

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
module Common

export Generator, Facet, Hull, nconic, nfacet, gtype
export get_canonical_winding
export except_index, indexof, isevenperm, isoddperm

export dominates, antidominates
export Face, facesof
export nbof, nbsof, generatorof, generatorsof, set_nb!

export add_facet!, del_facet!, mark_facet!, ismarked

abstract Generator
abstract Facet
abstract Hull

nconic(::None) = error("Unimplemented")

nfacet(x) = nconic(x) - 1

dot(f::None, g::None) = error("Unimplemented")
dominates(    g::Generator, f::Facet) = dot(f, g) > 0
antidominates(g::Generator, f::Facet) = dot(f, g) < 0

gtype(::None)                 = error("Unimplemented")
get_canonical_winding(::None) = error("Unimplemented")

except_index{T}(x::Vector{T}, k::Int) = [x[1:(k-1)]..., x[(k+1):end]...]

function indexof(v::Vector, x)
    for (k, y) in enumerate(v); if y == x;  return k;  end; end
    error("indexof: x = $x not found in v=$v")
end

isoddperm(p) = !isevenperm(p)
function isevenperm{T}(p::Vector{T})
    @assert isperm(p)
    n = length(p)
    used = falses(n)
    even = true
    for k = 1:n
        if used[k]; continue; end
        # Each even cycle flips even (an odd number of times)
        used[k] = true
        j = p[k]
        while !used[j]
            used[j] = true
            j = p[j]
            even = !even
        end
    end
    even
end

nbof(f::None, k::Int) = error("Unimplemented")
nbof(f::Facet, k::Union(Facet,Generator)) = nbof(f, indexof(f, k))

set_nb!(f::None, new_nb::None, k::Int) = error("Unimplemented")
set_nb!(f::Facet, new_nb::Facet, k::Union(Facet,Generator)) = set_nb!(f, new_nb, indexof(f, k))

generatorof(f::None, k::Int) = error("Unimplemented")
generatorof(f::Facet, k::Union(Facet,Generator)) = generatorof(f, indexof(f, k))

nbsof(f::None) = error("Unimplemented")
generatorsof(f::None) = error("Unimplemented")

# ---------------------------------- Face ------------------------------------

immutable Face{F<:Facet}
    parent::F
    k::Int
end
Face{F<:Facet}(parent::F, nb::F) = Face{F}(parent, indexof(parent, nb))
Face{F<:Facet}(parent::F, g::Generator) = Face{F}(parent, indexof(parent, g))

parentof(f::Face)    = f.parent
nbof(f::Face)        = nbof(f.parent, f.k)
generatorof(f::Face) = generatorof(f.parent, f.k)
flip(f::Face)        = (nb = nbof(f); Face(nb, indexof(nb, f)))

set_nb!{F}(f::Face{F}, newnb::F) = set_nb!(f.parent, newnb, f.k)

add_facet!(h::None, face::None, g::None)      = error("Unimplemented")
del_facet!(h::None, facet::None)              = error("Unimplemented")
mark_facet!(h::None, facet::None, mark::Bool) = error("Unimplemented")
mark_facet!(h, facet) = mark_facet!(h, facet, true)
ismarked(h::None, facet::None)                = error("Unimplemented")

# ---------------------------------- Faces -----------------------------------

immutable Faces{F<:Facet}; parent::F; end

Base.length(faces::Faces) = nfacet(faces.parent)
Base.start(faces::Faces) = 1
Base.next{F}(faces::Faces{F}, k) = (Face{F}(faces.parent, k), k+1)
Base.done(faces::Faces, k) = k == (nfacet(faces.parent)+1)
Base.eltype{F}(faces::Faces{F}) = Face{F}

facesof(f::Facet) = Faces(f)

end # module


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

# `ntuple`, for constructing tuples of a given length

"""
    ntuple(f::Function, n::Integer)

Create a tuple of length `n`, computing each element as `f(i)`,
where `i` is the index of the element.

# Examples
```jldoctest
julia> ntuple(i -> 2*i, 4)
(2, 4, 6, 8)
```
"""
@inline function ntuple(f::F, n::Integer) where F
    # marked inline since this benefits from constant propagation of `n`
    t = n == 0  ? () :
        n == 1  ? (f(1),) :
        n == 2  ? (f(1), f(2)) :
        n == 3  ? (f(1), f(2), f(3)) :
        n == 4  ? (f(1), f(2), f(3), f(4)) :
        n == 5  ? (f(1), f(2), f(3), f(4), f(5)) :
        n == 6  ? (f(1), f(2), f(3), f(4), f(5), f(6)) :
        n == 7  ? (f(1), f(2), f(3), f(4), f(5), f(6), f(7)) :
        n == 8  ? (f(1), f(2), f(3), f(4), f(5), f(6), f(7), f(8)) :
        n == 9  ? (f(1), f(2), f(3), f(4), f(5), f(6), f(7), f(8), f(9)) :
        n == 10 ? (f(1), f(2), f(3), f(4), f(5), f(6), f(7), f(8), f(9), f(10)) :
        _ntuple(f, n)
    return t
end

function _ntuple(f::F, n) where F
    @noinline
    (n >= 0) || throw(ArgumentError(string("tuple length should be ≥ 0, got ", n)))
    ([f(i) for i = 1:n]...,)
end

function ntupleany(f, n)
    @noinline
    (n >= 0) || throw(ArgumentError(string("tuple length should be ≥ 0, got ", n)))
    (Any[f(i) for i = 1:n]...,)
end

# inferrable ntuple (enough for bootstrapping)
ntuple(f, ::Val{0}) = ()
ntuple(f, ::Val{1}) = (@inline; (f(1),))
ntuple(f, ::Val{2}) = (@inline; (f(1), f(2)))
ntuple(f, ::Val{3}) = (@inline; (f(1), f(2), f(3)))

"""
    ntuple(f, ::Val{N})

Create a tuple of length `N`, computing each element as `f(i)`,
where `i` is the index of the element. By taking a `Val(N)`
argument, it is possible that this version of ntuple may
generate more efficient code than the version taking the
length as an integer. But `ntuple(f, N)` is preferable to
`ntuple(f, Val(N))` in cases where `N` cannot be determined
at compile time.

# Examples
```jldoctest
julia> ntuple(i -> 2*i, Val(4))
(2, 4, 6, 8)
```
"""
@inline function ntuple(f::F, ::Val{N}) where {F,N}
    N::Int
    (N >= 0) || throw(ArgumentError(string("tuple length should be ≥ 0, got ", N)))
    if @generated
        quote
            @nexprs $N i -> t_i = f(i)
            @ncall $N tuple t
        end
    else
        Tuple(f(i) for i = 1:N)
    end
end

@inline function fill_to_length(t::Tuple, val, ::Val{_N}) where {_N}
    M = length(t)
    N = _N::Int
    M > N && throw(ArgumentError("input tuple of length $M, requested $N"))
    if @generated
        quote
            (t..., $(fill(:val, (_N::Int) - length(t.parameters))...))
        end
    else
        (t..., fill(val, N-M)...)
    end
end

max_stars_repo_name: jonas-schulze/julia

max_issues_repo_name: jonas-schulze/julia

max_forks_repo_name: jonas-schulze/julia


text:
using Compat

if isa(Pkg.installed("IJulia"), VersionNumber)
    using IJulia
    python = strip(readline(IJulia.jupyter), ['\n',' ', '#','!'])
    ipywver = readstring(`$python -c 'import ipywidgets; print ipywidgets.__version__'`) |> strip |> VersionNumber
    info("ipywidgets version found: $ipywver")
    if ipywver < v"5.0.0"
        warn("""This version of Interact requires ipywidgets > 5.0 to work correctly.
                If you have ipywidgets version 4.x, run Pkg.checkout("Interact", "ipywidgets-4").""")
    else
        info("A compatible version of ipywidgets was found. All good.")
    end
else
    warn("""IJulia is not installed. run Pkg.build("Interact")
            once you install IJulia to use Interact inside IJulia""")
end


text:
@deprecate seeded_region_growing(
    img::AbstractArray,
    seeds::AbstractVector,
    kernel_dim::Vector{Int},
    diff_fn::Function = default_diff_fn) seeded_region_growing(img, seeds, (kernel_dim...,), diff_fn)

@deprecate unseeded_region_growing(
    img::AbstractArray,
    threshold::Real,
    kernel_dim::Vector{Int},
    diff_fn::Function = default_diff_fn) unseeded_region_growing(img, threshold, (kernel_dim...,), diff_fn)


text:
filterclass(x::Matrix, y::Vector{Int}, c::Int) = x[:,y .== c]
filterclass(x::Vector, y::Vector{Int}, c::Int) = x[y.== c]
filterclass(x::Tuple{A,B}, c::Int) where {A,B} = filterclass(x..., c)

"""
    detacc_at_fp(n, p, α = [0.005,0.01,0.05])

    detection accuracy at particular false positive rates

"""
detacc_at_fp(n::Vector{T}, p::Vector{T}, α = [0.005,0.01,0.05]) where {T<:Real} = mean(p' .> quantile(n,1 .- α), dims = 2)
detacc_at_fp(x::Vector{T}, y::Vector{Int}, α = [0.005,0.01,0.05]) where {T<:Real} = detacc_at_fp(x[y .== 1], x[y .== 2], α)

"""
    runtest(fit, ps, predicts, prnames,  dataset, anomaly_type, polution, variation, idir, odir, repetition = 1, steps = 50000)

    fit ---
"""
function runtest(fit, ps, predicts, prnames,  dataset, anomaly_type, polution, variation, idir, odir, name, repetition = 1, steps = 50000)
  println("processing knn ",dataset, anomaly_type, polution, variation)
  train, test, clusterdness = makeset(loaddataset(dataset,anomaly_type,idir)..., 0.75, variation)
  idim = size(train[1],1)
  data = RandomBatches((subsampleanomalous(train,polution)[1],), 100, steps)

  results = mapreduce(vcat, ps) do p
    m,info = fit(data,p...)

    #calculate areas under the curve on training and testing data for different prediction functions
    aucs = mapreduce(vcat,zip(predicts, prnames)) do i
      f, prname = i
      df = DataFrame(
        prname = prname,
      	auc_test = auc(f, m, test),
      	auc_train = auc(f, m, train),
      	auc_train005 = auc(f, m, subsampleanomalous(train,0.005,repetition)),
      	auc_train01 = auc(f, m, subsampleanomalous(train,0.01,repetition))
    	)
      hcat(df,fprstats(f, m, train,test,[0.005,0.01,0.05]))
    end
    aucs = join(info, aucs, kind = :cross)
  end

  results[:repetition] = repetition
  results[:clusterdness] = clusterdness
  ofname = joinpath(odir,dataset,@sprintf("%s_%s_%g_%s",name,anomaly_type,polution,variation))
  append2file(ofname*".jld2","auc",results)
  append2file(ofname*".csv",results)
end

"""
    anomalyexperiment(fit, trainstats, teststats, dataset; aparam, dataparts, repetition)

    trainstats --- statistics to be extracted on training data
    teststats --- statistics to be extracted on testing data
    aparam --- (type = "east", polution = "0.0", variation = "low") , variation parameters of generator of anomalies
    dataparts --- (0.75, 0.25) fraction of samples used for training and for testing
    repetition --- used as a seed for random number generator to initialize the data
"""
function anomalyexperiment(fit, trainstats, teststats, dataset; aparam = (type = "easy", polution = 0.0, variation = "low"), dataparts = (0.75, 0.25), repetition = 1)
  println("parameters of anomalies: ", aparam)
  trndata, tstdata, clusterdness = makeset(loaddataset(dataset, a.type ,idir)..., dataparts[1], aparam.variation)

  model = fit(subsampleanomalous(trndata, aparam.polution)[1])
end

"""
  fprstats(f, m, train,test,α)

  puts threshold on false positive rate on training data (assuming most of them are normal)
  and calculates on this threshold false positive rate and detection accuracy

"""
function fprstats(f, m, train,test,α::T) where {T<:Real}
  τ = quantile(f(m,train[1]), 1 - α)
  o = f(m, test[1])
  tstfp = mean( o[test[2][:] .== 1] .- τ .> 0)
  dacc = mean(o[test[2][:] .== 2] .- τ .> 0)
  npscore = max(0,tstfp - α)/α + mean(o[test[2][:] .== 2] .- τ .<= 0)
  names = map(s -> Symbol(replace(@sprintf("%s_%g",s,α),"." => "")),["fpr", "dacc", "dacctst", "npscore"])
  DataFrame([tstfp dacc detacc_at_fp(o,test[2],[α])[1] npscore] , names)
end

fprstats(f, m, train,test,α::Vector) = mapreduce(i -> fprstats(f, m, train, test, i), hcat, α)

auc(predict::Function, data) = auc(roccurve(predict(data[1]), data[2] .- 1)...)
auc(predict::Function, m, data) = auc(roccurve(predict(m, data[1]), data[2] .- 1)...)


text:
using Zygote, Test

D(f, x) = pushforward(f, x)(1)

@test D(x -> sin(cos(x)), 0.5) == -cos(cos(0.5))*sin(0.5)

@test D(x -> D(cos, x), 0.5) == -cos(0.5)

@test D(x -> x*D(y -> x*y, 1), 4) == 8

function pow(x, n)
  r = 1
  while n > 0
    n -= 1
    r *= x
  end
  return r
end

@test D(x -> pow(x, 3), 2) == 12

@test D(1) do x
  f(y) = x = x*y
  D(f, 1)
  D(f, 1)
end == 1

@test D(x -> D(y -> x = y, x)*x, 1) == 1

@test D(1) do x
  D(2) do y
    D(3) do z
      x = z * y
    end
  end
  x
end == 0

@test D(x -> abs(x+2im), 1) == gradient(x -> abs(x+2im), 1)[1]

using LinearAlgebra

@test D(3) do x
  A = zeros(5, 5)
  B = zeros(5, 5)
  A[1, 1] = x
  mul!(B, A, A)
  sum(B)
end == 6


text:
module TestDoctest

using Documenter
using Kaleido
using Test

@testset "doctest" begin
    doctest(Kaleido; manual = true)
end

end  # module


text:
using Test

using POMDPModels
# using POMDPSimulators
using POMDPTesting
using POMDPs
using POMDPModelTools
using BeliefUpdaters
using Random

let
    problem = BabyPOMDP()

    # starve policy
    # when the baby is never fed, the reward for starting in the hungry state should be -100
    sim = RolloutSimulator(eps=0.0001)
    ib = nothing
    policy = Starve()
    r = simulate(sim, problem, policy, updater(policy), ib, true)
    @test r ≈ -100.0 atol=0.01

    # test gen(::o,...)
    o = gen(DDNNode(:o), problem, true, MersenneTwister(1))
    @test o == 1
    # test vec
    ov = convert_s(Array{Float64}, true, problem)
    @test ov == [1.]
    o = convert_s(Bool, ov, problem)
    @test o == true

    POMDPTesting.probability_check(problem)

    bu = DiscreteUpdater(problem)
    bp =  update(bu,
                 initialize_belief(bu, BoolDistribution(0.0)),
                 false,
                 true)

    @test pdf(bp, true) ≈ 0.47058823529411764 atol=0.0001
    r = simulate(sim, problem, policy, DiscreteUpdater(problem), BoolDistribution(1.0))
    @test r ≈ -100.0 atol=0.01
end


text:
# Usage: `julia <script>.jl {gpu,cpu} {strong,weak}`
# e.g. `julia <script>.jl gpu strong
if !@isdefined USE_GPU
    const USE_GPU = length(ARGS) > 0 && ARGS[1] == "gpu"
end

using DataFrames
using CSV
using ImplicitGlobalGrid
import MPI
using ParallelStencil
using ParallelStencil.FiniteDifferences3D
if !ParallelStencil.is_initialized()
    @static if USE_GPU
        @init_parallel_stencil(CUDA, Float64, 3)
    else
        @init_parallel_stencil(Threads, Float64, 3)
    end
end
include("part1_kernel_programming.jl")

MPI.Init()

# Tolerance experiments reveal that 1e-6 gives sufficiently accurate results for this problem.
const SOLVER_TOL = 1e-6

for strong_scaling in [true, false]
    # We need 7 buffers for the solver
    # memory(nx) = nx^3 * sizeof(Float64) * 7
    # => with nx=ny=nz=2^9 we have memory(nx) ≈ 7.5GB => Fits on Titan X (12 GB)
    # After experimenting with 2^9 we found it very hard to converge.
    # Empirically, nx=ny=nz=2^7 gives a good trade-off between convergence speed and precision.
    nx, ny, nz =
        if strong_scaling
            dims_dict = Dict(
                             1 => (1, 1, 1),
                             2 => (2, 1, 1),
                             4 => (2, 2, 1),
                             8 => (2, 2, 2),
                            )
            2^7 .÷ dims_dict[MPI.Comm_size(MPI.COMM_WORLD)]  # global, gets divided
        else
            2^7 .* (1, 1, 1)  # local, gets solved on every node
        end

    filename = joinpath("benchmark-results", "bench_diffusion_scaling_$(USE_GPU ? "gpu" : "cpu").csv")

    for use_shared_memory in [true, false]
        _, _, bench_res = diffusion_3D_kernel_programming(nx=nx, ny=ny, nz=nz, ttot=2.0, tol=SOLVER_TOL, use_shared_memory=use_shared_memory, init_and_finalize_MPI=false, verbose=true, scale_physical_size=(strong_scaling ? false : true))
        bench_row = (delta_t=bench_res.Δt,
                     Work=bench_res.Work,
                     Performance=bench_res.Performance,
                     Memory=bench_res.Memory,
                     Intensity=bench_res.Intensity,
                     Throughput=bench_res.Throughput,
                     use_shared_memory=use_shared_memory,
                     use_gpu=USE_GPU,
                     strong_scaling=strong_scaling,
                     n_threads=Threads.nthreads(),
                     n_mpi_ranks=MPI.Comm_size(MPI.COMM_WORLD))

        # Note that this is potentially susceptible to a race condition
        if MPI.Comm_rank(MPI.COMM_WORLD) == 0
            bench_df =
                if isfile(filename)
                    DataFrame(CSV.File(filename))
                else
                    DataFrame()
                end
            push!(bench_df, bench_row)

            open(filename, "w") do io
                CSV.write(io, bench_df)
            end
        end
    end
end

MPI.Finalize()


text:
"""
    VARFilter(x,A,z0)

Create z Txn matrix from VAR model where z[t,:] = A1*z[t-1,:] +...+Ap*z[t-p,:]+ x[t,:].
A is an nxnxp array with cat(A1,A2,...,dims=3)
z0 is pxn initial values of z (for [t=-2;t=-1;t=0] for a VAR(3))
"""
function VARFilter(x,A,z0)

    p     = size(A,3)            #lag order
    (T,n) = (size(x,1),size(x,2))
    T0    = size(z0,1)
    (p != T0) && error("z0 must be $p x $n")

    z     = [z0;zeros(T,n)]
    x     = [zeros(p,n);x]
    for t = p+1:size(z,1)        #loop over t
        for s = 1:p              #loop over lags, z[t-1],z[t-2],...,z[t-p]
            z[t,:] = z[t,:] + A[:,:,s]*z[t-s,:]  #works since z[t,:] is a vector
        end
        z[t,:] = z[t,:] + x[t,:]
    end
    #printmat(z)
    z = z[p+1:end,:]                #cut padded part

    return z

end

max_stars_repo_name: duzaichuan/FinancialEconometrics

max_issues_repo_name: duzaichuan/FinancialEconometrics

max_forks_repo_name: duzaichuan/FinancialEconometrics


text:
# =================================================================
# Laub-Loomis model
# See https://easychair.org/publications/paper/gjfh
# =================================================================

using Reachability: Options, SafeStatesProperty
using MathematicalSystems, LazySets
using DynamicPolynomials, SemialgebraicSets

# ==============================
# Load model
# ==============================

"""
    laubloomis(; [T], [X0], [W])

Construct the Laub-Loomis model.

### Input

- `T`  --  (optional, default: `20.0`) the time horizon for the initial-value
            problem
- `W`  --  (optional, default: `0.01`) width of the initial states
- `X0` --  (optional, default: an axis-aligned box centered at
           `(1.2, 1.05, 1.5, 2.4, 1.0, 0.1, 0.45)` and radius `W`) set of initial states
- `unsafe_bound` -- (optional, default: `4.5`) bound for the variable `x₄` that
                    defines the unsafe region, of the form `x₄ ≥ unsafe_bound`
- `variables`    -- (optional, default: `PolyVar`) the set of variables used to
                    describe the polynomial ODE

### Output

The tuple `(𝑃, 𝑂)` where `𝑃` is an initial-value problem and `𝑂` are the options.
"""
function laubloomis(; T=20.0,
                      W=0.01,
                      X0=BallInf([1.2, 1.05, 1.5, 2.4, 1.0, 0.1, 0.45], W),
                      unsafe_bound=4.5,
                      variables=@polyvar x₁ x₂ x₃ x₄ x₅ x₆ x₇)

    𝑂 = Options()

    # unrwap the variables
    x₁, x₂, x₃, x₄, x₅, x₆, x₇ = variables
    𝑂[:variables] = variables
    𝑂[:vars] = [1:7;]

    # instantiate the polynomial system
    f = [1.4x₃ - 0.9x₁,
         2.5x₅ - 1.5x₂,
         0.6x₇ - 0.8x₂*x₃,
         2 - 1.3x₃*x₄,
         0.7x₁ - x₄*x₅,
         0.3x₁ - 3.1x₆,
         1.8x₆ - 1.6x₂*x₇]

    𝐹 = PolynomialContinuousSystem(f)

    # instantiate the IVP
    𝑃 = InitialValueProblem(𝐹, X0)

    # time horizon
    𝑂[:T] = T

    # variables to plot
    𝑂[:plot_vars] = [0, 4]

    # safety property
    𝑂[:property] = SafeStatesProperty(HalfSpace([0, 0, 0, -1., 0, 0, 0], -unsafe_bound))
    # @set x₄ ≥ 0.01, vars=(x₁, x₂, x₃, x₄, x₅, x₆, x₇)
    return (𝑃, 𝑂)
end

max_stars_repo_name: Aadesha/ReachabilityBenchmarks

max_issues_repo_name: Aadesha/ReachabilityBenchmarks

max_forks_repo_name: Aadesha/ReachabilityBenchmarks


text:
"""
    Initializeclassy(w0waCDMCosmology::w0waCDMCosmology)

This function, given a [`w0waCDMCosmology`](@ref), returns the
[`classyParams`](@ref) correctly initialized.
"""
function Initializeclassy(w0waCDMCosmology::w0waCDMCosmology)
    classyParamsDict = Dict("output" => "mPk",
    "non linear"=> "halofit",
    "Omega_b"=> w0waCDMCosmology.ΩB,
    "Omega_cdm"=> w0waCDMCosmology.ΩM-w0waCDMCosmology.ΩB-
    w0waCDMCosmology.Mν/(93.14*(w0waCDMCosmology.H0/100)^2),
    "N_ur"=> 2.0328,
    "h"=> w0waCDMCosmology.H0/100.,
    "sigma8" => w0waCDMCosmology.σ8,
    "n_s" => w0waCDMCosmology.ns,
    "m_ncdm" => w0waCDMCosmology.Mν,
    "P_k_max_1/Mpc" => 50.,
    "z_max_pk" =>  5.,
    "use_ppf" =>  "yes",
    "w0_fld" =>  w0waCDMCosmology.w0,
    "Omega_k" =>  w0waCDMCosmology.Ωk,
    "Omega_fld" =>  w0waCDMCosmology.ΩDE,
    "wa_fld" =>  w0waCDMCosmology.wa,
    "cs2_fld" =>  1.,
    "N_ncdm" =>  1,
    "tau_reio" =>  0.058)
    classyparams = classyParams(classyParamsDict = classyParamsDict)
    return classyparams
end

function Initializeclassy(Flatw0waCDMCosmology::Flatw0waCDMCosmology)
    classyParamsDict = Dict("output" => "mPk",
    "non linear"=> "halofit",
    "Omega_b"=> Flatw0waCDMCosmology.ΩB,
    "Omega_cdm"=> Flatw0waCDMCosmology.ΩM-Flatw0waCDMCosmology.ΩB-
    Flatw0waCDMCosmology.Mν/(93.14*(Flatw0waCDMCosmology.H0/100)^2),
    "N_ur"=> 2.0328,
    "h"=> Flatw0waCDMCosmology.H0/100.,
    "sigma8" => Flatw0waCDMCosmology.σ8,
    "n_s" => Flatw0waCDMCosmology.ns,
    "m_ncdm" => Flatw0waCDMCosmology.Mν,
    "P_k_max_1/Mpc" => 50.,
    "z_max_pk" =>  5.,
    "use_ppf" =>  "yes",
    "w0_fld" =>  Flatw0waCDMCosmology.w0,
    "Omega_k" =>  0,
    "Omega_fld" =>  1-Flatw0waCDMCosmology.ΩM,
    "wa_fld" =>  Flatw0waCDMCosmology.wa,
    "cs2_fld" =>  1.,
    "N_ncdm" =>  1,
    "tau_reio" =>  0.058)
    classyparams = classyParams(classyParamsDict = classyParamsDict)
    return classyparams
end

"""
    EvaluatePowerSpectrum!(classyParams:: classyParams, CosmologicalGrid::CosmologicalGrid,
    PowerSpectrum::PowerSpectrum)

This function runs classy to evaluate the Matter Power Spectrum over the ``k-z``
grid specified in [`CosmologicalGrid`](@ref)
"""
function EvaluatePowerSpectrum!(classyParams:: classyParams,
    CosmologicalGrid::CosmologicalGrid, PowerSpectrum::PowerSpectrum)
    cosmo = classy.Class()
    cosmo.set(classyParams.classyParamsDict)
    cosmo.compute()
    for (idxz, myz) in enumerate(CosmologicalGrid.ZArray)
        for (idxk, myk) in enumerate(CosmologicalGrid.KArray)
            PowerSpectrum.PowerSpectrumLinArray[idxk, idxz] =
            cosmo.pk_lin(myk, myz)
            PowerSpectrum.PowerSpectrumNonlinArray[idxk, idxz] =
            cosmo.pk(myk, myz)
        end
    end
end


text:
"""
    rename!(tier, name)
    rename!(tg, num, name)
    rename!(tg, name_orig, name)

Rename a `Tier` in place. The Tier can be passed as a `Tier` object or identified
by the number or the unique name of the tier.
"""
function rename!(tier::Tier, name::AbstractString)
    tier.name = name
    tier
end

rename!(tg::TextGrid, num::Int, name::AbstractString) = rename!(tg[num], name)

function rename!(tg::TextGrid, name_orig::AbstractString, name::AbstractString)
    tier_num = findfirst(x -> x.name == name_orig, tg)
    tg[tier_num].name = name
    tg
end

"""
    relabel!(interval, label)
    relabel!(tier, index, label)
    relabel!(tg, num, index, label)

Relabel an interval (or a point) in place. The Interval can be passed as an
`Interval` object or identified by its `index` from a `Tier` object. The `Tier`
itself can also be identified by its number in a `TextGrid`.
"""
function relabel!(interval::Interval, label::AbstractString)
    interval.label = label
    interval
end

function relabel!(tier::Tier, index::Int, label::AbstractString)
    interval_num = findfirst(x -> x.index == index, tier.contents)
    tier.contents[interval_num].label = label
    tier
end

relabel!(tg::TextGrid, num::Int, index::Int, label::AbstractString) = relabel!(tg[num], index, label)

"""
    relabel!(tier, indices, labels)
    relabel!(tier, labels)

Relabel multiple annotations in a Tier in place.
"""
function relabel!(tier::Tier, indices, labels)
    for (index, label) in zip(indices, labels)
        relabel!(tier, index, label)
    end
end

function relabel!(tier::Tier, labels)
    for (index, label) in enumerate(labels)
        relabel!(tier, index, label)
    end
    tier
end

"""
    insert_tier!(tg, name, num = 0; class)

Insert a Tier of `class` (interval or point) called `name` into TextGrid `tg` at
position `num`. Tier is inserted at the end of `tg` if `num` is unspecified or
outside the number of tiers in `tg`.

`insert_interval_tier!` and `insert_point_tier!` are class-specific wrappers.
"""
function insert_tier!(tg::TextGrid, name::AbstractString, num::Int = 0; class::AbstractString)
    class ∈ ("interval", "point") || error("Tier class not recognised.")
    num in 1:length(tg) || (num = length(tg) + 1)

    insert!(tg, num, Tier(num, class, name, 1, extrema(tg)..., Vector{Interval}([Interval(1, extrema(tg)..., "")])))
    num < length(tg) && renumber_tiers!(tg, from = num + 1)
    return tg
end

insert_interval_tier!(tg::TextGrid, name::AbstractString, num::Int = 0) = insert_tier!(tg, name, num, class = "interval")
insert_point_tier!(tg::TextGrid, name::AbstractString, num::Int = 0) = insert_tier!(tg, name, num, class = "point")

"""
    extract_tier(tg, tiers...)

Extract one or more Tiers from the input TextGrid `tg` and output another
`TextGrid`. Each `Tier` can be identified by number or by name, and a mix of
numbers and names can be used. Any argument that cannot be used to locate a
Tier is simply discarded. When multiple Tiers are extracted, they are
automatically renumbered.

This function can also be used to duplicate Tiers, by specifying the same Tier
multiple times.
"""
function extract_tier(tg::TextGrid, tiers...)
    tg_new = TextGrid()
    for tier in tiers
        if tier isa Integer
            tier_num = num
        elseif tier isa String
            tier_num = findfirst(x -> x.name == tier, tg)
        else
            tier_num = nothing
        end

        if !isnothing(tier_num)
            push!(tg_new, Tier(length(tg_new) + 1, tg[tier_num].class, tg[tier_num].name, tg[tier_num].size, tg[tier_num].xmin, tg[tier_num].xmax, tg[tier_num].contents))
        end
    end
    tg_new
end

"""
    remove_tier!(tg, tiers...)

Remove one or more Tiers from the TextGrid `tg`. Each `Tier` can be identified
by number or by name, and a mix of numbers and names can be used. Any argument
that cannot be used to locate a Tier is simply discarded.
"""
function remove_tier!(tg::TextGrid, tiers...)
    for tier in tiers
        if tier isa Integer
            deleteat!(tg, findfirst(x -> x.num == tier, tg))
        elseif tier isa String
            deleteat!(tg, findfirst(x -> x.name == tier, tg))
        end
    end
    renumber_tiers!(tg)
end

"""
    remove_empty_tiers!(tg; empty_labels = true)

Remove all empty Tiers from the TextGrid `tg`.

The keyword argument `empty_labels` (default `true`) determines if intervals or
points with empty labels still count as intervals or points. If `true`, only a
Point Tier with no points, or an Interval Tier with only a single empty Interval
will count as empty.
"""
function remove_empty_tiers!(tg::TextGrid; empty_labels::Bool = true)
    tiers = findall(x -> is_empty(x, empty_labels = empty_labels), tg)
    remove_tier!(tg, tiers...)
end

"""
is_empty(tier; empty_labels = true)

Check if `tier` is empty.

The keyword argument `empty_labels` (default `true`) determines if intervals or
points with empty labels still count as intervals or points. If `true`, only a
Point Tier with no points, or an Interval Tier with only a single empty Interval
will count as empty.
"""
function is_empty(tier::Tier; empty_labels::Bool = true)
    if tier.class == "point"
        tier.size == 0 && return true
    else
        if tier.size == 1
            length(tier.contents) == 0 && return true
            strip(tier.contents[1].label) == "" && return true
        end
    end

    if !empty_labels
        for interval in tier.contents
            strip(interval.label) == "" || return false
        end
        return true
    end
    return false
end

function resize_tier!(tier::Tier)
    if tier.class == "point"
        return length(tier.contents)
    end

    t = 0.0
    size = 0
    for interval in tier.contents
        if t < interval.xmin
            size += 1
        end
        size += 1
        interval.index = size
        t = interval.xmax
    end
    if t < tier.xmax
        size += 1
    end
    tier.size = size
    tier
end

resize_tier(tier::Tier) = resize_tier!(deepcopy(tier))

function renumber_tiers!(tg::TextGrid; from::Int = 1)
    for (i, tier) in zip(Iterators.countfrom(from), tg[from:end])
        tier.num = i
    end
    return tg
end

function reindex_intervals!(tier::Tier; from::Int = 1)
    for (i, interval) in zip(Iterators.countfrom(from), tier.contents[from:end])
        interval.index = i
    end
    return tier
end

reindex_intervals!(tg::TextGrid, num::Int) = reindex_intervals!(tg[num])


text:
include("arb_hypgeom.jl")
import .arb_hypgeom

using ArbNumerics
using QuadGK

"""
This file checks the binomial expansion of the denominators of the integrals:

    ∫_0^1 cosh(zt) / (4a^2/β^2 - t^2)^(n+3/2) dt
    ∫_0^1 sinh(zt) / (4a^2/β^2 - t^2)^(n+3/2) dt

which should give:

    ∑_{m=0}^∞ C(-n-3/2, m) (-1)^m (β/2a)^{2n+2m+3} ∫_0^1 cosh(zt) t^{2m} dt
    ∑_{m=0}^∞ C(-n-3/2, m) (-1)^m (β/2a)^{2n+2m+3} ∫_0^1 sinh(zt) t^{2m} dt

These are terms in hyperbolic_integral_five() from check_hyperbolic_integral.jl that depend on the iterator m of the second binomial expansion in check_hyperbolic_integral.jl (the positive integer n is the iterator from the first binomial expansion). This is reproduced here as we then want to check if substituting:

    ∫_0^1 cosh(zt) t^{2m} dt = 1F2(m+1/2; {1/2, m+3/2}; z^2/4)/(2m+1)
                             = ∑_{t=0}^∞ z^{2t} / ((2t+2m+1)⋅(2t)!)

    ∫_0^1 sinh(zt) t^{2m} dt = z ⋅ 1F2(m+1; {3/2, m+2}; z^2/4)/(2m+2)
                             = ∑_{t=0}^∞ z^{2t+1} / ((2t+2m+2)⋅(2t+1)!)

i.e. one_f_two() or one_f_two_fast() from arb_hypgeom.jl, accurately numerically reproduces the top two integrals. The reason for doing this over the expansions already done in check_hyperbolic_integral.jl (which is only as accurate as Float64) is that we want to ensure these integrals are evaluated accurately to precisions greater than Float64 (beyond the first ~16 digits) using the ArbReal summation algorithm first used in bessel_minus_struve.jl, and to hopefully produce cleaner, more digestible code by separating out all these nested sums rather than writing them out explicitly as one massive confusing functions.

Also, one_f_two() versus one_f_two_fast() are compared for efficiency (@time) and accuracy to the top two integrals.
"""

"""
This function uses quadgk to evalulate the integrals:

    ∫_0^1 cosh(zt) / (4a^2/β^2 - t^2)^(n+3/2) dt
    ∫_0^1 sinh(zt) / (4a^2/β^2 - t^2)^(n+3/2) dt.

In this function, h=0 gives the cosh version and h=1 gives the sinh version.
"""
function hyperbolic_integral(z, n, β, a, h)
    if h == 0
        integrand = function (t)
            cosh(z * t) / (4 * a^2 / β^2 - t^2)^(n + 3 // 2)
        end
    elseif h == 1
        integrand = function (t)
            sinh(z * t) / (4 * a^2 / β^2 - t^2)^(n + 3 // 2)
        end
    end
    integral = quadgk(t -> integrand(t), 0, 1)[1]
end

"""
This function then implements the binomial expansion of the denominator of the integrals:

    ∫_0^1 cosh(zt) / (4a^2/β^2 - t^2)^(n+3/2) dt
    ∫_0^1 sinh(zt) / (4a^2/β^2 - t^2)^(n+3/2) dt

i.e.:

    ∑_{m=0}^∞ C(-n-3/2, m) (-1)^m (β/2a)^{2n+2m+3} ∫_0^1 cosh(zt) t^{2m} dt
    ∑_{m=0}^∞ C(-n-3/2, m) (-1)^m (β/2a)^{2n+2m+3} ∫_0^1 sinh(zt) t^{2m} dt.

Here h=0 is the cosh version, h=1 is the sinh version.
"""
function hyperbolic_integral_mexpansion(z, n, β, a, h; prec = 64)
# h = 0 gives cosh, h = 1 gives sinh.

    # Initialise precision of ArbReal to prec.
    setextrabits(0)
    setprecision(ArbReal, prec + 8)

    z = ArbReal("$z")
    n = ArbReal("$n")
    β = ArbReal("$β")
    a = ArbReal("$a")

    m = ArbReal("0")
    result = ArbReal("0.0")
    err = eps(result)  # Machine accuracy of specified precision prec.

    m_binomial_coeff(m) = ArbReal(ArbNumerics.gamma(-n - 1//2) / (ArbNumerics.gamma(m + 1) * ArbNumerics.gamma(-n - m - 1//2)))

    if h == 0
        integrand = function (m, t)
            cosh(z * t) * t^(2 * m)
        end
    elseif h == 1
        integrand = function (m, t)
            sinh(z * t) * t^(2 * m)
        end
    end

    integral(m) = ArbReal(quadgk(t -> integrand(m, t), 0, 1)[1])

    while true
        term = ArbReal(m_binomial_coeff(m) * (-1)^m * (β / (2 * a))^(2 * n + 2 * m + 3) * integral(m))

        # Break loop if term smaller than accuracy of result.
        if abs(term) < err
            break
        end

        result += term
        m += ArbReal("1")

        # Double precision if rounding error in result exceeds accuracy specified by prec.
        if ball(result)[2] > err
            setprecision(ArbReal, precision(result) * 2)

            n = ArbReal("$n")
            z = ArbReal("$z")
            β = ArbReal("$β")
            a = ArbReal("$a")

            m = ArbReal("0")
            result = ArbReal("0.0")

            m_binomial_coeff(m) = ArbReal(ArbNumerics.gamma(-n - 1//2) / (ArbNumerics.gamma(m + 1) * ArbNumerics.gamma(-n - m - 1//2)))

            if h == 0
                integrand = function (m, t)
                    cosh(z * t) * t^(2 * m)
                end
            elseif h == 1
                integrand = function (m, t)
                    sinh(z * t) * t^(2 * m)
                end
            end
        end
    end
    ArbReal(result, bits = prec + 8)
end

"""
This function then implements the binomial expansion of the denominator of the integrals:

    ∫_0^1 cosh(zt) / (4a^2/β^2 - t^2)^(n+3/2) dt
    ∫_0^1 sinh(zt) / (4a^2/β^2 - t^2)^(n+3/2) dt

i.e.:

    ∑_{m=0}^∞ C(-n-3/2, m) (-1)^m (β/2a)^{2n+2m+3} ∫_0^1 cosh(zt) t^{2m} dt
    ∑_{m=0}^∞ C(-n-3/2, m) (-1)^m (β/2a)^{2n+2m+3} ∫_0^1 sinh(zt) t^{2m} dt.

but also substitutes the integrals:

    ∫_0^1 cosh(zt) t^{2m} dt = ∑_{t=0}^∞ z^{2t} / ((2t+2m+1)⋅(2t)!)
    ∫_0^1 sinh(zt) t^{2m} dt = ∑_{t=0}^∞ z^{2t+1} / ((2t+2m+2)⋅(2t+1)!)

by replacing those integrals in the previous function with one_F_two_fast() from arb_hypgeom.jl. h=0 gives cosh version, h=1 gives sinh version.
"""
function hyperbolic_hypergeom_mexpansion(z, n, β, a, h; prec = 64)
# h = 0 gives cosh, h = 1 gives sinh.

    # Initialise precision of ArbReal to prec.
    setextrabits(0)
    setprecision(ArbReal, prec + 8)

    z = ArbReal("$z")
    n = ArbReal("$n")
    β = ArbReal("$β")
    a = ArbReal("$a")

    m_binomial_coeff(m) = ArbReal(ArbNumerics.gamma(-n - 1//2) / (ArbNumerics.gamma(m + 1) * ArbNumerics.gamma(-n - m - 1//2)))

    S(m) = ArbReal(arb_hypgeom.one_f_two_fast(z, m, h; prec = prec))

    m = ArbReal("0")
    result = ArbReal("0.0")
    err = eps(result)  # Machine accuracy of specified precision prec.

    while true
        term = ArbReal(m_binomial_coeff(m) * (-1)^m * (β / (2 * a))^(2 * n + 2 * m + 3) * S(m))

        # Break loop if term smaller than accuracy of result.
        if abs(term) < err
            break
        end
        result += term
        m += ArbReal("1")

        # Double precision if rounding error in result exceeds accuracy specified by prec.
        if ball(result)[2] > err
            setprecision(ArbReal, precision(result) * 2)
            n = ArbReal("$n")
            z = ArbReal("$z")
            β = ArbReal("$β")
            a = ArbReal("$a")

            m = ArbReal("0")
            result = ArbReal("0.0")

            m_binomial_coeff(m) = ArbReal(ArbNumerics.gamma(-n - 1//2) / (ArbNumerics.gamma(m + 1) * ArbNumerics.gamma(-n - m - 1//2)))

            S(m) = ArbReal(arb_hypgeom.one_f_two_fast(z, m, h; prec = precision(result)))
        end
    end
    ArbReal(result, bits = prec + 8)
end

"""
This function then implements the binomial expansion of the denominator of the integrals:

    ∫_0^1 cosh(zt) / (4a^2/β^2 - t^2)^(n+3/2) dt
    ∫_0^1 sinh(zt) / (4a^2/β^2 - t^2)^(n+3/2) dt

i.e.:

    ∑_{m=0}^∞ C(-n-3/2, m) (-1)^m (β/2a)^{2n+2m+3} ∫_0^1 cosh(zt) t^{2m} dt
    ∑_{m=0}^∞ C(-n-3/2, m) (-1)^m (β/2a)^{2n+2m+3} ∫_0^1 sinh(zt) t^{2m} dt.

but also substitutes the integrals:

    ∫_0^1 cosh(zt) t^{2m} dt = 1F2(m+1/2; {1/2, m+3/2}; z^2/4)/(2m+1)
    ∫_0^1 sinh(zt) t^{2m} dt = z ⋅ 1F2(m+1; {3/2, m+2}; z^2/4)/(2m+2)

by replacing those integrals with one_F_two() from arb_hypgeom.jl. h=0 gives cosh version, h=1 gives sinh version.
"""
function hyperbolic_hypergeom_mexpansion2(z, n, β, a, h; prec = 64)
# h = 0 gives cosh, h = 1 gives sinh.

    # Initialise precision of ArbReal to prec.
    setextrabits(0)
    setprecision(ArbReal, prec + 8)

    z = ArbReal("$z")
    n = ArbReal("$n")
    β = ArbReal("$β")
    a = ArbReal("$a")

    m = ArbReal("0")
    result = ArbReal("0.0")
    err = eps(result)  # Machine accuracy of specified precision prec.

    m_binomial_coeff(m) = ArbReal(ArbNumerics.gamma(-n - 1//2) / (ArbNumerics.gamma(m + 1) * ArbNumerics.gamma(-n - m - 1//2)))

    if h == 0
        one_f_two = function (m)
            ArbReal(arb_hypgeom.one_f_two(m + 1/2, (1/2, m + 3/2), (z / 2)^2; prec = prec) / (2 * m + 1))
        end
    elseif h == 1
        one_f_two = function (m)
            ArbReal(z * arb_hypgeom.one_f_two(m + 1, (3/2, m + 2), (z / 2)^2; prec = prec) / (2 * m + 2))
        end
    end

    while true
        term = ArbReal(m_binomial_coeff(m) * (-1)^m * (β / (2 * a))^(2 * n + 2 * m + 3) * one_f_two(m))

        # Break loop if term smaller than accuracy of result.
        if abs(term) < err
            break
        end

        result += term
        m += ArbReal("1")

        # Double precision if rounding error in result exceeds accuracy specified by prec.
        if ball(result)[2] > err
            setprecision(ArbReal, precision(result) * 2)

            n = ArbReal("$n")
            z = ArbReal("$z")
            β = ArbReal("$β")
            a = ArbReal("$a")

            m = ArbReal("0")
            result = ArbReal("0.0")

            m_binomial_coeff(m) = ArbReal(ArbNumerics.gamma(-n - 1//2) / (ArbNumerics.gamma(m + 1) * ArbNumerics.gamma(-n - m - 1//2)))

            if h == 0
                one_f_two = function (m)
                    ArbReal(arb_hypgeom.one_f_two(m + 1/2, (1/2, m + 3/2), (z / 2)^2; prec = precision(result)) / (2 * m + 1))
                end
            elseif h == 1
                one_f_two = function (m)
                    ArbReal(z * arb_hypgeom.one_f_two(m + 1, (3/2, m + 2), (z / 2)^2; prec = precision(result)) / (2 * m + 2))
                end
            end
        end
    end
    ArbReal(result, bits = prec + 8)
end

"""
Plot and time each of these stages/implementations of the integrals and compare.
"""
using Plots
plotly()
Plots.PlotlyBackend()

n = 3
a = 6.23
β = 3.2
z_range = 0.01:1:200 # z=0 diverges
h = 1

# For n=3, a=6.23, β=3.2, z_range=0.01:1:200, h=1 I get:
# 0.161909 seconds (564.27 k allocations: 23.958 MiB, 5.98% gc time) for me.
@time original_integral = [hyperbolic_integral(z, n, β, a, h) for z in z_range]
# 15.588500 seconds (69.64 M allocations: 3.085 GiB, 13.72% gc time) for me.
@time m_expansion = [hyperbolic_integral_mexpansion(z, n, β, a, h) for z in z_range]
# Below uses one_F_two_fast()
# 150.126938 seconds (652.59 M allocations: 28.175 GiB, 15.40% gc time) for me.
@time hypgeom_expansion = [hyperbolic_hypergeom_mexpansion(z, n, β, a, h; prec = 64) for z in z_range]
# Below uses one_F_two()
# 322.089672 seconds (988.47 M allocations: 66.587 GiB, 7.72% gc time) for me.
@time hypgeom_expansion2 = [hyperbolic_hypergeom_mexpansion2(z, n, β, a, h; prec = 64) for z in z_range]
p = plot(z_range, hypgeom_expansion, yaxis=:log, label="1f2 expansion")
plot!(z_range, m_expansion, label="m expansion")
plot!(z_range, original_integral, label="original")
display(p) # plots all look identical.
@show(original_integral, m_expansion, hypgeom_expansion, hypgeom_expansion2)
# original_integral: [2.839772599386279e-8 ... 4.443248342101801e78]
# m_expansion: ArbReal{72}[2.8397725993862784615e-8 ... 4.44324834210182985933e+78]
# hypgeom_expansion = ArbReal{72}[2.8397725993862784742e-8 ... 4.44324834210182176148e+78]
# hypgeom_expansion2 = ArbReal{72}[2.8397725993862784741e-8 ... 4.4432483421018217615e+78]


text:
## PARTE 1
cd(@__DIR__)
conv(c) = c == '.' ? 0 : c == '#' ? 1 : missing
#init = conv.(reduce(hcat, collect.(eachline("input.txt"))))
init = conv.([
    [ '.' '#' '.' ];
    [ '.' '.' '#' ];
    [ '#' '#' '#' ]
])
sz = 3 * max(size(init)...)
cells = zeros(Int, sz, sz, sz)
cells[1:size(init, 1),1:size(init, 2),1] .= init;
new_cells = similar(cells);
start = deepcopy(cells);

##
wrap(x, y, z; sz=size(cells)) = CartesianIndex((mod(i, 1:si) for (i,si) in zip((x, y, z), sz))...)
function neighbours(c; sz=size(cells))
    x, y, z = c.I .- 1
    Δs = (-1, 0, 1)
    [wrap(x+a+1, y+b+1, z+c+1; sz) for a in Δs, b in Δs, c in Δs if (a,b,c) != (0,0,0)] |> vec
end
isactive(c) = c == 1
function conway!(new_cells, cells)
    @assert size(new_cells) == size(cells)
    _isactive(c) = isactive(cells[c])
    foreach(CartesianIndices(cells)) do c
        state = cells[c]
        ns = neighbours(c; sz=size(cells))
        num = count(_isactive, ns)
        cells[c] = if _isactive(c)
            num != 2 && num != 3 ? 0 : state
        else
            num == 3 ? 1 : state
        end
    end
end

##
cells = deepcopy(start)
new_cells .= 0
println("Epoch (0): active cells $(count(isactive, cells))")
foreach(1:6) do epoch
    global cells
    global new_cells
    conway!(new_cells, cells)
    new_cells, cells = cells, new_cells
    println("Epoch ($epoch): active cells $(count(isactive, cells))")
end


text:
type Dot <: Op; Dot(;o...)=new(); end

ninputs(::Dot)=2
canoverwrite(::Dot)=false
back_reads_x(::Dot)=true
back_reads_y(::Dot)=false

function forw(::Dot, x1, x2, y; o...)
    (y === x1 || y === x2) && error("No array sharing in dot.")
    (x1 == nothing || x2 == nothing) && return nothing
    A_mul_B!(y, x1, x2)
end

function back(::Dot, dy, dx1, dx2; x=nothing, o...)
    if dx1 != nothing
        (dx1 === dy || dx1 === x[2]) && error("No array sharing in dot.")
        (x[2] != nothing ? A_mul_Bt!(dx1, dy, x[2]) : fillsync!(dx1, 0))
    end
    if dx2 != nothing
        (dx2 === dy || dx2 === x[1]) && error("No array sharing in dot.")
        (x[1] != nothing ? At_mul_B!(dx2, x[1], dy) : fillsync!(dx2, 0))
    end
end

function infersize(d::Dot,a,b,c)
    a==b==c==nothing && return nothing
    # a,b may have more than 2 dims (convolution tensors etc.), in which case we group the first n-1 into a super-column
    # or they could have 1 dim (bias), in which case we assume a column vector
    # a,b can also be nothing
    dotsize(x)=(length(x)==1 ? (x[1],1) : (prod(x[1:end-1]),x[end]))
    a==nothing || (a=[a...]; (a1,a2)=dotsize(a))
    b==nothing || (b=[b...]; (b1,b2)=dotsize(b))
    (c1,c2) = (c == nothing ? (0,0) :
               length(c)==2 ? c :
               throw(DimensionMismatch()))
    if a != nothing && b != nothing
        a2 == b1 ? nothing :
        a2 == 0  ? a2=b1 :
        b1 == 0  ? b1=a2 :
        throw(DimensionMismatch())
    end
    if a != nothing
        a1 == c1 ? nothing :
        a1 == 0  ? a1=c1 :
        c1 == 0  ? c1=a1 :
        throw(DimensionMismatch())
        length(a) <= 2 && (a[1] = a1)
        length(a) >= 2 && (a[end] = a2)
    end
    if b != nothing
        b2 == c2 ? nothing :
        b2 == 0  ? b2=c2 :
        c2 == 0  ? c2=b2 :
        throw(DimensionMismatch())
        length(b) <= 2 && (b[1] = b1)
        length(b) >= 2 && (b[end] = b2)
    end
    (a==nothing ? a : tuple(a...),
     b==nothing ? b : tuple(b...),
     (c1,c2))
end

### DEAD CODE:

# import Base: dot

# TODO: Who implements averaging?  -- net should
# TODO: handle or don't use nothings? -- don't use

# dot(w,x,y)=(Dot(),w,x,y)


text:
using HierarchicalTemporalMemory, BenchmarkTools, CSV, Setfield, DataFrames
import Random.seed!
seed!(0)

include("setup.jl")

SUITE = BenchmarkGroup()
SUITE["SP"] = BenchmarkGroup()
SUITE["TM"] = BenchmarkGroup()

sp, refTM, data, encode, decoder= setupHTMexperiment()
z,_= encode(data,1)

seed!(0)
let grp= SUITE["SP"]
  # Run 1 iteration of SP adaptation repeatedly
  grp["1 iteration, global inhibit"]= @benchmarkable step!(sp,z)
  grp["1 iteration, local inhibit"]= @benchmarkable step!($(@set sp.params.enable_local_inhibit= true) ,z)
end

seed!(0)
let grp= SUITE["TM"]
  tN= size(data,1)
  a= Vector{sp(z)|>typeof}(undef,tN)
  for t= 1:tN
    z,_= encode(data,t)
    a[t]= step!(sp,z)
  end
  # The TM grows new dendrites/synapses at every step, and it will try to grow more
  # if the input is random. So try giving it reasonable input sequences.
  grp["first 100 steps"]= @benchmarkable for t= 1:100
      step!(tm,$a[t])
    end setup=( tm= deepcopy($refTM) )

  for t= 1 : tN-100
    step!(refTM,a[t])
  end
  grp["last 100 steps"]= @benchmarkable for t= $tN-99:$tN
      step!(tm,$a[t])
    end setup=( tm= deepcopy($refTM) )
end

#results= run(SUITE, verbose=true, seconds=20)


text:
## MIT license (C) 2019 by Andrew Lyasoff
##
## note the formula
## $\int_a^bf(t)dt={1\over2}(b-a)\int_{-1}^{+1}f\bigl({1\over2}(b-a)s+{1\over2}(a+b)\bigr)d s$
##(Used by FastGaussQuadrature)
##

begin
    using SpecialFunctions
    using Interpolations
    using Roots
    using FastGaussQuadrature

    function EUcall(S::Float64,t::Float64,K::Float64,σ::Float64,r::Float64,δ::Float64)
        local v1,v2
        v1=(δ-r+σ^2/2.0)*t+log(K/S)
        v2=(-δ + r + σ^2/2.0)*t - log(K/S)
        v3=σ*sqrt(2.0*t)
        return -K*(exp(-t*r)/2)+K*(exp(-t*r)/2)*erf(v1/v3)+S*(exp(-t*δ)/2)*erf(v2/v3)+S*(exp(-t*δ)/2)
    end

    function F(ϵ::Int64,t::Float64,u::Float64,v::Float64,r::Float64,δ::Float64,σ::Float64)
        v1=(r-δ+ϵ*σ^2/2)*(u-t)-log(v)
        v2=σ*sqrt(2*(u-t))
        return 1.0+erf(v1/v2)
    end

    function ah(t::Float64,z::Float64,r::Float64,δ::Float64,σ::Float64,f)
        return (exp(-δ*(z-t))*(δ/2)*F(1,t,z,f(z)/f(t),r,δ,σ))
    end

    function bh(t::Float64,z::Float64,r::Float64,δ::Float64,σ::Float64,K::Float64,f)
        return (exp(-r*(z-t))*(r*K/2)*F(-1,t,z,f(z)/f(t),r,δ,σ))
    end

    function make_grid0(step::Float64,size::Int64)
        return 0.0:step:(step+(size-1)*step)
    end
end

function mainF(start_iter::Int64,nmb_of_iter::Int64,conv_tol::Float64,K::Float64,σ::Float64,δ::Float64,r::Float64,T::Float64,Δ::Float64,nmb_grd_pnts::Int64,vls::Array{Float64,1},nds::Array{Float64,1},wghts::Array{Float64,1})
    local no_iter,conv_check,absc,absc0,valPrev,val,loc,f
    absc=range(0.0,length=nmb_grd_pnts+1,stop=T)
    absc0=range(0.0,length=nmb_grd_pnts,stop=(T-Δ));
    val=vls;
    f=CubicSplineInterpolation(absc,val,extrapolation_bc = Interpolations.Line())
    no_iter=start_iter;
    conv_check=100.0
    while no_iter<nmb_of_iter&&conv_check>conv_tol
        no_iter+=1
        loc=[max(K,K*(r/δ))]
        for ttt=Iterators.reverse(absc0)
            an=[(1/2)*(T-ttt)*ah(ttt,(1/2)*(T-ttt)*s+(1/2)*(ttt+T),r,δ,σ,f) for s in nodes]
            bn=[(1/2)*(T-ttt)*bh(ttt,(1/2)*(T-ttt)*s+(1/2)*(ttt+T),r,δ,σ,K,f) for s in nodes]
            aaa=weights'*an
            bbb=weights'*bn;
            LRT=find_zero(x->x-K-EUcall(x,T-ttt,K,σ,r,δ)-aaa*x+bbb,(K-10,K+20));
            pushfirst!(loc,LRT)
        end
        valPrev=val
        val=loc
        f=CubicSplineInterpolation(absc,val,extrapolation_bc = Interpolations.Line())
        conv_check=maximum(abs.(valPrev-val))
    end
    return absc,val,conv_check,no_iter
end

max_stars_repo_name: AndrewLyasoff/SMAP

max_issues_repo_name: lhyzh/SMAP

max_forks_repo_name: lhyzh/SMAP


text:
import LinearAlgebra.LAPACK: chkstride1, chkuplo, BlasInt, liblapack
import LinearAlgebra.BLAS: @blasfunc

struct ReinschQ{T} <: AbstractMatrix{T}
    h::Vector{T}
end

function Base.size(Q::ReinschQ)
    n = length(Q.h)
   (n+1, n-1)
end

function Base.getindex(Q::ReinschQ{T},i::Int,j::Int) where T
    h = Q.h
    if (i == j)
        1/h[i]
    elseif (i == j+1)
        -1/h[j] - 1/h[j+1]
    elseif (i == j+2)
        1/h[j+1]
    else
        zero(T)
    end
end

function LinearAlgebra.mul!(out::AbstractVector, Q::Transpose{T,ReinschQ{T}}, g::AbstractVector) where T
    n = length(out)
    h = Q.parent.h
    n == size(Q,1) || throw(DimensionMismatch())
    length(g)   == size(Q,2) || throw(DimensionMismatch())
    Δgp1 = (g[2] - g[1])/h[1]
    @inbounds for i=1:length(out)
        Δg = Δgp1
        Δgp1 = (g[i+2] - g[i+1])/h[i+1]
        out[i] = Δgp1 - Δg
    end
    out
end

function LinearAlgebra.mul!(out::AbstractVector, Q::ReinschQ, g::AbstractVector)
    n = length(out)
    n == size(Q,1) || throw(DimensionMismatch())
    length(g) == size(Q,2) || throw(DimensionMismatch())
    @inbounds for i = 1:n
        out[i] = 0
        for j = max(1,i-2):min(i,n-2)
            out[i] += g[j]*Q[i,j]
        end
    end
    out
end

struct ReinschR{T} <: AbstractMatrix{T}
    h::Vector{T}
end

function Base.size(R::ReinschR)
    n = length(R.h)
   (n-1, n-1)
end

function Base.getindex(R::ReinschR{T}, i::Int, j::Int) where T
    # TODO: add check bounds
    h = R.h
    if (i==j)
        (h[i] + h[i+1])/3
    elseif abs(i-j) == 1
        h[max(i,j)]/6
    else
        zero(T)
    end
end

function QtQpR(h::AbstractVector{T}, α::T, w::AbstractVector{T}=ones(T, length(h)+1)) where T<:Real
    # maybe has some redundant calculations but should be good enough for now
    n = length(h)-1
    Q = ReinschQ(h)
    R = ReinschR(h)
    out = zeros(T, 3, n)
    # start with main diagonal
    for i=1:n
        out[3, i] = α*((w[i+2]\Q[i+2,i] - w[i+1]\Q[i+1,i])/h[i+1] - (w[i+1]\Q[i+1,i] -w[i]\Q[i, i])/h[i]) +
                  R[i,i]
    end
    # 1st superdiagonal
    for i=1:(n-1)
        # idx (i, i+1)
        out[2,i+1] = α*((w[i+2]\Q[i+2, i+1] - w[i+1]\Q[i+1, i+1])/h[i+1] - w[i+1]\Q[i+1, i+1]/h[i]) +
                  R[i, i+1]
    end
    # 2nd superdiagonal
    for i=1:(n-2)
        # idx (i, i+2)
        out[1, i+2] = α*(w[i+2]\Q[i+2, i+2]/h[i+1]) + R[i,i+2]
    end
    out
end

for (pbtrf, pbtrs, elty) in
    ((:dpbtrf_,:dpbtrs_,:Float64),
     (:spbtrf_,:spbtrs_,:Float32),
     (:zpbtrf_,:zpbtrs_,:ComplexF64),
     (:cpbtrf_,:cpbtrs_,:ComplexF32))
    @eval begin
        # SUBROUTINE DPBTRF( UPLO, N, KD, AB, LDAB, INFO )
        # *     .. Scalar Arguments ..
        #       CHARACTER          UPLO
        #       INTEGER            INFO, KD, LDAB, N
        # *     ..
        # *     .. Array Arguments ..
        #       DOUBLE PRECISION   AB( LDAB, * )
        function pbtrf!(uplo::Char, kd::Integer, AB::StridedMatrix{$elty})
            chkuplo(uplo)
            chkstride1(AB)
            n    = size(AB, 2)
            info = Ref{BlasInt}()
            ccall((@blasfunc($pbtrf), liblapack), Cvoid,
                  (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},
                   Ptr{$elty}, Ref{BlasInt}, Ref{BlasInt}),
                   uplo, n, kd, AB, max(1,stride(AB,2)), info)
            AB
        end

        #SUBROUTINE DPBTRS( UPLO, N, KD, NRHS, AB, LDAB, B, LDB, INFO )
        # *     .. Scalar Arguments ..
        #     CHARACTER          UPLO
        #     INTEGER            INFO, KD, LDAB, LDB, N, NRHS
        # *     ..
        # *     .. Array Arguments ..
        #     DOUBLE PRECISION   AB( LDAB, * ), B( LDB, * )
        function pbtrs!(uplo::Char, kd::Integer,
                        AB::StridedMatrix{$elty}, B::StridedVecOrMat{$elty})
            chkuplo(uplo)
            chkstride1(AB, B)
            info = Ref{BlasInt}()
            n    = size(AB,2)
            if n != size(B,1)
                throw(DimensionMismatch("Matrix AB has dimensions $(size(AB)), but right hand side matrix B has dimensions $(size(B))"))
            end
            ccall((@blasfunc($pbtrs), liblapack), Cvoid,
                  (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt},  Ref{BlasInt},
                   Ptr{$elty}, Ref{BlasInt}, Ptr{$elty},   Ref{BlasInt},
                   Ref{BlasInt}),
                  uplo, n, kd, size(B,2), AB,  max(1,stride(AB,2)),
                  B, max(1,stride(B,2)), info)
            B
        end
    end
end


text:
#core.jl

function operators(n)
    rowb, colb, base, vacio = integer_digits(n)
    l = 2^n
    lb = n*2^(n-1)
    row = spzeros(lb+1)
    col = spzeros(lb+1)
    data = spzeros(lb+1)
    row[lb+1] = l
    col[lb+1] = 1
    data[lb+1] = 0
    for i in 1:lb
        j = floor(Int, rowb[i])-2^(n-floor(Int, colb[i]))
        sign = (-1)^(sum(base[floor(Int, rowb[i]),1:floor(Int, colb[i])])+1)
        row[i] = j
        col[i] = floor(Int, rowb[i])+(floor(Int, colb[i])-1)*l
        data[i] = sign
    end
    cm_tot = sparse(row, col, data)
    cd_tot = sparse(cm_tot')
    return cm_tot, cd_tot, l, base
end

#This function outputs elements for building
#the sparse matrix of basis
function integer_digits(n)
    rowb = spzeros(n*2^(n-1))
    colb = spzeros(n*2^(n-1))
    data = spzeros(n*2^(n-1))
    counter = 1
    for i in 0:(2^n-1)
        binary_base = bitstring(i)[65-n:64]
        bin_vector = splitter(binary_base)
        for j in 1:n
            if bin_vector[j] == 1
                rowb[counter] = (i+1)::Int
                colb[counter] = j::Int
                data[counter] = 1
                counter = counter + 1
            end
        end
    end
    base = sparse(rowb, colb, data)
    vacio = spzeros(2^n)
    vacio[1] = 1
    return rowb, colb, base, vacio
end

# The following function outputs the basis with fixed_state
# particle number together with the indices
function basis_m(n::Int64,m::Int64)
    len = binomial(n,m)
    rowb = spzeros(m*len)
    colb = spzeros(m*len)
    data = spzeros(m*len)
    #binarios will be the index for the basis element in
    #the complete basis
    binarios = spzeros(binomial(n,m))
    counter = 1
    counter2 = 1
    base = spzeros(len,n)
    for i in 0:(2^n-1)
        binary_base = bitstring(i)[65-n:64]
        bin_vector = splitter(binary_base)
        if sum(bin_vector) == m
            for j in 1:n
                if bin_vector[j] == 1
                    rowb[counter] = counter2::Int
                    colb[counter] = j::Int
                    data[counter] = 1
                    counter = counter + 1
                end
             end
            binarios[counter2] = i
            counter2 = counter2+1
        end
    end
    base = sparse(rowb, colb, data)
    # It could also output rowb and colb if necessary
    binarios = binarios .+1
    return base, binarios
end

function splitter(x)
    vectorized = []
    for i in x
        str_to_int = parse(Int, i)
        append!(vectorized, str_to_int)
    end
    return vectorized
end


text:
@doc raw"""
    cyclic_proximal_point(M, F, proxes, x)

perform a cyclic proximal point algorithm.

# Input
* `M` – a manifold $\mathcal M$
* `F` – a cost function $F\colon\mathcal M\to\mathbb R$ to minimize
* `proxes` – an Array of proximal maps (`Function`s) `(λ,x) -> y` for the summands of $F$
* `x` – an initial value $x ∈ \mathcal M$

# Optional
the default values are given in brackets
* `evaluationOrder` – ( [`LinearEvalOrder`](@ref) ) – whether
  to use a randomly permuted sequence ([`FixedRandomEvalOrder`](@ref)), a per
  cycle permuted sequence ([`RandomEvalOrder`](@ref)) or the default linear one.
* `λ` – ( `iter -> 1/iter` ) a function returning the (square summable but not
  summable) sequence of λi
* `stopping_criterion` – ([`StopWhenAny`](@ref)`(`[`StopAfterIteration`](@ref)`(5000),`[`StopWhenChangeLess`](@ref)`(10.0^-8))`) a [`StoppingCriterion`](@ref).
* `return_options` – (`false`) – if actiavated, the extended result, i.e. the
  complete [`Options`](@ref) are returned. This can be used to access recorded values.
  If set to false (default) just the optimal value `xOpt` if returned
and the ones that are passed to [`decorate_options`](@ref) for decorators.

# Output
* `xOpt` – the resulting (approximately critical) point of gradientDescent
OR
* `options` - the options returned by the solver (see `return_options`)
"""
function cyclic_proximal_point(M::MT,
  F::Function, proxes::Array{Function,N} where N, x0;
  evaluationOrder::EvalOrder = LinearEvalOrder(),
  stopping_criterion::StoppingCriterion = StopWhenAny( StopAfterIteration(5000), StopWhenChangeLess(10.0^-12)),
  λ = i -> 1/i,
  return_options=false,
  kwargs... #decorator options
  ) where {MT <: Manifold}
    p = ProximalProblem(M,F,proxes)
    o = CyclicProximalPointOptions(x0,stopping_criterion,λ,evaluationOrder)

    o = decorate_options(o; kwargs...)
    resultO = solve(p,o)
    if return_options
        return resultO
    else
        return get_solver_result(resultO)
    end
end
function initialize_solver!(p::ProximalProblem, o::CyclicProximalPointOptions)
    c = length(p.proxes)
    o.order = update_cpp_order(c,0,[1:c...],o.orderType)
end
function step_solver!(p::ProximalProblem, o::CyclicProximalPointOptions, iter)
    c = length(p.proxes)
    λi = o.λ(iter)
    for k=o.order
        o.x = getProximalMap(p,λi,o.x,k)
    end
    o.order = update_cpp_order(c,iter,o.order,o.orderType)
end
get_solver_result(o::CyclicProximalPointOptions) = o.x
update_cpp_order(n,i,o,::LinearEvalOrder) = o
update_cpp_order(n,i,o,::RandomEvalOrder) = collect(1:n)[randperm(length(o))]
update_cpp_order(n,i,o,::FixedRandomEvalOrder) = (i==0) ? collect(1:n)[randperm(length(o))] : o