text:
# # Building an MLJFlux regression model for the Boston house
# # price dataset

using Pkg
Pkg.activate("env")
Pkg.instantiate()

# **Julia version** is assumed to be 1.6.*

using MLJ
using MLJFlux
using Flux

MLJ.color_off()

using Plots

# This tutorial uses MLJ's `IteratedModel` wrapper to transform the
# MLJFlux `NeuralNetworkRegressor` into a model that **automatically
# selects the number of epochs** required to optimize an out-of-sample
# loss.

# We also show how to include the model in a **pipeline** to carry out
# standardization of the features and target.

# ## Loading data

data = OpenML.load(531); # Loads from https://www.openml.org/d/531

# The target `y` is `:MEDV` and everything else except `:CHAS` goes
# into the features `X`:

y, X = unpack(data, ==(:MEDV), !=(:CHAS); rng=123);

# We specified the seed `rng` to shuffle the observations. The Charles
# River dummy variable `:CHAS` is dropped, as not deemed to be
# relevant.

# Inspecting the scientific types:

scitype(y)

#-

schema(X)

# We'll regard `:RAD` (index of accessibility to radial highways) as
# `Continuous` as MLJFlux models don't handle ordered factors:

X = coerce(X, :RAD => Continuous);

# Let's split off a test set for final testing:

X, Xtest = partition(X, 0.7);
y, ytest = partition(y, 0.7);

# ## Defining a builder

# In the macro call below, `n_in` is expected to represent the number
# of inputs features and `rng` a RNG (builders are generic, ie can be
# applied to data with any number of input features):

builder = MLJFlux.@builder begin
    init=Flux.glorot_uniform(rng)
    Chain(Dense(n_in, 64, relu, init=init),
          Dense(64, 32, relu, init=init),
          Dense(32, 1, init=init))
end

# ## Defining a MLJFlux model:

NeuralNetworkRegressor = @load NeuralNetworkRegressor
    model = NeuralNetworkRegressor(builder=builder,
                                   rng=123,
                                   epochs=20)

# ## Standardization

# The following wraps our regressor in feature and target standardizations:

pipe = @pipeline Standardizer model target=Standardizer

# Notice that our original neural network model is now a
# hyper-parameter of the composite `pipe`, with the automatically
# generated name, `:neural_network_regressor`.

# ## Choosing a learning rate

# Let's see how the training losses look for the default optimiser. For
# MLJFlux models, `fit!` will print these losses if we bump the
# verbosity level (default is always 1):

mach = machine(pipe, X, y)
fit!(mach, verbosity=2)

# They are also extractable from the training report (which includes
# the pre-train loss):

report(mach).neural_network_regressor.training_losses

# Next, let's visually compare a few learning rates:

plt = plot()
rates = [5e-5, 1e-4, 0.005, 0.001, 0.05]

# By default, changing only the optimiser will not trigger a
# cold-restart when we `fit!` (to allow for adaptive learning rate
# control). So we call `fit!` with the `force=true`
# option. (Alternatively, one can change the hyper-parameter
# `pipe.neural_network_regressor.optimiser_changes_trigger_retraining`
# to `true`.)

# We'll skip the first few losses to get a better vertical scale in
# our plot.

foreach(rates) do η
    pipe.neural_network_regressor.optimiser.eta = η
    fit!(mach, force=true)
    losses = report(mach).neural_network_regressor.training_losses[3:end]
    plot!(1:length(losses), losses, label=η)
end
plt #!md

#-

savefig(joinpath("assets", "learning_rate.png"))

# ![](assets/learing_rate.png) #md

# We'll go with the second most conservative rate:

pipe.neural_network_regressor.optimiser.eta = 0.0001

# ## Wrapping in iteration control

# We want a model that trains until an out-of-sample loss satisfies
# the `NumberSinceBest(6)` stopping criterion. We'll add some fallback
# stopping criterion `InvalidValue` and `TimeLimit(1/60)`, and
# controls to print traces of the losses.

# For initializing or clearing the traces:

clear() = begin
    global losses = []
    global training_losses = []
    global epochs = []
    return nothing
end

# And to update the traces:

update_loss(loss) = push!(losses, loss)
update_training_loss(report) =
    push!(training_losses, report.neural_network_regressor.training_losses[end])
update_epochs(epoch) = push!(epochs, epoch)

# The controls to apply (see
# [here](https://alan-turing-institute.github.io/MLJ.jl/dev/controlling_iterative_models/#Controls-provided)
# for the complete list):

controls=[Step(1),
          NumberSinceBest(6),
          InvalidValue(),
          TimeLimit(1/60),
          WithLossDo(update_loss),
          WithReportDo(update_training_loss),
          WithIterationsDo(update_epochs)]

# Next we create a "self-iterating" version of the pipeline. Note
# that the iteration parameter is a nested hyperparameter:

iterated_pipe =
    IteratedModel(model=pipe,
                  controls=controls,
                  resampling=Holdout(fraction_train=0.8),
                  iteration_parameter=:(neural_network_regressor.epochs),
                  measure = l2)

# Training the wrapped model on all the train/validation data:

clear()
mach = machine(iterated_pipe, X, y)
fit!(mach)

# And plotting the traces:

plot(epochs, losses,
     xlab = "epoch",
     ylab = "mean sum of squares error",
     label="out-of-sample",
     legend = :topleft);
scatter!(twinx(), epochs, training_losses, label="training", color=:red) #!md

#-

savefig(joinpath("assets", "loss.png"))

# ![](assets/losses.png) #md

# **How `IteratedModel` works.** Training an `IteratedModel` means
# holding out some data (80% in this case) so an out-of-sample loss
# can be tracked and used in the specified stopping criterion,
# `NumberSinceBest(4)`. However, once the stop is triggered, the model
# wrapped by `IteratedModel` (our pipeline model) is retrained on all
# data for the same number of iterations. Calling `predict(mach,
# Xnew)` on new data uses the updated learned parameters.

# In other words, `iterated_model` is a "self-iterating" version of
# the original model, where `epochs` has been transformed from
# hyper-parameter to *learned* parameter.

# ## An evaluation the self-iterating model

# Here's an estimate of performance of our "self-iterating"
# model:

e = evaluate!(mach,
              resampling=CV(nfolds=8),
              measures=[l2, l1])

using Measurements
err = e.measurement[1] ± std(e.per_fold[1])/sqrt(7)
@show err

# which we can see has substantial uncertainty (not that CV estimates
# of the uncertainty are trustworthy; see e.g., [this
# survey](https://direct.mit.edu/neco/article-abstract/10/7/1895/6224/Approximate-Statistical-Tests-for-Comparing))

# ## Comparison with other models on the test set

# Although we cannot assign them statistical significance, here are
# comparisons, on the untouched test set, of the eror of our
# self-iterating neural network regressor with a couple of other
# models trained on the same data (using default hyperparameters):

function performance(model)
    mach = machine(model, X, y) |> fit!
    yhat = predict(mach, Xtest)
    l2(yhat, ytest) |> mean
end
performance(iterated_pipe)

three_models = [(@load EvoTreeRegressor)(), # tree boosting model
                (@load LinearRegressor pkg=MLJLinearModels)(),
                iterated_pipe]

errs = performance.(three_models)

(models=typeof.(three_models), mean_square_errors=errs) |> pretty


text:
using GraphMatFun, LinearAlgebra,Random
include("simulationtools.jl");

# Julia's base can change the random seed generators.
# We hardcode it to make the results reproducible over
# Julia versions
function hardcoded_rand(t)

    x=[  0.8872174233685589;  0.8377413644286786;  0.5608614162945871;  0.7142032679122071;  0.8217070528857874;  0.40399577722972235;  0.8659919887404998;  0.08029198288895634;  0.10052875305338604;  0.40627896679851794;  0.8342948272930204;  0.3820649650758666;  0.016757502407812153;  0.20137662172350335;  0.2555463313583042;  0.19248011642467033;  0.011618055692006912;  0.42854938782243956;  0.21254339113365306;  0.8654170264281245;  0.26087081576693616;  0.0330965543051257;  0.15171093545767034;  0.541356766187003;  0.7025653111402541;  0.3832283012090676;  0.24532744474494972;  0.8773921298355143;  0.28729058543087427;  0.02645861903949709;  0.10877066841853389;  0.25559799376084;  0.9817655834384857;  0.09350052884595439;  0.5129871676353079;  0.9371434490087653;  0.8916549976510515;  0.9426774471099136;  0.45772957349208;  0.9332990364480496;  0.8382243800948085;  0.40812659546896835;  0.653899520291725;  0.6110199290298604;  0.19036516891430988;  0.26967359818868164;  0.7360989887783765; 0.296931260360374];
    x=x[1:t[1]];
    return x;

end

# Construct Denman-Beavers iteration for the modified
# degree-optimal form.
function graph_generalized_denman_beavers()

    # Construct a Degopt and change :mult to :ldiv
    HA=[1.0  1.0   0.0  0.0  0.0
        0.5  0.0   0.5  0.0  0.0
        1.0  0.5   0.0  0.0  0.0
        0.25  0.0  0.25 0.0  0.5];

    HB=[1.0 0.0 0.0 0.0 0.0
        1.0 0.0 0.0 0.0 0.0
        1.0 0.0 0.0 0.0 0.0
        1.0 0.0 0.0 0.0 0.0];
    y=[1/4; 1/8;0;1/4;0;1/2];

    degopt=Degopt(HA,HB,y);
    (graph,cref)=graph_degopt(degopt);
    # Replace :mult with ldiv
    for k in keys(graph.operations)
        if (graph.operations[k] == :mult)
            graph.operations[k] = :ldiv
        end
    end

    return (graph,cref)
end
# Get all the coefficient valeus in the mod degopt stored in graph
function get_coeffs_mod_degopt(graph)
    graph=deepcopy(graph);
    for k in keys(graph.operations)
        if (graph.operations[k] == :ldiv)
            graph.operations[k] = :mult
        end
    end
    return get_degopt_coeffs(graph);
end
# Transform all rational functions such with a normalization
# of the maximum coefficient value in order to avoid
# large/small numbers
function normalize_mod_degopt!(graph)
    (HA,HB,y)=get_coeffs_mod_degopt(graph);

    for k=1:size(HA,1)
        alpha=maximum(abs.(HB[k,:]));
        HA[k,:] /= alpha
        HB[k,:] /= alpha
        beta=maximum(abs.(HA[k,:]));
        HA[k,:] /= beta;
        if (k+2 <= size(HA,2))
            HA[:,k+2] /= beta;
            HB[:,k+2] /= beta;
        end
        y[k+2] /= beta;
    end

    degopt=Degopt(HA,HB,y);
    (graph,cref)=graph_degopt(degopt);
    # Replace :mult with ldiv
    for k in keys(graph.operations)
        if (graph.operations[k] == :mult)
            graph.operations[k] = :ldiv
        end
    end

    return (graph,cref)

end

statelist=[];

f=s-> sqrt(s+1);

## Standard Denman-Beavers with shift *********************
(graph,cref)=graph_denman_beavers(2);
rename_node!(graph,:A,:A_shift,cref);
add_lincomb!(graph,:A_shift,1.0,:A,1.0,:I);
push!(cref,(:A_shift,1))
push!(cref,(:A_shift,2))

x=0.5;
f(x)-eval_graph(graph,x);
graph.outputs
o=graph.outputs[end];
empty!(graph.outputs)
add_output!(graph,o);
compress_graph!(graph,cref)

## Generalized Denman-Beavers (with shift) *********************
println("Optimizing Denman-Beavers (modified degree-optimal form) *********************");
T=BigFloat;
state=State(f,complex(T));
state.params[:rho]=0.5;
state.params[:target_n]=2000;
state.params[:n]=1000;
state.params[:graphname]="denman_beavers";
state.graph=graph;
showerr(state,showmeta=true)
push!(statelist,deepcopy(state))

dd=count(values(graph.operations) .== :ldiv)

println("Using $dd inverse multiplications");

x=get_coeffs(graph,cref);  x=x .+ 0.00001*hardcoded_rand(size(x)); set_coeffs!(graph,x,cref);

(graph,cref)= graph_generalized_denman_beavers();
state.graph=graph;
state.params[:graphname]="denman_beavers_opt";
discr=get_disc_discr(state);

opt_gauss_newton!(graph,f,discr,logger=1,cref=cref,
                  linlsqr=:real_svd,droptol=1e-8,
                  stoptol=1e-16,errtype=:relerr,maxit=5)
showerr(state,showmeta=true)
opt_gauss_newton!(graph,f,discr,logger=1,cref=cref,
                  linlsqr=:real_svd,droptol=1e-9,
                  stoptol=1e-16,errtype=:relerr,maxit=5)
showerr(state,showmeta=true)
opt_gauss_newton!(graph,f,discr,logger=1,cref=cref,
                  linlsqr=:real_svd,droptol=1e-11,
                  stoptol=1e-16,errtype=:relerr,maxit=5)
showerr(state,showmeta=true)
opt_gauss_newton!(graph,f,discr,logger=1,cref=cref,
                  linlsqr=:real_svd,droptol=1e-12,
                  stoptol=1e-16,errtype=:relerr,maxit=5)
showerr(state,showmeta=true)
opt_gauss_newton!(graph,f,discr,logger=1,cref=cref,
                  linlsqr=:real_svd,droptol=1e-13,
                  stoptol=1e-16,errtype=:relerr,maxit=5)
showerr(state,showmeta=true)
opt_gauss_newton!(graph,f,discr,logger=1,cref=cref,
                  linlsqr=:real_svd,droptol=1e-14,
                  stoptol=1e-16,errtype=:relerr,maxit=10)
count(values(graph.operations) .== :ldiv)
push!(statelist,deepcopy(state))

for state in deepcopy(statelist)
    #state.graph=Compgraph(Float64,state.graph)
    showerr(state,showmeta=true)
end

println("Optimizing PS  *********************");

# Computing the least squares match of coefficients on the disk
#(graph1,cref1)=graph_monomial(1:10.0) # If we want 4 multiplications
(graph1,cref1)=graph_monomial(1:13.0)
graph1=Compgraph(ComplexF64,graph1);
opt_linear_fit!(graph1,f,discr,cref1)
coeffs=get_coeffs(graph1,cref1);
(graph1,cref1)=graph_ps_degopt(real.(coeffs));

mm=count(values(graph1.operations) .== :mult)
println("Using $mm multiplications");
state1=deepcopy(state);
state1.graph=graph1;
state1.params[:n]=200;
state1.params[:graphname]="ps_m$(mm)";
showerr(state1,showmeta=true)
push!(statelist,deepcopy(state1))

## PS with degopt freedom optimization   *********************

x=get_coeffs(graph1,cref1);  x=x .+ 0.00001*hardcoded_rand(size(x,1)); set_coeffs!(graph1,x,cref1);

state1.params[:graphname]="ps_m$(mm)_opt";
discr=get_disc_discr(state1);
opt_gauss_newton!(graph1,f,discr,logger=1,cref=cref1,
                  linlsqr=:real_svd,droptol=1e-5,γ0=0.5,
                  stoptol=1e-14,errtype=:relerr,maxit=10)
showerr(state1,showmeta=true)
opt_gauss_newton!(graph1,f,discr,logger=1,cref=cref1,
                  linlsqr=:real_svd,droptol=1e-6,γ0=0.3,
                  stoptol=1e-14,errtype=:relerr,maxit=10)
showerr(state1,showmeta=true)
opt_gauss_newton!(graph1,f,discr,logger=1,cref=cref1,
                         linlsqr=:real_svd,droptol=1e-7,γ0=0.1,
                         stoptol=1e-14,errtype=:relerr,maxit=10)
showerr(state1,showmeta=true)
opt_gauss_newton!(graph1,f,discr,logger=1,cref=cref1,
                                linlsqr=:real_svd,droptol=1e-7,
                                stoptol=1e-14,errtype=:relerr,maxit=20)
showerr(state1,showmeta=true)
opt_gauss_newton!(graph1,f,discr,logger=1,cref=cref1,
                                linlsqr=:real_svd,droptol=1e-8,γ0=0.1,
                                stoptol=1e-14,errtype=:relerr,maxit=20)
showerr(state1,showmeta=true)
opt_gauss_newton!(graph1,f,discr,logger=1,cref=cref1,
                                linlsqr=:real_svd,droptol=1e-8,γ0=0.5,
                  stoptol=1e-14,errtype=:relerr,maxit=20)
showerr(state1,showmeta=true)
opt_gauss_newton!(graph1,f,discr,logger=1,cref=cref1,
                                       linlsqr=:real_svd,droptol=1e-8,γ0=1.0,
                                       stoptol=1e-14,errtype=:relerr,maxit=5)

showerr(state1,showmeta=true)
opt_gauss_newton!(graph1,f,discr,logger=1,cref=cref1,
                                linlsqr=:real_svd,droptol=1e-9,γ0=0.5,
                  stoptol=1e-14,errtype=:relerr,maxit=20)
showerr(state1,showmeta=true)
opt_gauss_newton!(graph1,f,discr,logger=1,cref=cref1,
                  linlsqr=:real_svd,droptol=1e-9,γ0=1.0,
                  stoptol=1e-14,errtype=:relerr,maxit=5)
showerr(state1,showmeta=true)

opt_gauss_newton!(graph1,f,discr,logger=1,cref=cref1,
                                       linlsqr=:real_svd,droptol=1e-10,γ0=0.5,
                         stoptol=1e-16,errtype=:relerr,maxit=30)
showerr(state1,showmeta=true)
opt_gauss_newton!(graph1,f,discr,logger=1,cref=cref1,
                  linlsqr=:real_svd,droptol=1e-10,γ0=1.0,
                  stoptol=1e-14,errtype=:relerr,maxit=5)

push!(statelist,deepcopy(state1))

Random.seed!(0);
A=randn(100,100);
A=0.5*A/opnorm(A);

for state in deepcopy(statelist)
    state.graph=Compgraph(Float64,state.graph)
    showerr(state,showmeta=true)

    graph=state.graph;
    if (:ldiv in values(graph.operations))

        try
            (graph,_)=normalize_mod_degopt!(graph)
            (HA,HB,y)=get_coeffs_mod_degopt(graph)

        catch
            HA=[];
            HB=[];
            y=[];
        end

    else
        (HA,HB,y)=get_degopt_coeffs(graph)
    end

    #@show latexify([HA HB])
    #@show latexify(y')

end

@show norm(eval_graph(graph,A)-sqrt(A+I))
@show norm(eval_graph(graph1,A)-sqrt(A+I))

saveall=false;
if (saveall)
    for state in deepcopy(statelist)
        graph=state.graph;
        if (:ldiv in values(graph.operations))

            try
                (graph,_)=normalize_mod_degopt!(graph)
                (HA,HB,y)=get_coeffs_mod_degopt(graph)

            catch
            end
        end
        graphname=state.params[:graphname]
        basename=joinpath("..","..","data","sqrt_plus_one",graphname);
        @show "$(basename).cgr"
        export_compgraph(graph,"$(basename).cgr")
    end
end

ext: jl

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
function mu()

end


text:
"""
    module FilterTestUtils

This module contains utilities for evaluating and debugging HDF5 Filters.
"""
module FilterTestUtils

import HDF5.API
import H5Zlz4: H5Z_filter_lz4
import H5Zzstd: H5Z_filter_zstd
import H5Zbzip2: H5Z_filter_bzip2
using Test

export test_filter

function test_filter_init(; cd_values = Cuint[], data = ones(UInt8, 1024))
    flags = Cuint(0)
    nbytes = sizeof(data)
    buf_size = Ref(Csize_t(sizeof(data)))
    databuf = Libc.malloc(sizeof(data))
    data = reinterpret(UInt8, data)
    unsafe_copyto!(Ptr{UInt8}(databuf), pointer(data), sizeof(data))
    buf = Ref(Ptr{Cvoid}(databuf))
    return flags, cd_values, nbytes, buf_size, buf
end

function test_filter_compress!(filter_func, flags::Cuint, cd_values::Vector{Cuint}, nbytes::Integer, buf_size::Ref{Csize_t}, buf::Ref{Ptr{Cvoid}})
    nbytes = Csize_t(nbytes)
    cd_nelmts = Csize_t(length(cd_values))
    GC.@preserve flags cd_nelmts cd_values nbytes buf_size buf begin
        ret_code = filter_func(
            flags,
            cd_nelmts,
            pointer(cd_values),
            Csize_t(nbytes),
            Base.unsafe_convert(Ptr{Csize_t}, buf_size),
            Base.unsafe_convert(Ptr{Ptr{Cvoid}}, buf)
        )
        @debug "Compression:" ret_code buf_size[]
        if ret_code <= 0
            error("Test compression failed: $ret_code.")
        end
    end
    return ret_code
end

function test_filter_decompress!(filter_func, flags::Cuint, cd_values::Vector{Cuint}, nbytes::Integer, buf_size::Ref{Csize_t}, buf::Ref{Ptr{Cvoid}})
    nbytes = Csize_t(nbytes)
    cd_nelmts = Csize_t(length(cd_values))
    flags |= UInt32(API.H5Z_FLAG_REVERSE)
    GC.@preserve flags cd_nelmts cd_values nbytes buf_size buf begin
        ret_code = filter_func(
            flags,
            cd_nelmts,
            pointer(cd_values),
            Csize_t(nbytes),
            Base.unsafe_convert(Ptr{Csize_t},buf_size),
            Base.unsafe_convert(Ptr{Ptr{Cvoid}}, buf)
        )
        @debug "Decompression:" ret_code buf_size[]
    end
    return ret_code
end

function test_filter_cleanup!(buf::Ref{Ptr{Cvoid}})
    Libc.free(buf[])
end

function test_filter(filter_func; cd_values::Vector{Cuint} = Cuint[], data = ones(UInt8, 1024))
    flags, cd_values, nbytes, buf_size, buf = test_filter_init(; cd_values = cd_values, data = data)
    nbytes_compressed, nbytes_decompressed = 0, 0
    try
        nbytes_compressed = test_filter_compress!(filter_func, flags, cd_values, nbytes, buf_size, buf)
        nbytes_decompressed = test_filter_decompress!(filter_func, flags, cd_values, nbytes_compressed, buf_size, buf)
        if nbytes_decompressed > 0
            # ret_code is the number of bytes out
            round_trip_data = unsafe_wrap(Array,Ptr{UInt8}(buf[]), nbytes_decompressed)
            @debug "Is the data the same after a roundtrip?" data == round_trip_data
        end
    catch err
        rethrow(err)
    finally
        test_filter_cleanup!(buf)
    end
    @debug "Compression Ratio" nbytes_compressed / nbytes_decompressed
    return nbytes_compressed, nbytes_decompressed
end

function test_bzip2_filter(data = ones(UInt8, 1024))
    cd_values = Cuint[8]
    test_filter(H5Z_filter_bzip2; cd_values = cd_values, data = data)
end

function test_lz4_filter(data = ones(UInt8, 1024))
    cd_values = Cuint[1024]
    test_filter(H5Z_filter_lz4; cd_values = cd_values, data = data)
end

function test_zstd_filter(data = ones(UInt8, 1024))
    cd_values = Cuint[3] # aggression
    test_filter(H5Z_filter_zstd; cd_values = cd_values, data = data)
end

function __init__()
    @testset "Compression Filter Unit Tests" begin
        @test argmin(test_bzip2_filter()) == 1
        @test argmin(test_lz4_filter()) == 1
        @test argmin(test_zstd_filter()) == 1
        str = codeunits(repeat("foobar", 1000))
        @test argmin(test_bzip2_filter(str)) == 1
        @test argmin(test_lz4_filter(str)) == 1
        @test argmin(test_zstd_filter(str)) == 1
    end
end

end


text:
# # Julia & Lux for the Uninitiated

# This is a quick intro to [Lux](https://github.com/avik-pal/:ux.jl) loosely based on:
#
# 1. [PyTorch's tutorial](https://pytorch.org/tutorials/beginner/deep_learning_60min_blitz.html).
# 2. [Flux's tutorial](https://fluxml.ai/tutorials/2020/09/15/deep-learning-flux.html).
# 3. [Flax's tutorial](https://flax.readthedocs.io/en/latest/notebooks/jax_for_the_impatient.html).
#
# It introduces basic Julia programming, as well `Zygote`, a source-to-source automatic differentiation (AD) framework in Julia. We'll use these tools to build a very simple neural network. Let's start with importing `Lux.jl`

using Lux, Random
using Pkg #hide
Pkg.activate(joinpath(dirname(pathof(Lux)), "..", "examples")) #hide

# Now let us control the randomness in our code using proper Pseudo Random Number Generator (PRNG)
rng = Random.default_rng()
Random.seed!(rng, 0)

# ## Arrays

# The starting point for all of our models is the `Array` (sometimes referred to as a `Tensor` in other frameworks). This is really just a list of numbers, which might be arranged into a shape like a square. Let's write down an array with three elements.

x = [1, 2, 3]

# Here's a matrix – a square array with four elements.

x = [1 2; 3 4]

# We often work with arrays of thousands of elements, and don't usually write them down by hand. Here's how we can create an array of 5×3 = 15 elements, each a random number from zero to one.

x = rand(rng, 5, 3)

# There's a few functions like this; try replacing `rand` with `ones`, `zeros`, or `randn`.

# By default, Julia works stores numbers is a high-precision format called `Float64`. In ML we often don't need all those digits, and can ask Julia to work with `Float32` instead. We can even ask for more digits using `BigFloat`.

x = rand(BigFloat, 5, 3)
#-
x = rand(Float32, 5, 3)

# We can ask the array how many elements it has.

length(x)

# Or, more specifically, what size it has.

size(x)

# We sometimes want to see some elements of the array on their own.

x
#-
x[2, 3]

# This means get the second row and the third column. We can also get every row of the third column.

x[:, 3]

# We can add arrays, and subtract them, which adds or subtracts each element of the array.

x + x
#-
x - x

# Julia supports a feature called *broadcasting*, using the `.` syntax. This tiles small arrays (or single numbers) to fill bigger ones.

x .+ 1

# We can see Julia tile the column vector `1:5` across all rows of the larger array.

zeros(5, 5) .+ (1:5)

# The x' syntax is used to transpose a column `1:5` into an equivalent row, and Julia will tile that across columns.

zeros(5, 5) .+ (1:5)'

# We can use this to make a times table.

(1:5) .* (1:5)'

# Finally, and importantly for machine learning, we can conveniently do things like matrix multiply.

W = randn(5, 10)
x = rand(10)
W * x

# Julia's arrays are very powerful, and you can learn more about what they can do [here](https://docs.julialang.org/en/v1/manual/arrays/).

# ### CUDA Arrays

