text:
function prepareLockModel(initialize, conf, counts, workloadName)
	if !isdefined(:oldConfigForPrepareLockModel)
		global oldConfigForPrepareLockModel = ""
	else
		global oldConfigForPrepareLockModel
	end

	if !isdefined(:oldWorkloadNameForPrepareLockModel)
		global oldWorkloadNameForPrepareLockModel = ""
	else
		global oldWorkloadNameForPrepareLockModel
	end

	if !isdefined(:JForPrepareLockModel)
		global JForPrepareLockModel = 1
	else
		global JForPrepareLockModel
	end

	if !isdefined(:KForPrepareLockModel)
		global KForPrepareLockModel = []
	else
		global KForPrepareLockModel
	end

	if !isdefined(:regionsForPrepareLockModel)
		global regionsForPrepareLockModel = 1
	else
		global regionsForPrepareLockModel
	end

	if !isdefined(:DForPrepareLockModel)
		global DForPrepareLockModel = []
	else
		global DForPrepareLockModel
	end

	if !isdefined(:SForPrepareLockModel)
		global SForPrepareLockModel = []
	else
		global SForPrepareLockModel
	end

	if !isdefined(:gForPrepareLockModel)
		global gForPrepareLockModel = []
	else
		global gForPrepareLockModel
	end

	if !isdefined(:readRowsForPrepareLockModel)
		global readRowsForPrepareLockModel = []
	else
		global readRowsForPrepareLockModel
	end

	if !isdefined(:updatedRowsForPrepareLockModel)
		global updatedRowsForPrepareLockModel = []
	else
		global updatedRowsForPrepareLockModel
	end

	#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	if workloadName == "LOCK1"
		beginCost = conf[1];
		commitCost = conf[2];
	elseif workloadName == "TPCC"
		beginCost = conf[1];
		interLockInterval = conf[2];
		DomainMultiplier =  conf[3]; #%0.0000000100;
		costMultiplier= conf[4];
	else
		error(string("Unknown workloadName: ", workloadName));
	end
	nano = 1000000000.0;
	#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	if initialize==1 || !isequal(conf, oldConfigForPrepareLockModel) || (workloadName != oldWorkloadNameForPrepareLockModel)
		oldConfigForPrepareLockModel = conf;
		oldWorkloadNameForPrepareLockModel = workloadName;
		if workloadName == "LOCK1"
			JForPrepareLockModel = 1; #%  of transaction classes, e.g. in TPC-C, J=5
			f = [1.0]; #% frequency matrix for each class
			KForPrepareLockModel = [2]; #% the number of locks requested in each transaction class
			regionsForPrepareLockModel = 1; #% number of database regions, perhaps should be equal to the number of tables?!
			DForPrepareLockModel = [1024]; #% number of data items in the i'th DB region
			maxK = Max(vec(KForPrepareLockModel));
			SForPrepareLockModel = zeros(convert(Int,JForPrepareLockModel),convert(Int,maxK)); #% S_jn is the processing of the n'th step of a transaction of class C_j
			S0 = beginCost; #% Barzan assumes that everything that precedes the first lock is constant time for all transactions
			SForPrepareLockModel[1,1] = 0.001;
			SForPrepareLockModel[1,2] = 0.001 + commitCost;

			gForPrepareLockModel = zeros(convert(Int,JForPrepareLockModel),convert(Int,maxK),convert(Int,regionsForPrepareLockModel));
			gForPrepareLockModel[1,:,1] = 1.0;
		elseif workloadName == "TPCC"
			JForPrepareLockModel = 5;
			f = [0.2 0.2 0.2 0.2 0.2];
			KForPrepareLockModel = [136 8 6 10 153]; #% info for transaction type 4 is not accurate!
			regionsForPrepareLockModel = 9;
			customer=1; district=2; history=3; item=4; new_order=5; oorder=6; order_line=7; stock=8; warehouse=9;
			DForPrepareLockModel = zeros(1, convert(Int,regionsForPrepareLockModel));
			DForPrepareLockModel[customer] = 960000; DForPrepareLockModel[district] = 320; DForPrepareLockModel[history] = 4390440; #% specs says 960000
			DForPrepareLockModel[item]=100000; DForPrepareLockModel[new_order]=42;#%or 0? or 32*9000
			DForPrepareLockModel[oorder] = 5255647; #% according to specs this should be 960000!
			DForPrepareLockModel[order_line]=32*300000;#%not sure about this one!
			DForPrepareLockModel[stock]=3200000; DForPrepareLockModel[warehouse]=32;

			maxK = Max(vec(KForPrepareLockModel));
			SForPrepareLockModel = zeros(convert(Int,JForPrepareLockModel),convert(Int,maxK));
			S0 = beginCost + dot(vec([2096  0.1*6403+0.9*3608 1753 1543 4170]), vec(f)) / nano;
			SForPrepareLockModel[1,1:136] = [interLockInterval 727780 609558 564452 604045 666491 574574 736047 568746 730061 563094 722544 563341 715844 565713 734241 562731 727079 570391 727549 562944 723694 569319 720982 571238 729174 repmat([interLockInterval], 1, 99) 2305666 repmat([interLockInterval], 1, 9) 5241200];
			SForPrepareLockModel[2,1:8] = 0.1*[589150 623830 590536 644205 946331 647169 797044 3018040] + 0.9*[586701 618534 588888 643966 950617 0 700665 3515563];
			SForPrepareLockModel[3,1:6] = [952166 77273 repmat([interLockInterval], 1, 3) 7834434];
			SForPrepareLockModel[4,1:10] = [640781 626863 621894 617058 617356 615745 613847 614119 615625 3111715];
			SForPrepareLockModel[5,1:153] = [591071 repmat([interLockInterval], 1, 151) 8500096];

			S0 = beginCost + (S0-beginCost) * costMultiplier;
			SForPrepareLockModel = (SForPrepareLockModel/nano) * costMultiplier;

			readRowsForPrepareLockModel = zeros(convert(Int,JForPrepareLockModel), convert(Int,regionsForPrepareLockModel)); #%b_ji= the probability that
			updatedRowsForPrepareLockModel = zeros(convert(Int,JForPrepareLockModel), convert(Int,regionsForPrepareLockModel));

			gForPrepareLockModel = zeros(convert(Int,JForPrepareLockModel),convert(Int,maxK),convert(Int,regionsForPrepareLockModel));
			gForPrepareLockModel[1,1,customer] = 1.0;
			gForPrepareLockModel[1,2,warehouse] = 1.0;
			gForPrepareLockModel[1,3,district] = 1.0;
			gForPrepareLockModel[1,4,new_order] = 1.0;
			gForPrepareLockModel[1,5,district] = 1.0;
			gForPrepareLockModel[1,6,oorder] = 1.0;
			gForPrepareLockModel[1,[7  9 11 13 15 17 19 21 23 25],item] = 1.0;
			gForPrepareLockModel[1,[8 10 12 14 16 18 20 22 24 26],stock] = 1.0;
			gForPrepareLockModel[1,27:126,order_line] = 1.0;
			gForPrepareLockModel[1,127:136,stock] = 1.0;
			readRowsForPrepareLockModel[1,customer] = 1;
			readRowsForPrepareLockModel[1,warehouse] = 1;
			updatedRowsForPrepareLockModel[1,district] = 2;
			readRowsForPrepareLockModel[1,new_order] = 1; #% insert
			readRowsForPrepareLockModel[1,oorder] = 1; #%insert
			readRowsForPrepareLockModel[1,item] = 10;
			updatedRowsForPrepareLockModel[1,stock] = 10; #% I only count them when they get the locks with "select for update"
			readRowsForPrepareLockModel[1,order_line] = 100; #% inserts
			gForPrepareLockModel[2,1,warehouse] = 1.0;
			gForPrepareLockModel[2,2,warehouse] = 1.0;
			gForPrepareLockModel[2,3,district] = 1.0;
			gForPrepareLockModel[2,4,district] = 1.0;
			gForPrepareLockModel[2,5,customer] = 1.0;
			gForPrepareLockModel[2,6,customer] = 1.0;
			gForPrepareLockModel[2,7,customer] = 1.0;
			gForPrepareLockModel[2,8,history] = 1.0; #% the trace says table 't' which I guess means 'history'
			readRowsForPrepareLockModel[2,warehouse] = 1;
			updatedRowsForPrepareLockModel[2,warehouse] = 1;
			readRowsForPrepareLockModel[2,district] = 1;
			updatedRowsForPrepareLockModel[2,district] = 1;
			readRowsForPrepareLockModel[2,customer] = 2;
			updatedRowsForPrepareLockModel[2,customer] = 1;
			readRowsForPrepareLockModel[2,history] = 1; #% insert
			gForPrepareLockModel[3,1,customer] = 1.0;
			gForPrepareLockModel[3,2,oorder] = 1.0;
			gForPrepareLockModel[3,3:6,order_line] = 1.0;
			readRowsForPrepareLockModel[3,customer] = 1;
			readRowsForPrepareLockModel[3,oorder] = 1;
			readRowsForPrepareLockModel[3,order_line] = 4;

			#% For transaction type 4:
			#% with probability (659+665) / (659+665+10774) = 0.11:
			#%    w_new_order_1 -> r_oorder_1 ->  w_oorder_1 -> w_order_line_5 ->  r_order_line_5 -> w_customer_1 ->
			#% with prob 1:
			#%    r_new_order_10
			tempA = 0.11;
			pBranch = 1-(1-tempA)^0.1;
			gForPrepareLockModel[4,1:10,new_order] = 1.0 - 3*pBranch;
			readRowsForPrepareLockModel[4,new_order] = 10; #% I have to yse number of rows as a probability instead of a simple count! this way I can also account for other types of transaction 4
			updatedRowsForPrepareLockModel[4, new_order] = tempA/(1.0 - 3*pBranch); #% to make sure in the end we update the new_order w prob 0.11 too!
			gForPrepareLockModel[4,1:10,oorder] = pBranch; #% this way the overall prob that tran type 4 accesses oorder adds up to tempA !
			readRowsForPrepareLockModel[4,oorder] = 1;
			updatedRowsForPrepareLockModel[4,oorder] = 1;
			gForPrepareLockModel[4,1:10,order_line] = pBranch;
			updatedRowsForPrepareLockModel[4,order_line] = 5;
			readRowsForPrepareLockModel[4,order_line] = 5;
			gForPrepareLockModel[4,1:10,customer] = pBranch;
			updatedRowsForPrepareLockModel[4,customer] = 1;

			#%w_new_order_1 w_oorder_1 w_order_line_5 w_customer_1

			gForPrepareLockModel[5,1,district] = 1.0;
			gForPrepareLockModel[5,2:152,order_line] = 1.0;
			gForPrepareLockModel[5,153,stock] = 1.0;
			readRowsForPrepareLockModel[5,district] = 1;
			readRowsForPrepareLockModel[5,order_line] = 151;
			readRowsForPrepareLockModel[5,stock] = 1;

		else
			error(string("Unknown workloadName: ", workloadName));
		end
	end #% end of initialization
	#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

	#% This is the main part!

	f = counts;
	tps = sum(f,2);
	if length(tps) == 1
		tps = tps[1]
	end

	#Julia: change from / to ./
	f = f ./ tps;

	S0 = beginCost + (dot(vec([2096  0.1*6403+0.9*3608 1753 1543 4170]), vec(f)) / nano) * costMultiplier;

	b = zeros(1, convert(Int,regionsForPrepareLockModel));
	for i=1:regionsForPrepareLockModel
		sm = (readRowsForPrepareLockModel[:,i]+updatedRowsForPrepareLockModel[:,i]);
		for j=1:length(sm)
			if sm[j] == 0
				sm[j] = 123456
			end
		end
		b[i] = dot(vec(f), vec((updatedRowsForPrepareLockModel[:,i] ./ sm)'));
	end

	for j=1:length(b)
		if b[j] == 0
			b[j] = 0.0000000100 #%DomainMultiplier;
		end
	end
	DForPrepareLockModel = DForPrepareLockModel ./ (1- (1-b.*b)); #% when we do have exclusive locks we should not shrink all the tables!

	DForPrepareLockModel = DForPrepareLockModel * DomainMultiplier;
	#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	tps_ = tps;
	J_ = JForPrepareLockModel;
	f_ = f;
	K_ = KForPrepareLockModel;
	regions_ = regionsForPrepareLockModel;
	D_ = DForPrepareLockModel;
	S0_ = S0;
	S_ = SForPrepareLockModel;
	g_ = gForPrepareLockModel;
	return tps_,J_,f_,K_,regions_,D_,S0_,S_,g_
end

max_stars_repo_name: barzan/dbseer

max_issues_repo_name: barzan/dbseer

max_forks_repo_name: barzan/dbseer


text:
# Score-Based Generative Modeling
#
# Score-Based Generative Modeling through Stochastic Differential Equations.
# Yang Song, Jascha Sohl-Dickstein, Diederik P. Kingma,
# Abhishek Kumar, Stefano Ermon, and Ben Poole.
# https://arxiv.org/pdf/2011.13456.pdf

using MLDatasets
using Flux
using Flux: @functor, chunk, params
using Flux.Data: DataLoader
using Parameters: @with_kw
using BSON
using CUDA
using Images
using Logging: with_logger
using ProgressMeter: Progress, next!
using TensorBoardLogger: TBLogger, tb_overwrite
using Random
using Statistics

"""
Projection of Gaussian Noise onto a time vector.

# Notes
This layer will help embed our random times onto the frequency domain. \n
W is not trainable and is sampled once upon construction - see assertions below.

# References
paper-  https://arxiv.org/abs/2006.10739
"""
function GaussianFourierProjection(embed_dim, scale)
    # Instantiate W once
    W = randn(Float32, embed_dim ÷ 2) .* scale
    # Return a function that always references the same W
    function GaussFourierProject(t)
        t_proj = t' .* W * Float32(2π)
        [sin.(t_proj); cos.(t_proj)]
    end
end

"""
Helper function that computes the *standard deviation* of 𝒫₀ₜ(𝘹(𝘵)|𝘹(0)).

# Notes
Derived from the Stochastic Differential Equation (SDE):    \n
                𝘥𝘹 = σᵗ𝘥𝘸,      𝘵 ∈ [0, 1]                   \n

We use properties of SDEs to analytically solve for the stddev
at time t conditioned on the data distribution. \n

We will be using this all over the codebase for computing our model's loss,
scaling our network output, and even sampling new images!
"""
marginal_prob_std(t, sigma=25.0f0) = sqrt.((sigma .^ (2t) .- 1.0f0) ./ 2.0f0 ./ log(sigma))

"""
Create a UNet architecture as a backbone to a diffusion model. \n

# Notes
Images stored in WHCN (width, height, channels, batch) order. \n
In our case, MNIST comes in as (28, 28, 1, batch). \n

# References
paper-  https://arxiv.org/abs/1505.04597
"""
struct UNet
    layers::NamedTuple
end

"""
User Facing API for UNet architecture.
"""
function UNet(channels=[32, 64, 128, 256], embed_dim=256, scale=30.0f0)
    return UNet((
        gaussfourierproj=GaussianFourierProjection(embed_dim, scale),
        linear=Dense(embed_dim, embed_dim, swish),
        # Encoding
        conv1=Conv((3, 3), 1 => channels[1], stride=1, bias=false),
        dense1=Dense(embed_dim, channels[1]),
        gnorm1=GroupNorm(channels[1], 4, swish),
        conv2=Conv((3, 3), channels[1] => channels[2], stride=2, bias=false),
        dense2=Dense(embed_dim, channels[2]),
        gnorm2=GroupNorm(channels[2], 32, swish),
        conv3=Conv((3, 3), channels[2] => channels[3], stride=2, bias=false),
        dense3=Dense(embed_dim, channels[3]),
        gnorm3=GroupNorm(channels[3], 32, swish),
        conv4=Conv((3, 3), channels[3] => channels[4], stride=2, bias=false),
        dense4=Dense(embed_dim, channels[4]),
        gnorm4=GroupNorm(channels[4], 32, swish),
        # Decoding
        tconv4=ConvTranspose((3, 3), channels[4] => channels[3], stride=2, bias=false),
        dense5=Dense(embed_dim, channels[3]),
        tgnorm4=GroupNorm(channels[3], 32, swish),
        tconv3=ConvTranspose((3, 3), channels[3] + channels[3] => channels[2], pad=(0, -1, 0, -1), stride=2, bias=false),
        dense6=Dense(embed_dim, channels[2]),
        tgnorm3=GroupNorm(channels[2], 32, swish),
        tconv2=ConvTranspose((3, 3), channels[2] + channels[2] => channels[1], pad=(0, -1, 0, -1), stride=2, bias=false),
        dense7=Dense(embed_dim, channels[1]),
        tgnorm2=GroupNorm(channels[1], 32, swish),
        tconv1=ConvTranspose((3, 3), channels[1] + channels[1] => 1, stride=1, bias=false),
    ))
end

@functor UNet

"""
Helper function that adds `dims` dimensions to the front of a `AbstractVecOrMat`.
Similar in spirit to TensorFlow's `expand_dims` function.

# References:
https://www.tensorflow.org/api_docs/python/tf/expand_dims
"""
expand_dims(x::AbstractVecOrMat, dims::Int=2) = reshape(x, (ntuple(i -> 1, dims)..., size(x)...))

"""
Makes the UNet struct callable and shows an example of a "Functional" API for modeling in Flux. \n
"""
function (unet::UNet)(x, t)
    # Embedding
    embed = unet.layers.gaussfourierproj(t)
    embed = unet.layers.linear(embed)
    # Encoder
    h1 = unet.layers.conv1(x)
    h1 = h1 .+ expand_dims(unet.layers.dense1(embed), 2)
    h1 = unet.layers.gnorm1(h1)
    h2 = unet.layers.conv2(h1)
    h2 = h2 .+ expand_dims(unet.layers.dense2(embed), 2)
    h2 = unet.layers.gnorm2(h2)
    h3 = unet.layers.conv3(h2)
    h3 = h3 .+ expand_dims(unet.layers.dense3(embed), 2)
    h3 = unet.layers.gnorm3(h3)
    h4 = unet.layers.conv4(h3)
    h4 = h4 .+ expand_dims(unet.layers.dense4(embed), 2)
    h4 = unet.layers.gnorm4(h4)
    # Decoder
    h = unet.layers.tconv4(h4)
    h = h .+ expand_dims(unet.layers.dense5(embed), 2)
    h = unet.layers.tgnorm4(h)
    h = unet.layers.tconv3(cat(h, h3; dims=3))
    h = h .+ expand_dims(unet.layers.dense6(embed), 2)
    h = unet.layers.tgnorm3(h)
    h = unet.layers.tconv2(cat(h, h2, dims=3))
    h = h .+ expand_dims(unet.layers.dense7(embed), 2)
    h = unet.layers.tgnorm2(h)
    h = unet.layers.tconv1(cat(h, h1, dims=3))
    # Scaling Factor
    h ./ expand_dims(marginal_prob_std(t), 3)
end

"""
Model loss following the denoising score matching objectives:

# Notes
Denoising score matching objective:
```julia
min wrt. θ (
    𝔼 wrt. 𝘵 ∼ 𝒰(0, 𝘛)[
        λ(𝘵) * 𝔼 wrt. 𝘹(0) ∼ 𝒫₀(𝘹) [
            𝔼 wrt. 𝘹(t) ∼ 𝒫₀ₜ(𝘹(𝘵)|𝘹(0)) [
                (||𝘚₀(𝘹(𝘵), 𝘵) - ∇ log [𝒫₀ₜ(𝘹(𝘵) | 𝘹(0))] ||₂)²
            ]
        ]
    ]
)
```
Where 𝒫₀ₜ(𝘹(𝘵) | 𝘹(0)) and λ(𝘵), are available analytically and
𝘚₀(𝘹(𝘵), 𝘵) is estimated by a U-Net architecture.

# References:
http://www.iro.umontreal.ca/~vincentp/Publications/smdae_techreport.pdf \n
https://yang-song.github.io/blog/2021/score/#estimating-the-reverse-sde-with-score-based-models-and-score-matching \n
https://yang-song.github.io/blog/2019/ssm/
"""
function model_loss(model, x, ϵ=1.0f-5)
    batch_size = size(x)[end]
    # (batch) of random times to approximate 𝔼[⋅] wrt. 𝘪 ∼ 𝒰(0, 𝘛)
    random_t = rand!(similar(x, batch_size)) .* (1.0f0 - ϵ) .+ ϵ
    # (batch) of perturbations to approximate 𝔼[⋅] wrt. 𝘹(0) ∼ 𝒫₀(𝘹)
    z = randn!(similar(x))
    std = expand_dims(marginal_prob_std(random_t), 3)
    # (batch) of perturbed 𝘹(𝘵)'s to approximate 𝔼 wrt. 𝘹(t) ∼ 𝒫₀ₜ(𝘹(𝘵)|𝘹(0))
    perturbed_x = x + z .* std
    # 𝘚₀(𝘹(𝘵), 𝘵)
    score = model(perturbed_x, random_t)
    # mean over batches
    mean(
        # L₂ norm over WHC dimensions
        sum((score .* std + z) .^ 2; dims=1:(ndims(x) - 1))
    )
end

"""
Helper function that loads MNIST images and returns loader.
"""
function get_data(batch_size)
    xtrain, ytrain = MLDatasets.MNIST(:train)[:]
    xtrain = reshape(xtrain, 28, 28, 1, :)
    DataLoader((xtrain, ytrain), batchsize=batch_size, shuffle=true)
end

"""
Helper function from DrWatson.jl to convert a struct to a dict
"""
function struct2dict(::Type{DT}, s) where {DT<:AbstractDict}
    DT(x => getfield(s, x) for x in fieldnames(typeof(s)))
end
struct2dict(s) = struct2dict(Dict, s)

# arguments for the `train` function
@with_kw mutable struct Args
    η = 1e-4                                        # learning rate
    batch_size = 32                                 # batch size
    epochs = 50                                     # number of epochs
    seed = 1                                        # random seed
    cuda = false                                    # use CPU
    verbose_freq = 10                               # logging for every verbose_freq iterations
    tblogger = true                                 # log training with tensorboard
    save_path = "output"                            # results path
end

function train(; kws...)
    # load hyperparamters
    args = Args(; kws...)
    args.seed > 0 && Random.seed!(args.seed)

    # GPU config
    if args.cuda && CUDA.has_cuda()
        device = gpu
        @info "Training on GPU"
    else
        device = cpu
        @info "Training on CPU"
    end

    # load MNIST images
    loader = get_data(args.batch_size)

    # initialize UNet model
    unet = UNet() |> device

    # ADAM optimizer
    opt = ADAM(args.η)

    # parameters
    ps = Flux.params(unet)

    !ispath(args.save_path) && mkpath(args.save_path)

    # logging by TensorBoard.jl
    if args.tblogger
        tblogger = TBLogger(args.save_path, tb_overwrite)
    end

    # Training
    train_steps = 0
    @info "Start Training, total $(args.epochs) epochs"
    for epoch = 1:args.epochs
        @info "Epoch $(epoch)"
        progress = Progress(length(loader))

        for (x, _) in loader
            x = device(x)
            loss, grad = Flux.withgradient(ps) do
                model_loss(unet, x)
            end
            Flux.Optimise.update!(opt, ps, grad)
            # progress meter
            next!(progress; showvalues=[(:loss, loss)])

            # logging with TensorBoard
            if args.tblogger && train_steps % args.verbose_freq == 0
                with_logger(tblogger) do
                    @info "train" loss = loss
                end
            end
            train_steps += 1
        end
    end

    # save model
    model_path = joinpath(args.save_path, "model.bson")
    let unet = cpu(unet), args = struct2dict(args)
        BSON.@save model_path unet args
        @info "Model saved: $(model_path)"
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    train()
end

ext: jl


text:
@testset "Ideals" begin
   Qx, x = PolynomialRing(FlintQQ, "x")

   K1, a1 = NumberField(x^3 - 2, "a")
   O1 = Order(K1, Hecke.FakeFmpqMat(FlintZZ[1 0 0; 0 2 0; 0 0 4], one(FlintZZ)))

   K2, a2 = NumberField(x^2 - 4^2*7^2*5, "a")
   O2 = Order(K2, [K2(1), a2])

  @testset "Construction" begin
    I = @inferred ideal(O1, -17)
    @test order(I) == O1
    @test I.princ_gen_special[1] == 1
    @test I.princ_gen_special[2] == 17
    @test I.princ_gen == O1(-17)
    @test basis_matrix(I) == MatrixSpace(FlintZZ, 3, 3)(17)

    J = @inferred ideal(O1, FlintZZ(-17))
    @test order(J) == O1
    @test J.princ_gen_special[1] == 2
    @test J.princ_gen_special[3] == FlintZZ(17)
    @test J.princ_gen == O1(-17)
    @test basis_matrix(J) == MatrixSpace(FlintZZ, 3, 3)(17)

    K = @inferred ideal(O1, O1(-17))
    @test K.princ_gen == O1(-17)
    @test basis_matrix(K) == MatrixSpace(FlintZZ, 3, 3)(17)

    M = @inferred O1(-17)*O1
    L = @inferred O1*O1(-17)

    @test I == J && J == K && K == M && M == L && L == I

    I2 = @inferred ideal(O2, O2(1 + a2))

    Ib = basis(I2)
    II = ideal(O2, Ib)
    @test I2 == II
  end

  I = ideal(O1, -17)
  J = ideal(O1, FlintZZ(-17))
  K = ideal(O1, O1(-17))
  M = O1(-17)*O1
  I2 = ideal(O2, O2(1 + a2))

  @testset "Minimum" begin
    @test 17 == @inferred minimum(I)
    @test 17 == @inferred minimum(J)
    @test 17 == @inferred minimum(K)
    @test 17 == @inferred minimum(I)
    @test 17 == @inferred minimum(J)
    @test 17 == @inferred minimum(K)

    # Test where gens are weakly normal and second generator is zero
    @testset begin
      R, x = PolynomialRing(FlintQQ, "x")
      _K, _a = NumberField(x, "a")
      _O = maximal_order(_K)
      _I = fractional_ideal(_O, _K(1))
      _J = _I*_K(fmpq(-1, 5))
      @test minimum(numerator(_J)) == 1
    end
  end

  @testset "Norm" begin
    @test 4913 == @inferred norm(I)
    @test 4913 == @inferred norm(J)
    @test 4913 == @inferred norm(K)
    @test 4913 == @inferred norm(I)
    @test 4913 == @inferred norm(J)
    @test 4913 == @inferred norm(K)
  end

  @testset "Deepcopy" begin
    L = @inferred deepcopy(I)
    @test L == I
    @test order(L) == O1
  end

  @testset "Basis" begin
    b = @inferred basis(I)
    @test b == NfOrdElem[ O1(17), O1(34*a1), O1(68*a1^2) ]
  end

  @testset "Basismatrix" begin
    M = @inferred ideal(O1, O1(4*a1^2))

    b = @inferred basis_matrix(M)
    @test b == FlintZZ[16 0 0; 0 16 0; 0 0 1]

    b = @inferred basis_mat_inv(M)
    @test b == Hecke.FakeFmpqMat(FlintZZ[1 0 0; 0 1 0; 0 0 16], FlintZZ(16))

    b = @inferred basis_matrix(M)
    @test b == FlintZZ[16 0 0; 0 16 0; 0 0 1]

    b = @inferred basis_mat_inv(M)
    @test b == Hecke.FakeFmpqMat(FlintZZ[1 0 0; 0 1 0; 0 0 16], FlintZZ(16))
  end

  @testset "Inclusion" begin
    M = @inferred ideal(O1, O1(4*a1^2))
    @test @inferred in(O1(4*a1^2), M)
    @test @inferred !in(O1(2*a1), M)
  end

  @testset "Binary operations" begin
    K = @inferred I + J # I == J
    @test K == I
    K = @inferred ideal(O1, 15) + ideal(O1, 10)
    @test K == ideal(O1, 5)

    K = @inferred intersect(ideal(O1, 15), ideal(O1, 10))
    @test K == ideal(O1, 30)
    K = @inferred lcm(ideal(O1, 15), ideal(O1, 10))
    @test K == ideal(O1, 30)

    K = @inferred ideal(O1, 15) * ideal(O1, 10)
    @test K == ideal(O1, 150)
    K = @inferred ideal(O1, O1(2*a1)) * ideal(O1, O1(4*a1^2))
    @test K == ideal(O1, 16)
  end

  @testset "Ad hoc binary operations" begin
    I = ideal(O1, O1(2*a1))
    J = ideal(O1, 3)
    K = ideal(O1, fmpz(3))

    @test ideal(O1, O1(10*a1)) == @inferred 5 * I
    @test ideal(O1, O1(10*a1)) == @inferred I * 5

    @test ideal(O1, O1(10*a1)) == @inferred fmpz(5) * I
    @test ideal(O1, O1(10*a1)) == @inferred I * fmpz(5)

    @test ideal(O1, 30) == @inferred 10 * J
    @test ideal(O1, 30) == @inferred J * 10

    @test ideal(O1, 30) == @inferred fmpz(10) * J
    @test ideal(O1, 30) == @inferred J * fmpz(10)

    @test ideal(O1, 30) == @inferred 10 * K
    @test ideal(O1, 30) == @inferred K * 10

    @test ideal(O1, 30) == @inferred fmpz(10) * K
    @test ideal(O1, 30) == @inferred K * fmpz(10)

    @test ideal(O1, O1(4*a1^2)) == @inferred I * O1(2*a1)
    @test ideal(O1, O1(4*a1^2)) == @inferred O1(2*a1) * I
  end

  @testset "Idempotents" begin
    I = ideal(O1, 2)
    J = ideal(O1, 3)
    e, f = @inferred idempotents(I, J)
    @test in(e, I)
    @test in(f, J)
    @test e + f == 1

    I = ideal(O1, O1(2*a1))
    J = ideal(O1, O1(1 - 2*a1))
    e, f = @inferred idempotents(I, J)
    @test in(e, I)
    @test in(f, J)
    @test e + f == 1

    @test_throws ErrorException idempotents(I, I)
  end

  @testset "Modular reduction" begin
    I = ideal(O1, 10)
    a = O1([11, 11, 11])
    b = @inferred mod(a, I)
    @test b == O1([1, 1, 1])

    I = ideal(O1, O1(2*a1 + 4*a1^2))
    a = O1([544, 1, 1])
    b = @inferred mod(a, I)
    @test iszero(b)

    c = [ Hecke.fmpz_preinvn_struct(fmpz(544)),
          Hecke.fmpz_preinvn_struct(fmpz(1)),
          Hecke.fmpz_preinvn_struct(fmpz(1)) ]

    b = @inferred mod(a, I, c)
    @test iszero(b)
  end

  @testset "p-Radical" begin
    I = @inferred pradical(O1, 2)
    @test I == ideal(O1, FlintZZ[2 0 0; 0 1 0; 0 0 1])
  end

  @testset "Prime Decomposition" begin
    L = NumberField(x^30-x^29+x^28-x^27+x^26+743*x^25-1363*x^24-3597*x^23-22009*x^22+458737*x^21+2608403*x^20+6374653*x^19-1890565*x^18-112632611*x^17-467834081*x^16-1365580319*x^15-1188283908*x^14+3831279180*x^13+28661663584*x^12+89106335984*x^11+226912479680*x^10+443487548480*x^9+719797891328*x^8+946994403328*x^7+1015828094976*x^6+878645952512*x^5+555353440256*x^4+124983967744*x^3+67515711488*x^2-5234491392*x+400505700352)[1]
    OL = maximal_order(L)
    @test length(prime_decomposition(OL, 2)) == 30
    Lns, gLns = number_field([x^2-5, x^2-13])
    OLns = maximal_order(Lns)
    lP = prime_decomposition(OLns, 5)
    @test length(lP) == 1
  end

  @testset "Frobenius automorphism" begin
    K = number_field(x^2+1)[1]
    OK = maximal_order(K)
    lp = prime_decomposition(OK, 5)
    P = lp[1][1]
    @test Hecke.frobenius_automorphism(P) == id_hom(K)
    lp = prime_decomposition(OK, 7)
    P = lp[1][1]
    @test Hecke.frobenius_automorphism(P) != id_hom(K)
  end

  include("Ideal/Prime.jl")
end


text:
mutable struct DMDc{K, B, Q, P} <: abstractKoopmanOperator
    koopman::K
    forcing::B

    Qₖ::Q # TODO implement for update
    Pₖ::P
end

function DMDc(X::AbstractArray, Y::AbstractArray, U::AbstractArray; B::AbstractArray = [], dt::Float64 = 0.0)
    @assert all(size(X) .== size(Y))
    @assert size(X)[2] .== size(U)[2]

    nₓ = size(X)[1]
    nᵤ = size(U)[1]

    if isempty(B)
        Ω = vcat(X, U)
        G = Y * pinv(Ω)

        Ã = G[:, 1:nₓ]
        B̃ = G[:, nₓ+1:end]
    else
        Ã = (Y - B*U)*pinv(X)
        B̃ = B
    end

    # Eigen Decomposition for solution
    Λ, W = eigen(Ã)

    # Right now now really useful, since we do not know
    # how to deal with B
    if dt > 0.0
        # Casting Complex enforces results
        ω = log.(Complex.(Λ)) / dt
    else
        ω = []
    end

    koopman = ExactDMD(Ã, Λ, ω, W, nothing, nothing)
    # TODO implement update method
    return DMDc(koopman, B̃, nothing, nothing)
end

function DMDc(X::AbstractArray, U::AbstractArray; B::AbstractArray = [], dt::Float64 = 0.0)
    @assert size(X)[2]-1 == size(U)[2] "Provide consistent input data."
    return DMDc(X[:, 1:end-1], X[:, 2:end], U, B = B, dt = dt)
end

# Some nice functions
LinearAlgebra.eigen(m::DMDc) = eigen(m.koopman)
LinearAlgebra.eigvals(m::DMDc) = eigvals(m.koopman)
LinearAlgebra.eigvecs(m::DMDc) = eigvecs(m.koopman)

DataDrivenDiffEq.isstable(m::DMDc) = isstable(m.koopman)

get_dynamics(m::DMDc) = m.koopman
get_input_map(m::DMDc) = m.forcing

# TODO this can be done better, maybe use macros
function dynamics(m::DMDc; control = nothing)
    control = isnothing(control) ? zero_callback(m) : control
    @inline function dudt_(u, p, t; y = control)
        m.koopman.Ã * u + m.forcing * y(u, p, t)
    end
    return dudt_
end

function zero_callback(m::DMDc)
    return length(size(m.forcing)) <= 1 ? (u, p, t) -> zero(eltype(m.forcing)) : (u, p, t) -> zeros(eltype(m.forcing), size(m.forcing)[2])
end


text:
using Test
using DelimitedFiles, LinearAlgebra, Statistics #, MLDatasets

using StableRNGs
#rng = StableRNG(123)
using BetaML

TESTRNG = FIXEDRNG # This could change...
#TESTRNG = StableRNG(123)

println("*** Testing Neural Network...")

# ==================================
# New test
# ==================================
println("Testing Learnable structure...")

L1a = Learnable(([1.0 2; 3 4.0], [1.1,2.2], Float64[], [1.1,2.2], [1.1]))
L1b = Learnable(([1.0 2; 3 4.0], [1.1,2.2], Float64[], [1.1,2.2], [1.1]))
L1c = Learnable(([1.0 2; 3 4.0], [1.1,2.2], Float64[], [1.1,2.2], [1.1]))
foo = (((((sum(L1a,L1b,L1c) - (L1a + L1b + L1c)) + 10) - 4) * 2 ) / 3) *
      (((((sum(L1a,L1b,L1c) - (L1a + L1b + L1c)) + 10) - 4) * 2 ) / 3)

@test foo.data[1][2] == 16.0 && foo.data[5][1] == 16.00
@test (L1a -1).data[1] ≈ (-1 * (1 - L1a)).data[1]
@test (2 / (L1a / 2)).data[2] ≈ (4/L1a).data[2]
@test sqrt(L1a).data[1][2,2] == 2.0

# ==================================
# New test
# ==================================
println("Testing basic NN behaviour...")
x = [0.1,1]
y = [1,0]
l1   = DenseNoBiasLayer(2,2,w=[2 1;1 1],f=identity,rng=copy(TESTRNG))
l2   = VectorFunctionLayer(2,f=softmax)
mynn = buildNetwork([l1,l2],squaredCost,name="Simple Multinomial logistic regression")
o1 = forward(l1,x)
@test o1 == [1.2,1.1]
#@code_warntype forward(l1,x)
o2 = forward(l2,o1)
@test o2 ≈ [0.5249791874789399, 0.47502081252106] ≈ softmax([1.2,1.1])
orig = predict(mynn,x')[1,:]
@test orig == o2
#@code_warntype Nn.predict(mynn,x')[1,:]
ϵ = squaredCost(o2,y)
#@code_warntype  squaredCost(o2,y)
lossOrig = loss(mynn,x',y')
@test ϵ == lossOrig
#@code_warntype  loss(mynn,x',y')
dϵ_do2 = dSquaredCost(o2,y)
@test dϵ_do2 == [-0.4750208125210601,0.47502081252106]
#@code_warntype dSquaredCost(o2,y)
dϵ_do1 = backward(l2,o1,dϵ_do2)
@test dϵ_do1 ≈ [-0.23691761847142412, 0.23691761847142412]
#@code_warntype backward(l2,o1,dϵ_do2)
dϵ_dX = backward(l1,x,dϵ_do1)
@test dϵ_dX ≈ [-0.23691761847142412, 0.0]
#@code_warntype backward(l1,x,dϵ_do1)
l1w = getParams(l1)
#@code_warntype
l2w = getParams(l2)
#@code_warntype
w = getParams(mynn)
#@code_warntype getParams(mynn)
#typeof(w)2-element Vector{Float64}:
origW = deepcopy(w)
l2dw = getGradient(l2,o1,dϵ_do2)
@test length(l2dw.data) == 0
#@code_warntype getGradient(l2,o1,dϵ_do2)
l1dw = getGradient(l1,x,dϵ_do1)
@test l1dw.data[1] ≈ [-0.023691761847142414 -0.23691761847142412; 0.023691761847142414 0.23691761847142412]
#@code_warntype
dw = getGradient(mynn,x,y)
#@code_warntype getGradient(mynn,x,y)
y_deltax1 = predict(mynn,[x[1]+0.001 x[2]])[1,:]
#@code_warntype Nn.predict(mynn,[x[1]+0.001 x[2]])[1,:]
lossDeltax1 = loss(mynn,[x[1]+0.001 x[2]],y')
#@code_warntype
deltaloss = dot(dϵ_dX,[0.001,0])
#@code_warntype
@test isapprox(lossDeltax1-lossOrig,deltaloss,atol=0.0000001)
l1wNew = l1w
l1wNew.data[1][1,1] += 0.001
setParams!(l1,l1wNew)
lossDelPar = loss(mynn,x',y')
#@code_warntype
deltaLossPar = 0.001*l1dw.data[1][1,1]
lossDelPar - lossOrig
@test isapprox(lossDelPar - lossOrig,deltaLossPar,atol=0.00000001)
η = 0.01
#w = gradientDescentSingleUpdate(w,dw,η)
#w = w - dw * η
w = w - dw * η
#@code_warntype gradSub.(w, gradMul.(dw,η))
#@code_warntype
setParams!(mynn,w)
loss2 = loss(mynn,x',y')
#@code_warntype
@test loss2 < lossOrig
for i in 1:10000
    local w  = getParams(mynn)
    local dw = getGradient(mynn,x,y)
    w  = w - dw * η
    setParams!(mynn,w)
end
lossFinal = loss(mynn,x',y')
@test predict(mynn,x')[1,1]>0.96
setParams!(mynn,origW)
train!(mynn,x',y',epochs=10000,batchSize=1,sequential=true,verbosity=NONE,optAlg=SGD(η=t->η,λ=1),rng=copy(TESTRNG))
#@code_warntype train!(mynn,x',y',epochs=10000,batchSize=1,sequential=true,verbosity=NONE,optAlg=SGD(η=t->η,λ=1))
lossTraining = loss(mynn,x',y')
#@code_warntype
@test isapprox(lossFinal,lossTraining,atol=0.00001)

li   = DenseLayer(2,2,w=[2 1;1 1],f=identity,rng=copy(TESTRNG))
@test getNParams(li) == 6

# ==================================
# NEW Test
println("Testing regression if it just works with manual derivatives...")

xtrain = [0.1 0.2; 0.3 0.5; 0.4 0.1; 0.5 0.4; 0.7 0.9; 0.2 0.1]
ytrain = [0.3; 0.8; 0.5; 0.9; 1.6; 0.3]
xtest = [0.5 0.6; 0.14 0.2; 0.3 0.7; 2.0 4.0]
ytest = [1.1; 0.36; 1.0; 6.0]
l1 = DenseLayer(2,3,w=[1 1; 1 1; 1 1], wb=[0 0 0], f=tanh, df=dtanh,rng=copy(TESTRNG))
l2 = DenseNoBiasLayer(3,2, w=[1 1 1; 1 1 1], f=relu, df=drelu,rng=copy(TESTRNG))
l3 = DenseLayer(2,1, w=[1 1], wb=[0], f=identity,df=didentity,rng=copy(TESTRNG))
mynn = buildNetwork([l1,l2,l3],squaredCost,name="Feed-forward Neural Network Model 1",dcf=dSquaredCost)
train!(mynn,xtrain,ytrain,batchSize=1,sequential=true,epochs=100,verbosity=NONE,optAlg=SGD(η=t -> 1/(1+t),λ=1),rng=copy(TESTRNG))
#@benchmark train!(mynn,xtrain,ytrain,batchSize=1,sequential=true,epochs=100,verbosity=NONE,optAlg=SGD(η=t -> 1/(1+t),λ=1))
avgLoss = loss(mynn,xtest,ytest)
@test  avgLoss ≈ 1.599729991966362
expectedOutput = [0.7360644412052633, 0.7360644412052633, 0.7360644412052633, 2.47093434438514]
predicted = dropdims(predict(mynn,xtest),dims=2)
@test any(isapprox(expectedOutput,predicted))

# With the ADAM optimizer...
l1 = DenseLayer(2,3,w=[1 1; 1 1; 1 1], wb=[0 0 0], f=tanh, df=dtanh,rng=copy(TESTRNG))
l2 = DenseNoBiasLayer(3,2, w=[1 1 1; 1 1 1], f=relu, df=drelu,rng=copy(TESTRNG))
l3 = DenseLayer(2,1, w=[1 1], wb=[0], f=identity,df=didentity,rng=copy(TESTRNG))
mynn = buildNetwork([l1,l2,l3],squaredCost,name="Feed-forward Neural Network with ADAM",dcf=dSquaredCost)
train!(mynn,xtrain,ytrain,batchSize=1,sequential=true,epochs=100,verbosity=NONE,optAlg=ADAM(),rng=copy(TESTRNG))
avgLoss = loss(mynn,xtest,ytest)
@test  avgLoss ≈ 0.9497779759064725
expectedOutput = [1.7020525792404175, -0.1074729043392682, 1.4998367847079956, 3.3985794704732717]
predicted = dropdims(predict(mynn,xtest),dims=2)
@test any(isapprox(expectedOutput,predicted))

# ==================================
# NEW TEST
# ==================================
println("Testing using AD...")

ϵtrain = [1.023,1.08,0.961,0.919,0.933,0.993,1.011,0.923,1.084,1.037,1.012]
ϵtest  = [1.056,0.902,0.998,0.977]
xtrain = [0.1 0.2; 0.3 0.5; 0.4 0.1; 0.5 0.4; 0.7 0.9; 0.2 0.1; 0.4 0.2; 0.3 0.3; 0.6 0.9; 0.3 0.4; 0.9 0.8]
ytrain = [(0.1*x[1]+0.2*x[2]+0.3)*ϵtrain[i] for (i,x) in enumerate(eachrow(xtrain))]
xtest  = [0.5 0.6; 0.14 0.2; 0.3 0.7; 20.0 40.0;]
ytest  = [(0.1*x[1]+0.2*x[2]+0.3)*ϵtest[i] for (i,x) in enumerate(eachrow(xtest))]

l1   = DenseLayer(2,3,w=ones(3,2), wb=zeros(3),rng=copy(TESTRNG))
l2   = DenseLayer(3,1, w=ones(1,3), wb=zeros(1),rng=copy(TESTRNG))
mynn = buildNetwork([l1,l2],squaredCost,name="Feed-forward Neural Network Model 1")
train!(mynn,xtrain,ytrain,epochs=1000,sequential=true,batchSize=1,verbosity=NONE,optAlg=SGD(η=t->0.01,λ=1),rng=copy(TESTRNG))
avgLoss = loss(mynn,xtest,ytest)
@test  avgLoss ≈ 0.0032018998005211886
ŷtestExpected = [0.4676699631752518,0.3448383593117405,0.4500863419692639,9.908883999376018]
ŷtrain = dropdims(predict(mynn,xtrain),dims=2)
ŷtest = dropdims(predict(mynn,xtest),dims=2)
@test any(isapprox(ŷtest,ŷtestExpected))
mreTrain = meanRelError(ŷtrain,ytrain)
@test mreTrain <= 0.06
mreTest  = meanRelError(ŷtest,ytest)
@test mreTest <= 0.05

#predicted = dropdims(Nn.predict(mynn,xtrain),dims=2)
#ytrain

# ==================================
# NEW TEST
# ==================================
println("Going through Multinomial logistic regression (using softmax)...")
#=
using RDatasets
using Random
using DataFrames: DataFrame
using CSV
Random.seed!(123);
iris = dataset("datasets", "iris")
iris = iris[shuffle(axes(iris, 1)), :]
CSV.write(joinpath(@__DIR__,"data","iris_shuffled.csv"),iris)
=#

iris     = readdlm(joinpath(@__DIR__,"data","iris_shuffled.csv"),',',skipstart=1)
x = convert(Array{Float64,2}, iris[:,1:4])
y = map(x->Dict("setosa" => 1, "versicolor" => 2, "virginica" =>3)[x],iris[:, 5])
y_oh = oneHotEncoder(y)

ntrain = Int64(round(size(x,1)*0.8))
xtrain = x[1:ntrain,:]
ytrain = y[1:ntrain]
ytrain_oh = y_oh[1:ntrain,:]
xtest = x[ntrain+1:end,:]
ytest = y[ntrain+1:end]

xtrain = scale(xtrain)
#pcaOut = pca(xtrain,error=0.01)
#(xtrain,P) = pcaOut.X, pcaOut.P
xtest = scale(xtest)
#xtest = xtest*P

l1   = DenseLayer(4,10, w=ones(10,4), wb=zeros(10),f=celu,rng=copy(TESTRNG))
l2   = DenseLayer(10,3, w=ones(3,10), wb=zeros(3),rng=copy(TESTRNG))
l3   = VectorFunctionLayer(3,f=softmax)
mynn = buildNetwork([l1,l2,l3],squaredCost,name="Multinomial logistic regression Model Sepal")
train!(mynn,xtrain,ytrain_oh,epochs=254,batchSize=8,sequential=true,verbosity=NONE,optAlg=SGD(η=t->0.001,λ=1),rng=copy(TESTRNG))
ŷtrain = predict(mynn,xtrain)
ŷtest  = predict(mynn,xtest)
trainAccuracy = accuracy(ŷtrain,ytrain,tol=1)
testAccuracy  = accuracy(ŷtest,ytest,tol=1)
@test testAccuracy >= 0.8 # set to random initialisation/training to have much better accuracy

# With ADAM
l1   = DenseLayer(4,10, w=ones(10,4), wb=zeros(10),f=celu,rng=copy(TESTRNG))
l2   = DenseLayer(10,3, w=ones(3,10), wb=zeros(3),rng=copy(TESTRNG))
l3   = VectorFunctionLayer(3,f=softmax)
mynn = buildNetwork([l1,l2,l3],squaredCost,name="Multinomial logistic regression Model Sepal")
train!(mynn,xtrain,ytrain_oh,epochs=10,batchSize=8,sequential=true,verbosity=NONE,optAlg=ADAM(η=t -> 1/(1+t), λ=0.5),rng=copy(TESTRNG))
ŷtrain = predict(mynn,xtrain)
ŷtest  = predict(mynn,xtest)
trainAccuracy = accuracy(ŷtrain,ytrain,tol=1)
testAccuracy  = accuracy(ŷtest,ytest,tol=1)
@test testAccuracy >= 1

#=
if "all" in ARGS
    # ==================================
    # NEW TEST
    # ==================================
    println("Testing colvolution layer with MINST data...")
    train_x, train_y = MNIST.traindata()
    test_x,  test_y  = MNIST.testdata()

    test = train_x[:,:,1]

    eltype(test)

    test .+ 1
end
=#

# ==================================
# NEW Test
if VERSION >= v"1.6"
    println("Testing VectorFunctionLayer with pool1d function...")
    println("Attention: this test requires at least Julia 1.6")

    x        = rand(copy(TESTRNG),300,5)
    y        = [norm(r[1:3])+2*norm(r[4:5],2) for r in eachrow(x) ]
    (N,D)    = size(x)
    l1       = DenseLayer(D,8, f=relu,rng=copy(TESTRNG))
    l2       = VectorFunctionLayer(size(l1)[2],f=(x->pool1d(x,2,f=mean)))
    l3       = DenseLayer(size(l2)[2],1,f=relu, rng=copy(TESTRNG))
    mynn     = buildNetwork([l1,l2,l3],squaredCost,name="Regression with a pooled layer")
    train!(mynn,x,y,epochs=50,verbosity=STD,rng=copy(TESTRNG))
    ŷ        = predict(mynn,x)
    mreTrain = meanRelError(ŷ,y,normRec=false)
    @test mreTrain  < 0.14
end


text:
@enum(CUjit_option, MAX_REGISTERS = Cint(0),
                    THREADS_PER_BLOCK,
                    WALL_TIME,
                    INFO_LOG_BUFFER,
                    INFO_LOG_BUFFER_SIZE_BYTES,
                    ERROR_LOG_BUFFER,
                    ERROR_LOG_BUFFER_SIZE_BYTES,
                    OPTIMIZATION_LEVEL,
                    TARGET_FROM_CUCONTEXT,
                    TARGET,
                    FALLBACK_STRATEGY,
                    GENERATE_DEBUG_INFO,
                    LOG_VERBOSE,
                    GENERATE_LINE_INFO,
                    CACHE_MODE)

@enum(CUjit_input, CUBIN = Cint(0),
                   PTX,
                   FATBINARY,
                   OBJECT,
                   LIBRARY)

function convert_bits{T}(::Type{T}, data::UInt)
    if sizeof(data) == sizeof(T)
        return reinterpret(T, data)
    elseif sizeof(data) == 8 && sizeof(T) == 4
        relevant = UInt32(data & 0x00000000ffffffff)
        return reinterpret(T, relevant)
    else
        error("don't know how to extract $(sizeof(T)) bytes out of $(sizeof(data))")
    end
end

function encode(options::Dict{CUjit_option,Any})
    keys = Array(CUjit_option, 0)
    vals = Array(Ptr{Void}, 0)

    for (opt, val) in options
        push!(keys, opt)
        if opt == GENERATE_LINE_INFO ||
           opt == GENERATE_DEBUG_INFO ||
           opt == LOG_VERBOSE
            push!(vals, convert(Ptr{Void}, convert(Int, val::Bool)))
        elseif opt == INFO_LOG_BUFFER
            buf = val::Vector{UInt8}
            push!(vals, pointer(buf))
            push!(keys, INFO_LOG_BUFFER_SIZE_BYTES)
            push!(vals, sizeof(buf))
        elseif opt == ERROR_LOG_BUFFER
            buf = val::Vector{UInt8}
            push!(vals, pointer(buf))
            push!(keys, ERROR_LOG_BUFFER_SIZE_BYTES)
            push!(vals, sizeof(buf))
        elseif isa(val, Ptr{Void})
            push!(vals, val)
        else
            error("cannot handle option $opt")
        end
    end

    @assert length(keys) == length(vals)
    return keys, vals
end

function decode(keys::Array{CUjit_option,1}, vals::Array{Ptr{Void}, 1})
    @assert length(keys) == length(vals)
    options = Dict{CUjit_option,Any}()

    # decode the raw option value bits to their proper type
    for (opt, val) = zip(keys, vals)
        data = reinterpret(UInt, val)
        if opt == WALL_TIME
            options[opt] = convert_bits(Cfloat, data)
        elseif opt == GENERATE_LINE_INFO ||
           opt == GENERATE_DEBUG_INFO ||
           opt == LOG_VERBOSE
            options[opt] = Bool(convert_bits(Cint, data))
        elseif opt == INFO_LOG_BUFFER_SIZE_BYTES ||
           opt == ERROR_LOG_BUFFER_SIZE_BYTES
            options[opt] = convert_bits(Cuint, data)
        elseif opt == INFO_LOG_BUFFER ||
           opt == ERROR_LOG_BUFFER
            options[opt] = reinterpret(Ptr{UInt8}, data)
        else
            options[opt] = data
        end
    end

    # convert some values to easier-to-handle types
    if haskey(options, INFO_LOG_BUFFER)
        buf = options[INFO_LOG_BUFFER]
        size = options[INFO_LOG_BUFFER_SIZE_BYTES]
        delete!(options, INFO_LOG_BUFFER_SIZE_BYTES)
        options[INFO_LOG_BUFFER] = unsafe_string(buf, Int(size / sizeof(UInt8)))
    end
    if haskey(options, ERROR_LOG_BUFFER)
        buf = options[ERROR_LOG_BUFFER]
        size = options[ERROR_LOG_BUFFER_SIZE_BYTES]
        delete!(options, ERROR_LOG_BUFFER_SIZE_BYTES)
        options[ERROR_LOG_BUFFER] = unsafe_string(buf, Int(size / sizeof(UInt8)))
    end

    return options
end


text:
mangling_char(node::Ast.Mangling) = begin
    if isa(node, Ast.ELFMangling)
        return 'e'
    elseif isa(node, Ast.MIPSMangling)
        return 'm'
    elseif isa(node, Ast.MachOMangling)
        return 'o'
    elseif isa(node, Ast.WindowsCOFFMangling)
        return 'w'
    else
        error("unknown type $(typeof(node))")
    end
end

alignment_char(node::Ast.AlignType) = begin
    if isa(node, Ast.IntegerAlign)
        return 'i'
    elseif isa(node, Ast.VectorAlign)
        return 'v'
    elseif isa(node, Ast.FloatAlign)
        return 'f'
    elseif isa(node, Ast.AggregateAlign)
        return 'a'
    else
        error("unknown type $(typeof(node))")
    end
end

endianness_str(node::Ast.BigEndian)    = "E"
endianness_str(node::Ast.LittleEndian) = "e"

mangling_str(node::Ast.Mangling) = "m:$(mangling_char(node))"

nativesizes_str(sizes::Set{Int}) = "n" * join(sort([s for s in sizes]), ":")

write_triple!(io::IO, t::(Int, Ast.AlignmentInfo)) = begin
    s, ai = t
    write(io, string(s))
    write(io, ':')
    write(io, string(ai.abi))
    if ai.perferred != nothing
        write(io, ':')
        write(io, string(ai.perferred))
    end
end

pointerlayouts_str(layouts::Ast.PointerLayoutMap) = begin
    io = IOBuffer()
    write(io, 'p')
    for (a, t) in layouts
        a.val > 0 && write(io, string(a.val))
        write(io, ':')
        write_triple!(io, t)
    end
    return bytestring(io)
end

typelayouts_str(layouts::Ast.TypeLayoutMap) = begin
    strs = {}
    for (t, ai) in layouts
        io = IOBuffer()
        at, s = t
        write(io, alignment_char(at))
        write_triple!(io, (s, ai))
        push!(strs, bytestring(io))
    end
    return strs
end

Base.print(io::IO, layout::Ast.DataLayout) = begin
    strs = {}
    if layout.endianness != nothing
        push!(strs, endianness_str(layout.endianness))
    end
    if layout.mangling != nothing
        push!(strs, mangling_str(layout.mangling))
    end
    if layout.stackalignment != nothing
        push!(strs, string("S", layout.stackalignment))
    end
    if !isempty(layout.pointerlayouts)
        push!(strs, pointerlayouts_str(layout.pointerlayouts))
    end
    if !isempty(layout.typelayouts)
        append!(strs, typelayouts_str(layout.typelayouts))
    end
    if layout.nativesizes != nothing
        push!(strs, nativesizes_str(layout.nativesizes))
    end
    write(io, join(strs, "-"))
end

parse_num(io::IO) = begin
    num = Char[]
    while !eof(io)
        c = read(io, Char)
        if !isdigit(c)
            skip(io, -1)
            break
        end
        push!(num, c)
    end
    return int(utf32(num))
end

parse_triple(io::IO) = begin
    s = parse_num(io)
    read(io, Char) === ':' || error("ill formed data layout")
    abi = parse_num(io)
    if eof(io)
        return (s, Ast.AlignmentInfo(abi, nothing))
    end
    c = read(io, Char)
    if c === ':'
        pref = parse_num(io)
        return (s, Ast.AlignmentInfo(abi, pref))
    elseif c !== '-'
        error("ill formed data layout")
    end
    skip(io, -1)
    return (s, Ast.AlignmentInfo(abi, nothing))
end

parse_spec(io::IO, dl::Ast.DataLayout) = begin
    eof(io) && error("cannot parse empty data layout spec")
    while !eof(io)
        c = read(io, Char)
        if c === 'e'
            dl.endianness == nothing || error("ill formed data layout string")
            dl.endianness = Ast.LittleEndian()
            (eof(io) || read(io, Char) === '-') || error("ill formed data layout")
            continue
        elseif c === 'E'
            dl.endianness == nothing || error("ill formed data layout string")
            dl.endianness = Ast.BigEndian()
            (eof(io) || read(io, Char) === '-') || error("ill formed data layout")
            continue
        elseif c === 'm'
            dl.mangling == nothing || error("ill formed data layout string")
            read(io, Char) === ':' || error("ill formed data layout string")
            c = read(io, Char)
            if c === 'e'
                dl.mangling = Ast.ELFMangling()
            elseif c === 'm'
                dl.mangling = Ast.MIPSMangling()
            elseif c === 'o'
                dl.mangling = Ast.MachOMangling()
            elseif c === 'w'
                dl.mangling = Ast.WindowsCOFFMangling()
            end
            (eof(io) || read(io, Char) === '-') || error("ill formed data layout string")
            continue
        elseif c === 'S'
            dl.stackalignment == nothing || error("ill formed data layout string")
            dl.stackalignment = parse_num(io)
            (eof(io) || read(io, Char) === '-') || error("ill formed data layout string")
            continue
        elseif c === 'p'
            c = read(io, Char)
            as = 0
            if isdigit(c)
                as = int(string(c))
                c = read(io, Char)
            end
            c === ':' || error("ill formed data layout string")
            t = parse_triple(io)
            dl.pointerlayouts[Ast.AddrSpace(as)] = t
            (eof(io) || read(io, Char) === '-') || error("ill formed data layout string")
            continue
        elseif c === 'i'
            at = Ast.IntegerAlign()
            (sz, ai) = parse_triple(io)
            (eof(io) || read(io, Char) === '-') || error("ill formed data layout string")
            dl.typelayouts[(at, sz)] = ai
            continue
        elseif c === 'v'
            at = Ast.VectorAlign()
            (sz, ai) = parse_triple(io)
            (eof(io) || read(io, Char) === '-') || error("ill formed data layout string")
            dl.typelayouts[(at, sz)] = ai
            continue
        elseif c === 'f'
            at = Ast.FloatAlign()
            (sz, ai) = parse_triple(io)
            (eof(io) || read(io, Char) === '-') || error("ill formed data layout string")
            dl.typelayouts[(at, sz)] = ai
            continue
        elseif c === 'a'
            at = Ast.AggregateAlign()
            (sz, ai) = parse_triple(io)
            (eof(io) || read(io, Char) === '-') || error("ill formed data layout string")
            dl.typelayouts[(at, sz)] = ai
            continue
        elseif c === 'n'
            dl.nativesizes == nothing || error("ill formed data layout string")
            ns = Set{Int}()
            while true
                push!(ns, parse_num(io))
                eof(io) && break
                read(io, Char) === ':' || error("ill formed data layout string")
            end
            dl.nativesizes = ns
            eof(io) || error("ill formed data layout string")
        end
    end
    return dl
end

parse_datalayout(io::IO) = parse_spec(io, Ast.DataLayout())
parse_datalayout(s::String) = parse_datalayout(IOBuffer(s))


text:
n = 5
println("The n = ", n, " spiral matrix:")
a = zeros(Int, (n, n))
for (i, s) in enumerate(spiral(n))
    a[s] = i-1
end
println(pretty(a))

m = 3
println()
println("Generalize to a non-square matrix (", m, "x", n, "):")
a = zeros(Int, (m, n))
for (i, s) in enumerate(spiral(m, n))
    a[s] = i-1
end
println(pretty(a))

p = primes(10^3)
n = 7
println()
println("An n = ", n, " prime spiral matrix:")
a = zeros(Int, (n, n))
for (i, s) in enumerate(spiral(n))
    a[s] = p[i]
end
println(pretty(a))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function create_response_file(pm::PredictionModel, test::Bool, target_dir::String, filename::String)
    response_file = joinpath(target_dir, filename)
    f = open(response_file, "w")
    @printf(f, "CellLine")
    drugs = collect(keys(pm.data.results))
    for drug in drugs
        @printf(f, "\t%s", drug.id)
    end
    @printf(f, "\n")

    for cl in collect(values(pm.data.cell_lines))
        # if these two are equal
        if test == cl.in_test_set
            @printf(f, "%s", cl.id)
            for drug in drugs
                if haskey(pm.data.results[drug].outcome_values, cl)
                    @printf(f, "\t%.9f", pm.data.results[drug].outcome_values[cl])
                else
                    @printf(f, "\tNA")
                end
            end
            @printf(f, "\n")
        end
    end
    close(f)
    response_file
end

function write_results(pm::PredictionModel, p::PredictionModelParameters,
        parent_dir::String, filename::String,
        predictions::Dict{Drug, Vector{Float64}}, ranks::Dict{Drug, Vector{Int64}})
    # write predictions and true values
    mkpath(joinpath(parent_dir, "prediction"))
    prediction_file = joinpath(parent_dir, "prediction", filename)
    f = open(prediction_file, "w")
    @printf(f, "DrugAnonID\tType")
    for d in collect(keys(pm.data.results))
        @printf(f, "\t%s (prediction)\t%s (measured)", d.id, d.id)
    end
    @printf(f, "\n")
    for (cl_id, cl) in enumerate(collect(values(pm.data.cell_lines)))
        @printf(f, "%s", cl.id)
        if cl.in_test_set @printf(f, "\tTest") else @printf(f, "\tTraining") end
        for drug in collect(keys(pm.data.results))
            # drug = pm.data.drugs["Drug$d_id"]
            target = 0.
            if haskey(pm.data.results, drug) && haskey(pm.data.results[drug].outcome_values, cl)
                target = pm.data.results[drug].outcome_values[cl]
            # elseif haskey(pm.data.test_results, drug) && haskey(pm.data.test_results[drug].outcome_values, cl)
            #     target = pm.data.test_results[drug].outcome_values[cl]
            end
            @printf(f, "\t%.5f\t%.5f", predictions[drug][cl_id], target)
        end
        @printf(f, "\n")
    end
    close(f)
    # write rankings
    mkpath(joinpath(parent_dir, "ranking"))
    ranking_file = joinpath(parent_dir, "ranking", filename)
    # ranks = Dict{Drug, Vector{Int64}}()
    # we predict neg log values, i.e. higher value means lower concentration
    # which means higher susceptibility of the cell line to the drug and thus higher rank (lower number)
    # hence, we sort in reverse order
    for d in collect(keys(predictions))
        ranks[d] = zeros(Int64, length(predictions[d]))
        ranks[d][sortperm(predictions[d], rev=true)] = collect(1:length(predictions[d]))
    end
    f = open(ranking_file, "w")
    @printf(f, "DrugAnonID")
    for drug in collect(keys(pm.data.results))
        @printf(f, ",%s", drug.id)
    end
    @printf(f, "\n")
    for (cl_id, cl) in enumerate(collect(values(pm.data.cell_lines)))
        @printf(f, "%s", cl.id)
        for drug in collect(keys(pm.data.results))
           @printf(f, ",%d", ranks[drug][cl_id])
       end
       @printf(f, "\n")
    end
    close(f)

    # write model
    mkpath(joinpath(parent_dir, "model"))
    model_file = joinpath(parent_dir, "model", filename)
    f = open(model_file, "w")
    # T, K, N
    @printf(f, "T\tK\n")
    @printf(f, "%d\t%d\n", p.T, p.K)
    @printf(f, "N\n")
    N = collect(values(p.N))
    [@printf(f, "%d\t", n) for n in N[1:end-1]]
    @printf(f, "%d\n", N[end])
    # gamma
    @printf(f, "ɣ\n")
    [@printf(f, "%.5f\t", expected_value(val)) for val in p.ɣ[1:end-1]]
    @printf(f, "%.5f\n", expected_value(p.ɣ[end]))
    # b
    @printf(f, "b\n")
    [@printf(f, "%.5f\t", expected_value(val)) for val in p.b[1:end-1]]
    @printf(f, "%.5f\n", expected_value(p.b[end]))
    # a
    @printf(f, "a\n")
    for a_t in p.a
        @printf(f, "%s\n", a_t.var_name)
        a_t_exp = expected_value(a_t)
        [@printf(f, "%.5f\t", val) for val in a_t_exp[1:end-1]]
        @printf(f, "%.5f\n", a_t_exp[end])
    end
    # lambda
    @printf(f, "λ\n")
    for (t, λ_t) in enumerate(p.λ)
        λ_t_exp = expected_value.(λ_t)
        @printf(f, "λ[%d]:\t", t)
        for val in λ_t_exp[1:end-1]
            @printf(f, "%.5f\t", val)
        end
        @printf(f, "%.5f\n", λ_t_exp[end])
    end
    # epsilon
    @printf(f, "ε\n")
    [@printf(f, "%.5f\t", expected_value(val)) for val in p.ε[1:end-1]]
    @printf(f, "%.5f\n", expected_value(p.ε[end]))
    # nu
    @printf(f, "ν\n")
    [@printf(f, "%.5f\t", expected_value(val)) for val in p.ν[1:end-1]]
    @printf(f, "%.5f\n", expected_value(p.ν[end]))
    # G
    @printf(f, "G\n")
    for t in 1:p.T, k in 1:p.K
        exp_g = expected_value(p.G[t,k])
        @printf(f, "%s:\t", p.G[t,k].var_name)
        [@printf(f, "%.5f\t", val) for val in exp_g[1:end-1]]
        @printf(f, "%.5f\n", exp_g[end])
    end
    # omega
    @printf(f, "⍵\n")
    [@printf(f, "%.5f\t", expected_value(val)) for val in p.⍵[1:end-1]]
    @printf(f, "%.5f\n", expected_value(p.⍵[end]))
    # e
    @printf(f, "e\n")
    [@printf(f, "%.5f\t", expected_value(val)) for val in p.e[1:end-1]]
    @printf(f, "%.5f\n", expected_value(p.e[end]))

    close(f)

end


text:
function totime(seconds)
  h,m = divrem(seconds ÷ 60, 60)
  "$h hour(s) and $m minute(s)"
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas