text:
### F0 estimation ###

const supported_otypes = ["pitch", "f0", "logf0"]

function swipe(x::StridedVector{Cdouble}, fs, hopsize=80;
               min::Float64=60.0,
               max::Float64=240.0,
               threshold::Float64=0.3,
               otype::Union{Int,AbstractString}="pitch")
    otype_int = (-1)::Int
    if isa(otype, AbstractString)
        otype_int = indexin([otype], supported_otypes)[1] - 1
    else
        otype_int = otype
    end

    if otype_int ∉ 0:2
        throw(ArgumentError("""
unsupported otype: $otype, must be (0) pitch, (1) f0, or (2) log(f0)"""))
    end

    if min >= max || max >= fs/2
        throw(ArgumentError("invalid min/max frequency parameters"))
    end

    expectedlen = ceil(Int, length(x) / hopsize)
    f0 = Array{Cdouble}(expectedlen)
    ccall((:swipe, libSPTK), Void,
          (Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint, Cint, Cdouble, Cdouble,
           Cdouble, Cint),
          x, f0, length(x), fs, hopsize, min, max, threshold, otype_int)
    f0
end

function rapt(x::StridedVector{Cfloat}, fs, hopsize=80;
              min::Float64=60.0,
              max::Float64=240.0,
              voice_bias::Float64=0.0,
              otype::Union{Int,AbstractString}="pitch")
    otype_int = (-1)::Int
    if isa(otype, AbstractString)
        otype_int = indexin([otype], supported_otypes)[1] - 1
    else
        otype_int = otype
    end

    if otype_int ∉ 0:2
        throw(ArgumentError("""
unsupported otype: $otype, must be (0) pitch, (1) f0, or (2) log(f0)"""))
    end

    if min >= max || max >= fs/2 || min <= fs/10000.0
        throw(ArgumentError("invalid min/max frequency parameters"))
    end

    frame_period = Float64(hopsize / fs)
    frame_period = trunc(0.5 + fs * frame_period) / fs
    if frame_period > 0.1 || frame_period < 1/fs
        throw(ArgumentError("frame period must be between [1/fs, 0.1]"))
    end

    expectedlen = ceil(Int, length(x) / hopsize)
    f0 = Array{Cfloat}(expectedlen)
    ret = ccall((:rapt, libSPTK), Int,
                (Ptr{Cfloat}, Ptr{Cfloat}, Cint, Cdouble, Cint, Cdouble, Cdouble,
                 Cdouble, Cint),
                x, f0, length(x), fs, hopsize, min, max, voice_bias, otype_int)

    if ret == 2
        error("input range too small for analysis by get_f0")
    elseif ret == 3
        error("problem in init_dp_f0()")
    end

    @assert ret == 0

    f0
end


text:
struct Multiplicative{f} <: LocalOperator
    field::f
end

kernelvals(biop::Multiplicative, x) = Nothing
integrand(op::Multiplicative, kernel, x, g, f) = dot(f[1], op.field(cartesian(x))*g[1])
scalartype(op::Multiplicative) = Float64

ext: jl


text:
using Feather, Test, CategoricalArrays
using DataFrames, Arrow, Tables
using Random, Dates

const ≅ = isequal

# whether or not to do python round-trip tests
# currently does not work in travis
const DO_PYTHON_ROUNDTRIP = false
const PYTHON_EXEC = "python3"

const testdir = joinpath(dirname(pathof(Feather)), "../test/data")
files = map(x -> joinpath(testdir, x), readdir(testdir))

temps = []

@testset "ReadWrite" for f in files
    @info("testing $f...")
    source = Feather.Source(f)
    df = DataFrame(source)
    temp = tempname()
    push!(temps, temp)
    sink = Feather.write(temp, df)
    df2 = Feather.read(temp)

    for (c1,c2) ∈ zip(getfield(df, :columns), getfield(df, :columns))
        for i = 1:length(c1)
            @test isequal(c1[i], c2[i])
        end
    end

    source2 = Feather.Source(sink)
    @test source.ctable.description == source2.ctable.description
    @test source.ctable.num_rows == source2.ctable.num_rows
    @test source.ctable.metadata == source2.ctable.metadata
    for (col1,col2) in zip(source.ctable.columns,source2.ctable.columns)
        @test col1.name == col2.name
        @test col1.metadata_type == col2.metadata_type
        @test typeof(col1.metadata) == typeof(col2.metadata)
        @test col1.user_metadata == col2.user_metadata

        v1 = col1.values; v2 = col2.values
        @test v1.dtype == v2.dtype
        @test v1.encoding == v2.encoding
        # @test v1.offset == v2.offset # currently not python/R compatible due to wesm/feather#182
        @test v1.length == v2.length
        @test v1.null_count == v2.null_count
        # @test v1.total_bytes == v2.total_bytes
    end
end

GC.gc(); GC.gc()
for t in temps
    rm(t)
end

const Nrows = 10

@testset "Buffer" begin
    io = IOBuffer()
    df = DataFrame(A=rand(Nrows), B=map(randstring, 1:Nrows))
    Feather.write(io, df)
    seekstart(io)
    df2 = Feather.read(io)
    for col ∈ names(df)
        @test df[col] == df2[col]
    end
end

include("arrowtests.jl")

GC.gc(); GC.gc()

@testset "issue#34" begin
    data = DataFrame(A=Union{Missing, String}[randstring(10) for i ∈ 1:100], B=rand(100))
    data[2, :A] = missing
    Feather.write("testfile.feather", data)
    dfo = Feather.read("testfile.feather")
    @test size(dfo) == (100, 2)
    GC.gc();
    rm("testfile.feather")
end

@testset "issue#124" begin
    df = DataFrame(A=rand(5), B=missings(5))
    @test_throws ArgumentError Feather.write("test124.feather", df)
    isfile("test124.feather") && rm("test124.feather")
end

if DO_PYTHON_ROUNDTRIP

@testset "PythonRoundtrip" begin
    try
        @info("Generating test.feather file from python...")
        run(`$PYTHON_EXEC runtests.py`)

        @info("Reading test.feather into Julia...")
        df = Feather.read("test.feather")

        dts = [Dates.DateTime(2016,1,1), Dates.DateTime(2016,1,2), Dates.DateTime(2016,1,3)]

        @test df[:Autf8][:] == ["hey","there","sailor"]
        @test df[:Abool][:] == [true, true, false]
        @test df[:Acat][:] == categorical(["a","b","c"])  # these violate Arrow standard by using Int8!!
        @test df[:Acatordered][:] == categorical(["d","e","f"])  # these violate Arrow standard by using Int8!!
        @test convert(Vector{Dates.DateTime}, df[:Adatetime][:]) == dts
        @test isequal(df[:Afloat32][:], [1.0, missing, 0.0])
        @test df[:Afloat64][:] == [Inf,1.0,0.0]

        df_ = Feather.read("test.feather"; use_mmap=false)

        @info("Writing test2.feather from Julia...")
        Feather.write("test2.feather", df)
        df2 = Feather.read("test2.feather")

        @test df2[:Autf8][:] == ["hey","there","sailor"]
        @test df2[:Abool][:] == [true, true, false]
        @test df2[:Acat][:] == categorical(["a","b","c"])  # these violate Arrow standard by using Int8!!
        @test df2[:Acatordered][:] == categorical(["d","e","f"])  # these violate Arrow standard by using Int8!!
        @test convert(Vector{Dates.DateTime}, df2[:Adatetime][:]) == dts
        @test isequal(df2[:Afloat32][:], [1.0, missing, 0.0])
        @test df2[:Afloat64][:] == [Inf,1.0,0.0]

        @info("Read test2.feather into Python...")
        @test (run(`$PYTHON_EXEC runtests2.py`); true)
    finally
        isfile("test.feather") && rm("test.feather")
        isfile("test2.feather") && rm("test2.feather")
    end
end

end


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
using GenerativeAD.Models
using BSON
using StatsBase: fit!, predict, sample

using Flux
using MLDataPattern

s = ArgParseSettings()
@add_arg_table! s begin
   "max_seed"
		required = true
		arg_type = Int
		help = "seed"
	"dataset"
		required = true
		arg_type = String
		help = "dataset"
    "contamination"
    	arg_type = Float64
    	help = "contamination rate of training data"
    	default = 0.0
end
parsed_args = parse_args(ARGS, s)
@unpack dataset, max_seed, contamination = parsed_args

modelname ="adVAE"

function sample_params()
	argnames = (
		:hdim,
		:zdim,
		:nlayers,
		:activation,
		:gamma,
		:lambda,
		:mx,
		:mz,
		:lr,
		:decay,
		:batch_size,
		:iters,
		:check_every,
		:patience,
		:init_seed,
	)
	par_vec = (
		2 .^(4:9),
		2 .^(1:8),
		3:4,
		["relu", "swish", "tanh"],
		[0.0005, 0.001, 0.005], # γ
		[0.005, 0.01, 0.05], # λ  # better params
		[1, 1.5],  #mx #0.5:0.5:2.5
 		[40, 50, 60], #mz 10:10:100
		10f0 .^ (-4:-3),
		0f0:0.1:0.5,
		2 .^ (5:6),
		[10000],
		[30],
		[10],
		1:Int(1e8),
	)

	return NamedTuple{argnames}(map(x->sample(x,1)[1], par_vec))
end

function fit(data, parameters)
	# define models (Generator, Discriminator)
	advae = GenerativeAD.Models.adVAE(;parameters...)

	# define optimiser
	try
		global info, fit_t, _, _, _ = @timed fit!(advae |> gpu, data, parameters)
	catch e
		println("Error caught.")
		return (fit_t = NaN, model = nothing, history = nothing, n_parameters = NaN), []
	end

	training_info = (
		fit_t = fit_t,
		model = info[2]|>cpu,
		history = info[1], # losses through time
		npars = info[3], # number of parameters
		iters = info[4] # optim iterations of model
		)

	return training_info,
	[(x -> GenerativeAD.Models.anomaly_score(advae|>cpu, x; dims=1, L=100)[:], merge(parameters, (L=100, )))]
	# L = samples for one x in anomaly_score computation
end

#_________________________________________________________________________________________________

try_counter = 0
max_tries = 10*max_seed
cont_string = (contamination == 0.0) ? "" : "_contamination-$contamination"
while try_counter < max_tries
    parameters = sample_params()

    for seed in 1:max_seed
		savepath = datadir("experiments/tabular$cont_string/$(modelname)/$(dataset)/seed=$(seed)")
		mkpath(savepath)

		# get data
		data = GenerativeAD.load_data(dataset, seed=seed, contamination=contamination)

		# update parameter
		parameters = merge(parameters, (idim=size(data[1][1],1), ))
		# here, check if a model with the same parameters was already tested
		if GenerativeAD.check_params(savepath, parameters)
			#(X_train,_), (X_val, y_val), (X_test, y_test) = data
			training_info, results = fit(data, parameters)
			# saving model separately
			if training_info.model !== nothing
				tagsave(joinpath(savepath, savename("model", parameters, "bson", digits=5)), Dict("model"=>training_info.model), safe = true)
				training_info = merge(training_info, (model = nothing,))
			end
			save_entries = merge(training_info, (modelname = modelname, seed = seed, dataset = dataset, contamination = contamination))
			# now loop over all anomaly score funs
			for result in results
				GenerativeAD.experiment(result..., data, savepath; save_entries...)
			end
			global try_counter = max_tries + 1
		else
			@info "Model already present, sampling new hyperparameters..."
			global try_counter += 1
		end
	end
end
(try_counter == max_tries) ? (@info "Reached $(max_tries) tries, giving up.") : nothing


text:
export Parameter, Parameters
export write_parameterfile
export write_inputfiles

@doc """
`ALPS.Parameter`, which is a typealias of `Dict{AbstractString, Any}`,
represents a set of mapping from parameter name to parameter value.
One `Parameter` corresponds to one ALPS task.
""" typealias Parameter Dict{AbstractString, Any}

@doc """
`ALPS.Parameters` is a typealias of `Vector{ALPS.Parameter}`.
One `Parameters` corresponds to one (whole) ALPS simulation.
""" typealias Parameters Vector{Parameter}

Parameters() = Parameter[]

@doc """
Generates a plain text parameter file from `parameter`,
and returns parameter filename.
""" ->
function write_parameterfile(filename :: AbstractString, parameter::Parameter)
  open(filename, "w")do io
    for (k,v) in parameter
      println(io, k, " = ", v)
    end
  end
  return filename
end

@doc """
Generates XML input files for ALPS simulation from `parameters`
and returns master XML file, `$(filename_prefix).in.xml`.
""" ->
function write_inputfiles(filename_prefix :: AbstractString, parameters::Parameters; baseseed :: Union(Unsigned, ()) = () )

  masterfilename = "$(filename_prefix).in.xml"

  if baseseed == ()
    baseseed = generate_seed()
  end

  bits = 31 - ndigits(length(parameters),2)

  io = open(masterfilename, "w")
  println(io, """<?xml version="1.0" encoding="UTF-8"?>""")
  println(io, """<?xml-stylesheet type="text/xsl" href="ALPS.xsl"?>""")
  print(io, """<JOB """)
  print(io, """xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" """)
  println(io, """xsi:noNamespaceSchemaLocation="http://xml.comp-phys.org/2003/8/job.xsd">""")
  println(io, """  <OUTPUT file="$filename_prefix.out.xml"/>""")

  for (i, p) in enumerate(parameters)
    println(io, """  <TASK status="new">""")
    println(io, """    <INPUT file="$filename_prefix.task$i.in.xml"/>""")
    println(io, """    <OUTPUT file="$filename_prefix.task$i.out.xml"/>""")
    println(io, """  </TASK>""")

    if !haskey(p, "SEED")
      seed = baseseed
      for j in 0:(div(32,bits)+1)
        seed $= ((i-1) << (j*bits))
        seed &= ((1<<30) | ((1<<30)-1))
      end
      p["SEED"] = seed
    end
    write_taskfile("$(filename_prefix).task$(i).in.xml",p)
  end
  println(io, """</JOB>""")
  close(io)
  copy_stylesheet()
  return masterfilename
end

function write_taskfile(filename::AbstractString, parameter::Parameter)
  io = open(filename, "w")
  println(io, """<?xml version="1.0" encoding="UTF-8"?>""")
  println(io, """<?xml-stylesheet type="text/xsl" href="ALPS.xsl"?>""")
  print(io, """<SIMULATION """)
  print(io, """xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" """)
  println(io, """xsi:noNamespaceSchemaLocation="http://xml.comp-phys.org/2003/8/job.xsd">""")
  println(io, """  <PARAMETERS>""")
  for (k,v) in parameter
    print(io, """    <PARAMETER name="$k">""")
    print(io, v)
    println(io, """</PARAMETER>""")
  end
  println(io, """  </PARAMETERS>""")
  println(io, """</SIMULATION>""")
  close(io)
end

function generate_seed()
  seed = uint(time()*1.0e9)
  (seed << 10) | (seed >> 22)
end

function copy_stylesheet(dst :: AbstractString = "./ALPS.xsl")
  if ispath(dst)
    return nothing
  end
  src = joinpath(dirname(@__FILE__()), "ALPS.xsl")
  cp(src, dst)
  return nothing
end


text:
##############
#    CSV     #
##############

mutable struct CSV
  file::IOStream
  file_timings::IOStream
  delimiter::Char
end

function CSV(file_name::String, delimiter::Char = ';')
  labels = reshape(["algorithm", "Julia"], (1, :))

  f = open(string(file_name, ".txt"), "w")
  writedlm(f, labels, delimiter)

  ft = open(string(file_name, "_timings.txt"), "w")
  return CSV(f, ft, delimiter)
end

function add_data(c::CSV, name::AbstractString, results::Results) where P
  data = [reshape([name], (1, :)) [results.min_time]]
  writedlm(c.file, data, c.delimiter)

  data = [reshape([name], (1, :)) transpose(results.timings)]
  writedlm(c.file_timings, data, c.delimiter)
end

function close(c::CSV)
  close(c.file)
  close(c.file_timings)
end


text:
export getFGMRESmem,FGMRES,FGMRES_relaxation

type FGMRESmem
	v_prec              ::ArrayTypes ## memory for the result of the preconditioner
	Az					::ArrayTypes ## memory for the result of A*z.
	V    				::ArrayTypes
	Z					::ArrayTypes
end

function resetMem(mem::FGMRESmem)
mem.v_prec[:]=0.0;
mem.Az[:]=0.0;
mem.V[:]=0.0;
mem.Z[:]=0.0;
end

function getFGMRESmem(n::Int64,flexible::Bool,T::Type,k::Int64,m::Int64=1)
v_prec = zeros(T,0);
Az = zeros(T,0);
if m==1
	Az = zeros(T,n);
	v_prec = zeros(T,n);
else
	Az = zeros(T,n,m);
	v_prec = zeros(T,n,m);
end

if flexible
	return FGMRESmem(v_prec,Az,zeros(T,n,k*m),zeros(T,n,k*m));
else
	return FGMRESmem(v_prec,Az,zeros(T,n,k*m),zeros(T,0));
end
end

function getEmptyFGMRESmem()
return FGMRESmem(zeros(0),zeros(0),zeros(0),zeros(0));
end
import Base.isempty
function isempty(mem::FGMRESmem)
return size(mem.V,1)==0;
end

function FGMRES(AT::SparseMatrixCSC,r0::Vector,x0::Vector,inner::Int64,prec::Function,TOL::Float64,
                 verbose::Bool,flexible::Bool,numCores::Int64, mem::FGMRESmem =  getEmptyFGMRESmem())
Az = zeros(eltype(r0),size(r0));
function Afun(z::Vector)
	SpMatMul(AT,z,Az,numCores);
	return Az;
end
return FGMRES(Afun,r0,x0,inner,prec,TOL,verbose,flexible,numCores,mem);
end

function FGMRES_relaxation(AT::SparseMatrixCSC,r0::Vector,x0::Vector,inner::Int64,prec::Function,TOL::Float64,
                 verbose::Bool,flexible::Bool,numCores::Int64, mem::FGMRESmem =  getEmptyFGMRESmem())
Az = zeros(eltype(r0),size(r0));
function Afun(z::Vector)
	SpMatMul(AT,z,Az,numCores);
	return Az;
end
return FGMRES_relaxation(Afun,r0,x0,inner,prec,TOL,verbose,flexible,numCores,mem);
end

function FGMRES(Afun::Function,r0::Vector,x0::Vector,inner::Int64,prec::Function,TOL::Float64,
                 verbose::Bool,flexible::Bool,numCores::Int64, mem::FGMRESmem = getEmptyFGMRESmem())
n = length(r0);
if isempty(mem)
	# warn("Allocating memory in GMRES")
	mem = getFGMRESmem(n,flexible,eltype(r0),inner);
else
    resetMem(mem);
	if size(mem.V,2)!= inner
		error("FGMRES: size of Krylov subspace is different than inner");
	end
end

betta = norm(r0);
rnorm0 = betta;

TYPE = eltype(r0);

H = zeros(TYPE,inner+1,inner);
xi = zeros(TYPE,inner+1);
t = zeros(TYPE,inner);
y = copy(t);
xi[1] = betta;

constOne = one(TYPE);
constZero = zero(TYPE);

Z = mem.Z;
V = mem.V;

BLAS.scal!(n,(1/betta)*constOne,r0,1); # w = w./betta

rnorms = zeros(inner);
w = 0;
for j = 1:inner
	if j==1
		V[:,j] = r0;  # no memory problem with this line....
		z = prec(r0)
	else
		V[:,j] = w;  # no memory problem with this line....
		z = prec(w)
	end

	if flexible
		Z[:,j] = z; # no memory problem with this line....
	end
	w = Afun(z);
    # Gram Schmidt:
	BLAS.gemv!('C', constOne, V, w,constZero,t);# t = V'*w;
	H[1:inner,j] = t;
	BLAS.gemv!('N', -constOne, V, t,constOne,w); # w = w - V*t

	## modified Gram Schmidt:
    # for i=1:j
        # H[i,j] = dot(vec(V[:,i]),w);
        # w = w - H[i,j]*vec(V[:,i]);
    # end

    betta = norm(w);
    H[j+1,j] = betta;

	BLAS.scal!(n,(1/betta)*constOne,w,1); # w = w*(1/betta);

	# the following 2 lines are equivalent to the 2 next
    # y = H[1:j+1,1:j]\xi[1:j+1];
    # y = norm(H[1:j+1,1:j]*y - xi[1:j+1])
    Q = qr(H[1:j+1,1:j],thin=false)[1];
	rnorms[j] = abs(Q[1,end]*xi[1]);
    if verbose
		if j==1
			resprev = rnorm0;
		else
			resprev = rnorms[j-1];
		end
        # println(string("Inner iter ",j,": FGMRES: residual norm: ",rnorms[j],", relres: ",rnorms[j]/rnorm0,", factor: ",rnorms[j]/resprev));
	end
    if rnorms[j] < TOL
        rnorms = rnorms[1:j];
        break;
    end
end
y[:] = pinv(H)*xi;
# norm(H*y - xi)
if flexible
	BLAS.gemv!('N', constOne, Z, y,constZero,w); # w = Z*y  #This is the correction that corresponds to the residual.
else
	BLAS.gemv!('N', constOne, V, y, constZero, w);
	z = prec(w);
	w = z;
end
addVectors(constOne,w,x0);

return x0,rnorms; # w is the correction vector;
end

function FGMRES_relaxation(Afun::Function,r0::Vector,x0::Vector,inner::Int64,prec::Function,TOL::Float64,
                 verbose::Bool,flexible::Bool,numCores::Int64, mem::FGMRESmem = getEmptyFGMRESmem())

n = length(r0);
if isempty(mem)
	warn("Allocating memory in GMRES")
	mem = getFGMRESmem(n,true,eltype(r0),inner);
else
    resetMem(mem);
	if size(mem.V,2)!= inner
		error("FGMRES_relaxation: size of Krylov subspace is different than inner");
	end
end

rnorm0 = norm(r0);

TYPE = eltype(r0);

H = zeros(TYPE,inner,inner);
xi = zeros(TYPE,inner);
t = zeros(TYPE,inner);

constOne = one(TYPE);
constZero = zero(TYPE);

Z = mem.Z;
AZ = mem.V;

# w[:] = r0;
w = 0;
rnorms = zeros(inner);

for j = 1:inner
	if j==1
		z = prec(r0); # here unfortunately we have an allocation....
	else
		z = prec(w);
	end
	Z[:,j] = z; # no memory problem with this line....

	w = Afun(z); # w = A'*z;# w = SpMatMul(A,z,w,numCores);
	AZ[:,j] = w;  # no memory problem with this line....

	xi[j] = dot(w,r0);

	BLAS.gemv!('C', constOne, AZ, w,constZero,t);

	H[:,j] = t;
	H[j,:] = t';
	H = 0.5*H + 0.5*H';
	t[:] = pinv(H)*xi;

	rnorms[j] = sqrt(real(dot(t,H*t)-2.0*dot(t,xi)+rnorm0^2));

    if verbose
		if j==1
			resprev = rnorm0;
		else
			resprev = rnorms[j-1];
		end
        println(string("Inner iter ",j,": FGMRES: residual norm: ",rnorms[j],", relres: ",rnorms[j]/rnorm0,", factor: ",rnorms[j]/resprev));
    end
    if rnorms[j] < TOL
        rnorms = rnorms[1:j];
        break;
    end
end
BLAS.gemv!('N', constOne, Z, t,constZero,w); # w = Z*y  #This is the correction that corresponds to the residual.
addVectors(constOne,w,x0);

return x0,rnorms; # w is the correction vector;
end


text:
module TestAqua

using Recyclers
using Aqua

test() = Aqua.test_all(Recyclers)

end  # module


text:
using CartesianGP
using Base.Test

# Test function mu_lambda() from Evolution.jl
# Needs to be revised

function parameters_half_adder()
    numperlevel = 1
    numlevels = 2
    numlevelsback = numlevels
    numinputs = 2
    numoutputs = 2

    mu = 1
    lambda = 4
    mutrate = 0.15
    targetfitness = 1.0
    funcs = [AND, OR, XOR]
    fitfunc = hamming_max

    return Parameters(mu, lambda, mutrate, targetfitness, numinputs, numoutputs, numperlevel, numlevels, numlevelsback, funcs, fitfunc)
end

# This version has one node per level
function chromosome_half_adder(numperlevel=1, numlevels=2, numlevelsback=2)
    numinputs = 2
    numoutputs = 2

    #p = Parameters(numinputs, numoutputs, numperlevel, numlevels, numlevelsback)
    p = parameters_half_adder()
    c = Chromosome(p)

    c.inputs = [InputNode(1), InputNode(2)]
    c.interiors = Array(InteriorNode,2,1)
    c.interiors[1,1] = InteriorNode(XOR, [(0, 1), (0, 2)])
    c.interiors[2,1] = InteriorNode(AND, [(0, 1), (0, 2)])
    c.outputs = [OutputNode((1, 1)), OutputNode((2, 1))]

    return c
end

function goal_half_adder()
    g = Goal(2, (0b0110, 0b1000))
    #g = Goal(2, (0b1000, 0b0110))

    return g
end

p = parameters_half_adder()
goal = goal_half_adder()
maxgens = 100
srand(1)

println("results from test/Evolution.jl")
(ch,gens) = mu_lambda(p, goal, maxgens)
# The intention here is to compare ch with chromosome_half_adder(), but == does not work for comparison.
print_chromosome(ch)
ch_ha = chromosome_half_adder()
execute_chromosome(ch_ha)   # Sets active nodes to active
print_chromosome(ch_ha)
#@test ch == ch_ha


text:
"""
    srm0(t::Real, I, V; lastspike, eta)
    srm0!(t::Real, I::AbstractArray{<:Real}, V::AbstractArray{<:Real}; lastspike::AbstractArray{<:Real}, eta)
    srm0!(t::Real, I::CuVector{<:Real}, V::CuVector{<:Real}; lastspike::CuVector{<:Real}, eta)

Evaluate a SRM0 neuron.
Use `CuVector` instead of `Vector` to evaluate on GPU.

# Fields
- `t`: current time in seconds
- `I`: external current
- `V`: current membrane potential
- `lastspike`: time of last output spike in seconds
- `eta`: post-synaptic response function
"""
srm0(t::Real, I, V; lastspike, eta) = eta(t - lastspike) + I
function srm0!(t::Real, I::AbstractArray{<:Real}, V::AbstractArray{<:Real}; lastspike::AbstractArray{<:Real}, eta)
    V .= map.(eta, (t .- lastspike))
    V .+= I
end
srm0!(t::Real, I::CuVector{<:Real}, V::CuVector{<:Real}; lastspike::CuVector{<:Real}, eta) =
    V .= map.(eta, t .- lastspike) .+ I


text:
module GateRegistry

import GateOps
GateLabel = GateOps.GateLabel
Gate = GateOps.Gate

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #
#                        Map GateLabel to Gates
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #

export pauli_x, pauli_y, pauli_z, hadamard, u
export r_x, r_y, r_z, r_phase, swap
export c_pauli_x, c_pauli_y, c_pauli_z, c_u
export c_r_x, c_r_y, c_r_z, c_r_phase

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #

# For each unique GateLabel we require an assigned matrix value
GateMap = Dict{GateLabel, Gate}(
    GateLabel(:x)=>Gate([0 1; 1 0] .+ 0im),
    GateLabel(:y)=>Gate([0 -1im; 1im 0]),
    GateLabel(:z)=>Gate([1 0; 0 -1] .+ 0im),
    GateLabel(:h)=>Gate([1 1; 1 -1]*(1/sqrt(2)) .+ 0im)
)

# List available operators
# Generalise later to accept various types
Operators = Dict{String, Function}(
    "x"=>GateOps.pauli_x,
    "y"=>GateOps.pauli_y,
    "z"=>GateOps.pauli_z,
    "h"=>GateOps.hadamard,
    "u"=>GateOps.u,
    "r_x"=>GateOps.r_x,
    "r_y"=>GateOps.r_y,
    "r_z"=>GateOps.r_z,
    "r_phase"=>GateOps.r_phase,
    "c_x"=>GateOps.c_pauli_x,
    "c_y"=>GateOps.c_pauli_y,
    "c_z"=>GateOps.c_pauli_z,
    "c_u"=>GateOps.c_u,
    "c_r_x"=>GateOps.c_r_x,
    "c_r_y"=>GateOps.c_r_y,
    "c_r_z"=>GateOps.c_r_z,
    "c_r_phase"=>GateOps.c_r_phase
)

# List available modifications to gate
GateModifiers = Dict{String, Function}(
    "sqrt"=>GateOps.sqrt,
    "adj"=>GateOps.adjoint,
    "e2m"=>GateOps.euler_to_mat,
    "m2e"=>GateOps.mat_to_euler,
    "no_op"=>x->x
)

struct Registry
    map_label_gate::Dict{GateLabel, Gate}
    operators::Dict{String, Function}
    gate_mods::Dict{String, Function}
end
GateReg = Registry(GateMap, Operators, GateModifiers)

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #
#                               Get and cache
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #
@memoize function apply_gate(gate_func::Function, gate_modifier::Union{Nothing, Function}, params::Tuple)
    if gate_modifier == nothing
        return gate_func(params...)
    else
        gate_modifier(gate_func(params...)) # Construct new gate, symbol, label, and cache before returning
    end
end

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #
#                               Gate cache
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #

function create_reg_entity(label::String, mat::Matrix{<:Number})
    g = Gate(mat)
    l = GateLabel(label)
    return g=>l
end

function create_reg_operator(label::String, f::Function, params::Tuple)
    l = GateLabel(label)
    g = f(params...)
    add_gate!(l, g)
    push!(Operators, )
    return g=>l
end

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #
#                          Boilerplate functions
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #
function add_gate!(g_map::Dict{GateLabel, Gate}, label::GateLabel, gate::Gate)
    if haskey(g_map, label)
        return
    else
        push!(g_map, label=>gate)
    end
end
function add_gate!(label::GateLabel, gate::Gate)
    if haskey(GateMap, label)
        return
    else
        push!(GateMap, label=>gate)
    end
end
function add_gate!(label_gate::Pair{GateLabel, Gate})
    if haskey(GateMap, label_gate.first)
        return
    else
        push!(GateMap, label_gate)
    end
end

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #

function remove_gate!(g_map::Dict{GateLabel, Gate}, label::GateLabel)
    delete!(g_map, label)
end

function remove_gate!(label::GateLabel)
    delete!(GateMap, label)
end

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #

function get_gate(label::GateLabel)
    if ~haskey(GateMap, label)
        return GateMap[label]
    else
        println("Key $(label) does not exist")
    end
end

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #

function has_gate(label::GateLabel)
    return haskey(GateMap, label)
end

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #

function list_gates(;val=false)
    if ~val
        return keys(GateMap)
    else
        return values(GateMap)
    end
end

function list_gate_map(;val=false)
    return GateMap
end

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #

@memoize function Base.:sqrt(g::GateOps.GateLabelP)
    m = sqrt(GateMap[g])
    GateOps.GateLabelP
end

function mat_to_euler(gate::GateOps.GateLabel)
    return mat_to_euler(GateMap[gate])
end
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ApplicationGatewayBackendHealthServer <: SwaggerModel
    address::Any # spec type: Union{ Nothing, String } # spec name: address
    ipConfiguration::Any # spec type: Union{ Nothing, NetworkInterfaceIPConfiguration } # spec name: ipConfiguration
    health::Any # spec type: Union{ Nothing, String } # spec name: health
    healthProbeLog::Any # spec type: Union{ Nothing, String } # spec name: healthProbeLog

    function ApplicationGatewayBackendHealthServer(;address=nothing, ipConfiguration=nothing, health=nothing, healthProbeLog=nothing)
        o = new()
        validate_property(ApplicationGatewayBackendHealthServer, Symbol("address"), address)
        setfield!(o, Symbol("address"), address)
        validate_property(ApplicationGatewayBackendHealthServer, Symbol("ipConfiguration"), ipConfiguration)
        setfield!(o, Symbol("ipConfiguration"), ipConfiguration)
        validate_property(ApplicationGatewayBackendHealthServer, Symbol("health"), health)
        setfield!(o, Symbol("health"), health)
        validate_property(ApplicationGatewayBackendHealthServer, Symbol("healthProbeLog"), healthProbeLog)
        setfield!(o, Symbol("healthProbeLog"), healthProbeLog)
        o
    end
end # type ApplicationGatewayBackendHealthServer

const _property_map_ApplicationGatewayBackendHealthServer = Dict{Symbol,Symbol}(Symbol("address")=>Symbol("address"), Symbol("ipConfiguration")=>Symbol("ipConfiguration"), Symbol("health")=>Symbol("health"), Symbol("healthProbeLog")=>Symbol("healthProbeLog"))
const _property_types_ApplicationGatewayBackendHealthServer = Dict{Symbol,String}(Symbol("address")=>"String", Symbol("ipConfiguration")=>"NetworkInterfaceIPConfiguration", Symbol("health")=>"String", Symbol("healthProbeLog")=>"String")
Base.propertynames(::Type{ ApplicationGatewayBackendHealthServer }) = collect(keys(_property_map_ApplicationGatewayBackendHealthServer))
Swagger.property_type(::Type{ ApplicationGatewayBackendHealthServer }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayBackendHealthServer[name]))}
Swagger.field_name(::Type{ ApplicationGatewayBackendHealthServer }, property_name::Symbol) =  _property_map_ApplicationGatewayBackendHealthServer[property_name]

const _allowed_ApplicationGatewayBackendHealthServer_health = ["Unknown", "Up", "Down", "Partial", "Draining"]

function check_required(o::ApplicationGatewayBackendHealthServer)
    true
end

function validate_property(::Type{ ApplicationGatewayBackendHealthServer }, name::Symbol, val)
    if name === Symbol("health")
        Swagger.validate_param(name, "ApplicationGatewayBackendHealthServer", :enum, val, _allowed_ApplicationGatewayBackendHealthServer_health)
    end
end


text:
using SimpleDirectMediaLayer.LibSDL2
using Test

SDL2_pkg_dir = joinpath(@__DIR__, "..", "..")
audio_example_assets_dir = joinpath(SDL2_pkg_dir, "examples", "audio_example")

# check that an audio device if available
SDL_Init(SDL_INIT_AUDIO)
device = Mix_OpenAudio(22050, MIX_DEFAULT_FORMAT, 2, 1024)

if device == 0
    Mix_CloseAudio()
    SDL_Quit()

    @testset "Init+Quit" begin
        # Test that you can init and quit SDL_mixer multiple times correctly.
        @test SDL_Init(SDL_INIT_AUDIO) == 0
        @test Mix_OpenAudio(22050, MIX_DEFAULT_FORMAT, 2, 1024) == 0
        Mix_CloseAudio()
        SDL_Quit()

        @test SDL_Init(SDL_INIT_AUDIO) == 0
        @test Mix_OpenAudio(22050, MIX_DEFAULT_FORMAT, 2, 1024) == 0
        Mix_CloseAudio()
        SDL_Quit()
    end

    @testset "Sounds" begin
        SDL_Init(SDL_INIT_AUDIO)
        Mix_OpenAudio(22050, MIX_DEFAULT_FORMAT, 2, 1024)

        med = Mix_LoadWAV_RW(SDL_RWFromFile(joinpath(audio_example_assets_dir, "medium.wav"), "rb"), 1)
        @test med != C_NULL
        @test Mix_PlayChannelTimed(-1, med, 0, 500) != -1
        # Test that can play multiple times successfully
        # (Note that if a sound overlaps with a previous sound, it will play on a
        # different channel. The return value is which channel it plays on.)
        @test Mix_PlayChannelTimed(-1, med, 0, 500) != -1
        @test Mix_PlayChannelTimed(-1, med, 0, 500) != -1

        # Test different overlapping sounds
        scratch = Mix_LoadWAV_RW(SDL_RWFromFile(joinpath(audio_example_assets_dir, "scratch.wav"), "rb"), 1)
        @test scratch != C_NULL
        @test Mix_PlayChannelTimed(-1, scratch, 0, 500) != -1
        @test Mix_PlayChannelTimed(-1, med, 0, 500) != -1
        @test Mix_PlayChannelTimed(-1, scratch, 0, 500) != -1

        Mix_CloseAudio()
        SDL_Quit()
    end

    @testset "Music" begin
        SDL_Init(SDL_INIT_AUDIO)
        Mix_OpenAudio(22050, MIX_DEFAULT_FORMAT, 2, 1024)

        # Load the music
        music = Mix_LoadMUS(joinpath(audio_example_assets_dir, "beat.wav"))
        @test music != C_NULL

        # Test playing/pausing the music
        @test Mix_PlayMusic(music, -1) == 0
        Mix_PauseMusic()
        Mix_ResumeMusic()
        @test Mix_HaltMusic() == 0

        # Test noops if no music is playing.
        @test Mix_HaltMusic() == 0
        Mix_PauseMusic()
        Mix_ResumeMusic()

        # Test playing multiple times
        @test Mix_PlayMusic(music, -1) == 0
        @test Mix_PlayMusic(music, -1) == 0

        Mix_CloseAudio()
        SDL_Quit()
    end

end
SDL_Quit()


text:
### A Pluto.jl notebook ###
# v0.12.4

using Markdown
using InteractiveUtils

# ╔═╡ cd661db2-0fde-11eb-30b0-df84785200a4
md"""
## Introduction to abstraction
- Abstraction is the process of hiding the details. Here we play with matrices of various types of elements, and show an example of abstraction where the mechanics of arrays is hidden from the mechanics of elements stored in the arrays.
- https://www.youtube.com/watch?v=foN1_hAGfNg
"""

# ╔═╡ f2522c40-0fde-11eb-12f2-c16e2a50ee02
	md"""	Creating an Array
-Notice - Value assigned to Element decides whether its Float of Integeer
"""

# ╔═╡ f23ef262-0fde-11eb-3770-4fe3e2e3a22f
element = 1

# ╔═╡ f220e310-0fde-11eb-3e91-1f8844c3a09b
	md"""now use this element in an array of 3by4"""

# ╔═╡ f20c97c0-0fde-11eb-0d65-43d676f26b8f
fill(element,3,4)

# ╔═╡ f1fba7ce-0fde-11eb-0b31-c50aa56730fc
	md"""Run type of command to see it data type"""

# ╔═╡ f1d3fb90-0fde-11eb-2fd0-6b4de2c6f079
typeof(element)

# ╔═╡ f1bf8930-0fde-11eb-0ad9-e39de1fcb0f9
	md"""Repeating the whole steps above this time as Float...watch"""

# ╔═╡ f1ae7230-0fde-11eb-3eb5-c7940bd9af8b
element01 = 1.0

# ╔═╡ f794c600-0fe2-11eb-11e4-15daa9f0b60a
fill(element01,3,4)

# ╔═╡ ffcea390-0fe2-11eb-1a3b-191c7f621105
typeof(element01)

# ╔═╡ f19a4df2-0fde-11eb-0bcc-73bbc9df9ec4
	md"""Below is very Important - below 2 commands do the same thing, but putting a comma makes it a vector and without comma makes it a array"""

# ╔═╡ f184ca20-0fde-11eb-3438-e1e2cf1ea97a
How_you_Make_Array = [1 1]

# ╔═╡ 9eabc730-0fe4-11eb-32e7-35867c6fd643
typeof(How_you_Make_Array)

# ╔═╡ 270fa930-0fe4-11eb-2446-5102fab49ac4
How_you_make_Vector = [1, 1]

# ╔═╡ 9dad6a00-0fe4-11eb-3465-b1fbcdefa2c8
typeof(How_you_make_Vector)

# ╔═╡ f16e5bf0-0fde-11eb-0843-eb66ae5e4ab3
	md"""Going back to our Element examples now"""

# ╔═╡ 4b308310-0fe5-11eb-2d54-a3c7b7f673b7
	md"""Same experiment with Strings ----"""

# ╔═╡ f1529690-0fde-11eb-0a49-b76a8afd8d29
element02 = "Hello"

# ╔═╡ 37dbc810-0fe5-11eb-29f6-cbdb796e3ff7
fill(element02,3,4)

# ╔═╡ 43ef0f40-0fe5-11eb-04c6-6f464fea37ec
typeof(element02)

# ╔═╡ f13b16f0-0fde-11eb-10dd-89a346d0ebfa
	md"""Same experiment with Rational --"""

# ╔═╡ f13327b0-0fde-11eb-3be5-8d3161ee2526
element03 = 1//1

# ╔═╡ 88746930-0fe5-11eb-37ba-fde0567d3584
fill(element03,3,4)

# ╔═╡ 88521420-0fe5-11eb-0b63-1f00a7558ebf
typeof(element03)

# ╔═╡ f0e46b70-0fde-11eb-101a-a35aff79f234
	md"""Showing them all together"""

# ╔═╡ f0bff37e-0fde-11eb-31df-8d0b8f5e0234
KeepTrack = [typeof(element), typeof(element01), typeof(element02), typeof(element03)]

# ╔═╡ f0aa96c0-0fde-11eb-2c6b-0f0a05260f8d
	md"""Can even see the type of KeepTrack, which will reveal it is a vector holding DataTypes"""

# ╔═╡ f090cd32-0fde-11eb-27ee-11b56e69b4b0
typeof(KeepTrack)

# ╔═╡ f07a8610-0fde-11eb-3047-4de6b9727e21
	md"""Lets see , if we can do the same on an image, before that load the required package"""

# ╔═╡ f7a62e50-0fe6-11eb-0ccc-e55edd0613b3
begin
	import Pkg
	Pkg.activate(mktempdir())
end

# ╔═╡ f05af020-0fde-11eb-12a3-915d7357a483
begin
	Pkg.add(["Images", "ImageIO", "ImageMagick", "TestImages"])
	using Images
end

# ╔═╡ f05ac910-0fde-11eb-1df2-a104c15a70bd
	md"""Download the Image"""

# ╔═╡ f0297fe0-0fde-11eb-33ee-2b7820d58ab3
thumpsup = "https://cdn1.sportngin.com/attachments/text_block/1723/1679/thumbs-up_small.png"

# ╔═╡ 3d3752ee-0fe7-11eb-2a47-f5e6e1c3a6b7
thumps = download(thumpsup, "Thumps_Image45.png")

# ╔═╡ 3d2307a0-0fe7-11eb-0cec-f93b5d63d8f0
ThumpsUp_01 = load(thumps)

# ╔═╡ 3d0c7260-0fe7-11eb-30ed-e95e8d2f1b3b
	md"""Beack to same thing, let create an array of it"""

# ╔═╡ 3cf9fbd0-0fe7-11eb-0135-19b0d155013c
element04 = ThumpsUp_01

# ╔═╡ 484a24a0-0fe8-11eb-2c66-a5a1f436af7e
fill(element04,3,4)

# ╔═╡ 59bcb220-0fe8-11eb-3f6b-c19546af91b8
KeepTrack1 = [typeof(element), typeof(element01), typeof(element02), typeof(element03), typeof(element04)]

# ╔═╡ 3ce36690-0fe7-11eb-226f-030d85b25e99
	md"""Last experimt was element as an Array, i can create an array of an array, by saying element05 = 2,3 4,5 something like that"""

# ╔═╡ 3cd2769e-0fe7-11eb-0cba-0144726b8b1e

# ╔═╡ 3cbbba50-0fe7-11eb-163c-afcb64398fe2
	md"""-------------------------"""

# ╔═╡ 3ca76f00-0fe7-11eb-0b72-37ebb4860d54

# ╔═╡ 3c94351e-0fe7-11eb-3170-3ff3c844cd4f
	md"""-------------------------"""

# ╔═╡ 3c7d9fde-0fe7-11eb-0746-9b56023f22b4

# ╔═╡ 1f08c5b0-0fe8-11eb-0741-b3d568e15342
	md"""-------------------------"""

# ╔═╡ 1ef5b2e0-0fe8-11eb-1f12-09dbbfec87e2

# ╔═╡ 1ee0a440-0fe8-11eb-12f4-8164440fa6b7
	md"""-------------------------"""

# ╔═╡ 1ec88860-0fe8-11eb-13de-530273037934

# ╔═╡ 1ea548f0-0fe8-11eb-1512-d14aaf5b6e43
	md"""-------------------------"""

# ╔═╡ 1e920f10-0fe8-11eb-0a99-9be2826c2e20

# ╔═╡ 1e7fbf90-0fe8-11eb-2510-1fcb2efc537f

# ╔═╡ 1e69c690-0fe8-11eb-1c05-713de2e85f8f

# ╔═╡ 1e54b7f0-0fe8-11eb-17d7-7d838b63f38c

# ╔═╡ 1e4352ce-0fe8-11eb-2198-addfa7f7cf92

# ╔═╡ 1e29d760-0fe8-11eb-1dab-8ff26e1cb786

# ╔═╡ Cell order:
# ╟─cd661db2-0fde-11eb-30b0-df84785200a4
# ╟─f2522c40-0fde-11eb-12f2-c16e2a50ee02
# ╠═f23ef262-0fde-11eb-3770-4fe3e2e3a22f
# ╟─f220e310-0fde-11eb-3e91-1f8844c3a09b
# ╠═f20c97c0-0fde-11eb-0d65-43d676f26b8f
# ╟─f1fba7ce-0fde-11eb-0b31-c50aa56730fc
# ╠═f1d3fb90-0fde-11eb-2fd0-6b4de2c6f079
# ╟─f1bf8930-0fde-11eb-0ad9-e39de1fcb0f9
# ╠═f1ae7230-0fde-11eb-3eb5-c7940bd9af8b
# ╠═f794c600-0fe2-11eb-11e4-15daa9f0b60a
# ╠═ffcea390-0fe2-11eb-1a3b-191c7f621105
# ╟─f19a4df2-0fde-11eb-0bcc-73bbc9df9ec4
# ╟─f184ca20-0fde-11eb-3438-e1e2cf1ea97a
# ╟─9eabc730-0fe4-11eb-32e7-35867c6fd643
# ╟─270fa930-0fe4-11eb-2446-5102fab49ac4
# ╟─9dad6a00-0fe4-11eb-3465-b1fbcdefa2c8
# ╟─f16e5bf0-0fde-11eb-0843-eb66ae5e4ab3
# ╟─4b308310-0fe5-11eb-2d54-a3c7b7f673b7
# ╟─f1529690-0fde-11eb-0a49-b76a8afd8d29
# ╟─37dbc810-0fe5-11eb-29f6-cbdb796e3ff7
# ╠═43ef0f40-0fe5-11eb-04c6-6f464fea37ec
# ╠═f13b16f0-0fde-11eb-10dd-89a346d0ebfa
# ╟─f13327b0-0fde-11eb-3be5-8d3161ee2526
# ╟─88746930-0fe5-11eb-37ba-fde0567d3584
# ╟─88521420-0fe5-11eb-0b63-1f00a7558ebf
# ╟─f0e46b70-0fde-11eb-101a-a35aff79f234
# ╟─f0bff37e-0fde-11eb-31df-8d0b8f5e0234
# ╟─f0aa96c0-0fde-11eb-2c6b-0f0a05260f8d
# ╟─f090cd32-0fde-11eb-27ee-11b56e69b4b0
# ╟─f07a8610-0fde-11eb-3047-4de6b9727e21
# ╠═f7a62e50-0fe6-11eb-0ccc-e55edd0613b3
# ╠═f05af020-0fde-11eb-12a3-915d7357a483
# ╠═f05ac910-0fde-11eb-1df2-a104c15a70bd
# ╠═f0297fe0-0fde-11eb-33ee-2b7820d58ab3
# ╠═3d3752ee-0fe7-11eb-2a47-f5e6e1c3a6b7
# ╠═3d2307a0-0fe7-11eb-0cec-f93b5d63d8f0
# ╠═3d0c7260-0fe7-11eb-30ed-e95e8d2f1b3b
# ╠═3cf9fbd0-0fe7-11eb-0135-19b0d155013c
# ╠═484a24a0-0fe8-11eb-2c66-a5a1f436af7e
# ╟─59bcb220-0fe8-11eb-3f6b-c19546af91b8
# ╟─3ce36690-0fe7-11eb-226f-030d85b25e99
# ╠═3cd2769e-0fe7-11eb-0cba-0144726b8b1e
# ╠═3cbbba50-0fe7-11eb-163c-afcb64398fe2
# ╠═3ca76f00-0fe7-11eb-0b72-37ebb4860d54
# ╠═3c94351e-0fe7-11eb-3170-3ff3c844cd4f
# ╠═3c7d9fde-0fe7-11eb-0746-9b56023f22b4
# ╠═1f08c5b0-0fe8-11eb-0741-b3d568e15342
# ╠═1ef5b2e0-0fe8-11eb-1f12-09dbbfec87e2
# ╠═1ee0a440-0fe8-11eb-12f4-8164440fa6b7
# ╠═1ec88860-0fe8-11eb-13de-530273037934
# ╠═1ea548f0-0fe8-11eb-1512-d14aaf5b6e43
# ╠═1e920f10-0fe8-11eb-0a99-9be2826c2e20
# ╠═1e7fbf90-0fe8-11eb-2510-1fcb2efc537f
# ╠═1e69c690-0fe8-11eb-1c05-713de2e85f8f
# ╠═1e54b7f0-0fe8-11eb-17d7-7d838b63f38c
# ╠═1e4352ce-0fe8-11eb-2198-addfa7f7cf92
# ╠═1e29d760-0fe8-11eb-1dab-8ff26e1cb786


text:
struct Foo <: Recommender
    data::DataAccessor
end

function test_not_implemented_error()
    println("-- Testing recommender not implemented case")
    data = DataAccessor(sparse([1 2 3; 4 5 6]))
    recommender = Foo(data)
    @test_throws ErrorException build!(recommender)
    @test_throws ErrorException predict(recommender, 1, 1)
    @test_throws ErrorException ranking(recommender, 1, 1)
end

function test_not_build_error()
    println("-- Testing recommender not built case")

    # non-personalized (MostPopular) recommendation for 3 items
    data = DataAccessor(sparse([1 0 0; 4 5 0]))
    recommender = MostPopular(data)

    # build!(recommender) <- should be called before recommend()

    @test_throws ErrorException recommend(recommender, 1, 3, [1, 2, 3])
end

function test_recommend()
    println("-- Testing recommender execution")

    # non-personalized (MostPopular) recommendation for 3 items
    data = DataAccessor(sparse([1 0 0; 4 5 0]))
    recommender = MostPopular(data)
    build!(recommender)
    pairs = recommend(recommender, 1, 3, [1, 2, 3])

    @test first(pairs[1]) == 1
    @test last(pairs[1]) == 2
    @test first(pairs[2]) == 2
    @test last(pairs[2]) == 1
    @test first(pairs[3]) == 3
    @test last(pairs[3]) == 0
end

test_not_build_error()
test_not_implemented_error()
test_recommend()


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module XMLRPC

using LightXML, Requests
export XMLRPCProxy, XMLRPCMethodCall, XMLRPCCall, xmlrpc_parse

"""
An XML RPC Proxy wrapper type for the server URL.
"""
immutable XMLRPCProxy
    url::AbstractString
end

"""
An XMLRPC call used for dispatch.
"""
immutable XMLRPCMethodCall
    proxy::XMLRPCProxy
    name::AbstractString
end

XMLRPCArguments = Union{Int32, Bool, AbstractString, Float64, DateTime, Dict, Vector} # base64 ommited

"""
A fully determined XMLRPC call with parameters specified
"""
immutable XMLRPCCall
    method::XMLRPCMethodCall
    parameters::Tuple
end

function Base.getindex(proxy::XMLRPCProxy, s::AbstractString)
    function _(m...)
        res = post(XMLRPCCall(XMLRPCMethodCall(proxy, ASCIIString(s)), m))
        h = headers(res)
        h["Content-Type"] == "text/xml" || error("Content is not text/xml!\n"*readall(res))
        r = readall(res)
        xmlrpc_parse(r)
    end
end

"""
Execute an XMLRPC call via HTTP POST.
"""
function Requests.post(x::XMLRPCCall)
    xdoc = XMLDocument(x)
    post(url(x); headers=Dict("Content-type" => "text/xml"), data=string(xdoc))
end

url(x::XMLRPCCall) = x.method.proxy.url

"""
Convert a `XMLRPCCall` into XML.
"""
function LightXML.XMLDocument(x::XMLRPCCall)
    xdoc = XMLDocument()
    xroot = create_root(xdoc, "methodCall")
    xs1 = new_child(xroot, "methodName")
    add_text(xs1, x.method.name)
    params = new_child(xroot, "params")
    for p in x.parameters
        rpc_arg(new_child(params, "param"), p)
    end
    xdoc
end

"""
Convert a value, Dict, or Vector into an XMLRPC snippet.
"""
function rpc_arg(x::XMLElement, p::Int32)
    add_text(new_child(new_child(x, "value"), "int"), string(p))
end

function rpc_arg(x::XMLElement, p::Int64)
    add_text(new_child(new_child(x, "value"), "int"), string(Int32(p)))
end

function rpc_arg(x::XMLElement, p::Bool)
    add_text(new_child(new_child(x, "value"), "boolean"), p?"1":"0")
end

function rpc_arg(x::XMLElement, p::Float64)
    add_text(new_child(new_child(x, "value"), "double"), string(p))
end

function rpc_arg(x::XMLElement, p::AbstractString)
    add_text(new_child(new_child(x, "value"), "string"), p)
end

function rpc_arg(x::XMLElement, p::DateTime)
    add_text(new_child(new_child(x, "value"), "dateTime.iso8601"), string(p))
end

# TODO new_child for base64

function rpc_arg(x::XMLElement, p::Vector)
    d = new_child(new_child(x, "array"), "data")
    for e in p
        rpc_arg(d, e)
    end
end

function rpc_arg(x::XMLElement, d::Dict)
    s = new_child(x, "struct")
    for p in d
        rpc_arg(s, p)
    end
end

function rpc_arg(x::XMLElement, p::Pair)
    m = new_child(x, "member")
    n = new_child(m, "name")
    add_text(n, string(p.first))
    rpc_arg(m, p.second)
end

function xmlrpc_parse(s::AbstractString)
    x = parse_string(s)
    xroot = root(x)
    name(xroot) == "methodResponse" || error("malformed XMLRPC response")
    xmlrpc_parse(collect(child_elements(xroot))[1])
end

function xmlrpc_parse(x::XMLElement)
    if name(x) == "value"
        c = collect(child_elements(x))[1]
        if name(c) == "i4" || name(c) == "int"
            return parse(Int32, content(c))
        elseif name(c) == "dateTime.iso8601"
            return DateTime(content(c))
        elseif name(c) == "boolean"
            return content(c) == "true" || content(c) == "1"
        elseif name(c) == "double"
            return parse(Float64, content(c))
        elseif name(c) == "base64"
            return base64decode(content(c))
        elseif name(c) == "string"
            return content(c)
        elseif name(c) == "array"
            c = collect(child_elements(c))[1] # <data>
            arr = []
            for elt in child_elements(c)
                push!(arr, xmlrpc_parse(elt))
            end
            return arr
        elseif name(c) == "struct"
            d = Dict()
            for elt in child_elements(c)
                push!(d, xmlrpc_parse(elt))
            end
            return d
        end
    elseif name(x) == "member"
        c = collect(child_elements(x))
        n = content(c[1]) # name
        v = xmlrpc_parse(c[2]) # value
        return Pair(n,v)
    elseif name(x) == "params" || name(x) == "param" # always one param on return
        return xmlrpc_parse(collect(child_elements(x))[1])
    elseif name(x) == "fault"
        error("XMLRPC Fault:\n$x")
    end
end

end # module


text:
module MeshAdaptiveDirectSearch

using StaticArrays, Random, ElasticArrays, LinearAlgebra, Primes
export MADS, LtMADS, OrthoMADS, RobustMADS, RobustLtMADS, RobustOrthoMADS,
Silent, Progress, minimize, NegReduction, NoReduction

# TODO:
# * QrMADS
# * MADS(suc,neg) and variants

####
#### Meshes
####

# implements Eq. 2.1 from Audet & Dennis 2006
mutable struct Mesh
    τ::Float64
    w⁺::Int
    w⁻::Int
    Δᵐ::Float64
end
Mesh(; τ = 4., Δᵐ = 1., w⁺ = 1, w⁻ = -1) = Mesh(τ, w⁺, w⁻, Δᵐ)
function update!(m::Mesh, i)
    i == 0 && return m
    if i > 0
        m.Δᵐ == 1. && return m
        w = rand(0:m.w⁺)
    else
        w = rand(m.w⁻:-1)
    end
    m.Δᵐ *= m.τ^w
    m
end
Δ(m::Mesh) = m.Δᵐ
ℓ(m::Mesh) = round(Int, -log(m.τ, m.Δᵐ))
# implements Eq. 2.1 from Audet & Dennis 2006 for w⁺ = -w⁻ = 1
mutable struct LogMesh
    τ::Int
    neglogΔᵐ::Int
end
LogMesh(; τ = 4, Δᵐ = 1) = LogMesh(τ, Int(-log(τ, Δᵐ)))
function update!(m::LogMesh, i)
    i == 0 && return m
    if i > 0
        m.neglogΔᵐ -= 1
    else
        m.neglogΔᵐ += 1
    end
    m
end
Δ(m::LogMesh) = min(1, (1/m.τ)^m.neglogΔᵐ)
ℓ(m::LogMesh) = m.neglogΔᵐ

####
#### LT
####

# Implements generation of the direction b(l), box on p. 203 from Audet &  Dennis 2006
struct LTDirectionGenerator{N}
    b::Dict{Int, NamedTuple{(:i, :v), Tuple{Int, SVector{N,Int}}}}
end
LTDirectionGenerator(N) = LTDirectionGenerator{N}(Dict{Int, NamedTuple{(:i, :v), Tuple{Int, SVector{N,Int}}}}())
function b(d::LTDirectionGenerator{N}, l::Int) where N
    haskey(d.b, l) && return d.b[l]
    i = rand(1:N)
    d.b[l] = (i = i,
              v = SVector{N}([rand((-1,1)) * (k == i ? 2^l : (2^l - 1)) for k in 1:N]))
end
iterator(g::LTDirectionGenerator, l) = LTDirectionIterator(g, max(0, l))

# Implements generation of the positive bias, box on p. 204 from Audet & Dennis 2006
struct LTDirectionIterator{N,Np}
    l::Int
    i::Int
    b::SVector{N,Int}
    rowperm::SVector{Np,Int}
    colperm::SVector{Np,Int}
    sum::MVector{N,Int}
end
function LTDirectionIterator(generator::LTDirectionGenerator{N}, l) where N
    i, v = b(generator, l)
    LTDirectionIterator{N,N-1}(l, i, v,
                               SVector{N-1}(randperm(N-1)),
                               SVector{N-1}(randperm(N-1)),
                               zeros(Int, N))
end

function L(i, j, l)
    j > i && return 0
    j == i && return rand((-1, 1)) * 2^l
    return rand(-2^l + 1:2^l - 1)
end

import Base: iterate, length
length(::LTDirectionIterator{N,Np}) where {N, Np} = N + 1
function iterate(it::LTDirectionIterator{N,Np}, j = 1) where {N, Np}
    j == N + 2 && return nothing
    j == N + 1 && return -SVector(it.sum), j + 1
    j == 1 && (it.sum .*= 0)
    if j == it.i
        v = it.b
    else
        v = SVector{N}([i == it.i ? 0 : L(it.rowperm[i - (i > it.i)],
                               it.colperm[j - (j > it.i)], it.l)
                        for i in 1:N])
    end
    it.sum .+= v
    v, j + 1
end

####
#### MADS
####

# implements general MADS, p. 193 from Audet & Dennis, 2006
abstract type AbstractMADS end
"""
    struct MADS{Tmesh,Tsearch,Tpoll} <: AbstractMADS
        mesh::Tmesh
        search::Tsearch
        poll::Tpoll
"""
struct MADS{Tmesh,Tsearch,Tpoll} <: AbstractMADS
    mesh::Tmesh
    search::Tsearch
    poll::Tpoll
end
function MADS(N; search = NoSearch(),
                 poll = LTDirectionGenerator(N),
                 mesh = LogMesh())
    MADS(mesh, search, poll)
end
"""
     LtMADS(N; search = NoSearch(), mesh = LogMesh())

Returns a `MADS` object with `poll = LTDirectionGenerator(N)` where `N` is the dimenionality of the problem.
See Audet & Dennis (2006), section 4, LTMADS
"""
function LtMADS(N; search = NoSearch(), mesh = LogMesh())
    MADS(mesh, search, LTDirectionGenerator(N))
end
"""
    OrthoMADS(N; search = NoSearch(), mesh = LogMesh(), reduction = NegReduction(N))

Returns a `MADS` object with `poll = OrthoDirectionGenerator(N)` where `N` is
the dimenionality of the problem and `reduction` can be
[`NegReduction(N)`](@ref) or [`NoReduction()`](@ref).
See Abramson et al. (2009), ORTHOMADS and Audet et al. 2014 for NegReduction.
"""
function OrthoMADS(N; search = NoSearch(),
                      mesh = LogMesh(),
                      reduction = NegReduction(N))
    MADS(mesh, search, OrthoDirectionGenerator(N, reduction = reduction))
end

# implements robust mads, from Audet et al. 2018
struct Cache
    x::ElasticArray{Float64,2,1}
    y::Vector{Float64}
    incumbents::Vector{Int}
end
Cache(N) = Cache(ElasticArray{Float64}(undef, N, 0), Float64[], Int[])
import Base.push!
push!(c::Cache, x, y) = begin append!(c.x, x); push!(c.y, y) end
"""
    struct RobustMADS{Tmesh,Tsearch,Tpoll,Tkernel} <: AbstractMADS
        mesh::Tmesh
        search::Tsearch
        poll::Tpoll
        kernel::Tkernel
        cache::Cache
        f::Vector{Float64}
        P::Vector{Float64}
"""
struct RobustMADS{Tmesh,Tsearch,Tpoll,Tkernel} <: AbstractMADS
    mesh::Tmesh
    search::Tsearch
    poll::Tpoll
    kernel::Tkernel
    cache::Cache
    f::Vector{Float64}
    P::Vector{Float64}