# CUDA functionality is provided separately by the [CUDA.jl package](https://github.com/JuliaGPU/CUDA.jl). If you have a GPU and CUDA available, `Lux` will automatically build the required CUDA dependencies using `CUDA.jl`.

# You can manually add `CUDA`. Once CUDA is loaded you can move any array to the GPU with the `cu` function (or the `gpu` function exported by `Lux``), and it supports all of the above operations with the same syntax.

## using CUDA
## x = cu(rand(5, 3))

# ## (Im)mutability

# Lux as you might have read is [Immutable by convention](http://lux.csail.mit.edu/dev/introduction/overview/#Design-Principles) which means that the core library is built without any form of mutation and all functions are pure. However, we don't enfore it in any form. We do **strongly recommend** that users extending this framework for their respective applications don't mutate their arrays.

x = reshape(1:8, 2, 4)

# To update this array, we should first copy the array.

x_copy = copy(x)
view(x_copy, :, 1) .= 0

println("Original Array ", x)
println("Mutated Array ", x_copy)

# Note that our current default AD engine (Zygote) is unable to differentiate through this mutation, however, for these specialized cases it is quite trivial to write custom backward passes. (This problem will be fixed once we move towards Enzyme.jl)

# ## Managing Randomness

# We relu on the Julia StdLib `Random` for managing the randomness in our execution. First, we create an PRNG and seed it.
rng = Random.default_rng() # Creates a Xoshiro PRNG
Random.seed!(rng, 0)

# If we call any function that relies on `rng` and uses it via `randn`, `rand`, etc. `rng` will be mutated. As we have already established we care a lot about immutability, hence we should use `Lux.replicate` on PRNG before using them.

# First, let us run a random number generator 3 times with the `replicate`d rng

for i in 1:3
    println("Iteration $i ", rand(Lux.replicate(rng), 10))
end

# As expected we get the same output. We can remove the `replicate` call and we will get different outputs

for i in 1:3
    println("Iteration $i ", rand(rng, 10))
end

# ## Automatic Differentiation

# Julia has quite a few (maybe too many) AD tools. For the purpose of this tutorial, we will use [AbstractDifferentiation.jl](https://github.com/JuliaDiff/AbstractDifferentiation.jl) which provides an uniform API across multiple AD backends. For the backends we will use:
#
# 1. [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) -- For Jacobian-Vector Product (JVP)
# 2. [Zygote.jl](https://github.com/FluxML/Zygote.jl) -- For Vector-Jacobian Product (VJP)
#
# *Slight Detour*: We have had several questions regarding if we will be considering any other AD system for the reverse-diff backend. For now we will stick to Zygote.jl, however once [Enzyme.jl](https://github.com/EnzymeAD/Enzyme.jl) has support for custom rules and we have tested Lux extensively with it, we will make the switch.

# Even though, theoretically, a VJP (Vector-Jacobian product - reverse autodiff) and a JVP (Jacobian-Vector product - forward-mode autodiff) are similar—they compute a product of a Jacobian and a vector—they differ by the computational complexity of the operation. In short, when you have a large number of parameters (hence a wide matrix), a JVP is less efficient computationally than a VJP, and, conversely, a JVP is more efficient when the Jacobian matrix is a tall matrix.

using ForwardDiff, Zygote, AbstractDifferentiation

# ### Gradients

# For our first example, consider a simple function computing ``f(x) = \frac{1}{2}x^T x``, where ``\nabla f(x) = x``

f(x) = x' * x / 2
∇f(x) = x
v = randn(rng, Float32, 4)

# Let's use AbstractDifferentiation and Zygote to compute the gradients

println("Actual Gradient: ", ∇f(v))
println("Computed Gradient via Reverse Mode AD (Zygote): ",
        AD.gradient(AD.ZygoteBackend(), f, v)[1])
println("Computed Gradient via Forward Mode AD (ForwardDiff): ",
        AD.gradient(AD.ForwardDiffBackend(), f, v)[1])

# Note that `AD.gradient` will only work for scalar valued outputs

# ### Jacobian-Vector Product

# I will defer the discussion on forward-mode AD to https://book.sciml.ai/notes/08/. Here let us just look at a mini example on how to use it.

f(x) = x .* x ./ 2
x = randn(rng, Float32, 5)
v = ones(Float32, 5)

# Construct the pushforward function.

pf_f = AD.value_and_pushforward_function(AD.ForwardDiffBackend(), f, x)

# Compute the jvp

val, jvp = pf_f(v)
println("Computed Value: f(", x, ") = ", val)
println("JVP: ", jvp[1])

# ### Vector-Jacobian Product

# Using the same function and inputs, let us compute the VJP

pb_f = AD.value_and_pullback_function(AD.ZygoteBackend(), f, x)

# Compute the vjp

val, vjp = pb_f(v)
println("Computed Value: f(", x, ") = ", val)
println("VJP: ", vjp[1])

# ## Linear Regression

# Finally, now let us consider a linear regression problem. From a set of data-points
# $\left\{ (x_i, y_i), i \in \left\{ 1, \dots, k \right\}, x_i \in \R^n, y_i \in \R^m \right\}$,
# we try to find a set of parameters $W$ and $b$, s.t. $f_{W,b}(x) = Wx + b$ minimizes the mean squared error:

# $$L(W, b) \longrightarrow \sum_{i = 1}^{k} \frac{1}{2} \| y_i - f_{W,b}(x_i) \|_2^2$$

# We can write `f` from scratch, but to demonstrate `Lux` let us use the `Dense` layer.

model = Dense(10 => 5)

rng = Random.default_rng()
Random.seed!(rng, 0)

# Let us initialize the parameters and states (in this case it is empty) for the model
ps, st = Lux.setup(rng, model)
ps = ps |> ComponentArray

# Set problem dimensions.
n_samples = 20
x_dim = 10
y_dim = 5

# Generate random ground truth W and b.
W = randn(rng, Float32, y_dim, x_dim)
b = randn(rng, Float32, y_dim)

# Generate samples with additional noise.
x_samples = randn(rng, Float32, x_dim, n_samples)
y_samples = W * x_samples .+ b .+ 0.01f0 .* randn(rng, Float32, y_dim, n_samples)
println("x shape: ", size(x_samples), "; y shape: ", size(y_samples))

# For updating our parameters let's use [Optimisers.jl](https://github.com/FluxML/Optimisers.jl)

using Optimisers

opt = Optimisers.Descent(0.01f0)

# Initialize the initial state of the optimiser

opt_state = Optimisers.setup(opt, ps)

# Define the loss function
mse(model, ps, st, X, y) = sum(abs2, model(X, ps, st)[1] .- y)
mse(weight, bias, X, y) = sum(abs2, weight * X .+ bias .- y)
loss_function(ps, X, y) = mse(model, ps, st, X, y)

println("Loss Value with ground true W & b: ", mse(W, b, x_samples, y_samples))

for i in 1:100
    ## In actual code, don't use globals. But here I will simply for the sake of demonstration
    global ps, st, opt_state
    ## Compute the gradient
    gs = gradient(loss_function, ps, x_samples, y_samples)[1]
    ## Perform parameter update
    opt_state, ps = Optimisers.update(opt_state, ps, gs)
    if i % 10 == 1 || i == 100
        println("Loss Value after $i iterations: ",
                mse(model, ps, st, x_samples, y_samples))
    end
end


text:
module TestCWrapper

using Ipopt
using Test

function test_hs071()
  # hs071
  # min x1 * x4 * (x1 + x2 + x3) + x3
  # st  x1 * x2 * x3 * x4 >= 25
  #     x1^2 + x2^2 + x3^2 + x4^2 = 40
  #     1 <= x1, x2, x3, x4 <= 5
  # Start at (1,5,5,1)
  # End at (1.000..., 4.743..., 3.821..., 1.379...)

  function eval_f(x::Vector{Float64})
    return x[1] * x[4] * (x[1] + x[2] + x[3]) + x[3]
  end

  function eval_g(x::Vector{Float64}, g::Vector{Float64})
    # Bad: g    = zeros(2)  # Allocates new array
    # OK:  g[:] = zeros(2)  # Modifies 'in place'
    g[1] = x[1]   * x[2]   * x[3]   * x[4]
    g[2] = x[1]^2 + x[2]^2 + x[3]^2 + x[4]^2
  end

  function eval_grad_f(x::Vector{Float64}, grad_f::Vector{Float64})
    # Bad: grad_f    = zeros(4)  # Allocates new array
    # OK:  grad_f[:] = zeros(4)  # Modifies 'in place'
    grad_f[1] = x[1] * x[4] + x[4] * (x[1] + x[2] + x[3])
    grad_f[2] = x[1] * x[4]
    grad_f[3] = x[1] * x[4] + 1
    grad_f[4] = x[1] * (x[1] + x[2] + x[3])
  end

  function eval_jac_g(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, values::Vector{Float64})
    if mode == :Structure
      # Constraint (row) 1
      rows[1] = 1; cols[1] = 1
      rows[2] = 1; cols[2] = 2
      rows[3] = 1; cols[3] = 3
      rows[4] = 1; cols[4] = 4
      # Constraint (row) 2
      rows[5] = 2; cols[5] = 1
      rows[6] = 2; cols[6] = 2
      rows[7] = 2; cols[7] = 3
      rows[8] = 2; cols[8] = 4
    else
      # Constraint (row) 1
      values[1] = x[2]*x[3]*x[4]  # 1,1
      values[2] = x[1]*x[3]*x[4]  # 1,2
      values[3] = x[1]*x[2]*x[4]  # 1,3
      values[4] = x[1]*x[2]*x[3]  # 1,4
      # Constraint (row) 2
      values[5] = 2*x[1]  # 2,1
      values[6] = 2*x[2]  # 2,2
      values[7] = 2*x[3]  # 2,3
      values[8] = 2*x[4]  # 2,4
    end
  end

  function eval_h(x::Vector{Float64}, mode, rows::Vector{Int32}, cols::Vector{Int32}, obj_factor::Float64, lambda::Vector{Float64}, values::Vector{Float64})
    if mode == :Structure
      # Symmetric matrix, fill the lower left triangle only
      idx = 1
      for row = 1:4
        for col = 1:row
          rows[idx] = row
          cols[idx] = col
          idx += 1
        end
      end
    else
      # Again, only lower left triangle
      # Objective
      values[1] = obj_factor * (2*x[4])  # 1,1
      values[2] = obj_factor * (  x[4])  # 2,1
      values[3] = 0                      # 2,2
      values[4] = obj_factor * (  x[4])  # 3,1
      values[5] = 0                      # 3,2
      values[6] = 0                      # 3,3
      values[7] = obj_factor * (2*x[1] + x[2] + x[3])  # 4,1
      values[8] = obj_factor * (  x[1])  # 4,2
      values[9] = obj_factor * (  x[1])  # 4,3
      values[10] = 0                     # 4,4

      # First constraint
      values[2] += lambda[1] * (x[3] * x[4])  # 2,1
      values[4] += lambda[1] * (x[2] * x[4])  # 3,1
      values[5] += lambda[1] * (x[1] * x[4])  # 3,2
      values[7] += lambda[1] * (x[2] * x[3])  # 4,1
      values[8] += lambda[1] * (x[1] * x[3])  # 4,2
      values[9] += lambda[1] * (x[1] * x[2])  # 4,3

      # Second constraint
      values[1]  += lambda[2] * 2  # 1,1
      values[3]  += lambda[2] * 2  # 2,2
      values[6]  += lambda[2] * 2  # 3,3
      values[10] += lambda[2] * 2  # 4,4
    end
  end

  n = 4
  x_L = [1.0, 1.0, 1.0, 1.0]
  x_U = [5.0, 5.0, 5.0, 5.0]

  m = 2
  g_L = [25.0, 40.0]
  g_U = [2.0e19, 40.0]

  prob = createProblem(n, x_L, x_U, m, g_L, g_U, 8, 10,
  eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h)

  prob.x = [1.0, 5.0, 5.0, 1.0]
  solvestat = solveProblem(prob)

  @test Ipopt.ApplicationReturnStatus[solvestat] == :Solve_Succeeded
  @test prob.x[1] ≈ 1.0000000000000000 atol=1e-5
  @test prob.x[2] ≈ 4.7429996418092970 atol=1e-5
  @test prob.x[3] ≈ 3.8211499817883077 atol=1e-5
  @test prob.x[4] ≈ 1.3794082897556983 atol=1e-5
  @test prob.obj_val ≈ 17.014017145179164 atol=1e-5

  # This tests callbacks.
  function intermediate(alg_mod::Int, iter_count::Int,
    obj_value::Float64, inf_pr::Float64, inf_du::Float64, mu::Float64,
    d_norm::Float64, regularization_size::Float64, alpha_du::Float64, alpha_pr::Float64,
    ls_trials::Int)
    return iter_count < 1  # Interrupts after one iteration.
  end

  setIntermediateCallback(prob, intermediate)

  solvestat = solveProblem(prob)

  @test Ipopt.ApplicationReturnStatus[solvestat] == :User_Requested_Stop

  # Test setting some options
  # String option
  println("\nString option")
  addOption(prob, "hessian_approximation", "exact")
  @test_throws ErrorException addOption(prob, "hessian_approximation", "badoption")
  println("\nInt option")
  # Int option
  addOption(prob, "file_print_level", 3) == nothing
  @test_throws ErrorException addOption(prob, "file_print_level", -1)
  # Double option
  println("\nFloat option")
  addOption(prob, "derivative_test_tol", 0.5)
  @test_throws ErrorException addOption(prob, "derivative_test_tol", -1.0)

  # Test opening an output file
  openOutputFile(prob, "blah.txt", 5)
  Ipopt.freeProblem(prob) # Needed before the `rm` on Windows.
  # unlink the output file
  rm("blah.txt")
