text:
const CORS = [
    "Access-Control-Allow-Origin" => "*",
    "Access-Control-Allow-Headers" => "*",
    "Access-Control-Allow-Methods" => "*"
]
const HEADERS_JSON = [CORS...,"Content-Type" => "application/json"]
const HEADERS_TEXT = [CORS..., "Content-Type" => "text/html"]

const JSONTypes = Union{
    Dict, NamedTuple, AbstractVector{<:Real},
    AbstractVector{<:Dict}, AbstractVector{<:NamedTuple}
}
const HeaderType = Vector{Pair{String, String}}

Res(body::HTTP.Response) = body

Res(status::Integer) = HTTP.Response(status)
Res(headers::HeaderType) = HTTP.Response(200, headers)

Res(body) = HTTP.Response(200, HEADERS_TEXT, body=string(body))
Res(body::AbstractString) = HTTP.Response(200, HEADERS_TEXT, body=body)
Res(body::JSONTypes) = HTTP.Response(200, HEADERS_JSON, body=JSON3.write(body))

Res(status::Integer, headers::HeaderType) = HTTP.Response(status, headers)

Res(status::Integer, body) = HTTP.Response(status, HEADERS_TEXT, body=string(body))
Res(status::Integer, body::AbstractString) = HTTP.Response(status, HEADERS_TEXT, body=body)
Res(status::Integer, body::JSONTypes) = HTTP.Response(status, HEADERS_JSON, body=JSON3.write(body))

Res(headers::HeaderType, body) = HTTP.Response(200, headers, body=string(body))
Res(headers::HeaderType, body::AbstractString) = HTTP.Response(200, headers, body=body)
Res(headers::HeaderType, body::JSONTypes) = HTTP.Response(200, headers, body=JSON3.write(body))

function Res(content_type::String, body)
    HTTP.Response(200, [CORS..., "Content-Type" => content_type], body=string(body))
end
function Res(content_type::String, body::AbstractString)
    HTTP.Response(200, [CORS..., "Content-Type" => content_type], body=body)
end

function Res(status::Integer, content_type::String, body)
    HTTP.Response(status, [CORS..., "Content-Type" => content_type], body=string(body))
end
function Res(status::Integer, content_type::String, body::AbstractString)
    HTTP.Response(status, [CORS..., "Content-Type" => content_type], body=body)
end

function Res(status::Integer, headers::HeaderType, body)
    HTTP.Response(status, headers, body=string(body))
end
function Res(status::Integer, headers::HeaderType, body::AbstractString)
    HTTP.Response(status, headers, body=body)
end
function Res(status::Integer, headers::HeaderType, body::JSONTypes)
    HTTP.Response(status, headers, body=JSON3.write(body))
end


text:
"""
    Hyperparameters of our MANOVA-DDP model, following the paper notation.
"""
@with_kw struct Model{A, B}
    D::Int
    r0::Int = 1
    ν0::Int = D + 2
    u0::A = zeros(D)
    S0::B = cholesky(collect(I(D)))
    a0::Float64 = 1.0
    b0::Float64 = 1.0
    ζ0::Float64 = 1.0
end

ext: jl


text:
# load julia packages
using PyCall
using Pandas

# python import
sk       = pyimport("sklearn")
pd       = pyimport("pandas")
np       = pyimport("numpy")
dt       = pyimport("lale.datasets")
pickle   = pyimport("pickle")
lale     = pyimport("lale")
RF       = pyimport("lale.lib.sklearn").RandomForestRegressor
LinReg   = pyimport("lale.lib.sklearn").LinearRegression
PCA      = pyimport("lale.lib.sklearn").PCA
ICA      = pyimport("lale.lib.autogen").FastICA
Tree     = pyimport("lale.lib.sklearn").DecisionTreeRegressor
Pipeline = pyimport("lale.operators").Pipeline

# @pyimport lale.lib.lale as lb
lb             = pyimport("lale.lib.lale")
Hyperopt       = lb.Hyperopt
GridSearchCV   = lb.GridSearchCV
NoOp           = lb.NoOp
ConcatFeatures = lb.ConcatFeatures

# @pyimport lale.operators as lops
laleops       = pyimport("lale.operators")
make_pipeline = laleops.make_pipeline
make_choice   = laleops.make_choice
make_union    = laleops.make_union
make_union_nc = laleops.make_union_no_concat

# define lale operators
import Base: >>, |,&, +
(&)(a::PyObject,b::PyObject) = make_union_nc(a,b)
>>(a::PyObject,b::PyObject)  = make_pipeline(a,b)
+(a::PyObject,b::PyObject)   = make_union(a,b)
|(a::PyObject,b::PyObject)   = make_choice(a,b) # UTF \times

# load data
(train_X, train_y), (test_X, test_y) = dt.california_housing_df()
pd.concat([train_X.head(), train_y.head()], axis=1)

# pipeline 1
pca_tree_planned = Pipeline(steps=[("tfm", PCA()), ("estim", Tree())])
pca_tree_planned.fit(train_X, train_y)
predicted = pca_tree_planned.predict(test_X);
print("R2 score: ",sk.metrics.r2_score(test_y, predicted))

# pipeline 2
pipeline    = PCA >> Tree
trainedhopt = pipeline.auto_configure(train_X, train_y, optimizer=Hyperopt, cv=3, max_evals=10, verbose=true)
predicted   = trainedhopt.predict(test_X);
print("R2 score: ",sk.metrics.r2_score(test_y, predicted))

# pipeline 3
planned    = (PCA + NoOp)   >>  (LinReg | Tree)
trained    = planned.auto_configure(train_X, train_y, optimizer=Hyperopt, cv=3, max_evals=3, verbose=true)
predicted  = trained.predict(test_X);
print("R2 score: ",sk.metrics.r2_score(test_y, predicted))

# pipeline 4
planned    = (PCA(svd_solver="auto") + NoOp) >>  (LinReg | Tree)
trained    = planned.auto_configure(train_X, train_y, optimizer=Hyperopt, cv=3,  max_evals=10,verbose=true)
predicted  = trained.predict(test_X);
print("R2 score: ",sk.metrics.r2_score(test_y, predicted))

# pipeline 5
pipeline    = PCA(svd_solver="auto") >> Tree
trainedgrid = pipeline.auto_configure(train_X, train_y, optimizer=Hyperopt,max_evals=5, cv=3,verbose=true)
predicted   = trainedgrid.predict(test_X);
print("R2 score: ",sk.metrics.r2_score(test_y, predicted))

# still in experimental stage
#=
using Distributed
nprocs() == 1 && addprocs(5;exeflags="--project")

@everywhere using PyCall
@everywhere np            = pyimport("numpy")
@everywhere sk            = pyimport("sklearn")
@everywhere pd            = pyimport("pandas")
@everywhere dt            = pyimport("lale.datasets")
@everywhere lale          = pyimport("lale")
@everywhere make_pipeline = pyimport("lale.operators").make_pipeline
@everywhere make_union    = pyimport("lale.operators").make_union
@everywhere make_choice   = pyimport("lale.operators").make_choice
@everywhere Hyperopt      = pyimport("lale.lib.lale").Hyperopt
@everywhere NoOp          = pyimport("lale.lib.lale").NoOp
@everywhere RF            = pyimport("lale.lib.sklearn").RandomForestRegressor
@everywhere LinReg        = pyimport("lale.lib.sklearn").LinearRegression
@everywhere PCA           = pyimport("lale.lib.sklearn").PCA
@everywhere ICA           = pyimport("lale.lib.autogen").FastICA
@everywhere Tree          = pyimport("lale.lib.sklearn").DecisionTreeRegressor
@everywhere (train_X, train_y), (test_X, test_y) = dt.california_housing_df()
@everywhere import Base: >>, |,+
@everywhere >>(a::PyObject,b::PyObject) = make_pipeline(a,b)
@everywhere +(a::PyObject,b::PyObject)  = make_union(a,b)  # UTF \oplus
@everywhere |(a::PyObject,b::PyObject)  = make_choice(a,b) # UTF \times

@everywhere pipelines=[PCA(svd_solver="auto") >> Tree() , (PCA(svd_solver="auto") + NoOp()) >>  (LinReg() | Tree()) , (ICA() + NoOp()) >> LinReg() , PCA() >> LinReg() , ICA() >> RF()]

res = @sync  @distributed (vcat) for pipeline in pipelines
   trainedpipe = pipeline.auto_configure(train_X, train_y,optimizer=Hyperopt, cv=3, max_evals=3,verbose=true)
   predicted = trainedpipe.predict(test_X)
   sk.metrics.r2_score(test_y, predicted)
end
=#


text:
using Documenter, ExtendedExtremes, Distributions, Random

makedocs(modules = [ExtendedExtremes, Distributions, Random],
        doctest = false,
        sitename="ExtendedExtremes.jl",
        pages = [
        "index.md",
        "starting.md",
        "distributions.md",
        "fit.md",
        "diagnostics.md",
        "contributing.md",
        #"Tutorial" =>["Getting started" => "tutorial/index.md"],
        #"contributing.md",
        "functions.md"
        ]
)

deploydocs(
        repo = "github.com/JuliaExtremes/ExtendedExtremes.jl.git",
)


text:
using Random
using Test
using FFTW
using CUDA
using TransformPlayground

Random.seed!(0)

FT = Float64
ε = eps(FT)

for N in (4, 5)
    A_cpu = randn(FT, N)
    A_gpu = CuArray(A_cpu)

    B_cpu = FFTW.r2r(A_cpu, FFTW.REDFT10)
    B_gpu = TransformPlayground.dct_makhoul_1d(A_gpu)

    @test isapprox.(B_cpu, Vector(B_gpu), atol=2ε) |> all

    C_cpu = FFTW.r2r(B_cpu, FFTW.REDFT01)
    C_cpu = C_cpu ./ 2N

    C_gpu = TransformPlayground.idct_makhoul_1d(B_gpu)

    @test isapprox.(C_cpu, Vector(C_gpu), atol=ε) |> any

    @test isapprox.(A_cpu, C_cpu, atol=ε) |> all
    @test isapprox.(Vector(A_gpu), Vector(C_gpu), atol=ε) |> all
end


text:
import HoTRG_lemon.LatticeModule:
		isZeroTemperature,
		isClassical,
		getModelname,
		getStates,
		isSymmetricFactorization,
		isPottsModel,
		getHamiltonian,
		getTemperature,
		getExternalfield,
		getEnvParameters,
		getTensorT,
		getFactorW,
		getCoarserate,
		getMeasureOperator,
		setTemperature!,
		setExternalfield!,
		setTensorT!,
		setEnvParameters!

#=
Simulator
getExpectationValue
getFreeEnergy(simulator)
getDimM(simulator)
getWholeiteration(simulator)
setNormalizationfactor!(simulator, normalizationfactor)
getNormalizationfactors(simulator)
=#

abstract Simulator{T}

#---
# build simulator according to the type of the lattice
"""
	buildSimulator

Build Simulator from `Lattice`.

# arguments
* `lattice`:
  * `Classical2dSquareLattice` with  `dimM`, `wholeiteration`
  * `Classical2dFractalLattice` with `dimM`, `wholeiteration`
  * `Classicl3dSquareLattice` with `dimM`, `wholeiteration`
  * `Quantum2dSquareLattice` with `dimM`
  * `Quantum2dFractalLattice` with `dimM`
* `dimM::Int`: the maximum tensor size
* `wholeiteration::Int`: For `ClassicalLattice`. How many times to  iterate.
For `QuantumLattice`, it is determined by `trotteriteration`.
"""
function buildSimulator{T}(lattice::Classical2dSquareLattice{T}, dimM::Int, wholeiteration::Int)
	return Classical2dSquareSimulator{T}(lattice, dimM, wholeiteration)
end

function buildSimulator{T}(lattice::Classical2dFractalLattice{T}, dimM::Int, wholeiteration::Int)
	return Classical2dFractalSimulator{T}(lattice, dimM, wholeiteration)
end

function buildSimulator{T}(lattice::Classical3dSquareLattice{T}, dimM::Int, wholeiteration::Int)
	return Classical3dSquareSimulator{T}(lattice, dimM, wholeiteration)
end

function buildSimulator{T}(lattice::Quantum2dSquareLattice{T}, dimM::Int)
	return Quantum2dSquareSimulator{T}(lattice, dimM)
end

function buildSimulator{T}(lattice::Quantum2dFractalLattice{T}, dimM::Int)
	return Quantum2dFractalSimulator{T}(lattice, dimM, 0)
end

function buildSimulator{T}(lattice::Quantum2dFractalLattice{T}, dimM::Int, inititeration::Int)
	return Quantum2dFractalSimulator{T}(lattice, dimM, inititeration)
end
#---
# functions for the expectation values
function getExpectationValue(simulator::Simulator)
	tensorT, impurityTensorTtilde = getTensorT(simulator)
	expectationvalue = traceTensorTPeriodic(impurityTensorTtilde) / traceTensorTPeriodic(tensorT)
	# println("traceT is \t", traceTensorTPeriodic(impurityTensorTtilde))
	# println("tractTt is\t",traceTensorTPeriodic(tensorT))
	# println("exp value\t", expectationvalue)
	return expectationvalue
end

##### ------- under this to be modified
#---
# functions inside of simulator

function getFreeEnergy(simulator::Simulator)

	freeenergy = - getFirstTerm(simulator) *(
										  termFromNormFactor(simulator) +
					 termCorrection(simulator)
					)
	return freeenergy
end

function termCorrection(simulator::Simulator)
	tensorT = getTensorT(simulator)[1]
	termcorrection = log(traceTensorTPeriodic(tensorT))/getNumberOfSites(simulator)
	return termcorrection
end

function termFromNormFactor(simulator::Simulator)
	termNorm = 0.0
	data4energy = getData4Energy(simulator)
	for i in 1:getCount(simulator)
		termNorm += data4energy[i]
	end
	return termNorm
end

#---
# functions for simulator
function getDimM(simulator::Simulator)
	return simulator.dimM
end

function getWholeiteration(simulator::Simulator)
	return simulator.wholeiteration
end

function getNormalizationfactors(simulator::Simulator)
	return simulator.normalizationFactor
end

function getData4Energy(simulator::Simulator)
	return simulator.data4energy
end

function getCount(simulator::Simulator)
	return simulator.countiteration
end

function getNumberOfSites(simulator::Simulator)
	a = getCoarserate(simulator)
	b = 2
	n = getSpaceCount(simulator)
	l = getTrotterCount(simulator)
	return exp(n*log(a) + l*log(b))
end

function countUp!(simulator::Simulator)
	simulator.countiteration = simulator.countiteration + 1
end

function countDown!(simulator::Simulator)
	simulator.countiteration = simulator.countiteration - 1
end

function isDone(simulator::Simulator)
	if getCount(simulator) >= getWholeiteration(simulator)
		return true
	else
		return false
	end
end

function setNormalizationfactor!{T}(simulator::Simulator{T}, normalizationfactor::T)
	iteration = getCount(simulator)
	simulator.normalizationFactor[iteration] = normalizationfactor
	dummy = log(normalizationfactor)/getNumberOfSites(simulator)
	simulator.data4energy[iteration] = dummy
end