end
"""
    RobustMADS(N; search = NoSearch(), poll = LTDirectionGenerator(N),
                  mesh = LogMesh(), kernel = GaussKernel(1, 1), cache = Cache(N))

Returns a `RobustMADS` object where `N` is the dimenionality of the problem.
See Audet et al. 2018.
"""
function RobustMADS(N; search = NoSearch(), poll = LTDirectionGenerator(N),
                    mesh = LogMesh(), kernel = GaussKernel(1, 1), cache = Cache(N))
    RobustMADS(mesh, search, poll, kernel, cache, Float64[], Float64[])
end
"""
    RobustLtMADS(N; kwargs...)

Returns a `RobustMADS` object with `poll =  LTDirectionGenerator(N)` where `N` is the dimenionality of the problem.
"""
function RobustLtMADS(N; kwargs...)
    RobustMADS(N, poll = LTDirectionGenerator(N), kwargs...)
end
"""
    RobustOrthoMADS(N; kwargs...)

Returns a `RobustMADS` object with `poll =  OrthoDirectionGenerator(N)` where `N` is the dimenionality of the problem.
"""
function RobustOrthoMADS(N; kwargs...)
    RobustMADS(N, poll = OrthoDirectionGenerator(N), kwargs...)
end
mutable struct GaussKernel
    β::Float64
    σ²::Float64
end
(g::GaussKernel)(x, y) = g() * reshape(exp.(-sum((x .- y).^2, dims = 1)/(2*g.σ²)), :)
(g::GaussKernel)() = 1/sqrt(2*π*g.σ²)
update!kernel!(g::GaussKernel, mesh) = g.σ² = (g.β * Δ(mesh))^2
function isnewincumbent(m::RobustMADS, x, fx, oldfx)
    update!kernel!(m.kernel, m.mesh)
    if length(m.cache.x) > 0
        ψ = m.kernel(m.cache.x, x)
        m.f .*= m.P
        m.f .+= ψ * fx
        m.P .+= ψ
        m.f ./= m.P
        push!(m.P, sum(ψ) + m.kernel())
        push!(m.f, (dot(m.cache.y, ψ) + m.kernel() * fx)/m.P[end])
    else
        push!(m.P, m.kernel())
        push!(m.f, fx)
    end
    append!(m.cache.x, x)
    push!(m.cache.y, fx)
    i = argmin(m.f)
    success = i == length(m.f)
    cachesuccess = length(m.cache.incumbents) == 0 || i != m.cache.incumbents[end]
    success || cachesuccess && push!(m.cache.incumbents, i)
    success ? 1 : cachesuccess ? 0 : -1, m.cache.x[:, i], m.f[i]
end

# barrier
function isvalid(cs, x)
    for c in cs
        c(x) || return false
    end
    return true
end

# poll stage
update!(::Any, ::Any) = nothing
init!(::Any, ::Any) = nothing
poll(m, f, constraints, x, fx) = poll(m, iterator(m.poll, ℓ(m.mesh)), Δ(m.mesh), f, constraints, x, fx)
isnewincumbent(m::MADS, x, fx, oldfx) = fx < oldfx ? 1 : -1, x, fx
@inline function poll(m, it, Δᵐ, f, constraints, x, fx)
    for v in it
        newx = clamp!(x .+ Δᵐ * v, -1., 1.)
        isvalid(constraints, newx) || continue
        fnewx = f(newx)
#         @show fnewx
        success, newincumbent, fnewx = isnewincumbent(m, newx, fnewx, fx)
        if success >= 0
            update!(m.poll, newincumbent)
            return (incumbent = newincumbent, fx = fnewx, hasimproved = success)
        end
    end
    (incumbent = x, fx = fx, hasimproved = -1)
end

# search stage
struct NoSearch end
search(m, f, constraints, x, fx) = search(m.search, f, constraints, x, fx)
search(::NoSearch, f, constraints, x, fx) = (incumbent = x, fx = fx, hasimproved = 0)

# implements OthoMADS, Abramson et al. 2009
@inline function haltonnumber(base::Integer, index::Integer)::Float64
    res = 0.
    f = 1 / base
    i = index
    while i > 0
       res += f * (i % base)
       i = div(i, base)
       f = f / base
    end
    res
end
struct HaltonIterator{N} end
function iterate(it::HaltonIterator{N}, i = 1) where N
    [haltonnumber(prime(k), i) for k in 1:N], i + 1
end
function normalized_halton_direction(u, l)
    α = 2^(l/2)/sqrt(length(u)) - 1/2
    q = normalize(2 * u .- ones(length(u)))
    while norm(round.(α * q)) < 2^(l/2) # TODO: make more efficient
        α += .1
    end
    round.(Int, (α - .1) * q)
end
scaledhouseholder(q) = sum(q.^2) * I - 2 * q * q'
"""
    NoReduction()

Implements no reduction for OrthoMADS.
See Abramson et al. 2009
"""
struct NoReduction end
struct NegReduction
    oldincumbent::Vector{Float64}
    w::Vector{Float64}
end
"""
     NegReduction(N)

Implements (suc, neg) reduction of OrthoMADS.
See Audet et al. 2014.
"""
NegReduction(N) = NegReduction(zeros(N), fill(NaN, N))
mutable struct OrthoDirectionGenerator{N,R}
    reduction::R
    t₀::Int
    ℓmax::Int
    tmax::Int
end
function OrthoDirectionGenerator(N; t0 = 2*N,
                                    reduction = NegReduction(N))
    OrthoDirectionGenerator{N,typeof(reduction)}(reduction, t0, 0, 0)
end
iterator(g::OrthoDirectionGenerator, l) = OrthoDirectionIterator(g, l)
init!(g::OrthoDirectionGenerator{N,NegReduction}, x) where N = g.reduction.oldincumbent .= x
function update!(g::OrthoDirectionGenerator{N,NegReduction}, x) where N
    @. g.reduction.w = x - g.reduction.oldincumbent
    @. g.reduction.oldincumbent = x
    return nothing
end
struct OrthoDirectionIterator{N, R}
    reduction::R
    H::Matrix{Float64}
end
function determine_t!(g::OrthoDirectionGenerator, ℓ)
    if ℓ >= g.ℓmax
        g.ℓmax = ℓ
        ℓ  + g.t₀ > g.tmax && (g.tmax = ℓ + g.t₀)
        return ℓ + g.t₀
    else
        g.tmax += 1
        return g.tmax
    end
end
function OrthoDirectionIterator(g::OrthoDirectionGenerator{N,R}, ℓ) where {N,R}
    t = determine_t!(g, ℓ)
    u = first(iterate(HaltonIterator{N}(), t))
    q = normalized_halton_direction(u, ℓ)
    H = scaledhouseholder(q)
    if R === NegReduction && isnan(g.reduction.w[1])
        reduction = NoReduction()
    else
        reduction = g.reduction
    end
    OrthoDirectionIterator{N,typeof(reduction)}(reduction, H)
end
length(::OrthoDirectionIterator{N}) where N = 2N
length(::OrthoDirectionIterator{N,NegReduction}) where N = N + 1
function iterate(it::OrthoDirectionIterator{N}, i = 1) where N
    i > 2N && return nothing
    i <= N && return @view(it.H[:, i]), i + 1
    return -@view(it.H[:, i - N]), i + 1
end
function iterate(it::OrthoDirectionIterator{N,NegReduction}, state = (1, zeros(N))) where N
    i, sumd = state
    i > N + 1 && return nothing
    if i <= N
        d = @view(it.H[:, i])
        sign = dot(d, it.reduction.w) >= 0 ? 1 : -1
        d .*= sign
        sumd .-= d
        return d, (i + 1, sumd)
    else
        return sumd, (i + 1, sumd)
    end
end

function standardtransformation(lowerbound, upperbound)
    d = upperbound - lowerbound
    (to = x -> @.((x + 1)/2 * (upperbound - lowerbound) + lowerbound),
     from = x -> @.(2 * (x - lowerbound) / (upperbound - lowerbound) - 1))
end

@enum Verbosity Silent Progress
@enum StoppingReason MaxIterations MinMeshSize
"""
    minimize(m, f, x0 = zeros(length(x0));
                   lowerbound = -ones(length(x0)),
                   upperbound = ones(length(x0)),
                   max_iterations = 10^4,
                   min_mesh_size = eps(Float64)/2,
                   constraints = [],
                   verbosity = Silent)

Minimize function `f` with method `m`.
For possible methods `m` see [`LtMADS`](@ref), [`OrthoMADS`](@ref), [`MADS`](@ref),
[`RobustLtMADS`](@ref), [`RobustOrthoMADS`](@ref), [`RobustMADS`](@ref).
Constraints can be defined by boolean functions, e.g.
`constraints = [x -> sum(x) > 1, x -> x[1]^2 < 3]`.
Possible `verbosity`-levels are `Silent` or `Progress`.
"""
function minimize(m::AbstractMADS, f, x0 = zeros(length(x0));
                  lowerbound = -ones(length(x0)),
                  upperbound = ones(length(x0)),
                  max_iterations = 10^4,
                  min_mesh_size = eps(Float64)/2,
                  verbosity = Silent,
                  constraints = [])
    isvalid(constraints, x0) || error("x0 = $x0 doesn't satisfy all constraints.")
    to, from = standardtransformation(lowerbound, upperbound)
    finternal = x -> f(to(x))
    cinternal = [x -> c(to(x)) for c in constraints]
    incumbent = from(x0)
    init!(m.poll, incumbent)
    fincumbent = finternal(incumbent)
    for k in 1:max_iterations
        incumbent, fincumbent, i = search(m, finternal, cinternal, incumbent, fincumbent)
        if i != 1
            incumbent, fincumbent, i = poll(m, finternal, cinternal, incumbent, fincumbent)
        end
        update!(m.mesh, i)
        Δ(m.mesh) < min_mesh_size && return (f = fincumbent,
                                             x = to(incumbent),
                                             stopping_reason = MinMeshSize,
                                             iterations = k)
        if verbosity > Silent
            @show k to(incumbent) fincumbent i Δ(m.mesh)
        end
    end
    (f = fincumbent, x = to(incumbent), stopping_reason = MaxIterations,
     iterations = max_iterations)
end
end # module


text:
using AbstractTensors
using Test, DirectSum, LinearAlgebra