end

function test_binary()
  # Test that the ipopt binary works
  @test Ipopt.amplexefun("-v") == 0
end

end  # TestCWrapper

runtests(TestCWrapper)


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

@doc raw"""Spec to control the desired behavior of rolling update.

    IoK8sApiAppsV1beta2RollingUpdateDeployment(;
        maxSurge=nothing,
        maxUnavailable=nothing,
    )

    - maxSurge::IoK8sApimachineryPkgUtilIntstrIntOrString : The maximum number of pods that can be scheduled above the desired number of pods. Value can be an absolute number (ex: 5) or a percentage of desired pods (ex: 10%). This can not be 0 if MaxUnavailable is 0. Absolute number is calculated from percentage by rounding up. Defaults to 25%. Example: when this is set to 30%, the new ReplicaSet can be scaled up immediately when the rolling update starts, such that the total number of old and new pods do not exceed 130% of desired pods. Once old pods have been killed, new ReplicaSet can be scaled up further, ensuring that total number of pods running at any time during the update is at most 130% of desired pods.
    - maxUnavailable::IoK8sApimachineryPkgUtilIntstrIntOrString : The maximum number of pods that can be unavailable during the update. Value can be an absolute number (ex: 5) or a percentage of desired pods (ex: 10%). Absolute number is calculated from percentage by rounding down. This can not be 0 if MaxSurge is 0. Defaults to 25%. Example: when this is set to 30%, the old ReplicaSet can be scaled down to 70% of desired pods immediately when the rolling update starts. Once new pods are ready, old ReplicaSet can be scaled down further, followed by scaling up the new ReplicaSet, ensuring that the total number of pods available at all times during the update is at least 70% of desired pods.
"""
mutable struct IoK8sApiAppsV1beta2RollingUpdateDeployment <: SwaggerModel
    maxSurge::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgUtilIntstrIntOrString } # spec name: maxSurge
    maxUnavailable::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgUtilIntstrIntOrString } # spec name: maxUnavailable

    function IoK8sApiAppsV1beta2RollingUpdateDeployment(;maxSurge=nothing, maxUnavailable=nothing)
        o = new()
        validate_property(IoK8sApiAppsV1beta2RollingUpdateDeployment, Symbol("maxSurge"), maxSurge)
        setfield!(o, Symbol("maxSurge"), maxSurge)
        validate_property(IoK8sApiAppsV1beta2RollingUpdateDeployment, Symbol("maxUnavailable"), maxUnavailable)
        setfield!(o, Symbol("maxUnavailable"), maxUnavailable)
        o
    end
end # type IoK8sApiAppsV1beta2RollingUpdateDeployment

const _property_map_IoK8sApiAppsV1beta2RollingUpdateDeployment = Dict{Symbol,Symbol}(Symbol("maxSurge")=>Symbol("maxSurge"), Symbol("maxUnavailable")=>Symbol("maxUnavailable"))
const _property_types_IoK8sApiAppsV1beta2RollingUpdateDeployment = Dict{Symbol,String}(Symbol("maxSurge")=>"IoK8sApimachineryPkgUtilIntstrIntOrString", Symbol("maxUnavailable")=>"IoK8sApimachineryPkgUtilIntstrIntOrString")
Base.propertynames(::Type{ IoK8sApiAppsV1beta2RollingUpdateDeployment }) = collect(keys(_property_map_IoK8sApiAppsV1beta2RollingUpdateDeployment))
Swagger.property_type(::Type{ IoK8sApiAppsV1beta2RollingUpdateDeployment }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAppsV1beta2RollingUpdateDeployment[name]))}
Swagger.field_name(::Type{ IoK8sApiAppsV1beta2RollingUpdateDeployment }, property_name::Symbol) =  _property_map_IoK8sApiAppsV1beta2RollingUpdateDeployment[property_name]

function check_required(o::IoK8sApiAppsV1beta2RollingUpdateDeployment)
    true
end

function validate_property(::Type{ IoK8sApiAppsV1beta2RollingUpdateDeployment }, name::Symbol, val)
end


text:
module Wilsonloop
    export make_staple,Wilsonline, make_staple_and_loop,derive_U,make_Cμ,
            make_plaq_staple,make_plaq, loops_staple_prime,get_position,derive_Udag,
            make_loops_fromname,make_chair,get_rightlinks,get_leftlinks,get_direction,
            loops_plaq,loops_rect,check_plaqset,isdag
    using LaTeXStrings
    using LinearAlgebra
    import Base
    import Base.:(==)

    abstract type Gaugelink{Dim} end

    struct GLink{Dim} <: Gaugelink{Dim}
        direction::Int8
        position::NTuple{Dim,Int64}
        isdag::Bool

        function GLink{Dim}(direction,position,isdag=false) where Dim
            return new{Dim}(direction,position,isdag)
        end
    end

    struct Adjoint_GLink{Dim} <: Gaugelink{Dim}
        parent::GLink{Dim}
    end

    function LinearAlgebra.adjoint(glink::Adjoint_GLink{Dim}) where Dim
        return GLink{Dim}(glink.direction,glink.position,true)#, glink.parent
    end

    function LinearAlgebra.adjoint(glink::GLink{Dim}) where Dim
        return GLink{Dim}(glink.direction,glink.position,!(glink.isdag))