#---
# functions for local tensors

function normalizeTensor(tensorT::Array; smallnumber = 1.0e-12)
	normalizationfactor = maximum(abs(tensorT))
	normalizationfactor = normalizationfactor <= smallnumber ? 1.0 : normalizationfactor
	normalizedTensorT = tensorT ./ normalizationfactor
	return normalizedTensorT, normalizationfactor
end

function traceTensorTPeriodic{T}(tensorT::Array{T,4})
	# trace tensorT(x,xp,y,yp) with periodic boundary condition
	xDimension = size(tensorT,1) == size(tensorT,2) ? size(tensorT,1) : error("size does not match")
	size(tensorT,3) == size(tensorT,4) ? yDimension = size(tensorT,3) : error("size does not match")
	tracetensor = 0.0
	for i = 1:xDimension, j = 1:yDimension
		tracetensor += tensorT[i,i,j,j]
	end
	return tracetensor
end

function traceTensorTPeriodic{T}(tensorT::Array{T,6})
	# trace tensorT(x,xp,y,yp,z,zp) with periodic boundary condition
	xDimension = size(tensorT,1) == size(tensorT,2) ? size(tensorT,1) : error("size does not match")
	size(tensorT,3) == size(tensorT,4) ? yDimension = size(tensorT,3) : error("size does not match")
	size(tensorT,5) == size(tensorT,6) ? zDimension = size(tensorT,5) : error("size does not match")
	tracetensor = 0.0
	for i = 1:xDimension, j = 1:yDimension, k = 1:zDimension
		tracetensor += tensorT[i,i,j,j,k,k]
	end
	return tracetensor
end

function rotateTensorT90{T}(tensorT::Array{T,4})
	# rotate 90 degrees tensorT to the clockwise
	# tenT[x,xp,y,yp] -> tenT[yp,y,x,xp]
	rotatedTensorT = permutedims(tensorT,[4,3,1,2])
	return rotatedTensorT
end

#---
# functions inheritted from lattice type
single_args_non_preserving = [
							:isZeroTemperature,
							:isClassical,
							:getModelname,
							:getStates,
							:getCoarserate,
							:isSymmetricFactorization,
							:isPottsModel,
							:getHamiltonian,
							:getTemperature,
							:getExternalfield,
							:getEnvParameters,
							:getTensorT,
							:getFactorW,
							:getMeasureOperator
							]

for f in single_args_non_preserving
	eval(quote

	  function $(f)(simulator::Simulator)
	  return $(f)(simulator.lattice)
		end

	end)
end

function getTensorT(simulator::Simulator, whichIsTensorT::Int)
	return getTensorT(simulator.lattice, whichIsTensorT)
end

set_functions = [
				:setTemperature!,
				:setExternalfield!
				]

for f in set_functions
	eval(quote

	  function $(f){T}(simulator::Simulator{T}, parameter::T)
		  $(f)(simulator.lattice, parameter)
		end

	end)
end

function setEnvParameters!{T}(simulator::Simulator{T}, temp::T, exf::T)
	setEnvParameters!(simulator.lattice, temp, exf)
end

function setTensorT!(simulator::Simulator, tensorT::Array; tilde = false)
	setTensorT!(simulator.lattice, tensorT; tilde = tilde)
end

function writeVector{T}(vector::Array{T,1}, filename::AbstractString = "data.txt")
	open(filename, "w") do f
		for i = 1:size(vector,1)
			write(f, "$(vector[i]) \n")
		end
	end
end


text:
# Before running this, please make sure to activate and instantiate the
# tutorial-specific package environment, using this
# [`Project.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/ISL-lab-5/Project.toml) and
# [this `Manifest.toml`](https://raw.githubusercontent.com/juliaai/DataScienceTutorials.jl/gh-pages/__generated/ISL-lab-5/Manifest.toml), or by following
# [these](https://juliaai.github.io/DataScienceTutorials.jl/#learning_by_doing) detailed instructions.

# ## Getting started

using MLJ
import RDatasets: dataset
import DataFrames: DataFrame, select
auto = dataset("ISLR", "Auto")
y, X = unpack(auto, ==(:MPG), col->true)
train, test = partition(eachindex(y), 0.5, shuffle=true, rng=444);

# Note the use of `rng=` to seed the shuffling of indices so that the results are reproducible.

# ### Polynomial regression

LR = @load LinearRegressor pkg=MLJLinearModels

# In this part we only build models with the `Horsepower` feature.

using PyPlot

figure(figsize=(8,6))
plot(X.Horsepower, y, ls="none", marker="o")

xlabel("Horsepower", fontsize=14)
xticks(50:50:250, fontsize=12)
yticks(10:10:50, fontsize=12)
ylabel("MPG", fontsize=14)

# \figalt{MPG v Horsepower}{ISL-lab-5-g1.svg}

# Let's get a baseline:

lm = LR()
mlm = machine(lm, select(X, :Horsepower), y)
fit!(mlm, rows=train)
rms(MLJ.predict(mlm, rows=test), y[test])^2

# Note that we square the measure to  match the results obtained in the ISL labs where the mean squared error (here we use the `rms` which is the square root of that).

xx = (Horsepower=range(50, 225, length=100) |> collect, )
yy = MLJ.predict(mlm, xx)

figure(figsize=(8,6))
plot(X.Horsepower, y, ls="none", marker="o")
plot(xx.Horsepower, yy, lw=3)

xlabel("Horsepower", fontsize=14)
xticks(50:50:250, fontsize=12)
yticks(10:10:50, fontsize=12)
ylabel("MPG", fontsize=14)

# \figalt{1st order baseline}{ISL-lab-5-g2.svg}

# We now want to build three polynomial models of degree 1, 2 and 3 respectively; we start by forming the corresponding feature matrix:

hp = X.Horsepower
Xhp = DataFrame(hp1=hp, hp2=hp.^2, hp3=hp.^3);

# Now we  can write a simple pipeline where the first step selects the features we want (and with it the degree of the polynomial) and the second is the linear regressor:

LinMod = Pipeline(
    FeatureSelector(features=[:hp1]),
    LR()
);

# Then we can  instantiate and fit 3 models where we specify the features each time:

lr1 = machine(LinMod, Xhp, y) # poly of degree 1 (line)
fit!(lr1, rows=train)

LinMod.feature_selector.features = [:hp1, :hp2] # poly of degree 2
lr2 = machine(LinMod, Xhp, y)
fit!(lr2, rows=train)

LinMod.feature_selector.features = [:hp1, :hp2, :hp3] # poly of degree 3
lr3 = machine(LinMod, Xhp, y)
fit!(lr3, rows=train)

# Let's check the performances on the test set

get_mse(lr) = rms(MLJ.predict(lr, rows=test), y[test])^2

@show get_mse(lr1)
@show get_mse(lr2)
@show get_mse(lr3)

# Let's visualise the models

hpn  = xx.Horsepower
Xnew = DataFrame(hp1=hpn, hp2=hpn.^2, hp3=hpn.^3)

yy1 = MLJ.predict(lr1, Xnew)
yy2 = MLJ.predict(lr2, Xnew)
yy3 = MLJ.predict(lr3, Xnew)

figure(figsize=(8,6))
plot(X.Horsepower, y, ls="none", marker="o")
plot(xx.Horsepower, yy1, lw=3, label="Order 1")
plot(xx.Horsepower, yy2, lw=3, label="Order 2")
plot(xx.Horsepower, yy3, lw=3, label="Order 3")

legend(fontsize=14)

xlabel("Horsepower", fontsize=14)
xticks(50:50:250, fontsize=12)
yticks(10:10:50, fontsize=12)
ylabel("MPG", fontsize=14)

# \figalt{1st, 2nd and 3d order fit}{ISL-lab-5-g3.svg}

# ## K-Folds Cross Validation
#
# Let's crossvalidate over the degree of the  polynomial.
#
# **Note**: there's a  bit of gymnastics here because MLJ doesn't directly support a polynomial regression; see our tutorial on [tuning models](/getting-started/model-tuning/) for a gentler introduction to model tuning.
# The gist of the following code is to create a dataframe where each column is a power of the `Horsepower` feature from 1 to 10 and we build a series of regression models using incrementally more of those features (higher degree):

Xhp = DataFrame([hp.^i for i in 1:10], :auto)

cases = [[Symbol("x$j") for j in 1:i] for i in 1:10]
r = range(LinMod, :(feature_selector.features), values=cases)

tm = TunedModel(model=LinMod, ranges=r, resampling=CV(nfolds=10), measure=rms)

# Now we're left with fitting the tuned model

mtm = machine(tm, Xhp, y)
fit!(mtm)
rep = report(mtm)

res = rep.plotting

@show round.(res.measurements.^2, digits=2)
@show argmin(res.measurements)

# So the conclusion here is that the 5th order polynomial does quite well.
#
# In ISL they use a different seed so the results are a bit different but comparable.

Xnew = DataFrame([hpn.^i for i in 1:10], :auto)
yy5 = MLJ.predict(mtm, Xnew)

figure(figsize=(8,6))
plot(X.Horsepower, y, ls="none", marker="o")
plot(xx.Horsepower, yy5, lw=3)

xlabel("Horsepower", fontsize=14)
xticks(50:50:250, fontsize=12)
yticks(10:10:50, fontsize=12)
ylabel("MPG", fontsize=14)

# \figalt{5th order fit}{ISL-lab-5-g4.svg}

# ## The Bootstrap
#
# _Bootstrapping is not currently supported in MLJ._

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
using CSV
using DataFrames

struct RegressionData
    X
    Y
    x_label
    y_label
end

# 物件価格）説明変数:居住面積
function regression_data1()
    path = "../data/house-prices-advanced-regression-techniques/train.csv"
    data = CSV.read(path, DataFrame)
    X = Array(data[data.MSSubClass .== 60, [:GrLivArea]])
    Y = data[data.MSSubClass .== 60, :SalePrice]
    x_label = "居住面積x[平方フィート]"
    y_label = "物件価格y[ドル]"
    return RegressionData(X, Y, x_label, y_label)
end

# 物件価格）説明変数:居住面積,車庫面積,プール面積,ベッド部屋数,全部屋
function regression_data2()
    path = "../data/house-prices-advanced-regression-techniques/train.csv"
    data = CSV.read(path, DataFrame)
    X = Array(data[data.MSSubClass .== 60, [:GrLivArea, :GarageArea, :PoolArea, :BedroomAbvGr, :TotRmsAbvGrd]])
    Y = data[data.MSSubClass .== 60, :SalePrice]
    x_label = ""
    y_label = "物件価格y[ドル]"
    return RegressionData(X, Y, x_label, y_label)
end

# 物件価格）説明変数:居住面積、物件価格に外れ値を追加
function regression_data3()
    path = "../data/house-prices-advanced-regression-techniques/train.csv"
    data = CSV.read(path, DataFrame)
    X = Array(data[data.MSSubClass .== 60, [:GrLivArea]])
    Y = data[data.MSSubClass .== 60, :SalePrice]
    Y[Y .> 700000] .-= 700000  # 外れ値
    x_label = "居住面積x[平方フィート]"
    y_label = "物件価格y[ドル]"
    return RegressionData(X, Y, x_label, y_label)
end


text:
using CoolPkg
using Documenter

makedocs(;
    modules=[CoolPkg],
    authors="Daniel Winkler <danielw2904@disroot.org> and contributors",
    repo="https://github.com/danielw2904/CoolPkg.jl/blob/{commit}{path}#L{line}",
    sitename="CoolPkg.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://danielw2904.github.io/CoolPkg.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/danielw2904/CoolPkg.jl",
)


text:
#------------------------------
# creating source wavelets
#------------------------------
function myricker2(tvec,f0,delay,derivative_number)
    # Gaussian=exp(-sigma*(t-t0)^2), sigma=-pi^2*f0^2. Mittet (eq 25), 1996, Geophysics
    nt=length(tvec)
    y=zeros(1,nt)
    if (derivative_number==2)
        println("Myricker: 2nd derivative of Gaussian (Ricker with a negative sign)")
    #---2nd derivative Gaussian (Ricker with a negative sign)
    for it=1:nt
     y[it] = -2*pi^2*f0^2*(1.0 - 2.0*(pi^2)*(f0^2)*((tvec[it]-delay).^2)) .* exp(-(pi^2)*(f0^2)*((tvec[it]-delay).^2))
    end
    elseif (derivative_number==1)
        println("Myricker: 1st derivative of Gaussian")
    #---1st derivative Gaussian
    for it=1:nt
     y[it] = -2*pi^2*f0^2*(tvec[it]-delay) .* exp(-(pi^2)*(f0^2)*((tvec[it]-delay).^2))
    end
    elseif (derivative_number==0)
        println("Myricker: Gaussian")
        #---Gaussian
    for it=1:nt
        y[it] = exp(-(pi^2)*(f0^2)*((tvec[it]-delay).^2))
    end
    elseif (derivative_number==3)
        println("Myricker: 3rd derivative of Gaussian")
    #---3rd derivative of Gaussian
    for it=1:nt
        y[it] = 4.0*pi^4*f0^4*(-2.0*pi^2*f0^2*(tvec[it]-delay)^3+3.0*(tvec[it]-delay))*
        exp(-(pi^2)*(f0^2)*((tvec[it]-delay).^2))
    end
    else
        println("Myricker: Error. Derivative_number should be 0, 1, 2 or 3.")
    end

    return y
end

#------------------------------
# converting src geometry to
# index numbers of a pressure grid
#------------------------------
function get_srcindex_monopole(srcgeom,dr,dz)
    #src index and amplitudes
    #--point input
    izsrc_ref=Int(ceil(srcgeom[1,1]/dz)+1)
    irsrc_ref=Int(ceil(srcgeom[1,2]/dr)+1)
    src_index=[izsrc_ref;irsrc_ref]
    src_dn=[1.]
    return src_index,src_dn
end

#------------------------------
# applying src term into grids
#------------------------------
function srcapply!(paramarray,src_index,src_dn,srcamp)
    #--e.g. fz src: vz=vz+fz (same form as vz update, thus fz is as-is and not time-derivatives)
    #--input array can be vx,vz etc
    for isrcind=1:length(src_dn)
       isz=src_index[1,isrcind]
       isr=src_index[2,isrcind]
       paramarray[isz,isr]=paramarray[isz,isr]+srcamp*src_dn[isrcind]
    end
end