# example data
struct SpecialTensor{V} <: TensorAlgebra{V} end
a,b = SpecialTensor{ℝ}(), SpecialTensor{ℝ'}()
@test ndims(+(a)) == ndims(b)

## tensor operation (trivial test)
op(s::SpecialTensor{V},::SpecialTensor{V}) where V = s
op(a::TensorAlgebra{V},b::TensorAlgebra{W}) where {V,W} = interop(op,a,b)
(W::Signature)(s::SpecialTensor{V}) where V = SpecialTensor{W}()
@test vectorspace(op(a,b)) == ℝ⊕ℝ'
@test vectorspace(interop(op,a,b)) == ℝ⊕ℝ'
@test vectorspace(op(a,a)) == ℝ
@test vectorspace(interop(op,a,a)) == ℝ
Base.:+(s::SpecialTensor{V},::SpecialTensor{V}) where V = s
@test vectorspace(+(a,b)) == ℝ⊕ℝ'
@test vectorspace(interop(+,a,b)) == ℝ⊕ℝ'
@test vectorspace(+(a,a)) == ℝ
@test vectorspace(interop(+,a,a)) == ℝ

## tensor pseudoscalar (trivial test)
op(a::TensorAlgebra{V},b::UniformScaling) where V = op(a,V(b))
op(a::UniformScaling,b::TensorAlgebra{V}) where V = op(V(a),b)
(W::Signature)(s::UniformScaling) where V = SpecialTensor{W}()
@test vectorspace(op(a,LinearAlgebra.I)) == ℝ
@test vectorspace(op(LinearAlgebra.I,a)) == ℝ
@test vectorspace(+(a,LinearAlgebra.I)) == ℝ
@test vectorspace(+(LinearAlgebra.I,a)) == ℝ

## tensor evaluation (trivial test)
(a::SpecialTensor{V})(b::SpecialTensor{V}) where V = a
(a::SpecialTensor{W})(b::SpecialTensor{V}) where {V,W} = interform(a,b)
@test vectorspace(b(a)) == ℝ⊕ℝ'
@test vectorspace(interform(b,a)) == ℝ⊕ℝ'
@test vectorspace(a(a)) == ℝ
@test vectorspace(interform(a,a)) == ℝ

## algebraic tests
@test !I == |(I)
@test (2)⁻¹ == 1/2
@test (im)ǂ == -im
@test (sqrt(2))ˣ == sqrt(2)
@test (sqrt(2))₊ == sqrt(2)
@test (sqrt(2))₋ == 0


text:
macro c(ret_type, func, arg_types, lib)
    local args_in = Any[Symbol(string('a',x)) for x in 1:length(arg_types.args)]
    quote
      $(esc(func))($(args_in...)) = ccall( ($(string(func)), $(esc(lib))), $(esc(ret_type)), $(esc(arg_types)), $(args_in...) )
    end
end

macro ctypedef(fake_t,real_t)
    quote
      const $(esc(fake_t)) = $(esc(real_t))
    end
end

@ctypedef XML_Char UInt8
@ctypedef XML_LChar UInt8
@ctypedef XML_Index Int32
@ctypedef XML_Size UInt32
@ctypedef XML_Parser Ptr{Nothing}
@ctypedef XML_Bool UInt8
# enum XML_Status
const XML_STATUS_ERROR = 0
const XML_STATUS_OK = 1
const XML_STATUS_SUSPENDED = 2
# end
# enum XML_Error
const XML_ERROR_NONE = 0
const XML_ERROR_NO_MEMORY = 1
const XML_ERROR_SYNTAX = 2
const XML_ERROR_NO_ELEMENTS = 3
const XML_ERROR_INVALID_TOKEN = 4
const XML_ERROR_UNCLOSED_TOKEN = 5
const XML_ERROR_PARTIAL_CHAR = 6
const XML_ERROR_TAG_MISMATCH = 7
const XML_ERROR_DUPLICATE_ATTRIBUTE = 8
const XML_ERROR_JUNK_AFTER_DOC_ELEMENT = 9
const XML_ERROR_PARAM_ENTITY_REF = 10
const XML_ERROR_UNDEFINED_ENTITY = 11
const XML_ERROR_RECURSIVE_ENTITY_REF = 12
const XML_ERROR_ASYNC_ENTITY = 13
const XML_ERROR_BAD_CHAR_REF = 14
const XML_ERROR_BINARY_ENTITY_REF = 15
const XML_ERROR_ATTRIBUTE_EXTERNAL_ENTITY_REF = 16
const XML_ERROR_MISPLACED_XML_PI = 17
const XML_ERROR_UNKNOWN_ENCODING = 18
const XML_ERROR_INCORRECT_ENCODING = 19
const XML_ERROR_UNCLOSED_CDATA_SECTION = 20
const XML_ERROR_EXTERNAL_ENTITY_HANDLING = 21
const XML_ERROR_NOT_STANDALONE = 22
const XML_ERROR_UNEXPECTED_STATE = 23
const XML_ERROR_ENTITY_DECLARED_IN_PE = 24
const XML_ERROR_FEATURE_REQUIRES_XML_DTD = 25
const XML_ERROR_CANT_CHANGE_FEATURE_ONCE_PARSING = 26
const XML_ERROR_UNBOUND_PREFIX = 27
const XML_ERROR_UNDECLARING_PREFIX = 28
const XML_ERROR_INCOMPLETE_PE = 29
const XML_ERROR_XML_DECL = 30
const XML_ERROR_TEXT_DECL = 31
const XML_ERROR_PUBLICID = 32
const XML_ERROR_SUSPENDED = 33
const XML_ERROR_NOT_SUSPENDED = 34
const XML_ERROR_ABORTED = 35
const XML_ERROR_FINISHED = 36
const XML_ERROR_SUSPEND_PE = 37
const XML_ERROR_RESERVED_PREFIX_XML = 38
const XML_ERROR_RESERVED_PREFIX_XMLNS = 39
const XML_ERROR_RESERVED_NAMESPACE_URI = 40
# end
# enum XML_Content_Type
const XML_CTYPE_EMPTY = 1
const XML_CTYPE_ANY = 2
const XML_CTYPE_MIXED = 3
const XML_CTYPE_NAME = 4
const XML_CTYPE_CHOICE = 5
const XML_CTYPE_SEQ = 6
# end
# enum XML_Content_Quant
const XML_CQUANT_NONE = 0
const XML_CQUANT_OPT = 1
const XML_CQUANT_REP = 2
const XML_CQUANT_PLUS = 3
# end
@ctypedef XML_Content Nothing
@ctypedef XML_ElementDeclHandler Ptr{Nothing}
@ctypedef XML_AttlistDeclHandler Ptr{Nothing}
@ctypedef XML_XmlDeclHandler Ptr{Nothing}
mutable struct XML_Memory_Handling_Suite
    malloc_fcn::Ptr{Nothing}
    realloc_fcn::Ptr{Nothing}
    free_fcn::Ptr{Nothing}
end
@ctypedef XML_StartElementHandler Ptr{Nothing}
@ctypedef XML_EndElementHandler Ptr{Nothing}
@ctypedef XML_CharacterDataHandler Ptr{Nothing}
@ctypedef XML_ProcessingInstructionHandler Ptr{Nothing}
@ctypedef XML_CommentHandler Ptr{Nothing}
@ctypedef XML_StartCdataSectionHandler Ptr{Nothing}
@ctypedef XML_EndCdataSectionHandler Ptr{Nothing}
@ctypedef XML_DefaultHandler Ptr{Nothing}
@ctypedef XML_StartDoctypeDeclHandler Ptr{Nothing}
@ctypedef XML_EndDoctypeDeclHandler Ptr{Nothing}
@ctypedef XML_EntityDeclHandler Ptr{Nothing}
@ctypedef XML_UnparsedEntityDeclHandler Ptr{Nothing}
@ctypedef XML_NotationDeclHandler Ptr{Nothing}
@ctypedef XML_StartNamespaceDeclHandler Ptr{Nothing}
@ctypedef XML_EndNamespaceDeclHandler Ptr{Nothing}
@ctypedef XML_NotStandaloneHandler Ptr{Nothing}
@ctypedef XML_ExternalEntityRefHandler Ptr{Nothing}
@ctypedef XML_SkippedEntityHandler Ptr{Nothing}
mutable struct XML_Encoding
    map::Nothing
    data::Ptr{Nothing}
    convert::Ptr{Nothing}
    release::Ptr{Nothing}
end
@ctypedef XML_UnknownEncodingHandler Ptr{Nothing}
# enum XML_Parsing
const XML_INITIALIZED = 0
const XML_PARSING = 1
const XML_FINISHED = 2
const XML_SUSPENDED = 3
# end
mutable struct XML_ParsingStatus
  parsing::Int32
  finalBuffer::XML_Bool
end
# enum XML_ParamEntityParsing
const XML_PARAM_ENTITY_PARSING_NEVER = 0
const XML_PARAM_ENTITY_PARSING_UNLESS_STANDALONE = 1
const XML_PARAM_ENTITY_PARSING_ALWAYS = 2
# end
mutable struct XML_Expat_Version
    major::Int32
    minor::Int32
    micro::Int32
end
# enum XML_FeatureEnum
const XML_FEATURE_END = 0
const XML_FEATURE_UNICODE = 1
const XML_FEATURE_UNICODE_WCHAR_T = 2
const XML_FEATURE_DTD = 3
const XML_FEATURE_CONTEXT_BYTES = 4
const XML_FEATURE_MIN_SIZE = 5
const XML_FEATURE_SIZEOF_XML_CHAR = 6
const XML_FEATURE_SIZEOF_XML_LCHAR = 7
const XML_FEATURE_NS = 8
const XML_FEATURE_LARGE_SIZE = 9
const XML_FEATURE_ATTR_INFO = 10
# end
mutable struct XML_Feature
    feature::Int32
    name::Ptr{XML_LChar}
    value::Int32
end


text:
using Statistics

function ising(L::Int, T::Number, next_I, prev_I; MCS=230_000::Int)
    lattice = ones(Int8, L,L)

    boltzmann_factor = map(ΔE -> exp(-ΔE / T), -8:4:8)

    #initial arrays (returned)
    magnetizations = zeros(0)
    energies = zeros(0)

    for k=1 : MCS
        for i=1:L, j=1:L
        #metropolis algorithm
            @inbounds ΔE::Int = 2 * lattice[i,j] * (lattice[next_I[i],j] + lattice[prev_I[i],j] + lattice[i,next_I[j]] + lattice[i,prev_I[j]])
            if (ΔE ≤ 0) || (rand() ≤ boltzmann_factor[ΔE÷4 + 3])
                @inbounds lattice[i,j] *= -1
            end
        end

        if (30_000 < k && k%100 == 0)
            m_k::Float64 = abs(mean(lattice))
            push!(magnetizations, m_k)
            E::Float64 = 0.0
            for i=1:L, j=1:L
                # E += 0.5 * lattice[i,j] * (lattice[next_I[i],j] + lattice[prev_I[i],j] + lattice[i,next_I[j]] + lattice[i,prev_I[j]])
                @inbounds E += lattice[i,j] * (lattice[next_I[i],j] + lattice[i,next_I[j]])
            end
            push!(energies, E)
        end
    end
    return magnetizations, energies
end

function run_MC(Ls::Array{Int}, Ts)
    result = Dict{Int, Dict{Char, Dict{Float64,Float64}}}(l => Dict(c => Dict() for c in ['m','χ','C']) for l in Ls)
    for L ∈ Ls
        next_I = collect(2:L+1)
        prev_I = collect(0:L-1)
        next_I[end] = 1
        prev_I[1] = L
        Threads.@threads for T ∈ Ts
            magnetizations, energies = ising(L,T,next_I, prev_I)

            result[L]['m'][T] = mean(magnetizations)                                                    # mean magnetization        <m>
            result[L]['χ'][T] = L^2 / T * var(magnetizations, mean=result[L]['m'][T], corrected=false)  # magnetic susceptibility    χ
            result[L]['C'][T] = var(energies, corrected=false)   / L^2 / T^2                            # heat capacity              C
        end
    end
    return result
end

function main()
    Ls = [10]
    # Ts = cat(0.5:0.05:1.75, 1.8:0.005:2.7, 2.75:0.05:3.5, dims=1)
    Ts = 0.5:0.1:3.5
    result_ret = @time run_MC(Ls,Ts)
end
main()


text:
function gradcheckCPU(net;showgrad=false,returngrad=false)
    #ONLY WORKS ON THE CPU AT THE MOMENT
    # (c) David Barber, University College London 2015
    println("GRADCHECK CPU")
    value=net.value
    node=net.node
    net2=deepcopy(net) ## CAREFUL -- DEEPCOPY IS BUGGY FOR CUDAARRAY

    ADforward!(net)
    ADbackward!(net)
    g=net.gradient

    epsilon=1e-7;
    tol=1e-6
    gemp=deepcopy(g)
    #for n=1:length(net.node) # loop over the nodes in the network
    for n in net.validnodes # loop over the nodes in the network
        if node[n].returnderivative
            println("------------------------------------------------------------------------------")
            println("node $n:")
            for linearidx in 1:prod(size(value[n]))
                net2.value[n][linearidx]=value[n][linearidx]+epsilon
                fplus=ADforward!(net2,returnf=true)
                net2.value[n][linearidx]=value[n][linearidx]-epsilon
                fminus=ADforward!(net2,returnf=true)
                net2.value[n][linearidx]=value[n][linearidx]
                gemp[n][linearidx]=0.5*(fplus-fminus)/epsilon
            end
            diff=mean(abs(g[n]-gemp[n]))
            reldiff=mean(abs((g[n]-gemp[n])./(realmin()+g[n])))
            println("absolute difference between analytic and empirical gradient=$diff")
            println("relative difference between analytic and empirical gradient=$reldiff")

                #if showgrad
                #println("analytic gradient:"); println(g[n])
                #println("empirical gradient:"); println(gemp[n])
                #end

            if diff>tol
                if showgrad
                println("analytic gradient:"); println(g[n])
                println("empirical gradient:"); println(gemp[n])
                end
                print_with_color(:red,"failed: analytic and empiricial gradient mismatch more than $tol\n")
            else
                print_with_color(:green,"passed\n")
            end
        end
    end
    if returngrad
        return gemp
    end
end

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
export ScaledHybridSystem

"""
    struct ScaledHybridSystem{T, H <: Union{DiscreteSwitchedLinearSystem, ConstrainedDiscreteSwitchedLinearSystem}} <: HybridSystems.AbstractHybridSystem
        system::H
        γ::T
    end

Discrete-time system where each reset map is scaled by `γ`, that is, the reset map `x ↦ Ax` of `system` is replaced by `x ↦ Ax/γ`.
"""
struct ScaledHybridSystem{T, H <: Union{DiscreteSwitchedLinearSystem, ConstrainedDiscreteSwitchedLinearSystem}} <: HybridSystems.AbstractHybridSystem
    system::H
    γ::T
end

# The hybrid system both acts like an automaton and a system

# Automaton
for f in [:state_property, :transition_property, :states, :nstates, :transitiontype, :transitions, :ntransitions, :source, :event, :target, :has_transition, :rem_transition!, :in_transitions, :out_transitions]
    @eval begin
        HybridSystems.$f(hs::ScaledHybridSystem, args...) = HybridSystems.$f(hs.system, args...)
    end
end
for f in [:state_property_type, :transition_property_type]
    @eval begin
        HybridSystems.$f(::Type{<:ScaledHybridSystem{T, H}}, args...) where {T, H} = HybridSystems.$f(H, args...)
    end
end

# System
MathematicalSystems.statedim(hs::ScaledHybridSystem, u) = MathematicalSystems.statedim(hs.system, u)
MathematicalSystems.stateset(hs::ScaledHybridSystem, u) = MathematicalSystems.stateset(hs.system, u)
function scale(lm::LinearMap, γ)
    return LinearMap(lm.A / γ)
end
scale(s::ContinuousIdentitySystem, γ) = s

function HybridSystems.mode(hs::ScaledHybridSystem, u)
    scale(HybridSystems.mode(hs.system, u), hs.γ)
end
function HybridSystems.resetmap(hs::ScaledHybridSystem, t)
    scale(HybridSystems.resetmap(hs.system, t), hs.γ)
end

# TODO should be able to remove it with SetProg
MultivariatePolynomials.variables(s::ScaledHybridSystem, state::Int) = MultivariatePolynomials.variables(s.system, state)


text:
import .ParallelKernel: get_body, set_body!, add_return, remove_return, substitute, literaltypes, push_to_signature!, add_loop, add_threadids, RANGES_VARNAME, RANGES_TYPE, RANGELENGTHS_VARNAMES

const PARALLEL_DOC = """
    @parallel kernel

Declare the `kernel` parallel and containing stencil computations be performed with one of the submodules `ParallelStencil.FiniteDifferences{1D|2D|3D}` (or with a compatible custom module or set of macros).

See also: [`@init_parallel_stencil`](@ref)

--------------------------------------------------------------------------------
$(replace(ParallelKernel.PARALLEL_DOC, "@init_parallel_kernel" => "@init_parallel_stencil"))
"""
@doc PARALLEL_DOC
macro parallel(args...) check_initialized(); checkargs_parallel(args...); esc(parallel(__module__, args...)); end

## MACROS FORCING PACKAGE, IGNORING INITIALIZATION

macro parallel_cuda(args...)    check_initialized(); checkargs_parallel(args...); esc(parallel(__module__, args...; package=PKG_CUDA)); end
macro parallel_threads(args...) check_initialized(); checkargs_parallel(args...); esc(parallel(__module__, args...; package=PKG_THREADS)); end

## ARGUMENT CHECKS

function checkargs_parallel(args...)
    if isempty(args) @ArgumentError("arguments missing.") end
    if is_function(args[end])  # Case: @parallel kernel
        if (length(args) != 1) @ArgumentError("wrong number of arguments in @parallel kernel call.") end
    elseif is_call(args[end])  # Case: @parallel <args...> kernelcall
        ParallelKernel.checkargs_parallel(args...)
    else
        @ArgumentError("the last argument must be a function definition or a kernel call (obtained: $(args[end])).")
    end
end

## GATEWAY FUNCTIONS

function parallel(caller::Module, args::Union{Symbol,Expr}...; package::Symbol=get_package())
    if is_function(args[end])
        numbertype = get_numbertype()
        ndims      = get_ndims()
        parallel_kernel(caller, package, numbertype, ndims, args...)
    elseif is_call(args[end])
        ParallelKernel.parallel(args...)
    end
end

## @PARALLEL KERNEL FUNCTIONS

function parallel_kernel(caller::Module, package::Symbol, numbertype::DataType, ndims::Integer, kernel::Expr)
    indices = get_indices_expr(ndims).args
    body = get_body(kernel)
    body = remove_return(body)
    validate_body(body)
    if (package == PKG_CUDA)
        kernel = substitute(kernel, :(Data.Array), :(Data.DeviceArray))
    end
    kernel = push_to_signature!(kernel, :($RANGES_VARNAME::$RANGES_TYPE))
    if (package == PKG_CUDA)
        kernel = push_to_signature!(kernel, :($(RANGELENGTHS_VARNAMES[1])::$INT_CUDA))
        kernel = push_to_signature!(kernel, :($(RANGELENGTHS_VARNAMES[2])::$INT_CUDA))
        kernel = push_to_signature!(kernel, :($(RANGELENGTHS_VARNAMES[3])::$INT_CUDA))
    end
    ranges = [:($RANGES_VARNAME[1]), :($RANGES_VARNAME[2]), :($RANGES_VARNAME[3])]
    if (package == PKG_CUDA)
        body = add_threadids(indices, ranges, body)
        body = literaltypes(numbertype, body)
    elseif (package == PKG_THREADS)
        body = add_loop(indices, ranges, body)
        body = literaltypes(numbertype, body)
    else
        @ArgumentError("$ERRMSG_UNSUPPORTED_PACKAGE (obtained: $package).")
    end
    check_mask_macro(caller)
    body = apply_masks(body, indices)
    body = add_return(body)
    set_body!(kernel, body)
    return kernel
end

## FUNCTIONS TO DEAL WITH MASKS (@WITHIN) AND INDICES

function check_mask_macro(caller::Module)
    if !isdefined(caller, Symbol("@within")) @MethodPluginError("the macro @within is not defined in the caller. You need to load one of the submodules ParallelStencil.FiniteDifferences{1|2|3}D (or a compatible custom module or set of macros).") end
    methods_str = string(methods(getfield(caller, Symbol("@within"))))
    if !occursin(r"(var\"@within\"|@within)\(__source__::LineNumberNode, __module__::Module, .*::String, .*::Symbol\)", methods_str) @MethodPluginError("the signature of the macro @within is not compatible with ParallelStencil (detected signature: \"$methods_str\"). The signature must correspond to the description in ParallelStencil.WITHIN_DOC. See in ParallelStencil.FiniteDifferences{1|2|3}D for examples.") end
end

function apply_masks(expr::Expr, indices::Array{Any})
    args = expr.args
    for i=1:length(args)
        if typeof(args[i]) == Expr
            e = args[i]
            if e.head == :(=) && typeof(e.args[1]) == Expr && e.args[1].head == :macrocall
                lefthand_macro = e.args[1].args[1]
                lefthand_var   = e.args[1].args[3]
                macroname = string(lefthand_macro)
                args[i] = quote
                              if (@within($macroname, $lefthand_var))
                                  $e
                              end
                          end
            else
                args[i] = apply_masks(e, indices)
            end
        end
    end
    return expr
end

function get_indices_expr(ndims::Integer)
    if ndims == 1
        return :($(INDICES[1]),)
    elseif ndims == 2
        return :($(INDICES[1]), $(INDICES[2]))
    elseif ndims == 3
        return :($(INDICES[1]), $(INDICES[2]), $(INDICES[3]))
    else
        @ModuleInternalError("argument 'ndims' must be 1, 2 or 3.")
    end
end


text:
import RowEchelon

#=
Homework 2
Section 1.4
Problem 4
=#

# Create the first Matrix
A = [1 3 -4 ;
     3 2 1]

# Create the second matrix
B = [1 ; 2 ; 1]

# Multiply the A and B
matrixProduct = A*B

# Print the results
println()
println("Original Matrix, A:")
println(A)
println("Original Matrix, B")
println(B)
println("Product Matrix, A*B:")
println(matrixProduct)
println()


text:
@testset "Mutations" begin

    @testset "GA" begin
        lx = [0.0, 0.0]
        ux = [2.0, 1.0]
        @test PM(lx, ux, Inf)([1.0, 2.0])[1] ∈ [0.0, 2.0]
        @test PM(lx, ux, Inf)([1.0, 2.0])[2] == 1.0
        @test MIPM(lx, ux, Inf, 1.0)(Real[1.0, 2])[1] ∈ [0.0, 2.0]
        @test MIPM(lx, ux, 10.0, Inf)(Real[1.0, 2])[2] == 1.0

        mut = uniform(3.0)
        @test all(mut(ones(3)) .!= 1)
        @test all(map(sum ∘ mut, [ones(3) for i in 1:10]) .!= 3)
        mut = uniform(0.0)
        @test all(map(sum ∘ mut, [ones(3) for i in 1:10]) .== 3)

        mut = gaussian(3.0)
        @test all(mut(ones(3)) .!= 1)
        @test all(map(sum ∘ mut, [ones(3) for i in 1:10]) .!= 3)
        mut = gaussian(0.0)
        @test all(map(sum ∘ mut, [ones(3) for i in 1:10]) .== 3)

    end

    @testset "DE" begin
        rec = ones(2)
        mut = [[0.5, 1.0], [1.0, 0.5]]
        @test_throws AssertionError Evolutionary.differentiation(rec, mut[1:1])
        @test Evolutionary.differentiation(rec, mut) == [0.5, 1.5]
    end

    @testset "GP" begin
        Random.seed!(9874984737482)
        H = 2
        tr = TreeGP(maxdepth=H, initialization=:grow);
        ex = rand(tr, H)

        # subtree mutation does not produce offspring expression longer then the parent
        mut = subtree(tr)
        off = [mut(copy(ex)) for i in 1:10]
        @testset "Offspring Height" for i in 1:10
            @test all(o->Evolutionary.height(o) <= H, map(mut, off))
        end

        # subtree mutation does not produce offspring expression longer then the parent
        mut = subtree(tr; growth=0.6)
        off = [mut(copy(ex)) for i in 1:10]
        @testset "Offspring Height" for i in 1:5
            @test all(o->Evolutionary.height(o) <= 3i, map(mut, off))
        end

        # hoist
        mut = subtree(tr; growth=0.6)
        off = [mut(copy(ex)) for i in 1:10]
        for i in 1:5; map(mut, off); end
        mut = hoist(tr)
        @testset "Hoist Mutation" for i in 1:10
            n = length(off[i])
            m = length(mut(off[i]))
            @test n > m
        end

        # shrink
        mut = subtree(tr; growth=0.6)
        off = [mut(copy(ex)) for i in 1:10]
        for i in 1:5; map(mut, off); end
        mut = shrink(tr)
        @testset "Shrink Mutation" for i in 1:10
            n = length(off[i])
            m = length(mut(off[i]))
            @test n >= m
        end

        # point
        mut = point(tr)
        @testset "Point Mutation" for i in 1:10
            mex = mut(copy(ex))
            @test length(ex) == length(mex) == 3
            @test sum(x == y for (x,y) in zip(ex.args, mex.args)) >= 2
        end

    end
end


text:
module Writer

import Dates

include("./status_codes.jl")

"""
    write(status_code::Int16, headers::Dict{String, String}; body::String="") :: String

Return HTTP response string
- Status line
- Headers
- Empty line
- Body (optional)
"""
function write(status_code::Int16, headers::Dict{String, String}; body::String="") :: String
    status_code_string::String = string(status_code) * " " * status_codes[status_code]
    response::String = "HTTP/1.1 " * status_code_string * "\n"

    if !haskey(headers, "Date")
        timestamp::Dates.DateTime = Dates.now(Dates.UTC)
        dayname::String = Dates.dayabbr(timestamp)
        day::Int64 = Dates.day(timestamp)
        monthname::String = Dates.monthabbr(timestamp)
        year::Int64 = Dates.year(timestamp)
        hour::Int64 = Dates.hour(timestamp)
        minute::Int64 = Dates.minute(timestamp)
        second::Int64 = Dates.second(timestamp)
        date_string::String = "$dayname, $day $monthname $year $hour:$minute:$second UTC"
        response *= "Date: " * date_string * "\n"
    end

    if !haskey(headers, "Content-Length")
        response *= "Content-Length: " * string(length(body)) * "\n"
    end

    for (key, val) in headers
        response *= "$key: $val\n"
    end

    # TODO: Only if Http 1.1
    # if Utils.has_header(headers, "Connection", "close")
    response *= "Connection: Closed"

    if body != ""
        response *= "\n\n" * body
    end

    return response
end

end