#        return Adjoint_GLink{Dim}(glink)
    end

    function get_direction(glink::GLink)
        return glink.direction
    end

    function get_position(glink::GLink)
        return glink.position
    end

    function isdag(glink::GLink)
        return glink.isdag
    end

    function set_position(glink::GLink{Dim},position) where Dim
        return GLink{Dim}(glink.direction,position,glink.isdag)
    end

    mutable struct Wilsonline{Dim}
        glinks::Vector{GLink{Dim}}
        #glinks::Array{Union{GLink{Dim},Adjoint_GLink{Dim}},1}

        Wilsonline(;Dim=4) = new{Dim}([])
        Wilsonline(glinks;Dim=4) = new{Dim}(glinks)
        function Wilsonline(segments_in::Array{Tuple{T,T},1};Dim=4) where T<: Integer
            segments = make_links(segments_in)
            numline = length(segments)
            glinks = Array{GLink{Dim},1}(undef,numline)
            #glinks = Array{Union{GLink{Dim},Adjoint_GLink{Dim}},1}(undef,numline)
            position = zeros(Int64,Dim)
            for (i,segment) in enumerate(segments)
                dimension = segment[1]
                hoppingdirection = segment[2]

                if hoppingdirection == 1
                    glinks[i] =  GLink{Dim}(dimension ,Tuple(position))
                    position[dimension] += 1
                elseif hoppingdirection == -1
                    position[dimension] += -1
                    glinks[i] =  GLink{Dim}(dimension ,Tuple(position))'
                else
                    error("hoppingdirection in segment should be 1 or -1. But now $hoppingdirection")
                end
            end
            return new{Dim}(glinks)
        end
    end

    function ==(x::GLink{Dim},y::GLink{Dim}) where Dim
        if x.isdag != y.isdag
            return false
        end

        if x.direction == y.direction && x.position == y.position
            return true
        else
            return false
        end
    end

    function ==(x::Wilsonline{Dim},y::Wilsonline{Dim}) where Dim
        flag = true
        if length(x) != length(y)
            return false
        end

        for i=1:length(x)
            if x[i] != y[i]
                return false
            end
        end
        return true
    end

    struct DwDU{Dim}
        parent::Wilsonline{Dim}
        insertindex::Int64
        position::NTuple{Dim,Int64}
        leftlinks::Wilsonline{Dim}
        rightlinks::Wilsonline{Dim}
        μ::Int8
    end

    function get_leftlinks(dw::DwDU)
        return dw.leftlinks
    end

    function get_rightlinks(dw::DwDU)
        return dw.rightlinks
    end

    function get_position(dw::DwDU)
        return dw.position
    end

    function Base.push!(w::Wilsonline,link)
        push!(w.glinks,link)
    end

    function Base.append!(w::Wilsonline,a::Wilsonline)
        append!(w.glinks,a.glinks)
    end

    function Base.length(w::Wilsonline)
        return length(w.glinks)
    end

    function Base.getindex(w::Wilsonline,i)
        return w.glinks[i]
    end

    function Base.lastindex(w::Wilsonline)
        return length(w)
    end

    function LinearAlgebra.adjoint(w::Wilsonline{Dim}) where Dim
        wa = Wilsonline(;Dim=Dim)
        numlinks = length(w)
        for i=numlinks:-1:1
            push!(wa,w[i]')
        end
        return wa
    end

    function LinearAlgebra.adjoint(ws::Array{<: Wilsonline{Dim},1}) where Dim
        num = length(ws)
        wad = Array{eltype(ws),1}(undef,num)
        for i=1:num
            wad[i] = ws[i]'
        end
        return wad
    end

    function Base.show(io::IO,ws::Array{<: Wilsonline{Dim},1}) where Dim
        for i=1:length(ws)
            if i==1
                st ="st"
            elseif i==2
                st ="nd"
            elseif i==3
                st ="rd"
            else
                st ="th"
            end
            println(io,"$i-$st loop")
            show(io,ws[i])
            #display(io,ws[i])
        end
    end

    function Base.display(ws::Array{<: Wilsonline{Dim},1}) where Dim
        for i=1:length(ws)
            if i==1
                st ="st"
            elseif i==2
                st ="nd"
            elseif i==3
                st ="rd"
            else
                st ="th"
            end
            println("$i-$st loop")
            display(ws[i])
        end
    end

    function get_printstring_direction(glink::Gaugelink{Dim}) where Dim

        nstring = "n"

        position = get_position(glink)
        for μ=1:Dim
            m = position[μ]
            if m != 0
                if abs(m)==1
                    if m >0
                        nstring = nstring*"+e_{$(μ)}"
                    else
                        nstring = nstring*"-e_{$(μ)}"
                    end
                else
                    if m >0
                        nstring = nstring*"+$(m)e_{$(μ)}"
                    else
                        nstring = nstring*"-$(abs(m))e_{$(μ)}"
                    end
                end
            end
        end
        return nstring
    end

    function get_printstring(glink::Gaugelink{Dim}) where Dim
        direction  = get_direction(glink)
        dagornot = ifelse(typeof(glink) <: GLink,"","^{\\dagger}")
        nstring = get_printstring_direction(glink)
        return "U$(dagornot)_{$(direction)}($(nstring))"
    end

    function Base.show(io::IO,glink::Gaugelink{Dim}) where Dim
        outputstring =get_printstring(glink)
        show(io,latexstring(outputstring))
        return latexstring(outputstring)
    end

    function Base.show(io::IO,w::Wilsonline{Dim}) where Dim
        outputstring = ""
        for (i,glink) in enumerate(w.glinks)
            outputstring = outputstring*get_printstring(glink)
        end
        show(io,latexstring(outputstring))
        println("\t")
        #println(io,outputstring)
        #return latexstring(outputstring)
    end

    function make_staple(w::Wilsonline{Dim},μ) where Dim
        dwdUs = derive_U(w,μ)
        numstaples = length(dwdUs)
        staple = Array{typeof(w),1}(undef,numstaples)
        for i=1:numstaples
            wi =Wilsonline(Dim=Dim)
            append!(wi,get_rightlinks(dwdUs[i]))
            append!(wi,get_leftlinks(dwdUs[i]))
            staple[i] = wi
        end
        return staple
    end

    function make_Cμ(w::Wilsonline{Dim},μ) where Dim
        V1 = make_staple(w,μ)
        V2 = make_staple(w',μ)
        C = eltype(V1)[]
        for i=1:length(V1)
            push!(C,V1[i]')
        end
        for i=1:length(V2)
            push!(C,V2[i]')
        end
        return C
    end

    function make_staple_and_loop(w::Wilsonline{Dim},μ) where Dim
        C = make_staple(w,μ)
        append!(C,make_staple(w',μ))
        numstaple = length(C)
        CUdag = Array{typeof(w),1}(undef,numstaple)
        Udag = GLink{Dim}(μ,(0,0,0,0))'
        for i=1:numstaple
            CUdag[i] = deepcopy(C[i])'
            push!(CUdag[i],Udag)
            #CUdag[i] = CUdag[i]'
        end
        return CUdag
    end

    function check_link(w,μ)
        numlinks = length(w)
        linkindices=Int64[]
        for i=1:numlinks
            link = w[i]
            if link.isdag == false
                #typeof(link) <: GLink
                if link.direction == μ
                    append!(linkindices,i)
                end
            end
        end
        return linkindices
    end

    function check_link_dag(w,μ)
        numlinks = length(w)
        linkindices=Int64[]
        for i=1:numlinks
            link = w[i]
            if link.isdag
            #if typeof(link) <: Adjoint_GLink
                if get_direction(link) == μ
                    append!(linkindices,i)
                end
            end
        end
        return linkindices
    end

    """
        like U U U U -> U U otimes U
    """
    function derive_U(w::Wilsonline{Dim},μ) where Dim
        numlinks = length(w)
        linkindices = check_link(w,μ)
        numstaples = length(linkindices)
        dwdU = Array{DwDU{Dim},1}(undef,numstaples)

        for (i,ith) in enumerate(linkindices)
            #wi =Wilsonline(Dim=Dim)
            rightlinks = Wilsonline(Dim=Dim)
            leftlinks = Wilsonline(Dim=Dim)
            origin = w[ith].position
            position = zero(collect(origin))
            position[w[ith].direction] += 1

            for j=ith+1:numlinks
                link = w[j]

                if link.isdag == false
                #if typeof(link) <: GLink
                    link_rev = set_position(link,Tuple(position))
                    position[get_direction(link)] += 1
                else
                    position[get_direction(link)] += -1
                    link_rev = set_position(link,Tuple(position))
                end
                push!(rightlinks,link_rev)

                #push!(rightlinks,link)
            end

            for j=1:ith-1
                link = w[j]

                position = collect(get_position(link)) .- origin
                link_rev =  set_position(link,Tuple(position))
                push!(leftlinks,link_rev)

                #push!(leftlinks,link)
            end
            dwdU[i] = DwDU{Dim}(w,ith,origin,leftlinks,rightlinks,μ)
            #println("μ = ",μ)
            #display(wi)
        end
        return dwdU
    end

    """
    like U U U U -> U U otimes U
"""
    function derive_Udag(w::Wilsonline{Dim},μ) where Dim
        #error("not yet")
        numlinks = length(w)
        linkindices = check_link_dag(w,μ)
        numstaples = length(linkindices)
        dwdUdag = Array{DwDU{Dim},1}(undef,numstaples)

        for (i,ith) in enumerate(linkindices)
            #wi =Wilsonline(Dim=Dim)
            rightlinks = Wilsonline(Dim=Dim)
            leftlinks = Wilsonline(Dim=Dim)
            origin = get_position(w[ith]) #.position
            position = zero(collect(origin))
            #position[w[ith].direction] += 1

            for j=ith+1:numlinks
                link = w[j]

                push!(rightlinks,link)
            end

            for j=1:ith-1
                link = w[j]

                push!(leftlinks,link)
            end
            dwdUdag[i] = DwDU{Dim}(w,ith,origin,leftlinks,rightlinks,μ)
            #println("μ = ",μ)
            #display(wi)
        end
        return dwdUdag
    end

    function Base.display(dwdU::DwDU{Dim}) where {Dim}
        outputstring = ""
        if length(dwdU.leftlinks.glinks) == 0
            outputstring =outputstring*"I "
        else
            for glink in dwdU.leftlinks.glinks
                outputstring =outputstring*get_printstring(glink)
            end
        end

        outputstring = outputstring*" \\otimes "

        if length(dwdU.rightlinks.glinks) == 0
            outputstring =outputstring*"I "
        else
            for glink in dwdU.rightlinks.glinks
                outputstring =outputstring*get_printstring(glink)
            end
        end

        nstring = get_printstring_direction(dwdU.parent.glinks[dwdU.insertindex])

        outputstring = outputstring*"\\delta_{m,$(nstring)}"

        println(outputstring)
        return outputstring
    end

    function Base.show(io::IO,dwdU::Array{DwDU{Dim},1}) where {Dim}
        for i=1:length(dwdU)
            if i==1
                st ="st"
            elseif i==2
                st ="nd"
            elseif i==3
                st ="rd"
            else
                st ="th"
            end
            println("$i-$st loop")
            show(dwdU[i])

        end
    end

    function Base.show(io::IO,dwdU::DwDU{Dim}) where {Dim}
        outputstring = ""
        if length(dwdU.leftlinks.glinks) == 0
            outputstring =outputstring*"I "
        else
            for glink in dwdU.leftlinks.glinks
                outputstring =outputstring*get_printstring(glink)
            end
        end

        outputstring = outputstring*" \\otimes "

        if length(dwdU.rightlinks.glinks) == 0
            outputstring =outputstring*"I "
        else
            for glink in dwdU.rightlinks.glinks
                outputstring =outputstring*get_printstring(glink)
            end
        end

        nstring = get_printstring_direction(dwdU.parent.glinks[dwdU.insertindex])

        outputstring = outputstring*"\\delta_{m,$(nstring)}"

        show(io,latexstring(outputstring))
        println("\t")
        #println(outputstring)
        return outputstring
    end

    function make_links(segments::Array{Tuple{T,T},1}) where T <: Integer
        links = Tuple{Int8,Int8}[]
        for segment in segments
            s = sign(segment[2])
            if segment[2] == 0
                push!(links,(segment[1],0))
            else
            #@assert segment[2] != 0
                for i=1:abs(segment[2])
                    push!(links,(segment[1],s*1))
                end
            end
        end
        return links
    end

    function make_plaq(μ,ν;Dim=4)
        return Wilsonline([(μ,1),(ν,1),(μ,-1),(ν,-1)],Dim = Dim)
    end

    function make_plaq(;Dim=4)
        loops = Wilsonline{Dim}[]
        for μ=1:Dim
            #for ν=1:4
            for ν=μ:Dim
                if ν == μ
                    continue
                end

                plaq = make_plaq(μ,ν,Dim=Dim)
                push!(loops,plaq)
            end
        end
        return loops
    end

    function construct_plaq()
        loops_plaq = Dict{Tuple{Int8,Int8,Int8},Any}()
        for Dim=1:4
            for μ=1:Dim
                for ν=μ:Dim
                    if μ == ν
                        continue
                    end
                    loops_plaq[(Dim,μ,ν)] = make_plaq(μ,ν,Dim=Dim)
                end
            end
        end
        return loops_plaq
    end

    function make_rect(;Dim=4)
        loops = Wilsonline{Dim}[]
        for μ=1:Dim
            for ν=μ:Dim
                if ν == μ
                    continue
                end
                #loop = make_links([(μ,1),(ν,2),(μ,-1),(ν,-2)])

                loop1 = Wilsonline([(μ,1),(ν,2),(μ,-1),(ν,-2)],Dim = Dim)
                #loop1 = Wilson_loop([(μ,1),(ν,2),(μ,-1),(ν,-2)])
                #loop1 = Wilson_loop(loop,Tuple(origin))
                push!(loops,loop1)
                loop1 = Wilsonline([(μ,2),(ν,1),(μ,-2),(ν,-1)],Dim = Dim)
                #loop1 = Wilson_loop([(μ,2),(ν,1),(μ,-2),(ν,-1)])
                push!(loops,loop1)
            end
        end
        return loops
    end

    function construct_rect()
        loops_rect = Dict{Tuple{Int8,Int8,Int8,Int8},Any}()
        for Dim=1:4
            for μ=1:Dim
                for ν=μ:Dim
                    if μ == ν
                        continue
                    end
                    loops_rect[(Dim,μ,ν,1)] = Wilsonline([(μ,1),(ν,2),(μ,-1),(ν,-2)],Dim = Dim)
                    loops_rect[(Dim,μ,ν,2)] = Wilsonline([(μ,2),(ν,1),(μ,-2),(ν,-1)],Dim = Dim)
                end
            end
        end
        return loops_rect
    end

    function make_cloverloops(μ,ν;Dim=4)
        loops = Wilsonline{Dim}[]
        loop_righttop = Wilsonline([(μ,1),(ν,1),(μ,-1),(ν,-1)])
        loop_lefttop = Wilsonline([(ν,1),(μ,-1),(ν,-1),(μ,1)])
        loop_rightbottom = Wilsonline([(ν,-1),(μ,1),(ν,1),(μ,-1)])
        loop_leftbottom= Wilsonline([(μ,-1),(ν,-1),(μ,1),(ν,1)])
        push!(loops,loop_righttop)
        push!(loops,loop_lefttop)
        push!(loops,loop_rightbottom)
        push!(loops,loop_leftbottom)
        return loops
    end

    function make_polyakov(μ,Lμ;Dim=4)
        loops = Wilsonline{Dim}[]
        loop1 = Wilsonline([(μ,Lμ)],Dim = Dim)
        push!(loops,loop1)
        return loops
    end

    function make_polyakov_xyz(Lμ;Dim=4)
        loops = Wilsonline{Dim}[]
        for μ=1:3
            loop1 = Wilsonline([(μ,Lμ)],Dim = Dim)
            push!(loops,loop1)
        end
        return loops
    end

    function make_loopforactions(couplinglist,L)
        Dim = length(L)
        loops = Array{Array{Wilsonline{Dim},1},1}(undef,length(couplinglist))
        for (i,name) in enumerate(couplinglist)
            if  name == "plaquette"
                loops[i] = make_plaq(Dim = Dim)
            elseif name == "rectangular"
                loops[i] = make_rect(Dim = Dim)
            elseif name =="chair"
                loops[i] = make_chair(Dim = Dim)
            elseif name == "polyakov_t"
                μ= Dim
                loops[i] = make_polyakov(μ,L[μ],Dim = Dim)
            elseif name == "polyakov_z"
                @assert Dim > 3 "Dimension should be Dim > 3 but now Dim = $Dim"
                μ= 3
                loops[i] = make_polyakov(μ,L[μ],Dim = Dim)
            elseif name == "polyakov_y"
                @assert Dim > 2 "Dimension should be Dim > 2 but now Dim = $Dim"
                μ= 2
                loops[i] = make_polyakov(μ,L[μ],Dim = Dim)
            elseif name == "polyakov_x"
                μ= 1
                loops[i] = make_polyakov(μ,L[μ],Dim = Dim)
            else
                error("$name is not supported!")
            end
        end
        return loops
    end

    function make_loops_fromname(name;Dim=4,L=nothing)
        if L != nothing
            @assert Dim == length(L)
        end

        if  name == "plaquette"
            loops = make_plaq(Dim = Dim)
        elseif name == "rectangular"
            loops = make_rect(Dim = Dim)
        elseif name =="chair"
            loops = make_chair(Dim = Dim)
        elseif name == "polyakov_t"
            @assert L != nothing "system size should be given to obtain polyakov loops. please do like make_loops(\"polyakov_t\";Dim=4,L=[4,4,4,4])"
            μ= Dim
            loops = make_polyakov(μ,L[μ],Dim = Dim)
        elseif name == "polyakov_z"
            @assert L != nothing "system size should be given to obtain polyakov loops. please do like make_loops(\"polyakov_z\";Dim=4,L=[4,4,4,4])"
            @assert Dim > 3 "Dimension should be Dim > 3 but now Dim = $Dim"
            μ= 3
            loops = make_polyakov(μ,L[μ],Dim = Dim)
        elseif name == "polyakov_y"
            @assert L != nothing "system size should be given to obtain polyakov loops. please do like make_loops(\"polyakov_y\";Dim=4,L=[4,4,4,4])"
            @assert Dim > 2 "Dimension should be Dim > 2 but now Dim = $Dim"
            μ= 2
            loops = make_polyakov(μ,L[μ],Dim = Dim)
        elseif name == "polyakov_x"
            @assert L != nothing "system size should be given to obtain polyakov loops. please do like make_loops(\"polyakov_x\";Dim=4,L=[4,4,4,4])"
            μ= 1
            loops = make_polyakov(μ,L[μ],Dim = Dim)
        else
            error("$name is not supported!")
        end

        return loops
    end

    function make_chair(;Dim=4)
        @assert Dim == 4 "only Dim = 4 is supported now"
        #loopset = []
        loopset = Wilsonline{Dim}[]
        set1 = (1,2,3)
        set2 = (1,2,4)
        set3 = (2,3,4)
        set4 = (1,3,4)

        for set in (set1,set2,set3,set4)
            mu,nu,rho = set
            origin = zeros(Int8,4)
            loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
            loop1 =  Wilsonline(loop,Dim = Dim)
            push!(loopset,loop1)

            mu,rho,nu = set
            loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
            loop1 =  Wilsonline(loop,Dim = Dim)
            push!(loopset,loop1)

            nu,rho,mu = set
            loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
            loop1 =  Wilsonline(loop,Dim = Dim)
            push!(loopset,loop1)

            nu,mu,rho = set
            loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
            loop1 =  Wilsonline(loop,Dim = Dim)
            push!(loopset,loop1)

            rho,mu,nu = set
            loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
            loop1 =  Wilsonline(loop,Dim = Dim)
            push!(loopset,loop1)

            rho,nu,mu = set
            loop = [(mu,1),(nu,1),(rho,1),(mu,-1),(rho,-1),(nu,-1)]
            loop1 =  Wilsonline(loop,Dim = Dim)
            push!(loopset,loop1)

        end
        return loopset
    end

    function make_plaq_staple(μ;Dim=4)
        loops = Wilsonline{Dim}[]
        plaqs = make_plaq(Dim = Dim)
        numplaq = length(plaqs)
        for i=1:numplaq
            plaq = plaqs[i]
            staples = make_staple(plaq,μ)
            for j=1:length(staples)
                push!(loops,staples[j])
            end

            plaqdag = plaqs[i]'
            staples = make_staple(plaqdag,μ)
            for j=1:length(staples)
                push!(loops,staples[j])
            end
        end

        return loops
    end

    function construct_staple_prime()
        loops_staple_prime = Dict{Tuple{Int8,Int8},Any}()
        for Dim=1:4
            for μ=1:Dim
                loops_staple_prime[(Dim,μ)] = make_plaq_staple(μ,Dim=Dim)'
            end
        end
        return loops_staple_prime
    end

    function construct_staple()
        loops_staple = Dict{Tuple{Int8,Int8},Any}()
        for Dim=1:4
            for μ=1:Dim
                loops_staple[(Dim,μ)] = make_plaq_staple(μ,Dim=Dim)
            end
        end
        return loops_staple
    end

    const loops_staple_prime = construct_staple_prime()
    const loops_staple = construct_staple()
    const loops_plaq = construct_plaq()
    const loops_rect = construct_rect()

    function check_plaqset(wi::Wilsonline{Dim}) where Dim
        flag = false
        for μ=1:Dim
            for ν=μ:Dim
                if μ==ν
                    continue
                end
                loop= loops_plaq[(Dim,μ,ν)]
                if loop == wi
                    #println("match!")
                    flag = true
                    break
                end
            end
            if flag
                continue
            end
        end

        return flag
    end

    function check_plaqset(w::Vector{Wilsonline{Dim}}) where Dim
        flag = false
        for wi in w
            flag = check_plaqset(wi)
            if flag != true
                return false
            end
        end

        return flag
    end

    function check_rectset(wi::Wilsonline{Dim}) where Dim
        flag = false
        direction = (0,0)
        for μ=1:Dim
            for ν=μ:Dim
                if μ==ν
                    continue
                end
                loop= loops_rect[(Dim,μ,ν,1)]
                if loop == wi
                    #println("match!")
                    direction = (μ,ν)
                    flag = true
                    break
                end

                loop= loops_rect[(Dim,μ,ν,2)]
                if loop == wi
                    #println("match!")
                    direction = (μ,ν)
                    flag = true
                    break
                end
            end
            if flag
                continue
            end
        end

        return flag,direction
    end

    function check_rectset(w::Vector{Wilsonline{Dim}}) where Dim
        flag = false
        direction = (0,0)
        for wi in w
            flag,direction = check_rectset(wi)
            if flag != true
                return false,direction
            end
        end

        return flag,direction
    end

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    DatamineAngles(α, β, θ)

Datamine ZXY rotation convention following the left-hand rule.
All angles are in degrees and the signal convention is CW, CW, CW
positive. Y is the principal axis.
"""
DatamineAngles(α, β, θ) = EulerAngles(-deg2rad(α - 90), deg2rad(β), -deg2rad(θ), :ZYX)

"""
    GslibAngles(α, β, θ)

GSLIB ZXY rotation convention following the left-hand rule.
All angles are in degrees and the signal convention is CW, CCW, CCW
positive. Y is the principal axis.

## Reference
* Deutsch, 2015. [The Angle Specification for GSLIB Software]
(https://geostatisticslessons.com/lessons/anglespecification)
"""
GslibAngles(α, β, θ) = EulerAngles(-deg2rad(α - 90), -deg2rad(β), deg2rad(θ), :ZYX)

"""
    VulcanAngles(α, β, θ)

GSLIB ZYX rotation convention following the right-hand rule.
All angles are in degrees and the signal convention is CW, CCW, CW
positive. X is the principal axis.

## Reference
* Deutsch, 2015. [The Angle Specification for GSLIB Software]
(https://geostatisticslessons.com/lessons/anglespecification)
"""
VulcanAngles(α, β, θ) = EulerAngles(-deg2rad(α - 90), -deg2rad(β), -deg2rad(θ), :ZYX)


text:
using Nature
using Serialization

pols = readdir("checkpoints/")
for (i, pol) in enumerate(pols)
    println(i, ". ", pol)
end
print("Select a run: ")
pol_name = "checkpoints/"*pols[parse(Int, readline())]

pols = readdir("$pol_name/")
for (i, pol) in enumerate(pols)
    println(i, ". ", pol)
end
print("Select a checkpoint: ")
pol_name = "$pol_name/"*pols[parse(Int, readline())]

p = deserialize(pol_name)[:policy]

WORLD_SIZE=(2,2)
env = NatureEnv(num_starting_players=2,
                world_size=WORLD_SIZE,
                window = 3,
                num_frames=2,
                episode_len=1000,
                vocab_size=5,
                food_generators=[
                    FoodGen([1,1],[0.0001,0.0001],max_pos=WORLD_SIZE),
                    FoodGen([2,2],[0.0001,0.0001],max_pos=WORLD_SIZE),
                   ])
step_through_env(env, p)


text:
module CUSPARSE

importall Base.SparseArrays

using CUDArt
using CUDArt.rt.cudaStream_t

const SparseChar = Char
import Base.one
import Base.zero

export CudaSparseMatrixCSC, CudaSparseMatrixCSR,
       CudaSparseMatrixHYB, CudaSparseMatrixBSR,
       CudaSparseMatrix, AbstractCudaSparseMatrix,
       CudaSparseVector

include("util.jl")
include("libcusparse_types.jl")

function statusmessage( status )
    if status == CUSPARSE_STATUS_SUCCESS
        return "cusparse success"
    end
    if status == CUSPARSE_STATUS_NOT_INITIALIZED
        return "cusparse not initialized"
    end
    if status == CUSPARSE_STATUS_ALLOC_FAILED
        return "cusparse allocation failed"
    end
    if status == CUSPARSE_STATUS_INVALID_VALUE
        return "cusparse invalid value"
    end
    if status == CUSPARSE_STATUS_ARCH_MISMATCH
        return "cusparse architecture mismatch"
    end
    if status == CUSPARSE_STATUS_MAPPING_ERROR
        return "cusparse mapping error"
    end
    if status == CUSPARSE_STATUS_EXECUTION_FAILED
        return "cusparse execution failed"
    end
    if status == CUSPARSE_STATUS_INTERNAL_ERROR
        return "cusparse internal error"
    end
    if status == CUSPARSE_STATUS_MATRIX_TYPE_NOT_SUPPORTED
        return "cusparse matrix type not supported"
    end
end

"""
Determine if a CUSPARSE call was successful.
If not, find the appropriate message to throw.
"""
function statuscheck(status)
    if status == CUSPARSE_STATUS_SUCCESS
        return nothing
    end
    warn("CUSPARSE error triggered from:")
    Base.show_backtrace(STDOUT, backtrace())
    println()
    throw(statusmessage( status ))
end

cuda_versions = ["80", "75", "70", "65"]
lib_list = vcat(["libcusparse", "cusparse"], "cusparse64_" .* cuda_versions, "cusparse32_" .* cuda_versions)
const libcusparse = Libdl.find_library(lib_list, ["/usr/lib/", "/usr/local/cuda", ""])

if isempty(libcusparse)
    error("CUSPARSE library not found!")
end

include("libcusparse.jl")

#setup handler for cusparse

cusparsehandle = cusparseHandle_t[0]
cusparseCreate( cusparsehandle )

#clean up handle at exit
atexit( ()->cusparseDestroy(cusparsehandle[1]) )

include("sparse.jl")

end


text:
# This file is a part of RadiationSpectra.jl, licensed under the MIT License (MIT).
using Test

# ENV["JULIA_DEBUG"] = "BAT"
# using BAT
using RadiationSpectra
using IntervalSets

@testset "Package RadiationSpectra" begin
    T = Float64

    h_uncal = RadiationSpectra.get_example_spectrum()

    gamma_lines = T[609.312, 911.204, 1120.287, 1460.830, 1764.494, 2614.533]
    h_cal, h_deconv, peakPositions, threshold, c, c_precal = RadiationSpectra.calibrate_spectrum(h_uncal, gamma_lines, min_n_peaks = 30, σ = 2.0 )

    @info "Calibration constant: c = $(c)"
    c_true = 0.011259105696794384
    @testset "Auto calibration" begin
        @test abs(c - c_true) / c_true < 0.001
    end

    ff_lsq = FitFunction(T, :GaussPlusLinearBackground)

    set_fitranges!(ff_lsq, ((1461 - 20, 1461 + 20),))
    set_initial_parameters!(ff_lsq, T[10000, 1.4, 1461, 20, 0])
    set_parameter_bounds!(ff_lsq, [0.0..10^7, 0.2..2.0, 1455..1465, 0..10^4, -30..30])

    lsqfit!(ff_lsq, h_cal)

    fitted_pars = collect(get_fitted_parameters(ff_lsq))
    @testset "LSQ Fit" begin
        @test abs(fitted_pars[3] - 1460.830) <= 3
    end

    ff_llh = FitFunction(T, :GaussPlusLinearBackground)
    set_fitranges!(ff_llh, ((1461 - 20, 1461 + 20),))
    fitted_pars[1] = 0.9 * fitted_pars[1]

    set_initial_parameters!(ff_llh, fitted_pars)
    llhfit!(ff_llh, h_cal)

    fitted_pars = get_fitted_parameters(ff_llh)
    @show ff_llh
    @testset "LLH Fit" begin
        @test abs(fitted_pars[3] - 1460.830) <= 3
    end

    # if @isdefined BAT
    #     ff_bat = FitFunction(T, :GaussPlusLinearBackground)
    #     set_fitranges!(ff_bat, ((1461 - 20, 1461 + 20),))
    #     set_parameter_bounds!(ff_bat, [0.0..10^7, 0.2..2.0, 1455..1465, 0..10^4, -30..30])

    #     batfit!(ff_bat, h_cal)
    #     # get_histogram_pdf(ff_bat, nbins = 3)
    #     # get_standard_deviations(ff_bat)
    #     # get_samples_inds_by_chain_id(ff_bat, 1)

    #     fitted_pars = get_fitted_parameters(ff_bat)
    #     @show ff_bat
    #     @testset "BAT Fit" begin
    #         @test abs(fitted_pars[3] - 1460.830) <= 3
    #     end
    # end

    @testset "General model functions" begin
        ff = FitFunction(T, :Gauss)
        RadiationSpectra._set_fitted_parameters!(ff, [1, 1, 0])
        ff.model(T[0, 1], collect(get_fitted_parameters(ff)))
        println(ff)
        @test round.(ff.model(T[0, 1], collect(get_fitted_parameters(ff))), digits = 5) == round.(T[0.3989422804014327, 0.24197072451914337], digits = 5)

        ff = FitFunction(:PDF_Gauss)
        ff = FitFunction(T, :PDF_Gauss)
        RadiationSpectra._set_fitted_parameters!(ff, [1, 0])
        ff.model(T[0, 1], collect(get_fitted_parameters(ff)))
        println(ff)
        @test round.(ff.model(T[0, 1], collect(get_fitted_parameters(ff))), digits = 5) == round.(T[0.3989422804014327, 0.24197072451914337], digits = 5)

        ff = FitFunction(T, :GaussPlusConstantBackground)
        RadiationSpectra._set_fitted_parameters!(ff, [1, 1, 0, 0])
        ff.model(T[0, 1], collect(get_fitted_parameters(ff)))
        println(ff)
        @test round.(ff.model(T[0, 1], collect(get_fitted_parameters(ff))), digits = 5) == round.(T[0.3989422804014327, 0.24197072451914337], digits = 5)

        ff = FitFunction(T, :GaussPlusLinearBackground)
        RadiationSpectra._set_fitted_parameters!(ff, [1, 1, 0, 0, 0])
        ff.model(T[0, 1], collect(get_fitted_parameters(ff)))
        println(ff)
        @test round.(ff.model(T[0, 1], collect(get_fitted_parameters(ff))), digits = 5) == round.(T[0.3989422804014327, 0.24197072451914337], digits = 5)
    end

end # testset


text:
"""Julia program to implement Insertion sort
In Insertion Sort we maintain a sorted part and an unsorted part,
and continously shift elements from unsorted part to the sorted part
thus getting the whole array sorted in the process.
"""

function insertion_sort(arr, n)
    # Begin the iteration from the second element considering, the first element as a part of the sorted part
    for i in 2:n
        temp = arr[i]
        j = i - 1
        # Find the first element from the beginning of the array that is smaller than the temp element
        # And shift all the elements till then towards right.
        while j > 0 && temp < arr[j]
            arr[j + 1] = arr[j]
            j = j - 1
        end
        arr[j + 1] = temp
    end
end

print("How many numbers do you want to sort? ")
n = readline()
n = parse(Int64, n)
if (n <= 0)
    println("No numbers to sort!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
insertion_sort(arr, n)
print("The numbers in sorted order is ")
print(arr)

"""
Time Complexity - O(n^2), where 'n' is the size of the array
Space Complexity - O(n)

SAMPLE INPUT AND OUTPUT

SAMPLE I
How many numbers do you want to sort? 6
Enter the numbers: 76 243 89 25 8337 72
The numbers in sorted order is [25, 72, 76, 89, 243, 8337]

"""


text:
#=
surface_to_obj:
- Julia version: 1.1.0
- Author: Jirik
- Date: 2019-04-10
=#

# include("../src/LarSurf.jl")
tim = time()
# using Revise
using LarSurf
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using Plasm, SparseArrays
using Pandas
using Seaborn

tim_prev = tim
tim = time()
println("using done in: ", tim - tim_prev)
threshold = 10

xystep = 10
zstep = 5
xystep = 20
zstep = 10
pth = Pio3d.datasets_join_path("medical/orig/3Dircadb1.1/MASKS_DICOM/liver")

datap = Pio3d.read3d(pth);
#
data3d_full = datap["data3d"]
println("orig size: ", size(data3d_full))
data3d = data3d_full[1:zstep:end, 1:xystep:end, 1:xystep:end];

data_size = LarSurf.size_as_array(size(data3d))
println("data size: ", data_size)
segmentation = data3d .> threshold;

tim_prev = tim
tim = time()
println("data read complete in time: ", tim - tim_prev)

# Run once to force compilation
## Generate data
# segmentation = LarSurf.generate_slope([9,10,11])
block_size = [5,5,5]
reduced_model, Flin, (bigV, tmodel) = LarSurf.get_surface_grid_per_block(segmentation, block_size)
bigVV, bigEV, bigFV, bigCV = tmodel
filtered_bigFV = reduced_model[2][1]

V = bigV
FV = filtered_bigFV
nV = size(bigV)[2]
# First method for computation of EV
EV = []
for f in FV
	push!(EV, [[f[1],f[2]],[f[3],f[4]],  [f[1],f[3]],[f[2],f[4]]])
end
doubleedges = sort(cat(EV))
doubleedges = convert(Lar.Cells, doubleedges)
EV = [doubleedges[k] for k=1:2:length(doubleedges)]

# Plasm.view(V,EV)

# Computing copFE
# kEV = Lar.characteristicMatrix(EV);
aFV = LarSurf.ll2array(FV)
kFV = LarSurf.characteristicMatrix(aFV, nV);

aEV = LarSurf.ll2array(EV)
kEV = LarSurf.characteristicMatrix(aEV, nV);
kFE = kFV * kEV'
I,J,Value = SparseArrays.findnz(kFE)
triples = hcat([[i,j,1] for (i,j,v) in zip(I,J,Value) if v == 2 ]...)
I,J,Value = [triples[k,:] for k=1:size(triples,1)]
copFE = SparseArrays.sparse(I,J,Value)
copFE = convert(SparseMatrixCSC{Int8,Int64},copFE)
FE = [findnz(copFE[k,:])[1] for k=1:length(FV)]
copFE = Lar.coboundary_1(V,kFV,kEV)

# Assembling the chain complex cc
copCF = SparseArrays.ones(Int8,1,size(kFV,1))
copCF = convert(Lar.ChainOp, copCF)
cc = [kEV, copFE, copCF]::Lar.ChainComplex

# View a surface triangulation
V = convert(Lar.Points,V')
triangles = cat(Lar.triangulate(V::Lar.Points, cc[1:2]))
TV = convert(Lar.Cells,triangles)
V = convert(Lar.Points,V')
Plasm.view(Plasm.hpc_exploded( (V,[TV]) )(1.2,1.2,1.2))

# Export an OBJ file of the surface
V = convert(Lar.Points,V')
open("testfile.obj","w") do f
    print(f, Lar.lar2obj(V::Lar.Points, cc) )
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

_h5io_keys(df::HDF5.H5DataStore) = keys(df)
_h5io_keys(df::HDF5.Dataset) = nothing

_h5io_objtype(df::HDF5.H5DataStore) = Val(:datafile)
_h5io_objtype(df::HDF5.Dataset) = Val(:dataset)

function _h5io_open(body::Function, filename::AbstractString, mode::AbstractString)
    HDF5.h5open(filename, mode) do f
        body(f)
    end
end

function _h5io_open(body::Function, fn_with_subpath::Tuple{AbstractString, AbstractString}, mode::AbstractString)
    filename, subpath = fn_with_subpath
    HDF5.h5open(filename, mode) do f
        body((f, subpath))
    end
end

"""
    bat_read(src::HDF5.H5DataStore)
    bat_read(src_with_subpath::Tuple{HDF5.H5DataStore, AbstractString})

Read data from HDF5 file or group `src` (optionally from an HDF5-path
relative to `src`).
"""
function bat_read(dest)
    _h5io_read(dest)
end

"""
    bat_write(dest::HDF5.H5DataStore, data)
    bat_write(dest_with_subpath::Tuple{HDF5.H5DataStore, AbstractString}, data)

Write `data` to HDF5 file or group `dest` (optionally to an HDF5-path
relative to `dest`).

`data` must be a table (i.e. implement the Tables.jl API).
"""
function bat_write(dest, data)
    _h5io_write!(dest, data)
end


text:
export CuSubArray
export view

type CuSubArray{T,N} <: AbstractCuArray{T,N}
    parent::CuArray{T}
    indexes::Tuple
    offset::Int
    dims::NTuple{N,Int}
    strides::NTuple{N,Int}
end

typealias CuSubVector{T} CuArray{T,1}
typealias CuSubMatrix{T} CuArray{T,2}
typealias CuSubVecOrMat{T} Union{CuSubVector{T}, CuSubMatrix{T}}

Base.size(a::CuSubArray) = a.dims
Base.size(a::CuSubArray, dim::Int) = a.dims[dim]

Base.strides(a::CuSubArray) = a.strides
Base.strides(a::CuSubArray, dim::Int) = a.strides[dim]

Base.length(a::CuSubArray) = prod(a.dims)

Base.similar(a::CuSubArray) = similar(a, size(a))
Base.similar(a::CuSubArray, dims::Int...) = similar(a, dims)
Base.similar{T,N}(a::CuSubArray{T}, dims::NTuple{N,Int}) = CuArray(T, dims)

function Base.pointer(a::CuSubArray, index::Int=1)
    index == 1 && return pointer(a.parent, a.offset+index)
    throw("Not implemented yet.")
end

function cupointer(a::CuSubArray, index::Int=1)
    index == 1 && return cupointer(a.parent, a.offset+index)
    throw("Not implemented yet.")
end

function view{T,N}(a::CuArray{T,N}, indexes::Union{UnitRange{Int},Colon,Int}...)
    dims = Int[]
    strides = Int[]
    stride = 1
    offset = 0
    for i = 1:length(indexes)
        r = indexes[i]
        t = typeof(r)
        if t == Colon
            push!(dims, size(a,i))
            push!(strides, stride)
        elseif t == Int
            offset += stride * (r-1)
        else
            push!(dims, length(r))
            push!(strides, stride)
            offset += stride * (first(r)-1)
        end
        stride *= size(a,i)
    end
    CuSubArray(a, indexes, offset, tuple(dims...), tuple(strides...))
end

Base.Array(a::CuSubArray) = Array(a.parent)[a.indexes...]
Base.SubArray(a::CuSubArray) = view(Array(a.parent), a.indexes)


text:
"""
Piecewise linear functions, stored in terms of an array of kinks (i.e. events).
"""
module Pwl

include("event.jl")
include("queues.jl")
include("clip.jl")

end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
module TestDatasets

using FileIO
using ColorTypes: Gray
using Pkg.Artifacts
using StringDistances

const artifacts_toml = abspath(joinpath(@__DIR__, "..", "Artifacts.toml"))

export testdataset

const remotedatasets = [
    "cameraman512_10",
    "cameraman256_10",
    "cameraman256_5",
    "cameraman128_5",
    "faces_5",
    "faces_10",
    "matches_5",
    "matches_10",
    "mandrill",
    "peppers",
    "playing_cards_5"
]

function testdataset(datasetname; download_only::Bool = false)

    datasetfiles = dataset_path(full_datasetname(datasetname))
    download_only && return datasetfiles
    dataset = load_dataset(datasetfiles)
    return dataset

end

function full_datasetname(datasetname)
    idx = findfirst(remotedatasets) do x
        startswith(x, datasetname)
    end
    if idx === nothing
        warn_msg = "\"$datasetname\" not found in `TestDatasets.remotedatasets`."

        best_match = _findnearest(datasetname)
        if isnothing(best_match[2])
            similar_matches = remotedatasets[_findall(datasetname)]
            if !isempty(similar_matches)
                similar_matches_msg = "  * \"" * join(similar_matches, "\"\n  * \"") * "\""
                warn_msg = "$(warn_msg) Do you mean one of the following?\n$(similar_matches_msg)"
            end
            throw(ArgumentError(warn_msg))
        else
            idx = best_match[2]
            @warn "$(warn_msg) Load \"$(remotedatasets[idx])\" instead."
        end
    end
    return remotedatasets[idx]
end

function dataset_path(datasetname)
    ensure_artifact_installed("datasets", artifacts_toml)

    dataset_dir = artifact_path(artifact_hash("datasets", artifacts_toml))
    return joinpath(dataset_dir, "datasets" ,datasetname)
end

function load_dataset(datasetfiles)
    image_pairs = readlines(joinpath(datasetfiles,"filelist.txt"))
    M,N = size(load(joinpath(datasetfiles,split(image_pairs[1],",")[1])))
    true_images = zeros(M,N,length(image_pairs))
    data_images = zeros(M,N,length(image_pairs))
    for i = 1:length(image_pairs)
        pair = split(image_pairs[i],",")
        true_images[:,:,i] = load(joinpath(datasetfiles,pair[1]))
        data_images[:,:,i] = load(joinpath(datasetfiles,pair[2]))
    end
    return true_images,data_images
end

_findall(name; min_score=0.6) = findall(name, remotefiles,JaroWinkler(), min_score=min_score)
_findnearest(name; min_score=0.8) = findnearest(name, remotefiles, JaroWinkler(), min_score=min_score)

end

max_stars_repo_name: dvillacis/TestDatasets

max_issues_repo_name: dvillacis/TestDatasets

max_forks_repo_name: dvillacis/TestDatasets


text:
@testset "Bloom Filter" begin
    @testset "Add" begin
        bf = BloomFilter(10, 5, 99)
        item = UInt8.(collect("Hello World"))
        Bitcoin.add!(bf, item)
        expected = "0000000a080000000140"
        @test bytes2hex(Bitcoin.filter_bytes(bf)) == expected
        item = UInt8.(collect("Goodbye!"))
        Bitcoin.add!(bf, item)
        expected = "4000600a080000010940"
        @test bytes2hex(Bitcoin.filter_bytes(bf)) == expected
    end
    @testset "Load" begin
        bf = BloomFilter(10, 5, 99)
        item = UInt8.(collect("Hello World"))
        Bitcoin.add!(bf, item)
        item = UInt8.(collect("Goodbye!"))
        Bitcoin.add!(bf, item)
        expected = "0a4000600a080000010940050000006300000001"
        @test bytes2hex(Bitcoin.serialize(Bitcoin.FilterLoadMessage(bf))) == expected
    end
end


text:
using CuTextures
using CuArrays
using CUDAnative

using Test

@inline function calcpoint(blockIdx, blockDim, threadIdx, size)
    i = (blockIdx - 1) * blockDim + threadIdx
    return i, Float32(i)
end
function kernel_texture_warp_native(dst::CuDeviceArray{T,1}, texture::CuDeviceTexture{T,1}) where {T}
    i, u = calcpoint(blockIdx().x, blockDim().x, threadIdx().x, size(dst)[1])
    @inbounds dst[i] = texture[u];
    return nothing
end
function kernel_texture_warp_native(dst::CuDeviceArray{T,2}, texture::CuDeviceTexture{T,2}) where {T}
    i, u = calcpoint(blockIdx().x, blockDim().x, threadIdx().x, size(dst)[1])
    j, v = calcpoint(blockIdx().y, blockDim().y, threadIdx().y, size(dst)[2])
    @inbounds dst[i,j] = texture[u,v];
    return nothing
end
function kernel_texture_warp_native(dst::CuDeviceArray{T,3}, texture::CuDeviceTexture{T,3}) where {T}
    i, u = calcpoint(blockIdx().x, blockDim().x, threadIdx().x, size(dst)[1])
    j, v = calcpoint(blockIdx().y, blockDim().y, threadIdx().y, size(dst)[2])
    k, w = calcpoint(blockIdx().z, blockDim().z, threadIdx().z, size(dst)[3])
    @inbounds dst[i,j,k] = texture[u,v,w];
    return nothing
end

function fetch_all(texture)
    dims = size(texture)
    d_out = CuArray{eltype(texture)}(undef, dims...)
    # @device_code_warntype @cuda threads = dims kernel_texture_warp_native(d_out, texture)
    # @device_code_llvm @cuda threads = dims kernel_texture_warp_native(d_out, texture)
    @cuda threads = dims kernel_texture_warp_native(d_out, texture)
    d_out
end

@testset "Using CuTextureArray initialized from device" begin
    testheight, testwidth, testdepth = 16, 16, 4
    a1D = convert(Array{Float32}, 1:testheight)
    a2D = convert(Array{Float32}, repeat(1:testheight, 1, testwidth) + repeat(0.01 * (1:testwidth)', testheight, 1))
    a3D = convert(Array{Float32}, repeat(a2D, 1, 1, testdepth))
    for k = 1:testdepth; a3D[:,:,k] .+= 0.0001 * k; end
    d_a1D = CuArray(a1D)
    d_a2D = CuArray(a2D)
    d_a3D = CuArray(a3D)

    texarr1D = CuTextureArray{Float32}(testheight)
    copyto!(texarr1D, d_a1D)
    tex1D = CuTexture(texarr1D)
    @test fetch_all(tex1D) == d_a1D

    texarr2D = CuTextureArray{Float32}(testheight, testwidth)
    tex2D = CuTexture(texarr2D)
    copyto!(texarr2D, d_a2D)
    @test fetch_all(tex2D) == d_a2D

    texarr2D_dir = CuTextureArray(d_a2D)
    tex2D_dir = CuTexture(texarr2D_dir)
    @test fetch_all(tex2D_dir) == d_a2D

    texarr3D = CuTextureArray{Float32}(testheight, testwidth, testdepth)
    tex3D = CuTexture(texarr3D)
    copyto!(texarr3D, d_a3D)
    @test fetch_all(tex3D) == d_a3D

    tex3D_direct = CuTexture{Float32}(testheight, testwidth, testdepth)
    copyto!(tex3D_direct.mem, d_a3D)
    @test fetch_all(tex3D_direct) == d_a3D
end

@testset "Using CuTextureArray initialized from host" begin
    testheight, testwidth, testdepth = 16, 16, 4
    a1D = convert(Array{Float32}, 1:testheight)
    a2D = convert(Array{Float32}, repeat(1:testheight, 1, testwidth) + repeat(0.01 * (1:testwidth)', testheight, 1))
    a3D = convert(Array{Float32}, repeat(a2D, 1, 1, testdepth))
    for k = 1:testdepth; a3D[:,:,k] .+= 0.0001 * k; end

    texarr1D = CuTextureArray{Float32}(testheight)
    copyto!(texarr1D, a1D)
    tex1D = CuTexture(texarr1D)
    @test Array(fetch_all(tex1D)) == a1D

    texarr2D = CuTextureArray{Float32}(testheight, testwidth)
    copyto!(texarr2D, a2D)
    tex2D = CuTexture(texarr2D)
    @test Array(fetch_all(tex2D)) == a2D

    tex2D_dir = CuTexture(CuTextureArray(a2D))
    @test Array(fetch_all(tex2D_dir)) == a2D

    texarr3D = CuTextureArray{Float32}(testheight, testwidth, testdepth)
    copyto!(texarr3D, a3D)
    tex3D = CuTexture(texarr3D)
    @test Array(fetch_all(tex3D)) == a3D
end

@testset "Wrapping CuArray" begin
    testheight, testwidth, testdepth = 16, 16, 4
    a1D = convert(Array{Float32}, 1:testheight)
    a2D = convert(Array{Float32}, repeat(1:testheight, 1, testwidth) + repeat(0.01 * (1:testwidth)', testheight, 1))
    d_a1D = CuArray(a1D)
    d_a2D = CuArray(a2D)

    # Strangely, this is not working
    texwrap1D = CuTexture(d_a1D)
    @test_broken fetch_all(texwrap1D) == d_a1D

    # This works as long as d_a2D is well pitched
    texwrap2D = CuTexture(d_a2D)
    @test fetch_all(texwrap2D) == d_a2D
end

@testset "All CUDA types" begin

    for T in (Int32, UInt32, Int16, UInt16, Int8, UInt8, Float32, Float16)
        testheight, testwidth, testdepth = 32, 32, 4
        a2D = rand(T, testheight, testwidth)
        d_a2D = CuArray(a2D)

        # Using CuTextureArray
        tex_2D = CuTexture{T}(testheight, testwidth)
        copyto!(tex_2D.mem, d_a2D)
        @test fetch_all(tex_2D) == d_a2D

        # Wrapping CuArray
        # This works as long as d_a2D is well pitched
        texwrap_2D = CuTexture(d_a2D)
        @test fetch_all(texwrap_2D) == d_a2D
    end

end

@testset "Multiple channels" begin
    testheight, testwidth, testdepth = 16, 16, 4
    a2D = [(Int32(i), Int32(j)) for i = 1:testheight, j = 1:testwidth]
    d_a2D = CuArray(a2D)
    texarr2D = CuTextureArray{eltype(d_a2D)}(size(d_a2D)...)
    copyto!(texarr2D, d_a2D)
    tex2D = CuTexture(texarr2D)
    @test fetch_all(tex2D) == d_a2D

    testheight, testwidth, testdepth = 16, 16, 4
    a2D = [(Int16(i), Int16(j), Int16(i + j), Int16(i - j)) for i = 1:testheight, j = 1:testwidth]
    d_a2D = CuArray(a2D)
    texarr2D = CuTextureArray{eltype(d_a2D)}(size(d_a2D)...)
    copyto!(texarr2D, d_a2D)
    tex2D = CuTexture(texarr2D)
    @test fetch_all(tex2D) == d_a2D
end

@testset "Custom type" begin
    @testset "Auto cast" begin
        struct AKindOfRGBA
            r::UInt8
            g::UInt8
            b::UInt8
            a::UInt8
        end

        @test cuda_texture_alias_type(AKindOfRGBA) == NTuple{4,UInt8}

        testheight, testwidth, testdepth = 16, 16, 4
        a2D = [AKindOfRGBA(UInt8(i), UInt8(j), UInt8(j + i), UInt8(j > i)) for i = 1:testheight, j = 1:testwidth]
        d_a2D = CuArray(a2D)

        texarr2D = CuTextureArray{eltype(d_a2D)}(size(d_a2D)...)
        copyto!(texarr2D, d_a2D)
        tex2D = CuTexture(texarr2D)

        @test fetch_all(tex2D) == d_a2D
    end

    @testset "Manual cast" begin
        primitive type AKindOfFloat32 32 end

        # It is not enough to define `cuda_texture_alias_type(::Type{AKindOfFloat32}) = Float32`, one has to define the whole type.
        #   This is due to `cuda_texture_alias_type` being a *`@generated`* function.
        CuTextures.cuda_texture_alias_type(::Type{NTuple{2,AKindOfFloat32}}) = NTuple{2,Float32}
        @test cuda_texture_alias_type(NTuple{2,AKindOfFloat32}) == NTuple{2,Float32}

        testheight, testwidth, testdepth = 16, 16, 4
        a2D = [
            (reinterpret(AKindOfFloat32, Float32(i)), reinterpret(AKindOfFloat32, Float32(j)))
            for i = 1:testheight, j = 1:testwidth]
        d_a2D = CuArray(a2D)

        texarr2D = CuTextureArray{eltype(d_a2D)}(size(d_a2D)...)
        copyto!(texarr2D, d_a2D)
        tex2D = CuTexture(texarr2D)

        @test fetch_all(tex2D) == d_a2D
    end
end

@testset "Access by normalized coordinates" begin
    testheight, testwidth, testdepth = 16, 16, 4
    a2D = convert(Array{Float32}, repeat(1:testheight, 1, testwidth) + repeat(0.01 * (1:testwidth)', testheight, 1))
    d_a2D = CuArray(a2D)

    tex2D = CuTexture{Float32}(testheight, testwidth)
    copyto!(tex2D.mem, d_a2D)

    function kernel_texture_warp_native_normalized(dst::CuDeviceArray{T,2}, texture::CuDeviceTexture{T,2}) where {T}
        i, u = calcpoint(blockIdx().x, blockDim().x, threadIdx().x, size(dst)[1])
        j, v = calcpoint(blockIdx().y, blockDim().y, threadIdx().y, size(dst)[2])
        x = (u - 0.5f0) / Float32(size(dst)[1])
        y = (v - 0.5f0) / Float32(size(dst)[2])
        @inbounds dst[i,j] = texture(x,y); # Access by normalized coordinates
        return nothing
    end
    d_out = CuArray{eltype(tex2D)}(undef, size(tex2D))
    @cuda threads = size(d_out) kernel_texture_warp_native_normalized(d_out, tex2D)

    @test d_out == d_a2D
end


text:
#=
test_array:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-11
=#

module TestArray
    # packages
    using Test

    # external modules
    include("array_initialization.jl")
    include("array_operation.jl")

    # methods
    function test()
        @testset "Array" begin
            @testset "ArrayInit" begin
                @test_nowarn ArrayInit.main()
            end
            @testset "ArrayOp" begin
                @test_nowarn ArrayOp.main()
            end
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .TestArray
    TestArray.test()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
using DataDeps

register(DataDep(
  "chess-king-rook-vs-king",
  "http://archive.ics.uci.edu/ml/datasets/Chess+%28King-Rook+vs.+King%29",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/chess/king-rook-vs-king/krkopt.data",
  "497b61d20d3d9e2ce310f9f4d2ca92dbdb0691d9e449d983ec6f7baf387dc7a9",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=7,
        feature_indices=1:6,
        categoric_indices=1:6,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
using HorizonSideRobots
include("../basics.jl")

md"""
        # Найти прямоугольную перегородку, включая внешнюю рамку, с наибольшей средней температурой клеток периметра
"""

"""Решение задачи 25"""
function solve_25!(r)

end


text:
"""
An instance of the budgeted m-set problem. The budget is a minimum total
weight that must be taken for a solution to be feasible.

All weights are supposed to be positive integers, the budget too (this
hypothesis is crucial to get efficient algorithms). Values are not restricted.

It can be formalised as follows:

``\\max \\sum_i \\mathrm{values}_i x_i``
``\\mathrm{s.t.} \\sum_i x_i \\leq m, \\quad \\sum_i \\mathrm{weights}_i x_i \\geq \\mathrm{budget}, \\qquad x \\in \\{0, 1\\}^d``

`max_weight` is a parameter that may be used by the algorithm and corresponds
to the largest value a single weight may take. If not present, it is computed
as `max(weights)`.

`budget` can be ignored; in that case, it will be computed as the maximum
budget that may have a feasible solution, i.e. `d * max_weight`, `d` being
the dimension of the problem (the number of objects). If a `budget` is given
and the algorithm computes all solutions for lower budgets, then the algorithm
is allowed to stop at the given budget (i.e. it will not be possible to
retrieve a solution for a larger budget).
"""
struct BudgetedMSetInstance
  values::Vector{Float64}
  weights::Vector{Int}
  m::Int
  budget::Int # weights * solution >= budget
  max_weight::Int

  function BudgetedMSetInstance(values::Vector{Float64}, weights::Vector{Int}, m::Int;
                                budget::Union{Nothing, Int}=nothing,
                                max_weight::Union{Nothing, Int}=nothing)
    # Error checking.
    if m < 0
      error("m is less than zero: there is no solution.")
    end

    if m == 0
      error("m is zero: the only solution is to take no items, if the budget constraint is satisfied.")
    end

    if ! isnothing(budget) && budget < 0
      error("Budget is present and negative.")
    end

    if any(weights .< 0)
      error("At least a weight is negative.")
    end

    if length(values) != length(weights)
      error("Not the same number of values and weights; these two vectors must have the same size.")
    end

    # Complete the optional parameters.
    if isnothing(max_weight)
      max_weight = maximum(weights)
    end

    if isnothing(budget)
      d = length(values) # Dimension of the problem.
      budget = d * max_weight
    end

    # Return a new instance.
    new(values, weights, m, budget, max_weight)
  end
end

values(i::BudgetedMSetInstance) = i.values
weights(i::BudgetedMSetInstance) = i.weights
m(i::BudgetedMSetInstance) = i.m
budget(i::BudgetedMSetInstance) = i.budget
max_weight(i::BudgetedMSetInstance) = i.max_weight

value(i::BudgetedMSetInstance, o::Int) = values(i)[o]
values(i::BudgetedMSetInstance, o) = values(i)[o]
weight(i::BudgetedMSetInstance, o::Int) = weights(i)[o]
weights(i::BudgetedMSetInstance, o) = weights(i)[o]
dimension(i::BudgetedMSetInstance) = length(values(i))

struct BudgetedMSetSolution
  instance::BudgetedMSetInstance
  items::Vector{Int} # Indices to the chosen items.
  state::Array{Float64, 3} # Data structure built by the dynamic-programming recursion.
  solutions::Dict{Tuple{Int, Int, Int}, Vector{Int}} # From the indices of state to the corresponding solution.
end

function value(s::BudgetedMSetSolution)
  return sum(s.instance.values[i] for i in s.items)
end

function items(s::BudgetedMSetSolution, budget::Int)
  return s.solutions[s.instance.m, 0, budget]
end

function items_all_budgets(s::BudgetedMSetSolution, max_budget::Int)
  sol = Dict{Int, Vector{Int}}()
  m = s.instance.m
  for budget in 0:max_budget
    sol[budget] = s.solutions[m, 0, budget]
  end
  return sol
end

function value(s::BudgetedMSetSolution, budget::Int)
  its = items(s, budget)
  if -1 in its
    return -Inf
  end
  return sum(s.instance.values[i] for i in its)
end


text:
using Gaugefields

function heatbathtest_4D(NX,NY,NZ,NT,β,NC)
    Dim = 4
    Nwing = 1

    U = Initialize_Gaugefields(NC,Nwing,NX,NY,NZ,NT,condition = "cold")
    println(typeof(U))

    gauge_action = GaugeAction(U)
    plaqloop = make_loops_fromname("plaquette",Dim=Dim)
    append!(plaqloop,plaqloop')
    βinp = β/2
    push!(gauge_action,βinp,plaqloop)

    rectloop = make_loops_fromname("rectangular",Dim=Dim)
    append!(rectloop,rectloop')
    βinp = β/2
    push!(gauge_action,βinp,rectloop)
    hnew = Heatbath_update(U,gauge_action)

    show(gauge_action)

    temp1 = similar(U[1])
    temp2 = similar(U[1])
    temp3 = similar(U[1])

    comb = 6
    factor = 1/(comb*U[1].NV*U[1].NC)
    @time plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
    println("plaq_t = $plaq_t")
    poly = calculate_Polyakov_loop(U,temp1,temp2)
    println("polyakov loop = $(real(poly)) $(imag(poly))")

    numhb = 1000
    for itrj = 1:numhb

        heatbath!(U,hnew)

        plaq_t = calculate_Plaquette(U,temp1,temp2)*factor
        poly = calculate_Polyakov_loop(U,temp1,temp2)

        if itrj % 40 == 0
            println("$itrj plaq_t = $plaq_t")
            println("$itrj polyakov loop = $(real(poly)) $(imag(poly))")
        end
    end

    close(fp)
    return plaq_t

end

NX = 4
NY = 4
NZ = 4
NT = 4
NC = 3
β = 5.7
heatbathtest_4D(NX,NY,NZ,NT,β,NC)