#-----------------------------------
# converting src geometry to
# index numbers of a pressure grid
# Rubust point src by Gaussian shape
#-----------------------------------
function get_srcindex_pGauss(srcgeom,dr,dz,windowsize,sigma)
    #window size->odd
    #src index and amplitudes
    #Correct scaling factor (as a dirac delta in cylindrical coordinate) only when source is located around r=0
    #--point input
    izsrc_ref=Int(ceil(srcgeom[1,1]/dz)+1)
    irsrc_ref=Int(ceil(srcgeom[1,2]/dr)+1)

    Gauss_Wmat=Gaussian_filter(windowsize,sigma)

    halfw=Int(round((windowsize-1)/2))

    #check if gaussian window exceeds at left boundary
    if(irsrc_ref-halfw < 1)
        irvec_1row=collect(range(irsrc_ref-halfw,irsrc_ref+halfw,length=windowsize))
        ir1=findall(x->x==1,irvec_1row)

        izvec_1col=collect(range(izsrc_ref-halfw,izsrc_ref+halfw,length=windowsize))
        izvec_mcol=repeat(izvec_1col,outer=windowsize-ir1[1]+1)

        irvec_1row=collect(range(1,irsrc_ref+halfw,length=windowsize-ir1[1]+1))
        irvec_mcol=repeat(irvec_1row,inner=windowsize)

        src_index=zeros(2,windowsize*(windowsize-ir1[1]+1))
        src_index[1,:]=izvec_mcol[:]
        src_index[2,:]=irvec_mcol[:]

        #extracting symmetric components around r=0
        Wmat=Gauss_Wmat[:,ir1[1]:windowsize]
        #evaluating volume integral to get scaing factor
        V=Check_VolumeIntegral(Wmat,dr,dz)
        #saling-corrected Gaussian
        src_dn=Wmat/V

        src_dn=src_dn[:]
    else

        src_index=zeros(2,windowsize^2)

        izvec_1col=collect(range(izsrc_ref-halfw,izsrc_ref+halfw,length=windowsize))
        izvec_mcol=repeat(izvec_1col,outer=windowsize)

        irvec_1row=collect(range(irsrc_ref-halfw,irsrc_ref+halfw,length=windowsize))
        irvec_mcol=repeat(irvec_1row,inner=windowsize)

        src_index[1,:]=izvec_mcol[:]
        src_index[2,:]=irvec_mcol[:]

        #Note: scaling factor does not represent a dirac delta function
        src_dn=Gauss_Wmat[:]
    end

    src_index=map(x->Int(x),src_index)
    return src_index,src_dn
end

#Gaussian filter
function Gaussian_filter(wsize, sigma)
#wsize should be odd number
    gauss=zeros(wsize,wsize); #2D filter matrix
    for i=-round((wsize-1)/2):1:round((wsize-1)/2)
        for j=-round((wsize-1)/2):1:round((wsize-1)/2)
            i0=(wsize+1)/2;
            j0=(wsize+1)/2;
            i1=Int(round(i+i0));
            j1=Int(round(j+j0));
            gauss[j1,i1]=exp(-((i1-i0)^2+(j1-j0)^2)/(2*sigma^2));
        end
    end
gauss=gauss/sum(gauss[:]);
return gauss
end

function Check_VolumeIntegral(Wmat,dr,dz)
    #evaluate volume integral cylindrical coordinate
    #V=2pi int int r dr dz f(r,z)
    #assumes (Gaussian) matrices symmetrical around r=0
    # in case of a scalar at r=0, it results in pi*dz*dr^2*1/4
    wsize_z,wsize_r=size(Wmat)

    V=0.
    for iz=1:wsize_z
        tmpvec_r=Wmat[iz,:]
        for ir=1:wsize_r
            if (ir==1)
                V=V+pi*dz*tmpvec_r[ir]*dr^2/4
            else
                V=V+pi*dz*tmpvec_r[ir]*dr^2/4*(ir-1)*8
            end
        end
    end

    return V
end


text:
export cache_data

using JLD2

function cache_data(func::Function, filename::String, path::String, args...; kwargs...)
  local prepared = nothing
  if isfile(filename)
    jldopen(filename, "r") do file
      if haskey(file, path)
        prepared = file[path]
      end
    end
  end
  if isnothing(prepared)
    prepared = func(args...; kwargs...)
    jldopen(filename, "a+") do file
      file[path] = prepared
    end
  end
  return prepared
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
println("\nTesting cost terms\n")

# Simple Terms

x = Variable(randn(10))
X = Variable(randn(3,4))
A = randn(4,10)
b = randn(4)

cf = norm(x, 0)
@test cf.lambda == 1
@test cf.f(~x) == norm(~x,0)

cf = 3*norm(x, 0)
@test cf.lambda == 3
@test cf.f(~x) == norm(~x,0)

cf = norm(x, 0) <= 3
@test cf.lambda == 1
@test cf.f(~x) == (IndBallL0(3))(~x)

cf = norm(x, 1)
@test cf.lambda == 1
@test cf.f(~x) == (NormL1())(~x)

cf = norm(x, 1) <= 1.5
@test cf.lambda == 1
@test cf.f(~x) == (IndBallL1(1.5))(~x)

cf = 10*norm(x, 1) <= 1.5
@test cf.lambda == 1
@test cf.f(~x) == (IndBallL1(1.5/10))(~x)

cf = norm(x)
@test cf.lambda == 1
@test cf.f(~x) == norm(~x)

cf = pi*norm(x,2)
@test cf.lambda - pi == 0
@test cf.f(~x) == norm(~x)

cf = 3*norm(X,2,1)
@test cf.lambda - 3 == 0
@test cf.f(~X) == sum(  sqrt.(sum((~X).^2, dims=1 )) )

cf = 4*norm(X,2,1,2)
@test cf.lambda - 4 == 0
@test cf.f(~X) == sum(  sqrt.(sum((~X).^2, dims=2 )) )

@test_throws ErrorException 4*norm(X,1,2)

cf = norm(x, 2) <= 2.3
@test cf.lambda == 1
@test cf.f(~x) == (IndBallL2(2.3))(~x)

cf = norm(x, 2) == 2.3
@test cf.lambda == 1
@test cf.f(~x) == (IndSphereL2(2.3))(~x)

cf = norm(x, Inf)
@test cf.lambda == 1
@test cf.f(~x) == norm(~x,Inf)

cf = norm(x, Inf) <= 5.0
@test cf.lambda == 1
@test cf.f(~x) == (IndBallLinf(5.0))(~x)

cf = x <= 3.0
@test cf.lambda == 1
@test cf.f(~x) == (IndBox(-Inf, 3.0))(~x)

cf = 3.0 <= x
@test cf.lambda == 1
@test cf.f(~x) == (IndBox(3.0, Inf))(~x)

cf = x >= 1.0
@test cf.lambda == 1
@test cf.f(~x) == (IndBox(1.0, Inf))(~x)

cf = 1.0 >= x
@test cf.lambda == 1
@test cf.f(~x) == (IndBox(-Inf,1.0))(~x)

cf = x in [-5.0, 5.0]
@test cf.lambda == 1
@test cf.f(~x) == (IndBox(-5.0, 5.0))(~x)

cf = norm(x, 2)^2
@test cf.lambda == 1
@test cf.f(~x) == norm(~x)^2

cf = 0.5*norm(x, 2)^2
@test cf.lambda == 0.5
@test cf.f(~x) == norm(~x)^2

cf = 7*(0.5*norm(x, 2))^2
@test cf.lambda == 7*0.25
@test cf.f(~x) == norm(~x)^2

cf = 2*rank(X) <= 6
@test cf.lambda == 1
@test cf.f(~X) == (IndBallRank(3))(~X)

cf = rank(X)
@test_throws MethodError cf.f(~X)

cf = norm(X,*)
U, S, V = svd(~X)
@test cf.lambda == 1
@test cf.f(~X) == sum(S)

cf = rank(X)
@test_throws MethodError cf.f(~X)

y = randn(size(~x))
cf = hingeloss(x,y)
@test cf.lambda == 1
@test cf.f(~x) == (HingeLoss(y))(~x)

y = randn(size(~x))
cf = sqrhingeloss(x,y)
@test cf.lambda == 1
@test cf.f(~x) == (SqrHingeLoss(y))(~x)

y = randn(size(~x))
cf = logisticloss(x,y)
@test cf.lambda == 1
@test cf.f(~x) == (LogisticLoss(y))(~x)

xp = Variable(rand(10))
bp = rand(Float64, size(~xp))
cf = crossentropy(xp,bp)
@test cf.lambda == 1
@test cf.f(~xp) == (CrossEntropy(bp))(~xp)

cf = logbarrier(x)
@test cf.lambda == 1
@test cf.f(~x) == (LogBarrier(1.0))(~x)

cf = maximum(x)
@test cf.lambda == 1
@test cf.f(~x) == (Maximum(1.0))(~x)

cf = sumpositive(x)
@test cf.lambda == 1
@test cf.f(~x) == (SumPositive())(~x)

a = 1.
cf = huberloss(x,a)
@test cf.lambda == 1
@test cf.f(~x) == (HuberLoss(a))(~x)

a = randn(size(x))
cf = dot(a,x)
@test cf.lambda == 1
@test cf.f(~x) == (Linear(a))(~x)

#IndBinary
lu = (-1.0,randn(length(~x)))
cf = x == lu
@test cf.lambda == 1
@test cf.f(~x) == (IndBinary(lu...))(~x)

# IndAffine (not working in julia < 1.1)
if VERSION.major >= 1 && VERSION.minor >= 1
  cf = A*x-b == 0
  @test cf.lambda == 1
  @test cf.f(~x) == (IndAffine(A,b))(~x)

  cf = (A*x == b)
  @test cf.lambda == 1
  @test cf.f(~x) == (IndAffine(A,-b))(~x)
end

cf = 2*norm(x,1)
ccf = conj(cf)
@test ccf.A == cf.A
@test ccf.f == Conjugate(Postcompose(NormL1(),2))
@test_throws ErrorException conj(norm(randn(2,10)*x,1))

cf = 2*norm(x,1)
ccf = smooth(cf,2.0)
@test ccf.A == cf.A
@test ccf.f(~x) == MoreauEnvelope(Postcompose(NormL1(),2),2.0)(~x)

# Summing terms

x = Variable(10)
cf = ls(x) + 10*norm(x, 1)
@test cf[1].lambda == 1
@test cf[1].f(~x) == 0.5*norm(~x)^2
@test cf[2].lambda == 10
@test cf[2].f(~x) == norm(~x,1)

x = Variable(10)
cf = () #empty cost function
cf += 10*norm(x, 1)
@test length(cf) == 1
@test cf[1].lambda == 10
@test cf[1].f(~x) == 10*norm(~x,1)

x = Variable(10)
cf = () #empty cost function
cf += ls(x) + 10*norm(x, 1)
@test cf[1].lambda == 1
@test cf[1].f(~x) == 0.5*norm(~x)^2
@test cf[2].lambda == 10
@test cf[2].f(~x) == norm(~x,1)

# More complex situations

x = Variable(10)
A = randn(5, 10)
y = Variable(7)
B = randn(5, 7)
b = randn(5)

cf = ls(A*x - b) + norm(x, 1)
@test cf[1].lambda == 1
@test cf[1].f(~x) == 0.5*norm(~x)^2
@test norm(affine(cf[1])*(~x) - (A*(~x)-b)) < 1e-12
@test cf[2].lambda == 1
@test cf[2].f(~x) == norm(~x,1)

cf = ls(A*x - B*y + b) + norm(y, 1) + 5*norm(y, 2)
@test cf[1].lambda == 1
@test cf[1].f(~x) == 0.5*norm(~x)^2
@test cf[2].lambda == 1
@test cf[2].f(~x) == norm(~x,1)
@test cf[3].lambda == 5
@test cf[3].f(~x) == norm(~x,2)

cf = 10*(ls(A*x - B*y + b) + norm(y, 1) + 5*norm(y, 2))
@test cf[1].lambda == 10
@test cf[1].f(~x) == 0.5*norm(~x)^2
@test cf[2].lambda == 10
@test cf[2].f(~x) == norm(~x,1)
@test cf[3].lambda == 50
@test cf[3].f(~x) == norm(~x,2)

cf = 0.5*norm(A*x - B*y + b, 2)^2 + norm(x, 1) + norm(y, 2)
@test cf[1].lambda == 0.5
@test cf[1].f(~x) == norm(~x)^2
@test cf[2].lambda == 1
@test cf[2].f(~x) == norm(~x,1)
@test cf[3].lambda == 1
@test cf[3].f(~x) == norm(~x,2)

# Properties
A = randn(5, 10)
u = Variable(5)
w = Variable(5)
z = Variable(5)

cf = norm(A*x + z)
@test StructuredOptimization.is_smooth(cf) == false
@test StructuredOptimization.is_smooth(cf^2) == true

cf = norm(w + z)^2
@test StructuredOptimization.is_smooth(cf) == true
@test StructuredOptimization.is_AcA_diagonal(cf) == false

cf = norm(x, 1) + norm(y, 2)
@test StructuredOptimization.is_smooth.(cf) == (false,false)
@test StructuredOptimization.is_AcA_diagonal.(cf) == (true,true)


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

JLLWrappers.@generate_main_file_header("CmdStanBinaries")
JLLWrappers.@generate_main_file("CmdStanBinaries", UUID("eecc378a-c28e-5a81-8849-a983a4161162"))
end  # module CmdStanBinaries_jll


text:
__precompile__()

module CompressedVolumes

using Images, HDF5, JLD, FunctionalData

export CompressedVolume, encode, encode!, decode, compact!

type CompressedVolume
	typ::Type
    siz::Array{Int,2}
    blocksiz
	nblocks
	blocksmi
	blocksma
    encoded
	decoded
	isdirty
	isencoded
	defaultvalue
	qp
end

defaultblocksiz = [64;64;64]

clamprange(a, mi, ma) = clamp(first(a),mi,ma):clamp(last(a),mi,ma)

function CompressedVolume{T}(a::Array{T,3}; kargs...)
    siz = reshape([size(a)...], (3,1))
	r = CompressedVolume(eltype(a), siz; kargs...)
    for o = 1:r.nblocks[3], n = 1:r.nblocks[2], m = 1:r.nblocks[1]
        r.decoded[m,n,o] = a[clamprange(r.blocksiz[1]*(m-1)+(1:r.blocksiz[1]), 1, siz[1]),
            clamprange(r.blocksiz[2]*(n-1)+(1:r.blocksiz[2]), 1, siz[2]),
            clamprange(r.blocksiz[3]*(o-1)+(1:r.blocksiz[3]), 1, siz[3])]
		r.isdirty[m,n,o] = true
    end
    r
end

CompressedVolume(size::Tuple; kargs...) = CompressedVolume([size...]'; kargs...)

CompressedVolume{T}(siz::Array{T,2}; kargs...) = CompressedVolume(Float32, siz; kargs...)
function CompressedVolume{T}(typ::Type, siz::Array{T,2}; blocksiz = defaultblocksiz, defaultvalue = 0.0, useblocks = true,
	encoded = :undefined, decoded = :undefined, qp = 28)

	siz = reshape(siz, (length(siz), 1))
	if !useblocks
		blocksiz = [typemax(Int); typemax(Int); typemax(Int)]
	end
    nblocks = ceil(Int, siz ./ blocksiz)
	blocksmi = Array(Any, nblocks...)
	blocksma = Array(Any, nblocks...)
	for o = 1:nblocks[3], n = 1:nblocks[2], m = 1:nblocks[1]
		blocksmi[m,n,o] = [blocksiz[1]*(m-1)+1; blocksiz[2]*(n-1)+1; blocksiz[3]*(o-1)+1]
		blocksma[m,n,o] = [min(blocksiz[1]*m, siz[1]); min(blocksiz[2]*n, siz[2]); min(blocksiz[3]*o, siz[3])]
	end
    #@show nblocks
	if encoded == :undefined
		encoded = {nothing for m = 1:nblocks[1], n = 1:nblocks[2], o = 1:nblocks[3]}
	end
	if decoded == :undefined
		decoded = {nothing for m = 1:nblocks[1], n = 1:nblocks[2], o = 1:nblocks[3]}
	end
    isdirty = falses(nblocks...)
    isencoded = falses(nblocks...)
    CompressedVolume(typ, siz, blocksiz, nblocks, blocksmi, blocksma, encoded, decoded, isdirty, isencoded, defaultvalue, qp)
end

encode(a; kargs...) = encode!(CompressedVolume(a; kargs...))
function encode!(a::CompressedVolume)
	for o = 1:a.nblocks[3], n = 1:a.nblocks[2], m = 1:a.nblocks[1]
        encodeblock!(a, m, n, o)
	end
	a
end

function decode!(a::CompressedVolume)
	for o = 1:a.nblocks[3], n = 1:a.nblocks[2], m = 1:a.nblocks[1]
		decodeblock!(a, m, n, o)
	end
	a
end

decode(a::CompressedVolume) = getindex(a, 1:a.siz[1], 1:a.siz[2], 1:a.siz[3])
# function decode(a::CompressedVolume)
#     r = zeros(a.siz...)
#     @show a.siz
#     for o = 1:a.nblocks[3], n = 1:a.nblocks[2], m = 1:a.nblocks[1]
# 		decodeblock!(a, m, n, o)
#         @show size(dec)
#         r[ clamprange(a.blocksiz[1]*(m-1)+(1:a.blocksiz[1]),1,a.siz[1]),
#            clamprange(a.blocksiz[2]*(n-1)+(1:a.blocksiz[2]),1,a.siz[2]),
#            clamprange(a.blocksiz[3]*(o-1)+(1:a.blocksiz[3]),1,a.siz[3])] = dec
#     end
#     r
# end

function encodeblock!(a::CompressedVolume, m, n, o)
	dec = a.decoded[m,n,o]

	if a.isdirty[m,n,o]
		# @show "dec2" dec[1]
		mi = minimum(dec)
		ma = maximum(dec)
		# @show "encodeblocks" size(dec)
		v = map(Float32,copy(dec))
		v -= mi
        d = ma - mi
        if d > 0
            v /= d
        end
		# v = uint16(v*2^16)

 		d = tempname()*"/"
		mkpath(d)
		for i = 1:size(v,3)
			save("$d$i.png",v[:,:,i])
		end
		outfile = d*"movie.mp4"
		run(`ffmpeg -v 0 -i $d%d.png -y -vcodec libx264 -qp $(a.qp) $outfile`)

		#@show "encodeblock" size(v)
		# a.encoded[m,n,o] = Dict(("nslices", "data", "mi", "ma"), (size(v,3), open(readbytes, outfile), mi, ma))
		a.encoded[m,n,o] = Dict("nslices" => size(v,3), "data" => open(readbytes, outfile), "mi" => mi, "ma" => ma)
		a.isdirty[m,n,o] = false
		a.isencoded[m,n,o] = true
	end
	nothing
end

blocksize(a::CompressedVolume, m, n, o) = a.blocksma[m,n,o] - a.blocksmi[m,n,o] + 1
function decodeblock!(a::CompressedVolume, m, n, o)
	if a.isencoded[m,n,o]
		enc = a.encoded[m,n,o]
		d = tempname()*"/"
		mkpath(d)
		outfile = d*"movie.mp4"
		open(io -> write(io, enc["data"]), outfile,"w")
		run(`ffmpeg -i $outfile -v 0 -y $d%d.png`)
		#@profile rawimread("$(d)1.png")

		slices = pmap(i -> rawimread("$d$i.png"), 1:enc["nslices"])
		dec = zeros([size(slices[1])...,length(slices)]...)
		# @show dec
		for i = 1:size(dec,3)
			dec[:,:,i] = slices[i]
		end
		dec = dec * (enc["ma"]-enc["mi"]) + enc["mi"]
		a.decoded[m,n,o] = dec
		a.isdirty[m,n,o] = false
	else
		a.decoded[m,n,o] = ones(a.typ, blocksize(a,m,n,o)...) * a.defaultvalue
		a.isdirty[m,n,o] = false
	end
	a
end

getindex(a::CompressedVolume, ms::Integer, ns::Integer, os::Integer) = accessindex(a, ms, ns, os, 0, false)[1]
getindex(a::CompressedVolume, ms, ns, os) = accessindex(a, ms, ns, os, 0, false)
setindex!(a::CompressedVolume, v, ms, ns, os) = accessindex(a, ms, ns, os, v, true)

function accessindex(a::CompressedVolume, ms, ns, os, v, set)
	mblocks = unique(ceil(Int,ms/a.blocksiz[1]))
	nblocks = unique(ceil(Int,ns/a.blocksiz[2]))
	oblocks = unique(ceil(Int,os/a.blocksiz[3]))
	# @show mblocks nblocks oblocks

	r = set ? 0. : fill(NaN, maximum(ms)-minimum(ms)+1, maximum(ns)-minimum(ns)+1, maximum(os)-minimum(os)+1)
	mi = [minimum(ms); minimum(ns); minimum(os)] - 1

	# @show a.blocksmi a.blocksma

	for o = oblocks, n = nblocks, m = mblocks
		if a.decoded[m,n,o] == nothing
			decodeblock!(a, m, n, o)
		end
		bmi = a.blocksmi[m,n,o]
		rm = clamprange(ms, bmi[1], a.blocksma[m,n,o][1])
		rn = clamprange(ns, bmi[2], a.blocksma[m,n,o][2])
		ro = clamprange(os, bmi[3], a.blocksma[m,n,o][3])

		# @show size(r)
		# @show rm rn ro mi
		# @show size(a.decoded[m,n,o])
		# @show r[ rm - mi[1], rn - mi[2], ro - mi[3] ]

		if set
			if typeof(v)<:Number
				a.decoded[m,n,o][rm - bmi[1] + 1, rn - bmi[2] + 1, ro - bmi[3] + 1] = v
			else
				a.decoded[m,n,o][rm - bmi[1] + 1, rn - bmi[2] + 1, ro - bmi[3] + 1] =
					v[ rm - mi[1], rn - mi[2], ro - mi[3] ]
			end
			a.isdirty[m,n,o] = true
		else
			r[ rm - mi[1], rn - mi[2], ro - mi[3] ] =
				a.decoded[m,n,o][rm - bmi[1] + 1, rn - bmi[2] + 1, ro - bmi[3] + 1]
		end
	end
	set ? a : r
end

function compact!(a::CompressedVolume)
	a.decoded[find(!a.isdirty)] = nothing
	a
end

function rawimread(a)
    img = load(a)
    img = convert(Array,img)
    try
        float(img)
    catch
        #@show typeof(img[1])
        r = zeros(size(img))
        for i = 1:length(img)
            r[i] = img[i].r
        end
        r
    end
end

end # module


text:
using SnoopCompile
using Base.Test

if VERSION >= v"0.5.0-dev"
    # Function as argument
    str = "Base.any($(typeof(identity)), Array{Bool, 1})"
    keep, pcstring, fpath, args = SnoopCompile.parse_call(str)
    @test keep
    @test pcstring == "    precompile(Base.any, (typeof(Base.identity), Array{Bool, 1},))"
    # Anonymous function as argument
    str = "Base.any($(typeof(x->x>0)), Array{Float32, 1})"
    keep, pcstring, fpath, args = SnoopCompile.parse_call(str)
    @test !keep
    # Function as a type
    str = "Base.Sort.sort!(Base.#sort!, Array{Any, 1}, Base.Sort.MergeSortAlg, Base.Order.By{Base.#string})"
    keep, pcstring, fpath, args = SnoopCompile.parse_call(str)
    @test keep
    @test pcstring == "    precompile(Base.Sort.sort!, (Array{Any, 1}, Base.Sort.MergeSortAlg, Base.Order.By{typeof(Base.string)},))"
end

include("colortypes.jl")


text:
@testset "Search" begin
    T = Float64
    @testset "NullSearch" begin
        p = DSProblem{T}(3, poll=LTMADS{T}())
        s = DS.NullSearch()
        @test DS.GenerateSearchPoints(p) == []
        @test DS.GenerateSearchPoints(p, s) == []
    end
    @testset "RandomSearch" begin
        N = 3
        c = generate_rand_cache(N, 10)
        d = 1.0
        s = DS.RandomSearch(10)

        points = DS.RandomPointsFromCache(N, c, d, s)

        for point in points
            mesh = keys(c.costs)
            closeness = [norm(point - m_point) ≈ d for m_point in mesh]
            @test true in closeness
        end
    end
end


text:
#=##############################################################################
# DESCRIPTION
    Unsteady simulation of eVTOL transition maneuver on tilt-wing, tandem,
    distributed propulsion aircraft.

    REFERENCES
    * Vahana geometry: Droandi, G., Syal, M., and Bower, G., “Tiltwing Multi-Rotor
    Aerodynamic Modeling in Hover, Transition and Cruise Flight Conditions,” AHS
    International 74th Annual Forum & Technology Display, 2018, p. 2018.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Oct 2019
  * License   : MIT
=###############################################################################

# ------------ MODULES ---------------------------------------------------------
import FLOWUnsteady

uns = FLOWUnsteady
vpm = uns.vpm
vlm = uns.vlm
gt = GeometricTools

# ------------ GLOBAL VARIABLES ------------------------------------------------
# Default path where to save data
extdrive_path = "/media/edoalvar/Samsung_T5/simulationdata202102/"

# Default data path where to find rotor and airfoil data
data_path = uns.def_data_path

# ------------ HEADERS ---------------------------------------------------------
for header_name in ["geometry", "kinematics", "monitor"]
    include("vahana_"*header_name*".jl")
end

# ------------ DRIVERS ---------------------------------------------------------

function run_simulation_vahana(;    save_path=extdrive_path*"vahana_sim09",
                                    prompt=true,
                                    run_name="vahana",
                                    verbose=true, v_lvl=1)

    # ----------------- PARAMETERS ---------------------------------------------

    # Geometry options
    n_factor = 1                            # Refinement factor
    add_rotors = true                       # Whether to include rotors

    # # Maneuver to perform
    # Vcruise = 0.125 * 125*0.44704         # Cruise speed
    # RPMh_w = 600                          # RPM of main wing rotors in hover
    # telapsed = 60.0                       # Total time to perform maneuver
    # nsteps = 9000                         # Time steps

    # Maneuver to perform
    Vcruise = 0.25 * 125*0.44704            # Cruise speed
    Vinf(x,t) = 1e-5*[1,0,-1]               # (m/s) freestream velocity, if 0 the simulation might crash
    # RPMh_w = 200                          # RPM of main wing rotors in hover
    # RPMh_w = 20
    RPMh_w = 100
    telapsed = 30.0                         # Total time to perform maneuver
    nsteps = 1500                           # Time steps
    # nsteps = 100

    dt = telapsed/nsteps

    # Simulation options
    rho = 1.225                             # (kg/m^3) air density
    mu = 1.81e-5                            # Air dynamic viscosity

    # Solver options
    R = 0.75                                # (m) blade radius as a reference
    lambda = 4.0                            # Target minimum core overlap
    p_per_step = 4                          # Particle sheds per time step
    overwrite_sigma = lambda * (2*pi*RPMh_w/60*R + Vcruise)*dt / p_per_step
    # vlm_sigma = R/25                      # VLM regularization
    # vlm_sigma = R
    vlm_sigma = -1
    vlm_rlx = 0.75                          # VLM relaxation (deactivated with -1)
    surf_sigma = R/10                       # Surface regularization
    shed_unsteady = false                   # Shed unsteady-loading particles
    VehicleType = uns.QVLMVehicle           # Type of vehicle to generate

    # ----------------- MANEUVER DEFINITION ------------------------------------
    # Generate maneuver
    maneuver = generate_maneuver_vahana1(; add_rotors=add_rotors)

    # Plot maneuver path and controls
    uns.plot_maneuver(maneuver; tstages=[0.2, 0.3, 0.5, 0.6])

    # ----------------- GEOMETRY GENERATION ------------------------------------
    # Generate geometry
    (vehicle, grounds) = generate_geometry_vahana(; n_factor=n_factor,
                                                    xfoil=false,
                                                    data_path=data_path,
                                                    run_name=run_name,
                                                    add_rotors=add_rotors,
                                                    VehicleType=VehicleType)

    # Move landing pad to landing area
    gt.lintransform!(grounds[2], eye(3), Vcruise*telapsed*[-0.25, 0, -0.0025])

    # ----------------- SIMULATION SETUP ---------------------------------------
    Vref = Vcruise
    RPMref = RPMh_w
    ttot = telapsed
    max_particles = ceil(Int, (nsteps+2)*(2*vlm.get_m(vehicle.vlm_system)+1)*p_per_step)

    Vinit = Vref*maneuver.Vvehicle(0)       # (m/s) initial vehicle velocity
                                            # (rad/s) initial vehicle angular velocity
    Winit = pi/180 * (maneuver.anglevehicle(0+1e-12)-
                                          maneuver.anglevehicle(0))/(ttot*1e-12)

    simulation = uns.Simulation(vehicle, maneuver, Vref, RPMref, ttot;
                                                    Vinit=Vinit, Winit=Winit)

    # ----------------- SIMULATION MONITOR -------------------------------------
    monitor = generate_monitor_vahana(vehicle, rho, RPMref, nsteps, save_path, Vinf)

    # ----------------- RUN SIMULATION -----------------------------------------
    pfield = uns.run_simulation(simulation, nsteps;
                                      # SIMULATION OPTIONS
                                      Vinf=Vinf,
                                      rho=rho,
                                      mu=mu,
                                      # SOLVERS OPTIONS
                                      p_per_step=p_per_step,
                                      overwrite_sigma=overwrite_sigma,
                                      vlm_sigma=vlm_sigma,
                                      vlm_rlx=vlm_rlx,
                                      surf_sigma=surf_sigma,
                                      max_particles=max_particles,
                                      shed_unsteady=shed_unsteady,
                                      extra_runtime_function=monitor,
                                      # OUTPUT OPTIONS
                                      save_path=save_path,
                                      run_name=run_name,
                                      prompt=prompt,
                                      verbose=verbose, v_lvl=v_lvl,
                                      save_code=splitdir(@__FILE__)[1]
                                      )

    # Save ground
    for (i, ground) in enumerate(grounds)
        gt.save(ground, run_name*"_Ground$i"; path=save_path)
    end

    return simulation, pfield
end

function visualize_maneuver_vahana(; save_path=extdrive_path*"vahana_maneuver00/",
                                        prompt=true,
                                        run_name="vahana",
                                        verbose=true, v_lvl=0,
                                        paraview=true,
                                        optargs...)

    # Maneuver to perform
    Vcruise = 0.25 * 125*0.44704            # Cruise speed
    RPMh_w = 400.0                          # RPM of main wing rotors in hover
    telapsed = 30.0                         # Total time to perform maneuver
    nsteps = 100                            # Time steps

    # # Maneuver to perform
    # Vcruise = 0.25 * 125*0.44704            # Cruise speed
    # RPMh_w = 1200.0                         # RPM of main wing rotors in hover
    # telapsed = 30.0                         # Total time to perform maneuver
    # nsteps = 9000                           # Time steps

    # # Maneuver to perform
    # Vcruise = 0.25 * 125*0.44704            # Cruise speed
    # RPMh_w = 1200.0                         # RPM of main wing rotors in hover
    # telapsed = 30.0                         # Total time to perform maneuver
    # nsteps = 27000                           # Time steps

    # # Maneuver to perform
    # Vcruise = 0.125 * 125*0.44704            # Cruise speed
    # RPMh_w = 600.0                          # RPM of main wing rotors in hover
    # telapsed = 60.0                         # Total time to perform maneuver
    # nsteps = 9000                           # Time steps

    dt = telapsed/nsteps

    # Generate maneuver
    maneuver = generate_maneuver_vahana1(; optargs...)

    # Plot maneuver path and controls
    uns.plot_maneuver(maneuver; tstages=[0.2, 0.3, 0.5, 0.6])

    # Generate geometry
    (vehicle, grounds) = generate_geometry_vahana(; n_factor=1,
                                                    xfoil=false,
                                                    data_path=data_path,
                                                    run_name=run_name,
                                                    optargs...)

    # Simulation setup
    Vref = Vcruise
    RPMref = RPMh_w
    ttot = telapsed

    Vinit = Vref*maneuver.Vvehicle(0)       # (m/s) initial vehicle velocity
                                            # (rad/s) initial vehicle angular velocity
    Winit = pi/180 * (maneuver.anglevehicle(0+1e-12)-
                                          maneuver.anglevehicle(0))/(ttot*1e-12)

    simulation = uns.Simulation(vehicle, maneuver, Vref, RPMref, ttot;
                                                    Vinit=Vinit, Winit=Winit)

    save_vtk_optsargs = [(:save_horseshoes, false)]

    # Visualize maneuver
    strn = uns.visualize_kinematics(simulation, nsteps, save_path;
                                    run_name=run_name,
                                    save_vtk_optsargs=save_vtk_optsargs,
                                    prompt=prompt, verbose=verbose, v_lvl=v_lvl,
                                    paraview=false
                                    )

    # Move landing pad to landing area
    vlm.vtk.lintransform!(grounds[2], eye(3), Vcruise*telapsed*[-0.25, 0, -0.0025])

    # Save ground
    for (i, ground) in enumerate(grounds)
        gt.save(ground, run_name*"_Ground$i"; path=save_path)
        strn *= run_name*"_Ground$i.vtk;"
    end

    # Call paraview
    if paraview
        run(`paraview --data="$save_path/$strn"`)
    end

    return strn
end

"""
    Generates geometry of Vahana aircraft, saves it as vtk files, and calls
    Paraview visualizing the VTK geometry.
"""
function visualize_geometry_vahana(; save_path=extdrive_path*"vahana_geometry00/",
                                     prompt=true, run_name="vahana", optargs...)

    (vehicle, grounds) = generate_geometry_vahana(; n_factor=1,
                                                     xfoil=false,
                                                     data_path=data_path,
                                                     run_name=run_name,
                                                     optargs...)

    # Setup dummy freestream and RPMs for horseshoe visualization
    Vinf(x,t) = [1,0,0]
    vlm.setVinf(vehicle.system, Vinf)

    for rotor_system in vehicle.rotor_systems
        for rotor in rotor_system
            vlm.setRPM(rotor, 6000)
        end
    end

    # Create save path
    gt.create_path(save_path, prompt)

    # Save vehicle
    strn = uns.save_vtk(vehicle, run_name; path=save_path, save_horseshoes=false)

    # Save ground
    for (i, ground) in enumerate(grounds)
        gt.save(ground, run_name*"_Ground$i"; path=save_path)
        strn *= run_name*"_Ground$i.vtk;"
    end

    # Call Paraview
    run(`paraview --data="$save_path/$strn"`)
end


text:
using FactCheck
FactCheck.clear_results()
if !isdefined(:throw_on_failure)
  throw_on_failure = true
end

module PkgTests
include("tools.jl")

include("test_internal.jl")
include("test_math.jl")
include("test_statistics.jl")
include("test_bases.jl")
include("test_quadrature.jl")
include("test_gpc.jl")
end

if throw_on_failure
  FactCheck.exitstatus()
end


text:
function generateCln(BiasDict13, f)
    b1 = BiasDict13["δ3"]
    bη = BiasDict13["η3"]
    bK² = BiasDict13["2tr[KK2]"]
    bδη = BiasDict13["δη2"]
    bη² = BiasDict13["2ηη2"]
    bKKpara = BiasDict13["2KK2∥"]
    btd = BiasDict13["Otd"]
    bδΠ2para = BiasDict13["δΠ2∥"]
    bηΠ2para = BiasDict13["ηΠ2∥"]
    bΠ2Kpara = BiasDict13["Π2K∥"]
    bΠ2para = -BiasDict13["s^k∂kΠ2∥"]
    bΠ3para = BiasDict13["Π3∥"] - 2bΠ2para

    Cln = zeros(Float64, 7, 5)

    Cln[1,1] = (2*b1^2)/21. + (10*b1*bK²)/7. + (4*b1*btd)/7. - (8*f*b1*bη)/63. -
    (10*f*bK²*bη)/21. - (4*f*btd*bη)/21. + (2*f^2*bη^2)/35. + (6*f^3*bη^2)/245. +
    (4*f^2*b1*bη²)/35. - (12*f^3*bη*bη²)/245. + (10*b1*bKKpara)/21. -
    (4*f*bη*bKKpara)/21. - (2*f*b1*bηΠ2para)/21. + (2*f^2*bη*bηΠ2para)/49. +
    (4*b1*bΠ2para)/63. + (4*f*bη*bΠ2para)/105. - (f^2*bη*bΠ2para)/49. +
    (5*b1*bΠ2Kpara)/21. - (2*f*bη*bΠ2Kpara)/21. + (2*b1*bΠ3para)/63. + (3*f*bη*bΠ3para)/70.

    Cln[1,2] = b1^2/14. + (f*b1^2)/7. - (15*b1*bK²)/7. - (6*b1*btd)/7. + (f*b1*bη)/21. -
    (6*f^2*b1*bη)/35. + (5*f*bK²*bη)/7. + (2*f*btd*bη)/7. - (3*f^2*bη^2)/70. +
    (6*f^3*bη^2)/245. - (6*f^2*b1*bη²)/35. + (18*f^3*bη*bη²)/245. -
    (5*b1*bKKpara)/7. + (2*f*bη*bKKpara)/7. + (f*b1*bηΠ2para)/7. -
    (3*f^2*bη*bηΠ2para)/49. - (2*b1*bΠ2para)/21. - (2*f*bη*bΠ2para)/35. +
    (3*f^2*bη*bΠ2para)/98. - (5*b1*bΠ2Kpara)/14. + (f*bη*bΠ2Kpara)/7. -
    (b1*bΠ3para)/21. - (9*f*bη*bΠ3para)/140.

    Cln[1,3] = -b1^2/6. - (f*b1^2)/7. + (5*b1*bK²)/7. + (2*b1*btd)/7. + (5*f*b1*bη)/63. +
    (2*f^2*b1*bη)/35. - (5*f*bK²*bη)/21. - (2*f*btd*bη)/21. -
    (4*f^2*bδη*bη)/35. - (f^2*bη^2)/70. + (12*f^3*bη^2)/245. +
    (2*f^2*b1*bη²)/35. + (18*f^3*bη*bη²)/245. + (5*b1*bKKpara)/21. -
    (2*f*bη*bKKpara)/63. + (4*f*bη*bδΠ2para)/21. - (f*b1*bηΠ2para)/21. -
    (3*f^2*bη*bηΠ2para)/49. + (2*b1*bΠ2para)/63. + (6*f*bη*bΠ2para)/35. -
    (5*f^2*bη*bΠ2para)/98. + (5*b1*bΠ2Kpara)/42. - (f*bη*bΠ2Kpara)/63. +
    (b1*bΠ3para)/63. + (61*f*bη*bΠ3para)/420.

    Cln[1,4] = (4*f^2*b1*bη)/35. + (4*f^2*bδη*bη)/35. - (24*f^3*bη^2)/245. -
    (24*f^3*bη*bη²)/245. - (4*f*bη*bKKpara)/63. - (4*f*bη*bδΠ2para)/21. +
    (4*f^2*bη*bηΠ2para)/49. - (16*f*bη*bΠ2para)/105. + (2*f^2*bη*bΠ2para)/49. -
    (2*f*bη*bΠ2Kpara)/63. - (13*f*bη*bΠ3para)/105.

    Cln[3,1] = (-16*f*b1*bη)/63. - (9*f^2*b1*bη)/49. - (20*f*bK²*bη)/21. -
    (8*f*btd*bη)/21. + (8*f^2*bη^2)/49. + (f^3*bη^2)/7. + (4*f^2*b1*bη²)/49. -
    (4*f^3*bη*bη²)/49. + (5*b1*bKKpara)/21. - (65*f*bη*bKKpara)/147. -
    (10*f*b1*bηΠ2para)/147. + (10*f^2*bη*bηΠ2para)/147. - (4*b1*bΠ2para)/9. +
    (15*f*b1*bΠ2para)/98. + (4*f*bη*bΠ2para)/21. - (5*f^2*bη*bΠ2para)/42. +
    (5*b1*bΠ2Kpara)/42. - (65*f*bη*bΠ2Kpara)/294. - (101*b1*bΠ3para)/252. +
    (37*f*bη*bΠ3para)/196.

    Cln[3,2] = (2*f*b1^2)/7. + (2*f*b1*bη)/21. - (3*f^2*b1*bη)/14. + (10*f*bK²*bη)/7. +
    (4*f*btd*bη)/7. - (6*f^2*bη^2)/49. - (f^3*bη^2)/98. - (6*f^2*b1*bη²)/49. +
    (6*f^3*bη*bη²)/49. - (5*b1*bKKpara)/14. + (65*f*bη*bKKpara)/98. +
    (5*f*b1*bηΠ2para)/49. - (5*f^2*bη*bηΠ2para)/49. + (2*b1*bΠ2para)/3. -
    (45*f*b1*bΠ2para)/196. - (2*f*bη*bΠ2para)/7. + (5*f^2*bη*bΠ2para)/28. -
    (5*b1*bΠ2Kpara)/28. + (65*f*bη*bΠ2Kpara)/196. + (101*b1*bΠ3para)/168. -
    (111*f*bη*bΠ3para)/392.

    Cln[3,3] = (4*f*b1^2)/7. + (6*f*b1*bδη)/7. + (10*f*b1*bη)/63. - (11*f^2*b1*bη)/14. -
    (10*f*bK²*bη)/21. - (4*f*btd*bη)/21. - (22*f^2*bδη*bη)/49. - (2*f^2*bη^2)/49. +
    (11*f^3*bη^2)/98. - (34*f^2*b1*bη²)/49. + (10*f^3*bη*bη²)/49. - (5*b1*bKKpara)/14. +
    (25*f*bη*bKKpara)/882. - (10*b1*bδΠ2para)/7. + (110*f*bη*bδΠ2para)/147. +
    (85*f*b1*bηΠ2para)/147. - (25*f^2*bη*bηΠ2para)/147. - (86*b1*bΠ2para)/63. +
    (75*f*b1*bΠ2para)/196. + (34*f*bη*bΠ2para)/49. - (5*f^2*bη*bΠ2para)/196. -
    (5*b1*bΠ2Kpara)/28. + (25*f*bη*bΠ2Kpara)/1764. - (569*b1*bΠ3para)/504. +
    (683*f*bη*bΠ3para)/1176.

    Cln[3,4] = (-6*f*b1^2)/7. - (6*f*b1*bδη)/7. + (58*f^2*b1*bη)/49. +
    (22*f^2*bδη*bη)/49. + (4*f^3*bη^2)/49. + (36*f^2*b1*bη²)/49. +
    (4*f^3*bη*bη²)/49. + (10*b1*bKKpara)/21. - (110*f*bη*bKKpara)/441. +
    (10*b1*bδΠ2para)/7. - (110*f*bη*bδΠ2para)/147. - (30*f*b1*bηΠ2para)/49. -
    (10*f^2*bη*bηΠ2para)/147. + (8*b1*bΠ2para)/7. - (15*f*b1*bΠ2para)/49. -
    (88*f*bη*bΠ2para)/147. - (15*f^2*bη*bΠ2para)/49. + (5*b1*bΠ2Kpara)/21. -
    (55*f*bη*bΠ2Kpara)/441. + (13*b1*bΠ3para)/14. - (143*f*bη*bΠ3para)/294.

    Cln[3,5] = (-16*f^3*bη^2)/49. - (16*f^3*bη*bη²)/49. + (40*f^2*bη*bηΠ2para)/147. +
    (40*f^2*bη*bΠ2para)/147.

    Cln[5,1] = (-12*f^2*b1*bη)/49. + (16*f^2*bη^2)/245. + (72*f^3*bη^2)/385. -
    (48*f^2*b1*bη²)/245. + (192*f^3*bη*bη²)/2695. - (4*f*bη*bKKpara)/49. +
    (8*f*b1*bηΠ2para)/49. - (32*f^2*bη*bηΠ2para)/539. + (10*f*b1*bΠ2para)/49. +
    (16*f*bη*bΠ2para)/105. - (12*f^2*bη*bΠ2para)/77. - (2*f*bη*bΠ2Kpara)/49. +
    (101*f*bη*bΠ3para)/735.

    Cln[5,2] = (6*f^2*b1*bη)/35. - (12*f^2*bη^2)/245. - (36*f^3*bη^2)/245. +
    (72*f^2*b1*bη²)/245. - (288*f^3*bη*bη²)/2695. + (6*f*bη*bKKpara)/49. -
    (12*f*b1*bηΠ2para)/49. + (48*f^2*bη*bηΠ2para)/539. - (15*f*b1*bΠ2para)/49. -
    (8*f*bη*bΠ2para)/35. + (18*f^2*bη*bΠ2para)/77. + (3*f*bη*bΠ2Kpara)/49. -
    (101*f*bη*bΠ3para)/490.

    Cln[5,3] = (18*f^2*b1*bη)/35. - (72*f^2*bδη*bη)/245. - (4*f^2*bη^2)/245. -
    (432*f^3*bη^2)/2695. + (156*f^2*b1*bη²)/245. - (228*f^3*bη*bη²)/2695. +
    (6*f*bη*bKKpara)/49. + (24*f*bη*bδΠ2para)/49. - (26*f*b1*bηΠ2para)/49. +
    (38*f^2*bη*bηΠ2para)/539. - (45*f*b1*bΠ2para)/49. + (344*f*bη*bΠ2para)/735. +
    (180*f^2*bη*bΠ2para)/539. + (3*f*bη*bΠ2Kpara)/49. + (569*f*bη*bΠ3para)/1470.

    Cln[5,4] = (-528*f^2*b1*bη)/245. + (72*f^2*bδη*bη)/245. + (2664*f^3*bη^2)/2695. -
    (120*f^2*b1*bη²)/49. + (2664*f^3*bη*bη²)/2695. - (8*f*bη*bKKpara)/49. -
    (24*f*bη*bδΠ2para)/49. + (100*f*b1*bηΠ2para)/49. - (444*f^2*bη*bηΠ2para)/539. +
    (120*f*b1*bΠ2para)/49. - (96*f*bη*bΠ2para)/245. - (612*f^2*bη*bΠ2para)/539. -
    (4*f*bη*bΠ2Kpara)/49. - (78*f*bη*bΠ3para)/245.

    Cln[5,5] = (12*f^2*b1*bη)/7. - (468*f^3*bη^2)/539. + (12*f^2*b1*bη²)/7. -
    (468*f^3*bη*bη²)/539. - (10*f*b1*bηΠ2para)/7. + (390*f^2*bη*bηΠ2para)/539. -
    (10*f*b1*bΠ2para)/7. + (390*f^2*bη*bΠ2para)/539.

    Cln[7,1] = (40*f^3*bη^2)/539. + (32*f^3*bη*bη²)/539. - (80*f^2*bη*bηΠ2para)/1617. -
    (100*f^2*bη*bΠ2para)/1617.

    Cln[7,2] = (-4*f^3*bη^2)/49. - (48*f^3*bη*bη²)/539. + (40*f^2*bη*bηΠ2para)/539. +
    (50*f^2*bη*bΠ2para)/539.

    Cln[7,3] = (-116*f^3*bη^2)/539. - (104*f^3*bη*bη²)/539. + (260*f^2*bη*bηΠ2para)/1617. +
    (150*f^2*bη*bΠ2para)/539.

    Cln[7,4] = (400*f^3*bη^2)/539. + (400*f^3*bη*bη²)/539. - (1000*f^2*bη*bηΠ2para)/1617. -
    (400*f^2*bη*bΠ2para)/539.

    Cln[7,5] = (-40*f^3*bη^2)/77. - (40*f^3*bη*bη²)/77. + (100*f^2*bη*bηΠ2para)/231. +
    (100*f^2*bη*bΠ2para)/231.

    return Cln
end

max_stars_repo_name: JosephTomlinson/GeneralBiasPk

max_issues_repo_name: JosephTomlinson/GeneralBiasPk

max_forks_repo_name: JosephTomlinson/GeneralBiasPk


text:
export tbaa_make_child, tbaa_addrspace

function tbaa_make_child(name::String, ctx::LLVM.Context=JuliaContext(); constant::Bool=false)
    tbaa_root = MDNode([MDString("custom_tbaa", ctx)], ctx)
    tbaa_struct_type =
        MDNode([MDString("custom_tbaa_$name", ctx),
                tbaa_root,
                LLVM.ConstantInt(0, ctx)], ctx)
    tbaa_access_tag =
        MDNode([tbaa_struct_type,
                tbaa_struct_type,
                LLVM.ConstantInt(0, ctx),
                LLVM.ConstantInt(constant ? 1 : 0, ctx)], ctx)

    return tbaa_access_tag
end

tbaa_addrspace(as, ctx::LLVM.Context=JuliaContext()) = tbaa_make_child("addrspace($(as))", ctx)


text:
#==============================================================================#
# SNS/test/runtests.jl
#
# Copyright OC Technology Pty Ltd 2014 - All rights reserved
#==============================================================================#

using AWSCore
using AWSLambda
using AWSSNS
using Retry
using SymDict
using JSON
using Base.Test
using Primes

AWSCore.set_debug_level(1)

try
    AWSLambda.@lambda_eval run(`uname`)
catch e
    if ecode(e) == "404"
        AWSLambda.deploy_jl_lambda_eval()
    else
        rethrow(e)
    end
end

#-------------------------------------------------------------------------------
# Lambda tests using base jl_labda_eval function
#-------------------------------------------------------------------------------

f = AWSLambda.lambda_function(readstring)
@test chomp(f(`uname`)) == "Linux"

f = AWSLambda.lambda_function(sum)
@test f([1,2,3.5]) == 6.5

f = AWSLambda.lambda_function((a::String, b::Int) -> begin
    repeat(a, b)
end)
@test f("FOO", 2) == "FOOFOO"

module Foo

    import ..AWSLambda

    f = AWSLambda.lambda_function(eval(Main, :((a::String, b::Int) -> begin
        repeat(a, b)
    end)))

end
@test Foo.f("FOO", 2) == "FOOFOO"

AWSLambda.@lambda function foo(a::String, b::Int)
    repeat(a, b)
end
@test foo("FOO", 2) == "FOOFOO"

AWSLambda.@lambda function lreadstring(a)
    readstring(a)
end
@test chomp(lreadstring(`uname`)) == "Linux"

module Foo2
    using AWSLambda
    AWSLambda.@lambda function foo(a::Cmd)
        chomp(readstring(a))
    end
end
@test Foo2.foo(`uname`) == "Linux"

@test AWSLambda.lambda_eval(quote
        open("/proc/cpuinfo") do io
            return Dict(strip(k) => strip(v)
                        for (k, v) in (split(l, ":")
                        for l in Iterators.filter(x -> contains(x, ":"), eachline(io))))
        end
    end)["model name"] in ["Intel(R) Xeon(R) CPU E5-2666 v3 @ 2.90GHz",
                           "Intel(R) Xeon(R) CPU E5-2676 v3 @ 2.40GHz",
                           "Intel(R) Xeon(R) CPU E5-2680 v2 @ 2.80GHz",
                           "Intel(R) Xeon(R) CPU E5-2686 v4 @ 2.30GHz"]

@test AWSLambda.lambda_include_string("1 + 1") == 2

@test (AWSLambda.@lambda_eval begin
    HTTP.header(HTTP.get("http://httpbin.org/ip"), "Content-Type")
end) == "application/json"

r = AWSLambda.@lambda_eval begin

    module Foo

        export foo

        using HTTP
        using JSON

        const url = "http://httpbin.org/ip"

        function foo()
            JSON.parse(String(HTTP.get(url).body))
        end
    end

    using .Foo

    foo()
end

@test ismatch(r"^[0-9.]+$", r["origin"])

# Count primes in the cloud...

AWSLambda.@lambda function count_primes(low::Int, high::Int)

    function is_prime(n)
        if n ≤ 1
            return false
        elseif n ≤ 3
            return true
        elseif n % 2 == 0 || n % 3 == 0
            return false
        end
        i = 5
        while i * i ≤ n
            if n % i == 0 || n % (i + 2) == 0
                return false
            end
            i += 6
        end
        return true
    end

    c = count(is_prime, low:high)
    println("$c primes between $low and $high.")
    return c
end

@test count_primes(10, 100) == 21

                                 if !haskey(ENV, "AWS_LAMBDA_JL_TEST_SKIP_SLOW")
# Count primes in parallel...
function pcount_primes(low::Int, high::Int)
    w = 5000000
    counts = asyncmap(x->count_primes(x...),
                      [(i, min(high,i + w)) for i = low:w:high])
    count = sum(counts)
    println("$count primes between $low and $high.")
    return count
end

@test pcount_primes(10, 100000000) == 5761451
                                                                             end
AWSLambda.@lambda using Primes function count_primes2(low::Int, high::Int)

    c = length(Primes.primes(low, high))
    println("$c primes between $low and $high.")
    return c
end

@test count_primes2(10, 100) == 21

                                 if !haskey(ENV, "AWS_LAMBDA_JL_TEST_SKIP_SLOW")
# Count primes in parallel...
function pcount_primes2(low::Int, high::Int)
    w = 100000000
    counts = asyncmap(x->count_primes2(x...),
                      [(i, min(high,i + w)) for i = low:w:high])
    count = sum(counts)
    println("$count primes between $low and $high.")
    return count
end

@test pcount_primes2(10, 10^10) == 455052507
                                                                             end

#@test invoke_lambda(aws, "count_primes", low = 10, high = 100)[:jl_data] == "21"

mktempdir() do tmp
    cd(tmp) do

        # Create a test module under "tmp"...

        mkpath("TestModule")
        open(io->write(io, """
            __precompile__()

            module TestModule

            export test_function

            test_function(x) = x * x

            end
        """), "TestModule/TestModule.jl", "w")

        push!(LOAD_PATH, "TestModule")

        eval(Main, :(using TestModule))
        @test Base.invokelatest(test_function,5) == 25

        # Create a lambda that uses the TestModule...
        eval(:(AWSLambda.@lambda using TestModule function lambda_test(x)
            return test_function(x)
        end))

        @test Base.invokelatest(lambda_test,4) == 16

        # Create a lambda that uses the TestModule (with explicit aws config)...
        eval(:(AWSLambda.@lambda using TestModule function lambda_test2(x)
            return test_function(x)
        end AWSCore.default_aws_config()))

        @test Base.invokelatest(lambda_test2,5) == 25
    end
end

#-------------------------------------------------------------------------------
# Lambda tests using deployed lambda functions
#-------------------------------------------------------------------------------

AWSLambda.@deploy [
   :MemorySize => 1024,
   :Timeout => 60
] function count_primes_slow(low::Int, high::Int)

    function is_prime(n)
        if n ≤ 1
            return false
        elseif n ≤ 3
            return true
        elseif n % 2 == 0 || n % 3 == 0
            return false
        end
        i = 5
        while i * i ≤ n
            if n % i == 0 || n % (i + 2) == 0
                return false
            end
            i += 6
        end
        return true
    end

    c = count(is_prime, low:high)
    println("$c primes between $low and $high.")
    return c
end

@test AWSLambda.invoke_jl_lambda("count_primes_slow", 10, 100) == 21
r = AWSLambda.invoke_lambda("count_primes_slow", low=10, high=100)
@test r == 21

AWSLambda.@deploy using Primes function count_primes_fast(low::Int, high::Int)

    c = length(Primes.primes(low, high))
    println("$c primes between $low and $high.")
    return c
end

@test AWSLambda.invoke_jl_lambda("count_primes_fast", 10, 100) == 21

mktempdir() do tmp
    cd(tmp) do
        run(`aws lambda invoke --function-name count_primes_fast
                               --payload "{\"low\": 10, \"high\": 100}"
                               output.txt`)
        r = JSON.parse(readstring("output.txt"))

        @test r == 21
    end
end

#==============================================================================#
# End of file.
#==============================================================================#


text:
struct StaticObjectDetector{TDetection <: AbstractObjectDetection} <: AbstractObjectDetector{TDetection}
    detections::Vector{TDetection}
end

function detect(detector::StaticObjectDetector, image)
    return detector.detections
end


text:
@testset "238.product-of-array-except-self.jl" begin
    @test product_except_self([1, 2, 3, 4]) == [24, 12, 8, 6]
    @test product_except_self([1, 2, 3, 4, 0]) == [0, 0, 0, 0, 24]
    @test product_except_self([1, 2, 0, 4, 0]) == [0, 0, 0, 0, 0]
end


text:
include("./imports.jl")
include("cacliumscoring.jl")


text: include("layers.jl")


text:
export RandomFourier, FT

"""
    RandomFourier(phases = true)

A surrogate that randomizes the Fourier components
of the signal in some manner. If `phases==true`, the phases are randomized,
otherwise the amplitudes are randomized. `FT` is an alias for `RandomFourier`.

Random Fourier phase surrogates[^Theiler1991] preserve the
autocorrelation function, or power spectrum, of the original signal.
Random Fourier amplitude surrogates preserve the mean and autocorrelation
function but do not preserve the variance of the original. Random
amplitude surrogates are not common in the literature, but are provided
for convenience.

Random phase surrogates can be used to test the null hypothesis that
the original signal was produced by a linear Gaussian process [^Theiler1991].

[^Theiler1991]: J. Theiler, S. Eubank, A. Longtin, B. Galdrikian, J. Farmer, Testing for nonlinearity in time series: The method of surrogate data, Physica D 58 (1–4) (1992) 77–94.
"""
struct RandomFourier <: Surrogate
    phases::Bool
end
RandomFourier() = RandomFourier(true)
const FT = RandomFourier

function surrogenerator(x::AbstractVector, rf::RandomFourier, rng = Random.default_rng())
    forward = plan_rfft(x)
    inverse = plan_irfft(forward*x, length(x))
    m = mean(x)
    𝓕 = forward*(x .- m)
    shuffled𝓕 = zero(𝓕)
    s = similar(x)
    n = length(𝓕)
    r = abs.(𝓕)
    ϕ = angle.(𝓕)
    coeffs = zero(r)

    init = (inverse = inverse, m = m, coeffs = coeffs, n = n, r = r,
            ϕ = ϕ, shuffled𝓕 = shuffled𝓕)
    return SurrogateGenerator(rf, x, s, init, rng)
end

function (sg::SurrogateGenerator{<:RandomFourier})()
    inverse, m, coeffs, n, r, ϕ, shuffled𝓕 =
        getfield.(Ref(sg.init),
        (:inverse, :m, :coeffs, :n, :r, :ϕ, :shuffled𝓕))
    s, rng, phases = sg.s, sg.rng, sg.method.phases

    if phases
        coeffs .= rand(rng, Uniform(0, 2π), n)
        shuffled𝓕 .= r .* exp.(coeffs .* 1im)
    else
        coeffs .= r .* rand(rng, Uniform(0, 2π), n)
        shuffled𝓕 .= coeffs .* exp.(ϕ .* 1im)
    end
    s .= inverse * shuffled𝓕 .+ m
    return s
end


text:
module LatBoTests
using LatBo
using FactCheck: facts, context, @fact, not, roughly, exitstatus, exactly, @fact_throws

include("indices/cached.jl")
# include("units.jl")
include("playground.jl")
include("geometry.jl")
include("indices/indices.jl")
include("lb/lattice.jl")
include("lb/thermodynamics.jl")
include("lb/collision.jl")
include("lb/streaming.jl")
include("lb/iolet.jl")
include("lb/initialization.jl")
include("lb/local_kernel.jl")
exitstatus()
end


text:
"""
$(TYPEDEF)

Extended state definition used within integrators when there are
`VariableRateJump`s in a system. For detailed examples and usage information see
the
- [Tutorial](https://diffeq.sciml.ai/stable/tutorials/discrete_stochastic_example/)

### Fields

$(FIELDS)

## Examples
```julia
using DiffEqJump, OrdinaryDiffEq
f(du,u,p,t) = du .= 0
rate(u,p,t) = (1+t)*u[1]*u[2]

# suppose we wish to decrease each of the two variables by one
# when a jump occurs
function affect!(integrator)
   # Method 1, direct indexing works like normal
   integrator.u[1] -= 1
   integrator.u[2] -= 1

   # Method 2, if we want to broadcast or use array operations we need
   # to access integrator.u.u which is the actual state object.
   # So equivalently to above we could have said:
   # integrator.u.u .-= 1
end

u0 = [10.0, 10.0]
vrj = VariableRateJump(rate, affect!)
oprob = ODEProblem(f, u0, (0.0,2.0))
jprob = JumpProblem(oprob, Direct(), vrj)
sol = solve(jprob,Tsit5())
```

## Notes
- If `ueja isa ExtendedJumpArray` with `ueja.u` of size `N` and `ueja.jump_u` of
  size `num_variableratejumps` then
  ```julia
  # for 1 <= i <= N
  ueja[i] == ueja.u[i]

  # for N < i <= (N+num_variableratejumps)
  ueja[i] == ueja.jump_u[i]
  ```
- In a system with `VariableRateJump`s all callback, `ConstantRateJump`, and
  `VariableRateJump` `affect!` functions will receive integrators with
  `integrator.u` an `ExtendedJumpArray`.
- As such, `affect!` functions that wish to modify the state via vector
  operations should use `ueja.u.u` to obtain the aliased state object.
"""
struct ExtendedJumpArray{T3<:Number, T1, T<:AbstractArray{T3,T1},T2} <: AbstractArray{T3,1}
  """The current state."""
  u::T
  """The current rate (i.e. hazard, intensity, or propensity) values for the `VariableRateJump`s."""
  jump_u::T2
end

Base.length(A::ExtendedJumpArray) = length(A.u) + length(A.jump_u)
Base.size(A::ExtendedJumpArray) = (length(A),)
@inline function Base.getindex(A::ExtendedJumpArray,i::Int)
  i <= length(A.u) ? A.u[i] : A.jump_u[i-length(A.u)]
end
@inline function Base.getindex(A::ExtendedJumpArray,I::Int...)
  prod(I) <= length(A.u) ? A.u[I...] : A.jump_u[prod(I)-length(A.u)]
end
@inline function Base.getindex(A::ExtendedJumpArray,I::CartesianIndex{1})
  A[I[1]]
end
@inline Base.setindex!(A::ExtendedJumpArray,v,I...) = (A[CartesianIndices(A.u,I...)] = v)
@inline Base.setindex!(A::ExtendedJumpArray,v,I::CartesianIndex{1}) = (A[I[1]] = v)
@inline function Base.setindex!(A::ExtendedJumpArray,v,i::Int)
  i <= length(A.u) ? (A.u[i] = v) : (A.jump_u[i-length(A.u)] = v)
end

Base.IndexStyle(::Type{<:ExtendedJumpArray}) = IndexLinear()
Base.similar(A::ExtendedJumpArray) = ExtendedJumpArray(similar(A.u),similar(A.jump_u))
Base.similar(A::ExtendedJumpArray,::Type{S}) where {S} = ExtendedJumpArray(similar(A.u,S),similar(A.jump_u,S))
Base.zero(A::ExtendedJumpArray) = fill!(similar(A),0)

# Required for non-diagonal noise
LinearAlgebra.mul!(c::ExtendedJumpArray,A::AbstractVecOrMat,u::AbstractVector) = mul!(c.u,A,u)

# Ignore axes
Base.similar(A::ExtendedJumpArray,::Type{S},axes::Tuple{Base.OneTo{Int}}) where {S} = ExtendedJumpArray(similar(A.u,S),similar(A.jump_u,S))

# Stiff ODE solver
function ArrayInterfaceCore.zeromatrix(A::ExtendedJumpArray)
  u = [vec(A.u);vec(A.jump_u)]
  u .* u' .* false
end
LinearAlgebra.ldiv!(A,b::ExtendedJumpArray) = LinearAlgebra.ldiv!(A,[vec(b.u);vec(b.jump_u)])

function recursivecopy!(dest::T, src::T) where T<:ExtendedJumpArray
  recursivecopy!(dest.u,src.u)
  recursivecopy!(dest.jump_u,src.jump_u)
end
Base.show(io::IO,A::ExtendedJumpArray) = show(io,A.u)
TreeViews.hastreeview(x::ExtendedJumpArray) = true
plot_indices(A::ExtendedJumpArray) = eachindex(A.u)

## broadcasting

struct ExtendedJumpArrayStyle{Style <: Broadcast.BroadcastStyle} <: Broadcast.AbstractArrayStyle{Any} end
ExtendedJumpArrayStyle(::S) where {S} = ExtendedJumpArrayStyle{S}()
ExtendedJumpArrayStyle(::S, ::Val{N}) where {S,N} = ExtendedJumpArrayStyle(S(Val(N)))
ExtendedJumpArrayStyle(::Val{N}) where N = ExtendedJumpArrayStyle{Broadcast.DefaultArrayStyle{N}}()

# promotion rules
@inline function Broadcast.BroadcastStyle(::ExtendedJumpArrayStyle{AStyle}, ::ExtendedJumpArrayStyle{BStyle}) where {AStyle, BStyle}
    ExtendedJumpArrayStyle(Broadcast.BroadcastStyle(AStyle(), BStyle()))
end
Broadcast.BroadcastStyle(::ExtendedJumpArrayStyle{Style}, ::Broadcast.DefaultArrayStyle{0}) where Style<:Broadcast.BroadcastStyle = ExtendedJumpArrayStyle{Style}()
Broadcast.BroadcastStyle(::ExtendedJumpArrayStyle, ::Broadcast.DefaultArrayStyle{N}) where N = Broadcast.DefaultArrayStyle{N}()

combine_styles(args::Tuple{})         = Broadcast.DefaultArrayStyle{0}()
@inline combine_styles(args::Tuple{Any})      = Broadcast.result_style(Broadcast.BroadcastStyle(args[1]))
@inline combine_styles(args::Tuple{Any, Any}) = Broadcast.result_style(Broadcast.BroadcastStyle(args[1]), Broadcast.BroadcastStyle(args[2]))
@inline combine_styles(args::Tuple)   = Broadcast.result_style(Broadcast.BroadcastStyle(args[1]), combine_styles(Base.tail(args)))

function Broadcast.BroadcastStyle(::Type{ExtendedJumpArray{T,S}}) where {T, S}
    ExtendedJumpArrayStyle(Broadcast.result_style(Broadcast.BroadcastStyle(T)))
end

@inline function Base.copy(bc::Broadcast.Broadcasted{ExtendedJumpArrayStyle{Style}}) where Style
    ExtendedJumpArray(copy(unpack(bc, Val(:u))),copy(unpack(bc, Val(:jump_u))))
end

@inline function Base.copyto!(dest::ExtendedJumpArray, bc::Broadcast.Broadcasted{ExtendedJumpArrayStyle{Style}}) where Style
    copyto!(dest.u,unpack(bc, Val(:u)))
    copyto!(dest.jump_u,unpack(bc, Val(:jump_u)))
    dest
end

# drop axes because it is easier to recompute
@inline unpack(bc::Broadcast.Broadcasted{Style}, i) where Style = Broadcast.Broadcasted{Style}(bc.f, unpack_args(i, bc.args))
@inline unpack(bc::Broadcast.Broadcasted{ExtendedJumpArrayStyle{Style}}, i) where Style = Broadcast.Broadcasted{Style}(bc.f, unpack_args(i, bc.args))
unpack(x,::Any) = x
unpack(x::ExtendedJumpArray, ::Val{:u}) = x.u
unpack(x::ExtendedJumpArray, ::Val{:jump_u}) = x.jump_u

@inline unpack_args(i, args::Tuple) = (unpack(args[1], i), unpack_args(i, Base.tail(args))...)
unpack_args(i, args::Tuple{Any}) = (unpack(args[1], i),)
unpack_args(::Any, args::Tuple{}) = ()


text:
#####
##### High level training procedure
#####

"""
    Env

Type for an AlphZero environment.

The environment features the current neural network, the best neural network
seen so far that is used for data generation, a memory buffer
and an iteration counter.

# Constructor

    Env(game_spec, params, curnn, bestnn=copy(curnn), experience=[], itc=0)

Construct a new AlphaZero environment:
- `game_spec` specified the game being played
- `params` has type [`Params`](@ref)
- `curnn` is the current neural network and has type [`AbstractNetwork`](@ref)
- `bestnn` is the best neural network so far, which is used for data generation
- `experience` is the initial content of the memory buffer
   as a vector of [`TrainingSample`](@ref)
- `itc` is the value of the iteration counter (0 at the start of training)
"""
mutable struct Env{GameSpec, Network, State}
  gspec  :: GameSpec
  params :: Params
  curnn  :: Network
  bestnn :: Network
  memory :: MemoryBuffer{GameSpec, State}
  itc    :: Int
  function Env(
      gspec::AbstractGameSpec,
      params, curnn, bestnn=copy(curnn), experience=[], itc=0)
    msize = max(params.mem_buffer_size[itc], length(experience))
    memory = MemoryBuffer(gspec, msize, experience)
    return new{typeof(gspec), typeof(curnn), GI.state_type(gspec)}(
      gspec, params, curnn, bestnn, memory, itc)
  end
end

#####
##### Training handlers
#####

"""
    Handlers

Namespace for the callback functions that are used during training.
This enables logging, saving and plotting to be implemented separately.
An example handler object is `Session`.

All callback functions take a handler object `h` as their first argument
and sometimes a second argment `r` that consists in a report.

| Callback                    | Comment                                        |
|:----------------------------|:-----------------------------------------------|
| `iteration_started(h)`      | called at the beggining of an iteration        |
| `self_play_started(h)`      | called once per iter before self play starts   |
| `game_played(h)`            | called after each game of self play            |
| `self_play_finished(h, r)`  | sends report: [`Report.SelfPlay`](@ref)        |
| `memory_analyzed(h, r)`     | sends report: [`Report.Memory`](@ref)          |
| `learning_started(h)`       | called at the beginning of the learning phase  |
| `updates_started(h, r)`     | sends report: [`Report.LearningStatus`](@ref)  |
| `updates_finished(h, r)`    | sends report: [`Report.LearningStatus`](@ref)  |
| `checkpoint_started(h)`     | called before a checkpoint evaluation starts   |
| `checkpoint_game_played(h)` | called after each arena game                   |
| `checkpoint_finished(h, r)` | sends report: [`Report.Checkpoint`](@ref)      |
| `learning_finished(h, r)`   | sends report: [`Report.Learning`](@ref)        |
| `iteration_finished(h, r)`  | sends report: [`Report.Iteration`](@ref)       |
| `training_finished(h)`      | called once at the end of training             |
"""
module Handlers

  import ..Report

  function iteration_started(h)      return end
  function self_play_started(h)      return end
  function game_played(h)            return end
  function self_play_finished(h, r)  return end
  function memory_analyzed(h, r)     return end
  function learning_started(h)       return end
  function updates_started(h)        return end
  function updates_finished(h, r)    return end
  function checkpoint_started(h)     return end
  function checkpoint_game_played(h) return end
  function checkpoint_finished(h, r) return end
  function learning_finished(h, r)   return end
  function iteration_finished(h, r)  return end
  function training_finished(h)      return end

end

#####
##### Public utilities
#####

"""
    get_experience(env::Env)

Return the content of the agent's memory as a vector of [`TrainingSample`](@ref).
"""
get_experience(env::Env) = get_experience(env.memory)

"""
    initial_report(env::Env)

Return a report summarizing the configuration of agent before training starts,
as an object of type [`Report.Initial`](@ref).
"""
function initial_report(env::Env)
  num_network_parameters = Network.num_parameters(env.curnn)
  num_reg_params = Network.num_regularized_parameters(env.curnn)
  player = MctsPlayer(env.gspec, env.curnn, env.params.self_play.mcts)
  mcts_footprint_per_node = MCTS.memory_footprint_per_node(env.gspec)
  errs, warns = check_params(env.gspec, env.params)
  return Report.Initial(
    num_network_parameters, num_reg_params, mcts_footprint_per_node, errs, warns)
end

#####
##### Evaluating networks
#####

# Have a "contender" network play against a "baseline" network (params::ArenaParams)
# Return (rewards vector, redundancy)
# Version for two-player games
function pit_networks(gspec, contender, baseline, params, handler)
  make_oracles() = (
    Network.copy(contender, on_gpu=params.sim.use_gpu, test_mode=true),
    Network.copy(baseline, on_gpu=params.sim.use_gpu, test_mode=true))
  simulator = Simulator(make_oracles, record_trace) do oracles
    white = MctsPlayer(gspec, oracles[1], params.mcts)
    black = MctsPlayer(gspec, oracles[2], params.mcts)
    return TwoPlayers(white, black)
  end
  samples = simulate(
    simulator, gspec, params.sim,
    game_simulated=(() -> Handlers.checkpoint_game_played(handler)))
  return rewards_and_redundancy(samples, gamma=params.mcts.gamma)
end

# Evaluate a single neural network for a one-player game (params::ArenaParams)
function evaluate_network(gspec, net, params, handler)
  make_oracles() = Network.copy(net, on_gpu=params.sim.use_gpu, test_mode=true)
  simulator = Simulator(make_oracles, record_trace) do oracle
    MctsPlayer(gspec, oracle, params.mcts)
  end
  samples = simulate(
    simulator, gspec, params.sim,
    game_simulated=(() -> Handlers.checkpoint_game_played(handler)))
  return rewards_and_redundancy(samples, gamma=params.mcts.gamma)
end

# Compare two versions of a neural network (params::ArenaParams)
# Works for both two-player and single-player games
function compare_networks(gspec, contender, baseline, params, handler)
  legend = "Most recent NN versus best NN so far"
  if GI.two_players(gspec)
    (rewards_c, red), t =
      @timed pit_networks(gspec, contender, baseline, params, handler)
    avgr = mean(rewards_c)
    rewards_b = nothing
  else
    (rewards_c, red_c), tc = @timed evaluate_network(gspec, contender, params, handler)
    (rewards_b, red_b), tb = @timed evaluate_network(gspec, baseline, params, handler)
    avgr = mean(rewards_c) - mean(rewards_b)
    red = mean([red_c, red_b])
    t = tc + tb
  end
  return Report.Evaluation(legend, avgr, red, rewards_c, rewards_b, t)
end

#####
##### Training Loop
#####

function resize_memory!(env::Env, n)
  exp = get_experience(env.memory)
  env.memory = MemoryBuffer(env.gspec, n, exp)
  return
end

function dummy_learning_report()
  eps = 1e-16
  dummy_loss = Report.Loss(0, 0, 0, 0, 0)
  dummy_status = Report.LearningStatus(dummy_loss, 0, 0)
  return Report.Learning(eps, eps, eps, eps, dummy_status, [], [], false)
end

function learning_step!(env::Env, handler)
  ap = env.params.arena
  lp = env.params.learning
  checkpoints = Report.Checkpoint[]
  losses = Float32[]
  tloss, teval, ttrain = 0., 0., 0.
  experience = get_experience(env.memory)
  if env.params.use_symmetries
    experience = augment_with_symmetries(env.gspec, experience)
  end
  if isempty(experience)
    # Skipping the learning phase
    return dummy_learning_report()
  end
  trainer, tconvert = @timed Trainer(env.gspec, env.curnn, experience, lp)
  init_status = learning_status(trainer)
  status = init_status
  Handlers.learning_started(handler)
  # Compute the number of batches between each checkpoint
  nbatches = lp.max_batches_per_checkpoint
  if !iszero(lp.min_checkpoints_per_epoch)
    ntotal = num_batches_total(trainer)
    nbatches = min(nbatches, ntotal ÷ lp.min_checkpoints_per_epoch)
  end
  # Loop state variables
  best_evalr = isnothing(ap) ? nothing : ap.update_threshold
  nn_replaced = false

  for k in 1:lp.num_checkpoints
    # Execute a series of batch updates
    Handlers.updates_started(handler, status)
    dlosses, dttrain = @timed batch_updates!(trainer, nbatches)
    status, dtloss = @timed learning_status(trainer)
    Handlers.updates_finished(handler, status)
    tloss += dtloss
    ttrain += dttrain
    append!(losses, dlosses)
    # Run a checkpoint evaluation if the arena parameter is provided
    if isnothing(ap)
      env.curnn = get_trained_network(trainer)
      env.bestnn = copy(env.curnn)
      nn_replaced = true
    else
      Handlers.checkpoint_started(handler)
      env.curnn = get_trained_network(trainer)
      eval_report =
        compare_networks(env.gspec, env.curnn, env.bestnn, ap, handler)
      teval += eval_report.time
      # If eval is good enough, replace network
      success = (eval_report.avgr >= best_evalr)
      if success
        nn_replaced = true
        env.bestnn = copy(env.curnn)
        best_evalr = eval_report.avgr
      end
      checkpoint_report =
        Report.Checkpoint(k * nbatches, eval_report, status, success)
      push!(checkpoints, checkpoint_report)
      Handlers.checkpoint_finished(handler, checkpoint_report)
    end
  end
  report = Report.Learning(
    tconvert, tloss, ttrain, teval,
    init_status, losses, checkpoints, nn_replaced)
  Handlers.learning_finished(handler, report)
  return report
end

function simple_memory_stats(env::Env)
  mem = get_experience(env)
  nsamples = length(mem)
  ndistinct = length(merge_by_state(mem))
  return nsamples, ndistinct
end

# To be given as an argument to `Simulator`
function self_play_measurements(trace, _, player)
  mem = MCTS.approximate_memory_footprint(player.mcts)
  edepth = MCTS.average_exploration_depth(player.mcts)
  return (trace=trace, mem=mem, edepth=edepth)
end

function self_play_step!(env::Env, handler)
  params = env.params.self_play
  Handlers.self_play_started(handler)
  make_oracle() =
    Network.copy(env.bestnn, on_gpu=params.sim.use_gpu, test_mode=true)
  simulator = Simulator(make_oracle, self_play_measurements) do oracle
    return MctsPlayer(env.gspec, oracle, params.mcts)
  end
  # Run the simulations
  results, elapsed = @timed simulate_distributed(
    simulator, env.gspec, params.sim,
    game_simulated=()->Handlers.game_played(handler))
  # Add the collected samples in memory
  new_batch!(env.memory)
  for x in results
    push_trace!(env.memory, x.trace, params.mcts.gamma)
  end
  speed = cur_batch_size(env.memory) / elapsed
  edepth = mean([x.edepth for x in results])
  mem_footprint = maximum([x.mem for x in results])
  memsize, memdistinct = simple_memory_stats(env)
  report = Report.SelfPlay(
    speed, edepth, mem_footprint, memsize, memdistinct)
  Handlers.self_play_finished(handler, report)
  return report
end

function memory_report(env::Env, handler)
  if isnothing(env.params.memory_analysis)
    return nothing
  else
    report = memory_report(
      env.memory, env.curnn, env.params.learning, env.params.memory_analysis)
    Handlers.memory_analyzed(handler, report)
    return report
  end
end

"""
    train!(env::Env, handler=nothing)

Start or resume the training of an AlphaZero agent.

A `handler` object can be passed that implements a subset of the callback
functions defined in [`Handlers`](@ref).
"""
function train!(env::Env, handler=nothing)
  while env.itc < env.params.num_iters
    Handlers.iteration_started(handler)
    resize_memory!(env, env.params.mem_buffer_size[env.itc])
    sprep, spperfs = Report.@timed self_play_step!(env, handler)
    mrep, mperfs = Report.@timed memory_report(env, handler)
    lrep, lperfs = Report.@timed learning_step!(env, handler)
    rep = Report.Iteration(spperfs, mperfs, lperfs, sprep, mrep, lrep)
    env.itc += 1
    Handlers.iteration_finished(handler, rep)
  end
  Handlers.training_finished(handler)
end

#####
##### AlphaZero player
#####

function guess_mcts_arena_params(env::Env)
  p = env.params
  return isnothing(p.arena) ? p.self_play.mcts : p.arena.mcts
end

function guess_use_gpu(env::Env)
  p = env.params
  return isnothing(p.arena) ? p.self_play.sim.use_gpu : p.arena.sim.use_gpu
end

"""
    AlphaZeroPlayer(::Env; [timeout, mcts_params, use_gpu])

Create an AlphaZero player from the current [training environment](@ref environment).

Note that the returned player may be slow as it does not batch MCTS requests.
"""
function AlphaZeroPlayer(env::Env; timeout=2.0, mcts_params=nothing, use_gpu=nothing)
  isnothing(mcts_params) && (mcts_params = guess_mcts_arena_params(env))
  isnothing(use_gpu) && (use_gpu = guess_use_gpu(env))
  net = Network.copy(env.bestnn, on_gpu=use_gpu, test_mode=true)
  return MctsPlayer(env.gspec, net, mcts_params, timeout=timeout)
end


text:
using Islander
using WCS
using Test

@testset "Islander.jl" begin
    # write some test for what will go in Islander
end

@testset "Pixell.jl" begin

    a = Enmap(3.0 .* ones(2,3), WCSTransform(2))
    b = Enmap(2.0 .* ones(2,3), WCSTransform(2))

    # broadcasting
    @test a .+ b ≈ (5.0 .* ones(2,3))
    @test a .* b ≈ (6.0 .* ones(2,3))
    @test a .- b ≈ (1.0 .* ones(2,3))
    @test a ./ b ≈ (1.5 .* ones(2,3))

    # get / set
    c = Enmap(3.0 .* ones(2,3), WCSTransform(2))
    c[1,1] = 1.0
    @test c[1,1] ≈ 1.0
    @test c[2,1] ≈ 3.0

end


text:
export Grid3D

"""
    Grid3D(;xlims=(-1,1), ylim=(-1, 1), zlims=(-1, 1), res=(21, 21, 21))

Create a 3D grid.
"""
function Grid3D(;xlims=(-1, 1), ylims=(-1, 1), zlims=(-1, 1), res=(21, 21, 21))
    Grid3D(
        central_range(xlims..., res[1]),
        central_range(ylims..., res[2]),
        central_range(zlims..., res[3]))
end

function Grid3D(lims, res::Int)
    Grid3D(
        central_range(lims..., res),
        central_range(lims..., res),
        central_range(lims..., res))
end

Base.broadcastable(G::Grid3D) = Ref(G)
Base.size(G::Grid3D) = (length(G.xrange), length(G.yrange), length(G.zrange))
Base.size(G::Grid3D, i) = size(G)[i]
Base.length(G::Grid3D) = prod(size(G))
Base.ndims(G::Grid3D) = 3

@inline function Base.getindex(G::Grid3D, s::Integer)
    i, j, k = _ind2sub(size(G), s)
    G[i, j, k]
end
Base.getindex(G::Grid3D, I::CartesianIndex{3}) = G[Tuple(I)...]

@inline function Base.getindex(G::Grid3D, i::Integer, j::Integer, k::Integer)
    G.xrange[i], G.yrange[j], G.zrange[k]
end
Base.getindex(G::Grid3D, ijk::NTuple{3, Int}) = G[ijk[1], ijk[2], ijk[3]]
Base.getindex(G::Grid3D, xyz::NTuple{3, Float64}) = gridpoint(G, xyz[1], xyz[2], xyz[3])

xcoordinates(G::Grid3D) = G.xrange
ycoordinates(G::Grid3D) = G.yrange
zcoordinates(G::Grid3D) = G.zrange

xcoordinate(G::Grid3D, i) = G.xrange[i]
ycoordinate(G::Grid3D, i) = G.yrange[i]
zcoordinate(G::Grid3D, i) = G.urange[i]

xgridpoint(G::Grid3D, x::Float64) = rangeindex(G.xrange, x)
ygridpoint(G::Grid3D, y::Float64) = rangeindex(G.yrange, y)
zgridpoint(G::Grid3D, z::Float64) = rangeindex(G.zrange, z)
gridpoint(G::Grid3D, x, y, z) = (xgridpoint(G, x), ygridpoint(G, y), zgridpoint(G, z))


text:
# protein type that keeps centered MD data for a protein
# 2014 - murozturk@ku.edu.tr

immutable protein
	R::Matrix{Float64}
	index::Vector{Int}
	means::Vector{Float64}
	stdPerRes::Vector{Float64}
	covariance::Matrix{Float64}

	function protein(R::Matrix{Float64}, index=1:size(R,2))
		m = mean(R,1)
		R .-= m

		println ("Calculating covariance matrix...")
		gamma = cov(R, mean=0, corrected=false)

		println ("Calculating real fluctuations...")
		rF = realFluc(R, centered=true)

		new(R, index[:], m[:], std(rF,1)[:], gamma  )
	end
end

# Returns decentered coordinates
# recovers original md data

deCenter(p::protein) = p->p.R .+ p.means'

# Returns real fluctuations per residue

function realFluc(R::Matrix{Float64}; centered=false)
	!centered && (R .-= mean(R,1))
	sqrt(R[:,1:3:end].^2 + R[:,2:3:end].^2 + R[:,3:3:end].^2)
end
realFluc(p::protein) = realFluc(p.R, centered=true)


text:
include("get_observation.jl")
include("get_advisory.jl")

"""
	Update the dynamics of an Aircraft object when there is alert
	ac.advisory != :COC
"""
function update_dynamics_with_alert(ac::Aircraft)
	# get turn rate
	ac.dyn.turn_rate = deg2rad(ac.cas.Advisory_to_Action_dict[ac.advisory])
	# update heading
	ac.dyn.heading += ac.dyn.turn_rate * ac.dyn.dt
	# constrain heading
	if ac.dyn.heading > pi
		ac.dyn.heading -= 2 * pi
	elseif ac.dyn.heading < - pi
		ac.dyn.heading += 2 * pi
	end
	# update coordinates
	ac.dyn.x += cos(ac.dyn.heading) * ac.dyn.v * ac.dyn.dt
	ac.dyn.y += sin(ac.dyn.heading) * ac.dyn.v * ac.dyn.dt
	ac.actual_route_len += ac.dyn.v * ac.dyn.dt
end

"""
	Update the dynamics of an Aircraft object when there is NO alert
	ac.advisory == :COC || isempty(ac.sensor.obs)
"""
function go_to_dest(ac::Aircraft)
	# set threshold for confirming the correct heading towards destination
	threshold = deg2rad(10.)

	heading_vec = [cos(ac.dyn.heading); sin(ac.dyn.heading); 0.]
	desti_dir_angle = atan2(ac.y_dest - ac.dyn.y, ac.x_dest - ac.dyn.x)
	heading_error = abs(desti_dir_angle - ac.dyn.heading)

	if heading_error < threshold
		ac.dyn.turn_rate = 0.
		ac.dyn.heading = desti_dir_angle
	else
		desti_vec = [cos(desti_dir_angle); sin(desti_dir_angle); 0.]
		cross_result = cross(heading_vec, desti_vec)

		if cross_result[end] < 0
			side = :Right
		elseif cross_result[end] > 0
			side = :Left
		else
			side = :OnDir
		end

		# get turn rate
		if side == :Right
			ac.dyn.turn_rate = deg2rad(ac.cas.Advisory_to_Action_dict[:SR])
		elseif side == :Left
			ac.dyn.turn_rate = deg2rad(ac.cas.Advisory_to_Action_dict[:SL])
		else # dest on heading
			ac.dyn.turn_rate = 0. # keep heading
		end
	end
	# update heading
	ac.dyn.heading = ac.dyn.heading + ac.dyn.turn_rate * ac.dyn.dt
	# constrain heading
	if ac.dyn.heading > pi
		ac.dyn.heading = ac.dyn.heading - 2 * pi
	elseif ac.dyn.heading < - pi
		ac.dyn.heading = ac.dyn.heading + 2 * pi
	end
	# update coordinates
	ac.dyn.x = ac.dyn.x + cos(ac.dyn.heading) * ac.dyn.v * ac.dyn.dt
	ac.dyn.y = ac.dyn.y + sin(ac.dyn.heading) * ac.dyn.v * ac.dyn.dt
	ac.actual_route_len += ac.dyn.v * ac.dyn.dt
end

"""
	Warp up
	update dynamics for all the agents
"""
function update_all_dynamics(All_AC::Array{Aircraft, 1}; correction=false)
	for ac in All_AC
		if correction == false
			# update observation (state)
			get_observation!(ac, All_AC, ac.sensor)
			# update action (advisory)
			get_advisory(ac, ac.cas, ac.sensor)
		else
			full_state, lowfi_state = get_full_and_low_fi_state(correction_type, ac, All_AC, policy.env.num_sections)
			normalize_obs!(ac, policy.env.num_sections, full_state)
			get_advisory(ac, full_state, lowfi_state)
		end
	end

	for ac in All_AC
		if isempty(ac.sensor.obs) == true || ac.advisory == :COC
			ac.advisory = :COC
			go_to_dest(ac)
		else
			update_dynamics_with_alert(ac)
		end
	end
end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
using LabView0mqJl

try
    include(get_script_path(raw"/users/eben60/Desktop/Julia/LabView0mqJl.jl/src/LabVIEW/jl-scripts/09_combined_example.jl"))
    LabView0mqJl.setglobals(isOK=true)
catch excep
    LabView0mqJl.setglobals(isOK=false, excpn=(;excep, stack_trace=backtrace()))
end

server_0mq4lv(fns)


text:
const func_to_nl = Dict(
:+ => 0,
:- => 1,
:* => 2,
:/ => 3,
# :rem => 4,
:^ => 5,
# :less => 6,
:min => 11,
:max => 12,
# :floor => 13,
# :ceil => 14,
:abs => 15,
:neg => 16,
:|| => 20,
:&& => 21,
:< => 22,
:<= => 23,
:(==) => 24,
:>= => 28,
:> => 29,
:!= => 30,
:ifelse => 35,
# :not => 34,
:tanh => 37,
:tan => 38,
:sqrt => 39,
:sinh => 40,
:sin => 41,
:log10 => 42,
:log => 43,
:exp => 44,
:cosh => 45,
:cos => 46,
:atanh => 47,
# :atan2 => 48,
:atan => 49,
:asinh => 50,
:asin => 51,
:acosh => 52,
:acos => 53,
:sum => 54,
# :intdiv => 55,
# :precision => 56,
# :round => 57,
# :trunc => 58,
# :count => 59,
# :numberof => 60,
# :numberofs => 61,
# :ifs => 65,
# :and_n => 70,
# :or_n => 71,
# :implies => 72,
# :iff => 73,
# :alldiff => 74,
)

const sense_to_nl = Dict(
:Min => 0,
:Max => 1,
)

const relation_to_nl = Dict(
:multiple => 0,
:<=       => 1,
:>=       => 2,
:(==)     => 4,
)

const reverse_relation = Dict(
:0 => 0,
:1 => 2,
:2 => 1,
:4 => 4,
)

const nary_functions = Set{Symbol}([
:min,
:max,
:sum,
# :count,
# :numberof,
# :numberofs,
# :and_n,
# :or_n,
# :alldiff,
])