text:
# # Introduction
#
# *A gentle introduction to parallelization and GPU programming in Julia*
#
# [Julia](https://julialang.org/) has first-class support for GPU programming: you can use
# high-level abstractions or obtain fine-grained control, all without ever leaving your
# favorite programming language. The purpose of this tutorial is to help Julia users take
# their first step into GPU computing. In this tutorial, you'll compare CPU and GPU
# implementations of a simple calculation, and learn about a few of the factors that
# influence the performance you obtain.
#
# This tutorial is inspired partly by a blog post by Mark Harris, [An Even Easier
# Introduction to CUDA](https://devblogs.nvidia.com/even-easier-introduction-cuda/), which
# introduced CUDA using the C++ programming language. You do not need to read that
# tutorial, as this one starts from the beginning.

# ## A simple example on the CPU

# We'll consider the following demo, a simple calculation on the CPU.

N = 2^20
x = fill(1.0f0, N)  # a vector filled with 1.0 (Float32)
y = fill(2.0f0, N)  # a vector filled with 2.0

y .+= x             # increment each element of y with the corresponding element of x

# check that we got the right answer
using Test
@test all(y .== 3.0f0)

# From the `Test Passed` line we know everything is in order. We used `Float32` numbers in
# preparation for the switch to GPU computations: GPUs are faster (sometimes, much faster)
# when working with `Float32` than with `Float64`.

# A distinguishing feature of this calculation is that every element of `y` is being
# updated using the same operation. This suggests that we might be able to parallelize
# this.

# ### Parallelization on the CPU

# First let's do the parallelization on the CPU. We'll create a "kernel function" (the
# computational core of the algorithm) in two implementations, first a sequential version:

function sequential_add!(y, x)
    for i in eachindex(y, x)
        @inbounds y[i] += x[i]
    end
    return nothing
end

fill!(y, 2)
sequential_add!(y, x)
@test all(y .== 3.0f0)

# And now a parallel implementation:

function parallel_add!(y, x)
    Threads.@threads for i in eachindex(y, x)
        @inbounds y[i] += x[i]
    end
    return nothing
end

fill!(y, 2)
parallel_add!(y, x)
@test all(y .== 3.0f0)

# Now if I've started Julia with `JULIA_NUM_THREADS=4` on a machine with at least 4 cores,
# I get the following:

@assert Threads.nthreads() == 4     #src

# ```julia
# using BenchmarkTools
# @btime sequential_add!($y, $x)
# ```

# ```
#   487.303 μs (0 allocations: 0 bytes)
# ```

# versus

# ```julia
# @btime parallel_add!($y, $x)
# ```

# ```
#   259.587 μs (13 allocations: 1.48 KiB)
# ```

# You can see there's a performance benefit to parallelization, though not by a factor of 4
# due to the overhead for starting threads. With larger arrays, the overhead would be
# "diluted" by a larger amount of "real work"; these would demonstrate scaling that is
# closer to linear in the number of cores. Conversely, with small arrays, the parallel
# version might be slower than the serial version.

# ## Your first GPU computation

# ### Installation

# For most of this tutorial you need to have a computer with a compatible GPU and have
# installed [CUDA](https://developer.nvidia.com/cuda-downloads). You should also install
# the following packages using Julia's [package
# manager](https://docs.julialang.org/en/latest/stdlib/Pkg/):

# ```julia
# pkg> add CUDA
# ```

# If this is your first time, it's not a bad idea to test whether your GPU is working by
# testing the CUDA.jl package:

# ```julia
# pkg> add CUDA
# pkg> test CUDA
# ```

# ### Parallelization on the GPU

# We'll first demonstrate GPU computations at a high level using the `CuArray` type,
# without explicitly writing a kernel function:

using CUDA

x_d = CUDA.fill(1.0f0, N)  # a vector stored on the GPU filled with 1.0 (Float32)
y_d = CUDA.fill(2.0f0, N)  # a vector stored on the GPU filled with 2.0

# Here the `d` means "device," in contrast with "host". Now let's do the increment:

y_d .+= x_d
@test all(Array(y_d) .== 3.0f0)

# The statement `Array(y_d)` moves the data in `y_d` back to the host for testing. If we
# want to benchmark this, let's put it in a function:

function add_broadcast!(y, x)
    CUDA.@sync y .+= x
    return
end

# ```julia
# @btime add_broadcast!($y_d, $x_d)
# ```

# ```
#   67.047 μs (84 allocations: 2.66 KiB)
# ```

# The most interesting part of this is the call to `CUDA.@sync`. The CPU can assign
# jobs to the GPU and then go do other stuff (such as assigning *more* jobs to the GPU)
# while the GPU completes its tasks. Wrapping the execution in a `CUDA.@sync` block
# will make the CPU block until the queued GPU tasks are done, similar to how `Base.@sync`
# waits for distributed CPU tasks. Without such synchronization, you'd be measuring the
# time takes to launch the computation, not the time to perform the computation. But most
# of the time you don't need to synchronize explicitly: many operations, like copying
# memory from the GPU to the CPU, implicitly synchronize execution.

# For this particular computer and GPU, you can see the GPU computation was significantly
# faster than the single-threaded CPU computation, and that the use of multiple CPU threads makes
# the CPU implementation competitive. Depending on your hardware you may get different
# results.

# ### Writing your first GPU kernel

# Using the high-level GPU array functionality made it easy to perform this computation
# on the GPU. However, we didn't learn about what's going on under the hood, and that's the
# main goal of this tutorial. So let's implement the same functionality with a GPU kernel:

function gpu_add1!(y, x)
    for i = 1:length(y)
        @inbounds y[i] += x[i]
    end
    return nothing
end

fill!(y_d, 2)
@cuda gpu_add1!(y_d, x_d)
@test all(Array(y_d) .== 3.0f0)

# Aside from using the `CuArray`s `x_d` and `y_d`, the only GPU-specific part of this is the
# *kernel launch* via `@cuda`. The first time you issue this `@cuda` statement, it will
# compile the kernel (`gpu_add1!`) for execution on the GPU. Once compiled, future
# invocations are fast. You can see what `@cuda` expands to using `?@cuda` from the Julia
# prompt.

# Let's benchmark this:

function bench_gpu1!(y, x)
    CUDA.@sync begin
        @cuda gpu_add1!(y, x)
    end
end

# ```julia
# @btime bench_gpu1!($y_d, $x_d)
# ```

# ```
#   119.783 ms (47 allocations: 1.23 KiB)
# ```

# That's a *lot* slower than the version above based on broadcasting. What happened?

# ### Profiling

# When you don't get the performance you expect, usually your first step should be to
# profile the code and see where it's spending its time. For that, you'll need to be able to
# run NVIDIA's [`nvprof`
# tool](https://devblogs.nvidia.com/cuda-pro-tip-nvprof-your-handy-universal-gpu-profiler/).
# On Unix systems, launch Julia this way:
#
# ```sh
# $ nvprof --profile-from-start off /path/to/julia
# ```
#
# replacing the `/path/to/julia` with the path to your Julia binary. Note that we don't
# immediately start the profiler, but instead call into the CUDA APIs and manually start the
# profiler with `CUDA.@profile` (thus excluding the time to compile our kernel):

bench_gpu1!(y_d, x_d)  # run it once to force compilation
CUDA.@profile bench_gpu1!(y_d, x_d)

# When we quit the Julia REPL, the profiler process will print information about the
# executed kernels and API calls:

# ```
# ==2574== Profiling result:
#             Type  Time(%)      Time     Calls       Avg       Min       Max  Name
#  GPU activities:  100.00%  247.61ms         1  247.61ms  247.61ms  247.61ms  ptxcall_gpu_add1__1
#       API calls:   99.54%  247.83ms         1  247.83ms  247.83ms  247.83ms  cuEventSynchronize
#                     0.46%  1.1343ms         1  1.1343ms  1.1343ms  1.1343ms  cuLaunchKernel
#                     0.00%  4.9490us         1  4.9490us  4.9490us  4.9490us  cuEventRecord
#                     0.00%  4.4190us         1  4.4190us  4.4190us  4.4190us  cuEventCreate
#                     0.00%     960ns         2     480ns     358ns     602ns  cuCtxGetCurrent
# ```

# You can see that 100% of the time was spent in `ptxcall_gpu_add1__1`, the name of the
# kernel that CUDA.jl assigned when compiling `gpu_add1!` for these inputs. (Had you created
# arrays of multiple data types, e.g., `xu_d = CUDA.fill(0x01, N)`, you might have also seen
# `ptxcall_gpu_add1__2` and so on. Like the rest of Julia, you can define a single method
# and it will be specialized at compile time for the particular data types you're using.)

# For further insight, run the profiling with the option `--print-gpu-trace`. You can also
# invoke Julia with as argument the path to a file containing all commands you want to run
# (including a call to `CUDA.@profile`):
#
# ```sh
# $ nvprof --profile-from-start off --print-gpu-trace /path/to/julia /path/to/script.jl
#      Start  Duration   Grid Size   Block Size     Regs*    SSMem*    DSMem*           Device   Context    Stream  Name
#   13.3134s  245.04ms     (1 1 1)      (1 1 1)        20        0B        0B  GeForce GTX TIT         1         7  ptxcall_gpu_add1__1 [34]
# ```

# The key thing to note here is the `(1 1 1)` in the "Grid Size" and "Block Size" columns.
# These terms will be explained shortly, but for now, suffice it to say that this is an
# indication that this computation ran sequentially. Of note, sequential processing with
# GPUs is much slower than with CPUs; where GPUs shine is with large-scale parallelism.

# ### Writing a parallel GPU kernel

# To speed up the kernel, we want to parallelize it, which means assigning different tasks
# to different threads.  To facilitate the assignment of work, each CUDA thread gets access
# to variables that indicate its own unique identity, much as
# [`Threads.threadid()`](https://docs.julialang.org/en/latest/manual/parallel-computing/#Multi-Threading-(Experimental)-1)
# does for CPU threads. The CUDA analogs of `threadid` and `nthreads` are called
# `threadIdx` and `blockDim`, respectively; one difference is that these return a
# 3-dimensional structure with fields `x`, `y`, and `z` to simplify cartesian indexing for
# up to 3-dimensional arrays. Consequently we can assign unique work in the following way:

function gpu_add2!(y, x)
    index = threadIdx().x    # this example only requires linear indexing, so just use `x`
    stride = blockDim().x
    for i = index:stride:length(y)
        @inbounds y[i] += x[i]
    end
    return nothing
end

fill!(y_d, 2)
@cuda threads=256 gpu_add2!(y_d, x_d)
@test all(Array(y_d) .== 3.0f0)

# Note the `threads=256` here, which divides the work among 256 threads numbered in a
# linear pattern. (For a two-dimensional array, we might have used `threads=(16, 16)` and
# then both `x` and `y` would be relevant.)

# Now let's try benchmarking it:

function bench_gpu2!(y, x)
    CUDA.@sync begin
        @cuda threads=256 gpu_add2!(y, x)
    end
end

# ```julia
# @btime bench_gpu2!($y_d, $x_d)
# ```

# ```
#   1.873 ms (47 allocations: 1.23 KiB)
# ```

# Much better!

# But obviously we still have a ways to go to match the initial broadcasting result. To do
# even better, we need to parallelize more. GPUs have a limited number of threads they can
# run on a single *streaming multiprocessor* (SM), but they also have multiple SMs. To take
# advantage of them all, we need to run a kernel with multiple *blocks*. We'll divide up
# the work like this:
#
# ![block grid](intro1.png)
#
# This diagram was [borrowed from a description of the C/C++
# library](https://devblogs.nvidia.com/even-easier-introduction-cuda/); in Julia, threads
# and blocks begin numbering with 1 instead of 0. In this diagram, the 4096 blocks of 256
# threads (making 1048576 = 2^20 threads) ensures that each thread increments just a single
# entry; however, to ensure that arrays of arbitrary size can be handled, let's still use a
# loop:

function gpu_add3!(y, x)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = blockDim().x * gridDim().x
    for i = index:stride:length(y)
        @inbounds y[i] += x[i]
    end
    return
end

numblocks = ceil(Int, N/256)

fill!(y_d, 2)
@cuda threads=256 blocks=numblocks gpu_add3!(y_d, x_d)
@test all(Array(y_d) .== 3.0f0)

# The benchmark:

function bench_gpu3!(y, x)
    numblocks = ceil(Int, length(y)/256)
    CUDA.@sync begin
        @cuda threads=256 blocks=numblocks gpu_add3!(y, x)
    end
end

# ```julia
# @btime bench_gpu3!($y_d, $x_d)
# ```

# ```
#   67.268 μs (52 allocations: 1.31 KiB)
# ```

# Finally, we've achieved the similar performance to what we got with the broadcasted
# version. Let's run `nvprof` again to confirm this launch configuration:
#
# ```
# ==23972== Profiling result:
#    Start  Duration            Grid Size      Block Size     Regs*    SSMem*    DSMem*           Device   Context    Stream  Name
# 13.3526s  101.22us           (4096 1 1)       (256 1 1)        32        0B        0B  GeForce GTX TIT         1         7  ptxcall_gpu_add3__1 [34]
# ```

# ### Printing

# When debugging, it's not uncommon to want to print some values. This is achieved with
# `@cuprint`:

function gpu_add2_print!(y, x)
    index = threadIdx().x    # this example only requires linear indexing, so just use `x`
    stride = blockDim().x
    @cuprintln("thread $index, block $stride")
    for i = index:stride:length(y)
        @inbounds y[i] += x[i]
    end
    return nothing
end

@cuda threads=16 gpu_add2_print!(y_d, x_d)
synchronize()

# Note that the printed output is only generated when synchronizing the entire GPU with
# `synchronize()`. This is similar to `CUDA.@sync`, and is the counterpart of
# `cudaDeviceSynchronize` in CUDA C++.

# ### Error-handling

# The final topic of this intro concerns the handling of errors. Note that the kernels
# above used `@inbounds`, but did not check whether `y` and `x` have the same length. If
# your kernel does not respect these bounds, you will run into nasty errors:

# ```
# ERROR: CUDA error: an illegal memory access was encountered (code #700, ERROR_ILLEGAL_ADDRESS)
# Stacktrace:
#  [1] ...
# ```

# If you remove the `@inbounds` annotation, instead you get
#
# ```
# ERROR: a exception was thrown during kernel execution.
#        Run Julia on debug level 2 for device stack traces.
# ```

# As the error message mentions, a higher level of debug information will result in a more
# detailed report. Let's run the same code with with `-g2`:
#
# ```
# ERROR: a exception was thrown during kernel execution.
# Stacktrace:
#  [1] throw_boundserror at abstractarray.jl:484
#  [2] checkbounds at abstractarray.jl:449
#  [3] setindex! at /home/tbesard/Julia/CUDA/src/device/array.jl:79
#  [4] some_kernel at /tmp/tmpIMYANH:6
# ```

# !!! warning
#
#     On older GPUs (with a compute capability below `sm_70`) these errors are fatal,
#     and effectively kill the CUDA environment. On such GPUs, it's often a good idea to
#     perform your "sanity checks" using code that runs on the CPU and only turn over the
#     computation to the GPU once you've deemed it to be safe.

# ## Summary

# Keep in mind that the high-level functionality of CUDA often means that you don't
# need to worry about writing kernels at such a low level. However, there are many cases
# where computations can be optimized using clever low-level manipulations. Hopefully, you
# now feel comfortable taking the plunge.


text:
using PGFPlots
using Colors
using ColorBrewer

function plot_estimated_prob(tree)
    nbin = 200
    xmin = tree.lbs[1]
    xmax = tree.ubs[1]
    ymin = tree.lbs[2]
    ymax = tree.ubs[2]

    #safe_color = RGB(0.0, 0.0, 1.0) # blue
    safe_color = RGB(169.0 / 255.0, 169.0 / 255.0, 169.0 / 255.0) # gray
    #safe_color = RGB(105.0 / 255.0, 105.0 / 255.0, 105.0 / 255.0) # dark gray
    #unsafe_color = RGB(1.0, 0.0, 0.0) # red
    unsafe_color = RGB(139.0 / 255.0, 0.0, 0.0) # dark red

    colors = [safe_color, unsafe_color]

    ax = Axis(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, width="7cm", height="8cm",
    ylabel=L"$\theta$ (degrees)", xlabel=L"$x$ (meters)", title="Safe Cells")

    function get_heat(x, y)
        traj = sim_traj_true_next([x, y], tree)
        xvals = [traj[i][1] for i = 1:length(traj)]
        prob = any(xvals .< -10) || any(xvals .> 10)
        return prob
    end

    push!(ax, Plots.Image(get_heat, (xmin, xmax), (ymin, ymax), zmin = 0, zmax = 1,
        xbins = nbin, ybins = nbin, colormap = ColorMaps.RGBArrayMap(colors), colorbar=false))

    return ax
end

function plot_prob(tree)
    nbin = 200
    xmin = tree.lbs[1]
    xmax = tree.ubs[1]
    ymin = tree.lbs[2]
    ymax = tree.ubs[2]

    #safe_color = RGB(0.0, 0.0, 1.0) # blue
    safe_color = RGB(169.0 / 255.0, 169.0 / 255.0, 169.0 / 255.0) # gray
    #safe_color = RGB(105.0 / 255.0, 105.0 / 255.0, 105.0 / 255.0) # dark gray
    #unsafe_color = RGB(1.0, 0.0, 0.0) # red
    unsafe_color = RGB(139.0 / 255.0, 0.0, 0.0) # dark red

    colors = [safe_color, unsafe_color]

    ax = Axis(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, width="7cm", height="8cm",
    ylabel=L"$\theta$ (degrees)", xlabel=L"$x$ (meters)", title="Safe Cells")

    function get_heat(x, y)
        leaf = get_leaf(tree.root_node, [x, y])
        return leaf.prob
    end

    push!(ax, Plots.Image(get_heat, (xmin, xmax), (ymin, ymax), zmin = 0, zmax = 1,
        xbins = nbin, ybins = nbin, colormap = ColorMaps.RGBArrayMap(colors), colorbar=false))

    return ax
end

function plot_control_range(tree; nbin = 200)
    xmin = tree.lbs[1]
    xmax = tree.ubs[1]
    ymin = tree.lbs[2]
    ymax = tree.ubs[2]

    ax = Axis(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, width="7cm", height="8cm",
    ylabel=L"$\theta$ (degrees)", xlabel=L"$x$ (meters)", title="Control Range")

    function get_heat(x, y)
        leaf = get_leaf(tree.root_node, [x, y])
        return leaf.max_control - leaf.min_control
    end

    push!(ax, Plots.Image(get_heat, (xmin, xmax), (ymin, ymax),
        xbins = nbin, ybins = nbin, colormap = pasteljet, colorbar=false))

    return ax
end

function plot_reachable(tree)
    nbin = 200
    xmin = tree.lbs[1]
    xmax = tree.ubs[1]
    ymin = tree.lbs[2]
    ymax = tree.ubs[2]

    unreach_color = RGB(1.0, 1.0, 1.0) # white
    #reach_color = RGB(0.0, 0.0, 0.502) # navy blue
    reach_color = RGB(0.0, 0.5, 0.5) # teal

    colors = [unreach_color, reach_color]

    ax = Axis(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, width="7cm", height="8cm",
    ylabel=L"$\theta$ (degrees)", xlabel=L"$x$ (meters)", title="Reachable Cells")

    function get_heat(x, y)
        leaf = get_leaf(tree.root_node, [x, y])
        return leaf.prob
    end

    push!(ax, Plots.Image(get_heat, (xmin, xmax), (ymin, ymax), zmin = 0, zmax = 1,
        xbins = nbin, ybins = nbin, colormap = ColorMaps.RGBArrayMap(colors), colorbar=false))

    return ax
end

function plot_cells(tree; color = "blue", opacity = 1.0, lw = 1)
    xmin = tree.lbs[1]
    xmax = tree.ubs[1]
    ymin = tree.lbs[2]
    ymax = tree.ubs[2]

    #safe_color = RGB(0.0, 0.0, 1.0) # blue
    safe_color = RGB(169.0 / 255.0, 169.0 / 255.0, 169.0 / 255.0) # gray
    #safe_color = RGB(105.0 / 255.0, 105.0 / 255.0, 105.0 / 255.0) # dark gray
    #unsafe_color = RGB(1.0, 0.0, 0.0) # red
    unsafe_color = RGB(139.0 / 255.0, 0.0, 0.0) # dark red

    colors = [safe_color, unsafe_color]

    ax = Axis(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, width="7cm", height="8cm",
    ylabel=L"$\theta$ (degrees)", xlabel=L"$x$ (meters)")

    lbs, ubs, probs = get_bounds_and_probs(tree)

    for i = 1:length(lbs)
        if probs[i] == 1.0
            #push!(ax, Plots.Command(get_filled_rectangle_border(lbs[i], ubs[i], color,
                            #opacity = opacity, lw = lw)))
            push!(ax, Plots.Command(get_rectangle(lbs[i], ubs[i], color, lw)))
        else
            push!(ax, Plots.Command(get_rectangle(lbs[i], ubs[i], "black", 0.1)))
        end
    end

    return ax
end

function plot_reachable_xy_cells(xy_tree; color = "blue", opacity = 1.0)
    nbin = 200
    xmin = xy_tree.lbs[1]
    xmax = xy_tree.ubs[1]
    ymin = xy_tree.lbs[2]
    ymax = xy_tree.ubs[2]

    ax = Axis(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, width="7cm", height="8cm",
    ylabel=L"$y$ (meters)", xlabel=L"$x$ (meters)", title="Reachable Cells")

    lbs, ubs, probs = get_bounds_and_probs(xy_tree)

    for i = 1:length(lbs)
        if probs[i] == 1.0
            push!(ax, Plots.Command(get_filled_rectangle(lbs[i], ubs[i], color, opacity = opacity)))
        end
    end

    return ax
end

function get_road(;xmin=0.0, xmax=150.0, ymin=-12.0, ymax=12.0)
	ax = Axis(PGFPlots.Plots.Linear([xmin, xmax], [-10.0, -10.0], mark="none",
            style="black,solid,thick, name path=A"))
	push!(ax, PGFPlots.Plots.Linear([xmin, xmax], [10.0, 10.0], mark="none",
            style="black,solid,thick, name path=B"))
	push!(ax, PGFPlots.Plots.Linear([xmin, xmax], [0.0, 0.0], mark="none",
            style="white,dashed,thick"))
	push!(ax, PGFPlots.Plots.Command("\\path[name path=axis] (axis cs:$xmin,$ymin) -- (axis cs:$xmax,$ymin);"))
	push!(ax, PGFPlots.Plots.Command("\\path[name path=axisright] (axis cs:$xmin,$ymax) -- (axis cs:$xmax,$ymax);"))
	push!(ax, PGFPlots.Plots.Command("\\addplot[pastelGreen!40] fill between[of=A and axis];"))
	push!(ax, PGFPlots.Plots.Command("\\addplot[pastelGreen!40] fill between[of=B and axisright];"))
	push!(ax, PGFPlots.Plots.Command("\\addplot[black!40] fill between[of=A and B];"))
	ax.width = "15cm"
	ax.height = "4cm"
    ax.xmin = xmin
    ax.xmax = xmax
    ax.ymin = ymin
    ax.ymax = ymax
    ax.xlabel = "Downtrack Position (m)"
    ax.ylabel = "Crosstrack Error (m)"
	#ax.axisEqualImage = true
	return ax
end

function plot_reachable_cells_on_road(xy_tree::KDTREE; color = "blue", opacity = 1.0,
                                xmin=0.0, xmax=150.0, ymin=-12.0, ymax=12.0)
    ax = get_road(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax)

    lbs, ubs, probs = get_bounds_and_probs(xy_tree)

    for i = 1:length(lbs)
        if probs[i] == 1.0
            push!(ax, Plots.Command(get_filled_rectangle(reverse(lbs[i]), reverse(ubs[i]),
                                                            color, opacity = opacity)))
        end
    end

    return ax
end

function plot_reachable_cells_on_road(lbs, ubs; color = "blue", opacity = 1.0,
                                xmin=0.0, xmax=150.0, ymin=-12.0, ymax=12.0)
    ax = get_road(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax)

    for i = 1:length(lbs)
        push!(ax, Plots.Command(get_filled_rectangle(reverse(lbs[i]), reverse(ubs[i]),
                                                        color, opacity = opacity)))
    end

    return ax
end

function plot_reachable_cells_and_trajs_on_road(xy_tree, xres, yres; color = "blue", opacity = 1.0,
                                xmin=0.0, xmax=150.0, ymin=-12.0, ymax=12.0)
    ax = get_road(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax)
    for i = 1:length(xres)
        push!(ax, Plots.Linear(yres[i][1:5:end], xres[i][1:5:end], mark = "none", style = "solid, blue"))
    end

    lbs, ubs, probs = get_bounds_and_probs(xy_tree)

    for i = 1:length(lbs)
        if probs[i] == 1.0
            push!(ax, Plots.Command(get_filled_rectangle(reverse(lbs[i]), reverse(ubs[i]),
                                                            color, opacity = opacity)))
        end
    end

    return ax
end

""" Plotting Utils """

function get_filled_rectangle(lb, ub, color; opacity = 1.0)
    return "\\fill [$(color), opacity=$(opacity)] (axis cs:$(string(lb[1])),$(string(lb[2]))) rectangle (axis cs:$(string(ub[1])),$(string(ub[2])));"
end

function get_filled_rectangle_border(lb, ub, color; opacity = 1.0, lw = 1)
    return "\\draw [fill=$(color), $(color), opacity=$(opacity), line width = $(lw)] (axis cs:$(string(lb[1])),$(string(lb[2]))) rectangle (axis cs:$(string(ub[1])),$(string(ub[2])));"
end

function get_rectangle(lb, ub)
    return "\\draw ($(string(lb[1])),$(string(lb[2]))) rectangle ($(string(ub[1])),$(string(ub[2])));"
end

function get_rectangle(lb, ub, color, lw)
    return "\\draw [$(color), line width = $(lw)] ($(string(lb[1])),$(string(lb[2]))) rectangle ($(string(ub[1])),$(string(ub[2])));"
end

ext: jl

max_stars_repo_name: sisl/VerifyGAN

max_issues_repo_name: sisl/VerifyGAN

max_forks_repo_name: sisl/VerifyGAN


text:
"""Squared deviations from the array mean."""
function SqDeviation(X::Array{T,1}) where {T<:AbstractFloat}

    s = zero(T)
    μ = mean(X)

    @simd for xi in X
        s += (xi - μ)^2
    end

    return s
end

"""Linear deviations from the array mean - i.e. L1-Norm."""
function LinDeviation(X::Array{T,1}) where {T<:AbstractFloat}

    s = zero(T)
    μ = mean(X)

    @simd for xi in X
        s += abs(xi - μ)
    end

    return s
end

"""Squared deviations from the class means."""
function SDCM(X::Array{T,1},breaks::Array{Int,1}) where {T<:AbstractFloat}

    n = length(breaks)-1  # number of classes (aka intervals)

    @boundscheck breaks[end] == length(X)+1 || throw(BoundsError())
    @boundscheck breaks[1] == 1 || throw(BoundsError())

    s = zero(T)

    for iclass in 1:n
        s += SqDeviation(ClassValues(X,breaks,iclass))
    end

    return s
end

"""Goodness of variance fit for data array X (assumed to be sorted) and classes determined by break indices."""
function GVF(X::Array{T,1},breaks::Array{Int,1}) where {T<:AbstractFloat}
    sdam = SqDeviation(X)       # squared deviation from array mean
    sdcm = SDCM(X,breaks)       # sum  of squared deviations from class mean for each class
    return 1 - sdcm/sdam
end

"""Average rounding error for data array X (assumed to be sorted) and classes determind by break indices."""
function ARE(X::Array{T,1},breaks::Array{Int,1}) where {T<:AbstractFloat}
    s = zero(T)
    ndata = length(X)
    for iclass in 1:(length(breaks)-1)
        s += LinDeviation(ClassValues(X,breaks,iclass))
    end
    return s/ndata
end


text:
using PastaQ
using ITensors
using Test
using LinearAlgebra

@testset "Gate generation: 1-qubit gates" begin
  i = Index(2, tags = "Qubit")

  g = gate("I", i)
  @test plev(inds(g)[1]) == 1
  @test plev(inds(g)[2]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ Matrix{Int}(I, 2, 2)

  g = gate("X", i)
  @test plev(inds(g)[1]) == 1
  @test plev(inds(g)[2]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ Matrix{Int}(I, 2, 2)
  @test g ≈ gate("σx",i)
  @test g ≈ gate("σ1",i)

  g = gate("Y", i)
  @test plev(inds(g)[1]) == 1
  @test plev(inds(g)[2]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ Matrix{Int}(I, 2, 2)
  @test g ≈ gate("σy",i)
  @test g ≈ gate("σ2",i)

  g = gate("Z", i)
  @test plev(inds(g)[1]) == 1
  @test plev(inds(g)[2]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ Matrix{Int}(I, 2, 2)
  @test g ≈ gate("σz",i)
  @test g ≈ gate("σ3",i)

  g = gate("√X", i)
  @test plev(inds(g)[1]) == 1
  @test plev(inds(g)[2]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ Matrix{Int}(I, 2, 2)
  @test g ≈ gate("√NOT",i)

  g = gate("H", i)
  @test plev(inds(g)[1]) == 1
  @test plev(inds(g)[2]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ Matrix{Int}(I, 2, 2)

  g = gate("S", i)
  @test plev(inds(g)[1]) == 1
  @test plev(inds(g)[2]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ Matrix{Int}(I, 2, 2)
  @test g ≈ gate("P",i)

  g = gate("T", i)
  @test plev(inds(g)[1]) == 1
  @test plev(inds(g)[2]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ Matrix{Int}(I, 2, 2)

  g = state("X+", i)
  @test plev(ind(g, 1)) == 0

  g = state("X-", i)
  @test plev(ind(g, 1)) == 0

  g = state("Y+", i)
  @test plev(ind(g, 1)) == 0

  g = state("Y-", i)
  @test plev(ind(g, 1)) == 0

  g = state("Z+", i)
  @test plev(ind(g, 1)) == 0

  g = state("Z-", i)
  @test plev(ind(g, 1)) == 0

  θ = π * rand()
  g = gate("Rx", i; θ=θ)
  @test plev(inds(g)[1]) == 1
  @test plev(inds(g)[2]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ Matrix{Int}(I, 2, 2)

  θ = π * rand()
  g = gate("Ry", i; θ=θ)
  @test plev(inds(g)[1]) == 1
  @test plev(inds(g)[2]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ Matrix{Int}(I, 2, 2)

  ϕ = 2π * rand()
  g = gate("Rz", i; ϕ=ϕ)
  @test plev(inds(g)[1]) == 1
  @test plev(inds(g)[2]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ Matrix{Int}(I, 2, 2)

  angles = rand(3)
  θ = π * angles[1]
  ϕ = 2π * angles[2]
  λ = 2π * angles[3]
  g = gate("Rn", i; θ=θ, ϕ=ϕ, λ=λ)
  @test plev(inds(g)[1]) == 1
  @test plev(inds(g)[2]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ Matrix{Int}(I, 2, 2)
  @test g ≈ gate("Rn̂", i; θ=θ, ϕ=ϕ, λ=λ)
end

@testset "Gate generation: 2-qubit gates" begin
  i = Index(2)
  j = Index(2)
  angles = rand(3)
  θ = π * angles[1]
  ϕ = 2π * angles[2]
  λ = 2π * angles[3]

  g = gate("SWAP", i, j)
  @test plev(inds(g)[1]) == 1 && plev(inds(g)[2]) == 1
  @test plev(inds(g)[3]) == 0 && plev(inds(g)[4]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ reshape(Matrix{ComplexF64}(I, 4, 4), (2, 2, 2, 2))
  @test g ≈ gate("Swap", i, j)
  @test g ≈ gate("Sw", i, j)

  g = gate("√SWAP", i, j)
  @test plev(inds(g)[1]) == 1 && plev(inds(g)[2]) == 1
  @test plev(inds(g)[3]) == 0 && plev(inds(g)[4]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ reshape(Matrix{ComplexF64}(I, 4, 4), (2, 2, 2, 2))
  @test g ≈ gate("√Swap", i, j)
  @test g ≈ gate("√Sw", i, j)

  g = gate("iSWAP", i, j)
  @test plev(inds(g)[1]) == 1 && plev(inds(g)[2]) == 1
  @test plev(inds(g)[3]) == 0 && plev(inds(g)[4]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test g ≈ gate("iSwap", i, j)
  @test g ≈ gate("iSw", i, j)

  g = gate("CX", i, j)
  @test plev(inds(g)[1]) == 1 && plev(inds(g)[2]) == 1
  @test plev(inds(g)[3]) == 0 && plev(inds(g)[4]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ reshape(Matrix{ComplexF64}(I, 4, 4), (2, 2, 2, 2))

  g = gate("CY", i, j)
  @test plev(inds(g)[1]) == 1 && plev(inds(g)[2]) == 1
  @test plev(inds(g)[3]) == 0 && plev(inds(g)[4]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ reshape(Matrix{ComplexF64}(I, 4, 4), (2, 2, 2, 2))

  g = gate("CZ", i, j)
  @test plev(inds(g)[1]) == 1 && plev(inds(g)[2]) == 1
  @test plev(inds(g)[3]) == 0 && plev(inds(g)[4]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ reshape(Matrix{ComplexF64}(I, 4, 4), (2, 2, 2, 2))

  g = gate("CRz", i, j; ϕ = ϕ)
  @test plev(inds(g)[1]) == 1 && plev(inds(g)[2]) == 1
  @test plev(inds(g)[3]) == 0 && plev(inds(g)[4]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ reshape(Matrix{ComplexF64}(I, 4, 4), (2, 2, 2, 2))

  g = gate("CRn", i, j; θ=θ, ϕ=ϕ, λ=λ)
  @test plev(inds(g)[1]) == 1 && plev(inds(g)[2]) == 1
  @test plev(inds(g)[3]) == 0 && plev(inds(g)[4]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ reshape(Matrix{ComplexF64}(I, 4, 4), (2, 2, 2, 2))
  @test g ≈ gate("CRn̂", i, j; θ=θ, ϕ=ϕ, λ=λ)
end

@testset "Gate generation: 3-qubit gates" begin
  i = Index(2)
  j = Index(2)
  k = Index(2)

  g = gate("Toffoli", i, j, k)
  @test plev(inds(g)[1]) == 1 && plev(inds(g)[2]) == 1 && plev(inds(g)[3]) == 1
  @test plev(inds(g)[4]) == 0 && plev(inds(g)[5]) == 0 && plev(inds(g)[6]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ reshape(Matrix{ComplexF64}(I, 8, 8), (2, 2, 2, 2, 2, 2))
  @test g ≈ gate("CCNOT", i, j, k)
  @test g ≈ gate("CCX", i, j, k)
  @test g ≈ gate("TOFF", i, j, k)

  g = gate("Fredkin", i, j, k)
  @test plev(inds(g)[1]) == 1 && plev(inds(g)[2]) == 1 && plev(inds(g)[3]) == 1
  @test plev(inds(g)[4]) == 0 && plev(inds(g)[5]) == 0 && plev(inds(g)[6]) == 0
  ggdag = g * prime(dag(g), 1; plev=1)
  @test ITensors.array(ggdag) ≈ reshape(Matrix{ComplexF64}(I, 8, 8), (2, 2, 2, 2, 2, 2))
  @test g ≈ gate("CSWAP", i, j, k)
  @test g ≈ gate("CSwap", i, j, k)
  @test g ≈ gate("CSw", i, j, k)

end

@testset "Gate generation: Haar gates" begin
  i = Index(2)
  j = Index(2)

  g = gate("randU", i)
  @test hasinds(g, i', i)
  Id = g * swapprime(dag(g)', 2 => 1)
  for n in 1:dim(i), n′ in 1:dim(i)
    if n == n′
      @test Id[n, n′] ≈ 1
    else
      @test Id[n, n′] ≈ 0 atol = 1e-15
    end
  end

  g = gate("randU", i, j)
  @test hasinds(g, i', j', i, j)
  Id = g * swapprime(dag(g)', 2 => 1)
  for n in 1:dim(i), m in 1:dim(j), n′ in 1:dim(i), m′ in 1:dim(j)
    if (n, m) == (n′, m′)
      @test Id[n, m, n′, m′] ≈ 1
    else
      @test Id[n, m, n′, m′] ≈ 0 atol = 1e-14
    end
  end

@testset "Custom gate with long name" begin
    PastaQ.gate(::GateName"my_favorite_gate") = [0.11 0.12; 0.21 0.22]
    s = Index(2, "Qubit, Site")
    gate("my_favorite_gate", s)
    g = gate("my_favorite_gate", s)
    @test g[s' => 1, s => 1] == 0.11
    @test g[s' => 1, s => 2] == 0.12
    @test g[s' => 2, s => 1] == 0.21
    @test g[s' => 2, s => 2] == 0.22
  end
end

@testset "Bases rotations" begin

  zp = [1.0,0.0]
  zm = [0.0,1.0]
  xp = [1.0,1.0]/√2
  xm = [1.0,-1.0]/√2
  yp = [1.0,im]/√2
  ym = [1.0,-im]/√2

  Ux = gate("basisX"; dag = true)
  @test Ux * xp ≈ zp
  @test Ux * xm ≈ zm

  Ux = gate("basisY"; dag = true)
  @test Ux * yp ≈ zp
  @test Ux * ym ≈ zm

end

@testset "apply gate: Id" begin
  psi = productstate(2)
  gate_data = ("I", 1)
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi) ≈ [1.0, 0.0, 0.0, 0.0]
end

@testset "apply gate: X" begin
  # Build gate first, then apply using an ITensor
  psi = productstate(2)
  site = 1
  gate_data = ("X", 1)
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi) ≈ [0.0, 0.0, 1.0, 0.0]
end

@testset "apply gate: Y" begin
  psi = productstate(2)
  site = 1
  gate_data = ("Y", 1)
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi) ≈ [0.0, 0.0, im, 0.0]
end

@testset "apply gate: Z" begin
  # Build gate first, then apply using an ITensor
  psi = productstate(2)
  site = 1
  gate_data = ("Z", 1)
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi) ≈ [1.0, 0.0, 0.0, 0.0]
end

@testset "apply gate: H" begin
  # Build gate first, then apply using an ITensor
  psi = productstate(2)
  site = 1
  gate_data = ("H", 1)
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi) ≈ 1 / sqrt(2.0) * [1.0, 0.0, 1.0, 0.0]
end

@testset "apply gate: S" begin
  psi = productstate(1)
  gate_data = ("S", 1)
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi[1]) ≈ [1.0, 0.0]
  psi = productstate(1)
  gate_data = ("X", 1)
  psi = runcircuit(psi, gate_data)
  gate_data = ("S", 1)
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi[1]) ≈ [0.0, im]
end

@testset "apply gate: T" begin
  psi = productstate(1)
  gate_data = ("T", 1)
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi[1]) ≈ [1.0, 0.0]
  psi = productstate(1)
  gate_data = ("X", 1)
  psi = runcircuit(psi, gate_data)
  gate_data = ("T", 1)
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi[1]) ≈ [0.0, exp(im * π / 4)]
end

@testset "apply gate: Rx" begin
  θ = π * rand()
  psi = productstate(1)
  gate_data = ("Rx", 1, (θ=θ,))
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi[1]) ≈ [cos(θ / 2.0), -im * sin(θ / 2.0)]
  psi = productstate(1)
  psi = runcircuit(psi, ("X", 1))
  #PastaQ.applygate!(psi,"X",1)
  gate_data = ("Rx", 1, (θ=θ,))
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi[1]) ≈ [-im * sin(θ / 2.0), cos(θ / 2.0)]
end

@testset "apply gate: Ry" begin
  θ = π * rand()
  psi = productstate(1)
  gate_data = ("Ry", 1, (θ=θ,))
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi[1]) ≈ [cos(θ / 2.0), sin(θ / 2.0)]
  psi = productstate(1)
  psi = runcircuit(psi, ("X", 1))
  gate_data = ("Ry", 1, (θ=θ,))
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi[1]) ≈ [-sin(θ / 2.0), cos(θ / 2.0)]
end

@testset "apply gate: Rz" begin
  ϕ = 2π * rand()
  psi = productstate(1)
  gate_data = ("Rz", 1, (ϕ=ϕ,))
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi[1]) ≈ [1.0, 0.0]
  psi = productstate(1)
  psi = runcircuit(psi, ("X", 1))
  gate_data = ("Rz", 1, (ϕ=ϕ,))
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi[1]) ≈ [0.0, exp(im * ϕ)]
end

@testset "apply gate: Rn" begin
  θ = 1.0
  ϕ = 2.0
  λ = 3.0
  psi = productstate(1)
  gate_data = ("Rn", 1, (θ=θ, ϕ=ϕ, λ=λ))
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi[1]) ≈ [cos(θ / 2.0), exp(im * ϕ) * sin(θ / 2.0)]
  psi = productstate(1)
  psi = runcircuit(psi, ("X", 1))
  gate_data = ("Rn", 1, (θ=θ, ϕ=ϕ, λ=λ))
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi[1]) ≈
        [-exp(im * λ) * sin(θ / 2.0), exp(im * (ϕ + λ)) * cos(θ / 2.0)]
end

@testset "apply gate: prep X+/X-" begin
  psi = productstate(1, ["X+"])
  @test PastaQ.array(psi) ≈ 1 / sqrt(2.0) * [1.0, 1.0]

  psi = productstate(1, ["X-"])
  @test PastaQ.array(psi) ≈ 1 / sqrt(2.0) * [1.0, -1.0]
end

@testset "apply gate: prep Y+/Y-" begin
  psi = productstate(1, ["Y+"])
  @test PastaQ.array(psi) ≈ 1 / sqrt(2.0) * [1.0, im]

  psi = productstate(1, ["Y-"])
  @test PastaQ.array(psi) ≈ 1 / sqrt(2.0) * [1.0, -im]
end

@testset "apply gate: prep Z+/Z-" begin
  psi = productstate(1, ["Z+"])
  @test PastaQ.array(psi) ≈ [1.0, 0.0]

  psi = productstate(1, ["Z-"])
  @test PastaQ.array(psi) ≈ [0.0, 1.0]
end

@testset "apply gate: meas X" begin
  psi = productstate(1, ["X+"])
  gate_data = ("basisX", 1, (dag=true,))
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi) ≈ [1.0, 0.0]
  psi = productstate(1, ["X-"])
  gate_data = ("basisX", 1, (dag=true,))
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi) ≈ [0.0, 1.0]
end

@testset "apply gate: meas Y" begin
  psi = productstate(1, ["Y+"])
  gate_data = ("basisY", 1, (dag=true,))
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi) ≈ [1.0, 0.0]
  psi = productstate(1, ["Y-"])
  gate_data = ("basisY", 1, (dag=true,))
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi) ≈ [0.0, 1.0]
end

@testset "apply gate: meas Z" begin
  psi = productstate(1)
  gate_data = ("basisZ", 1, (dag=true,))
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi) ≈ [1.0, 0.0]
  psi = productstate(1, ["Z-"])
  gate_data = ("basisZ", 1, (dag=true,))
  psi = runcircuit(psi, gate_data)
  @test PastaQ.array(psi) ≈ [0.0, 1.0]
end

@testset "apply gate: CX" begin
  # CONTROL - TARGET
  psi = productstate(2)
  # |00> -> |00> = (1 0 0 0) (natural order)
  gate_data = ("CX", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [1.0, 0.0, 0.0, 0.0]

  psi = productstate(2)
  # |10> -> |11> = (0 0 0 1) (natural order)
  gate_data = ("X", 1)
  psi = runcircuit(psi, gate_data)
  gate_data = ("CX", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 0.0, 0.0, 1.0]

  psi = productstate(2)
  # |01> -> |01> = (0 1 0 0) (natural order)
  gate_data = ("X", 2)
  psi = runcircuit(psi, gate_data)
  gate_data = ("CX", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 1.0, 0.0, 0.0]

  psi = productstate(2)
  # |11> -> |10> = (0 0 1 0) (natural order)
  gate_data = ("X", 1)
  psi = runcircuit(psi, gate_data)
  gate_data = ("X", 2)
  psi = runcircuit(psi, gate_data)
  gate_data = ("CX", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 0.0, 1.0, 0.0]
end

@testset "apply gate: CY" begin
  # CONTROL - TARGET
  psi = productstate(2)
  # |00> -> |00> = (1 0 0 0) (natural order)
  gate_data = ("CY", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [1.0, 0.0, 0.0, 0.0]

  psi = productstate(2)
  # |10> -> i|11> = (0 0 0 i) (natural order)
  gate_data = ("X", 1)
  psi = runcircuit(psi, gate_data)
  gate_data = ("CY", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 0.0, 0.0, im]

  psi = productstate(2)
  # |01> -> |01> = (0 1 0 0) (natural order)
  gate_data = ("X", 2)
  psi = runcircuit(psi, gate_data)
  gate_data = ("CY", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 1.0, 0.0, 0.0]

  psi = productstate(2)
  # |11> -> -i|10> = (0 0 -i 0) (natural order)
  gate_data = ("X", 1)
  psi = runcircuit(psi, gate_data)
  gate_data = ("X", 2)
  psi = runcircuit(psi, gate_data)
  gate_data = ("CY", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 0.0, -im, 0.0]

  # TARGET - CONTROL
  psi = productstate(2)
  # |00> -> |00> = (1 0 0 0) (natural order)
  gate_data = ("CY", (2, 1))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [1.0, 0.0, 0.0, 0.0]

  psi = productstate(2)
  # |10> -> |10> = (0 0 1 0) (natural order)
  gate_data = ("X", 1)
  psi = runcircuit(psi, gate_data)
  gate_data = ("CY", (2, 1))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 0.0, 1.0, 0.0]

  psi = productstate(2)
  # |01> -> i|11> = (0 0 0 i) (natural order)
  gate_data = ("X", 2)
  psi = runcircuit(psi, gate_data)
  gate_data = ("CY", (2, 1))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 0.0, 0.0, im]

  psi = productstate(2)
  # |11> -> -i|01> = (0 -i 0 0) (natural order)
  gate_data = ("X", 1)
  psi = runcircuit(psi, gate_data)
  gate_data = ("X", 2)
  psi = runcircuit(psi, gate_data)
  gate_data = ("CY", (2, 1))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, -im, 0.0, 0.0]
end

@testset "apply gate: CZ" begin
  # CONTROL - TARGET
  psi = productstate(2)
  # |00> -> |00> = (1 0 0 0) (natural order)
  gate_data = ("CZ", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [1.0, 0.0, 0.0, 0.0]

  psi = productstate(2)
  # |10> -> |10> = (0 0 1 0) (natural order)
  gate_data = ("X", 1)
  psi = runcircuit(psi, gate_data)
  gate_data = ("CZ", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 0.0, 1.0, 0.0]

  psi = productstate(2)
  # |01> -> |01> = (0 1 0 0) (natural order)
  gate_data = ("X", 2)
  psi = runcircuit(psi, gate_data)
  gate_data = ("CZ", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 1.0, 0.0, 0.0]

  psi = productstate(2)
  # |11> -> -|11> = (0 0 0 -1) (natural order)
  gate_data = ("X", 1)
  psi = runcircuit(psi, gate_data)
  gate_data = ("X", 2)
  psi = runcircuit(psi, gate_data)
  gate_data = ("CZ", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 0.0, 0.0, -1.0]

  # CONTROL - TARGET
  psi = productstate(2)
  # |00> -> |00> = (1 0 0 0) (natural order)
  gate_data = ("CZ", (2, 1))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [1.0, 0.0, 0.0, 0.0]

  psi = productstate(2)
  # |10> -> |10> = (0 0 1 0) (natural order)
  gate_data = ("X", 1)
  psi = runcircuit(psi, gate_data)
  gate_data = ("CZ", (2, 1))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 0.0, 1.0, 0.0]

  psi = productstate(2)
  # |01> -> |01> = (0 1 0 0) (natural order)
  gate_data = ("X", 2)
  psi = runcircuit(psi, gate_data)
  gate_data = ("CZ", (2, 1))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 1.0, 0.0, 0.0]

  psi = productstate(2)
  # |11> -> -|11> = (0 0 0 -1) (natural order)
  gate_data = ("X", 1)
  psi = runcircuit(psi, gate_data)
  gate_data = ("X", 2)
  psi = runcircuit(psi, gate_data)
  gate_data = ("CZ", (2, 1))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 0.0, 0.0, -1.0]
end

@testset "apply gate: Sw" begin
  # CONTROL - TARGET
  psi = productstate(2)
  # |00> -> |00> = (1 0 0 0) (natural order)
  gate_data = ("Sw", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [1.0, 0.0, 0.0, 0.0]

  psi = productstate(2)
  # |10> -> |01> = (0 1 0 0) (natural order)
  gate_data = ("X", 1)
  psi = runcircuit(psi, gate_data)
  gate_data = ("Sw", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 1.0, 0.0, 0.0]

  psi = productstate(2)
  # |01> -> |10> = (0 0 1 0) (natural order)
  gate_data = ("X", 2)
  psi = runcircuit(psi, gate_data)
  gate_data = ("Sw", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 0.0, 1.0, 0.0]

  psi = productstate(2)
  # |11> -> |11> = (0 0 0 1) (natural order)
  gate_data = ("X", 1)
  psi = runcircuit(psi, gate_data)
  gate_data = ("X", 2)
  psi = runcircuit(psi, gate_data)
  gate_data = ("Sw", (1, 2))
  psi = runcircuit(psi, gate_data)
  psi_vec = PastaQ.array(psi)
  @test psi_vec ≈ [0.0, 0.0, 0.0, 1.0]
end

@testset "function applied to a gate" begin
  s = siteinds("Qubit", 2)

  θ = 0.1
  rx = gate("Rx"; θ = 0.1)
  exp_rx = exp(rx)
  gtest = gate("Rx",s[1]; θ = 0.1, f = x -> exp(x))
  @test PastaQ.array(gtest) ≈ exp_rx

  cx = 0.1*gate("CX")
  exp_cx = exp(cx)
  gtest = gate("CX",s[1],s[2]; f = x -> exp(0.1*x))
  @test PastaQ.array(gtest) ≈ exp_cx

end

@testset "gate built out of elementary gates" begin

  s = siteinds("Qubit", 2)

  G = gate("S") * gate("S")
  gtest = gate("S * S", s[1])
  @test PastaQ.array(gtest) ≈ G

  cx = gate("CX")

  G = cx * cx * cx
  gtest = gate("CX*CX*CX", s[1],s[2])
  @test PastaQ.array(gtest) ≈ G

  G = gate("S") + gate("T")
  gtest = gate("S + T", s[1])
  @test PastaQ.array(gtest) ≈ G

  G = gate("S") + gate("T") * gate("X")
  gtest = gate("S + T * X", s[1])
  @test PastaQ.array(gtest) ≈ G

  G = gate("S") + gate("T") * gate("X") * gate("Y")
  gtest = gate("S + T * X * Y", s[1])
  @test PastaQ.array(gtest) ≈ G

  G = gate("Z")*gate("S") + gate("T") * gate("X") * gate("Y")
  gtest = gate("Z*S + T * X * Y", s[1])
  @test PastaQ.array(gtest) ≈ G

  exp_cx = exp(0.1*cx)
  gtest = gate("CX",s[1],s[2]; f = x -> exp(0.1*x))
  @test PastaQ.array(gtest) ≈ exp_cx

  d = 3
  q = siteinds("Qudit", 4; dim = d)
  g1 = gate("a†a", q[1]; f = x -> exp(0.1*x))
  g2 = gate("a† * a", q[1]; f = x -> exp(0.1*x))
  @test g1 ≈ g2
end

@testset "qudit gates" begin
  dim = 3
  s = siteinds("Qudit", 4; dim = dim)
  dims = (dim,)
  @test gate("a†", dims)      ≈ [0 0 0; 1 0 0; 0 √2 0]
  @test gate("a", dims)       ≈ [0 1 0; 0 0 √2; 0 0 0]

  @test gate("a†", dims) * gate("a", dims) ≈ [0 0 0; 0 1 0; 0 0 2]

  @test PastaQ.array(gate("a†", s[1])) ≈ [0 0 0; 1 0 0; 0 √2 0]
  @test PastaQ.array(gate("a", s[1])) ≈ [0 1 0; 0 0 √2; 0 0 0]

  @test PastaQ.array(gate("a†a", s[1])) ≈ [0 0 0; 0 1 0; 0 0 2]
  @test PastaQ.array(gate("aa†", s[1])) ≈ [1 0 0; 0 2 0; 0 0 0]
  @test PastaQ.array(gate("a† * a", s[1])) ≈ [0 0 0; 0 1 0; 0 0 2]
  @test PastaQ.array(gate("a * a†", s[1])) ≈ [1 0 0; 0 2 0; 0 0 0]

  dim = 10
  s = siteinds("Qudit", 4; dim = dim)
  dims = (dim,)
  @test PastaQ.array(gate("a† * a† * a * a", s[1])) ≈ gate("a†", dims)  * gate("a†", dims)  * gate("a", dims) * gate("a", dims)

  @test PastaQ.array(gate("a†a", s[1], s[2])) ≈ kron(gate("a†", dims), gate("a", dims))
  @test PastaQ.array(gate("aa†", s[1], s[2])) ≈ kron(gate("a", dims), gate("a†", dims))

  @test PastaQ.array(gate("a†a+aa†", s[1], s[2])) ≈ kron(gate("a†", dims), gate("a", dims)) + kron(gate("a", dims), gate("a†", dims))
end


text:
module Models

export IncompressibleModel, NonDimensionalModel

include("IncompressibleModels/IncompressibleModels.jl")
include("ShallowWaterModels/ShallowWaterModels.jl")

using .IncompressibleModels: IncompressibleModel, NonDimensionalModel
using .ShallowWaterModels: ShallowWaterModel

end


text: @show gamma(1)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test
using QuantumLattices.Prerequisites.VectorSpaces
using QuantumLattices.Interfaces: dimension, rank, ⊕, ⊗

import QuantumLattices.Prerequisites.VectorSpaces: VectorSpaceStyle, shape
import QuantumLattices.Prerequisites.Traits: contentnames, getcontent

struct SimpleVectorSpace{B, N} <: VectorSpace{B}
    sorted::Bool
    table::NTuple{N, B}
    SimpleVectorSpace(sorted::Bool, table::NTuple{N, B}) where {B, N} = new{B, N}(sorted, table)
end
@inline VectorSpaceStyle(::Type{<:SimpleVectorSpace}) = VectorSpaceEnumerative()
@inline contentnames(::Type{<:SimpleVectorSpace}) = (:sorted, :table)
@inline Base.issorted(vs::SimpleVectorSpace) = vs.sorted

@testset "VectorSpaceEnumerative" begin
    id0, id4 = (1, 0), (1, 4)
    id1, id2, id3 = (1, 1), (1, 2), (1, 3)
    vs = SimpleVectorSpace(true, (id1, id2, id3))
    @test vs == deepcopy(vs)
    @test isequal(vs, deepcopy(vs))
    @test vs|>size == (3,)
    @test vs|>dimension == 3
    @test vs|>collect == [id1, id2, id3]
    @test vs[1]==id1 && vs[2]==id2 && vs[3]==id3
    @test searchsortedfirst(vs, id0)==1 && searchsortedfirst(vs, id4)==4
    @test searchsortedfirst(vs, id1)==1 && searchsortedfirst(vs, id2)==2 && searchsortedfirst(vs, id3)==3
    @test isnothing(findfirst(id0, vs)) && isnothing(findfirst(id4, vs))
    @test findfirst(id1, vs)==1 && findfirst(id2, vs)==2 && findfirst(id3, vs)==3
    @test (id0∉vs) && (id4 ∉ vs)
    @test (id1∈vs) && (id2∈vs) && (id3 ∈ vs)

    vs = SimpleVectorSpace(false, (id1, id2, id3))
    @test searchsortedfirst(vs, id0)==4 && searchsortedfirst(vs, id4)==4
    @test searchsortedfirst(vs, id1)==1 && searchsortedfirst(vs, id2)==2 && searchsortedfirst(vs, id3)==3
    @test isnothing(findfirst(id0, vs)) && isnothing(findfirst(id4, vs))
    @test findfirst(id1, vs)==1 && findfirst(id2, vs)==2 && findfirst(id3, vs)==3
end

struct SimpleIndices{N} <: VectorSpace{CartesianIndex{N}}
    shape::NTuple{N, UnitRange{Int}}
    SimpleIndices(shape::NTuple{N, UnitRange{Int}}) where N = new{N}(shape)
end
@inline SimpleIndices(shape::UnitRange{Int}...) = SimpleIndices(shape)
@inline VectorSpaceStyle(::Type{<:SimpleIndices}) = VectorSpaceCartesian()
@inline shape(vs::SimpleIndices) = vs.shape
@inline Base.CartesianIndex(basis::CartesianIndex{N}, ::SimpleIndices{N}) where N = basis

@testset "VectorSpaceCartesian" begin
    foi = SimpleIndices(2:3, 2:3, 2:3)
    @test dimension(foi) == 8
    @test issorted(foi) == true
    @test foi|>collect == CartesianIndex.([(2, 2, 2), (3, 2, 2), (2, 3, 2), (3, 3, 2), (2, 2, 3), (3, 2, 3), (2, 3, 3), (3, 3, 3)])
    for (i, index) in enumerate(CartesianIndices((2:3, 2:3, 2:3)))
        @test foi[i] == foi[index] == index
        @test findfirst(index, foi) == i
        @test searchsortedfirst(foi, index) == i
        @test index∈foi
    end
    i1, i2 = CartesianIndex(1, 1, 1), CartesianIndex(4, 4, 4)
    @test isnothing(findfirst(i1, foi)) && isnothing(findfirst(i2, foi))
    @test searchsortedfirst(foi, i1) == 1 && searchsortedfirst(foi, i2) == 9
    @test i1∉foi && i2∉foi
end

struct DirectSummedVectorSpace{T<:Tuple, B} <: VectorSpace{B}
    contents::T
    DirectSummedVectorSpace(contents::Tuple) = new{typeof(contents), mapreduce(eltype, typejoin, contents)}(contents)
end
@inline VectorSpaceStyle(::Type{<:DirectSummedVectorSpace}) = VectorSpaceDirectSummed()
@inline DirectSummedVectorSpace(contents...) = DirectSummedVectorSpace(contents)

@testset "VectorSpaceDirectSummed" begin
    a = SimpleIndices(1:3)
    b = SimpleIndices(3:4, 7:7)
    c = DirectSummedVectorSpace(a, b)
    @test dimension(c) == 5
    @test c[1]==CartesianIndex(1) && c[2]==CartesianIndex(2) && c[3]==CartesianIndex(3) && c[4]==CartesianIndex(3, 7) && c[5]==CartesianIndex(4, 7)
end

@testset "NamedVectorSpace" begin
    t = ParameterSpace{:t}(1:2)
    @test contentnames(typeof(t)) == (:content,)
    @test dimension(t) == 2
    @test t≠ParameterSpace{:μ}(1:2) && t==ParameterSpace{:t}(1:2)
    @test !isequal(t, ParameterSpace{:μ}(1:2)) && isequal(t, ParameterSpace{:t}(1:2))
    @test keys(t) == keys(typeof(t)) == (:t,)
    @test t[1]==1 && t[2]==2
    @test t[1]∈t && t[2]∈t
    ps = pairs(t)
    @test size(ps) == (2,)
    @test eltype(ps) == eltype(typeof(ps)) == NamedTuple{(:t,), Tuple{Int}}
    @test ps[1]==(t=1,) && ps[2]==(t=2,)

    U = ParameterSpace{:U}([8.0, 9.0])
    zps = ZippedNamedVectorSpace(t, U)
    @test VectorSpaceStyle(zps) == VectorSpaceZipped()
    @test eltype(zps) == eltype(typeof(zps)) == Tuple{Int, Float64}
    @test keys(zps) == keys(typeof(zps)) == (:t, :U)
    @test dimension(zps) == 2
    @test zps[1]==(1, 8.0) && zps[2]==(2, 9.0)
    ps = pairs(zps)
    @test size(ps) == (2,)
    @test eltype(ps) == eltype(typeof(ps)) == NamedTuple{(:t, :U), Tuple{Int, Float64}}
    @test ps[1]==(t=1, U=8.0) && ps[2]==(t=2, U=9.0)
    pps = DirectProductedNamedVectorSpace(t, U)
    @test VectorSpaceStyle(pps) == VectorSpaceDirectProducted()
    @test eltype(pps) == eltype(typeof(pps)) == Tuple{Int, Float64}
    @test keys(pps) == keys(typeof(pps)) == (:t, :U)
    @test dimension(pps) == 4
    @test pps[1]==(1, 8.0) && pps[2]==(2, 8.0) && pps[3]==(1, 9.0) && pps[4]==(2, 9.0)

    μ = ParameterSpace{:μ}([11, 12])
    @test t⊕U == zps
    @test μ⊕zps == ZippedNamedVectorSpace(μ, t, U)
    @test zps⊕μ == ZippedNamedVectorSpace(t, U, μ)
    @test zps⊕ZippedNamedVectorSpace(μ) == ZippedNamedVectorSpace(t, U, μ)
    @test t⊗U == pps
    @test μ⊗pps == DirectProductedNamedVectorSpace(μ, t, U)
    @test pps⊗μ == DirectProductedNamedVectorSpace(t, U, μ)
    @test pps⊗DirectProductedNamedVectorSpace(μ) == DirectProductedNamedVectorSpace(t, U, μ)
end


text:
using UnivariateFunctions
using Dates

tol = 10*eps()

# Testing constructors.
und = Undefined_Function()
before = PE_Function(0.0,0.0,0.0,1)
first_ = PE_Function(1.0,0.0,0.0,1)
second= PE_Function(20.0,0.0,0.0,2)
last_  = Sum_Of_Functions([first_, second])

f0 = Piecewise_Function([-1.0,3.0, 10.0], [before, first_, second])
!ismissing(f0(-0.75))
f0 = trim_piecewise_function(f0, -0.5, 9.0)
ismissing(f0(-0.75))
f0 = Piecewise_Function([-Inf,3.0, 10.0], [before, first_, second])
!ismissing(f0(-0.75))

f1 = Piecewise_Function([-Inf, -1.0,3.0, 10.0], [before, first_, second, last_])
abs(evaluate(f1, -10.0) - 0.0 ) < tol
abs(evaluate(f1,   1.0) - 1.0 ) < tol
abs(evaluate(f1,   5.0) - 500.0) < tol
abs(evaluate(f1,  20.0) - 20.0 - 20*20*20) < tol
f2 = Piecewise_Function([-0.1, 0.0,2.0, 40.0], [before, first_, second, last_])

f3 =  Piecewise_Function([-0.1, 0.0,2.0, 40.0], [before, f1, second, last_])
f3.starts_ == [-Inf, -0.100, 0.00, 2.00, 40.0]
isa(f3.functions_[1],UnivariateFunctions.Undefined_Function)
isa(f3.functions_[2],UnivariateFunctions.PE_Function)
isa(f3.functions_[3],UnivariateFunctions.PE_Function)
isa(f3.functions_[4],UnivariateFunctions.PE_Function)
isa(f3.functions_[5],UnivariateFunctions.Sum_Of_Functions)

f4 =  Piecewise_Function([-2.1, -1.1,4.0, 40.0], [second, f1, second, last_])
f4.starts_ == [-Inf, -2.10, -1.10, -1.00, 3.00, 4.00, 40.0]
isa(f4.functions_[1], UnivariateFunctions.Undefined_Function)
isa(f4.functions_[2], UnivariateFunctions.PE_Function)
isa(f4.functions_[3], UnivariateFunctions.PE_Function)
isa(f4.functions_[4], UnivariateFunctions.PE_Function)
isa(f4.functions_[5], UnivariateFunctions.PE_Function)
isa(f4.functions_[6], UnivariateFunctions.PE_Function)
isa(f4.functions_[7], UnivariateFunctions.Sum_Of_Functions)

# Testing the other Piecewise_Function constructors
f5 = Piecewise_Function(Date(2020,05,01), second)
ismissing(f5(Date(2019,05,01)))
f6 = Piecewise_Function([Date(2020,05,01)], [second])
ismissing(f6(Date(2019,05,01)))
f7 = Piecewise_Function([77.0, Inf], [second, second])
ismissing(f7(Date(2019,05,01)))

function test_result(func, eval0, eval5, len = 1)
    val_test0 = abs(evaluate(func, 0.0) - eval0) < 1e-09
    if (!val_test0)
        print("Failed Val Test at 0.0")
    end
    val_test5 = abs(evaluate(func, 5.0) - eval5) < 1e-09
    if (!val_test5)
        print("Failed Val Test at 5.0")
    end
    if isa(func, UnivariateFunctions.Piecewise_Function)
        len_test = length(func.functions_) == len
        if (!len_test)
            print("Failed Length Test")
        end
        return all([val_test0,val_test5,len_test])
    else
        return all([val_test0,val_test5])
    end
end

# Testing with Ints
test_result(f1 + 5, 5.0, 505.0, 4)
test_result(5 + f1, 5.0, 505.0, 4)
test_result(f1 - 5, -5.0, 495.0, 4)
test_result(5 - f1, 5.0, -495.0, 4)
test_result(f1 * 5, 0.0, 2500.0, 4)
test_result(5 * f1, 0.0, 2500.0, 4)
test_result(f1 / 5, 0.0, 100.0, 4)
# And with Float64s
test_result(f1 + 5.0, 5.0, 505.0, 4)
test_result(5.0 + f1, 5.0, 505.0, 4)
test_result(f1 - 5.0, -5.0, 495.0, 4)
test_result(5.0 - f1, 5.0, -495.0, 4)
test_result(f1 * 5.0, 0.0, 2500.0, 4)
test_result(5.0 * f1, 0.0, 2500.0, 4)
test_result(f1 / 5.0, 0.0, 100.0, 4)

# Testing with other functions
f5 = f4 + f1
f5.starts_[2] == -2.1
typeof(f5.functions_[1]) == UnivariateFunctions.Undefined_Function

typeof(f4 + und) == UnivariateFunctions.Undefined_Function
typeof(und + f4) == UnivariateFunctions.Undefined_Function
typeof(f4 - und) == UnivariateFunctions.Undefined_Function
typeof(und - f4) == UnivariateFunctions.Undefined_Function
typeof(f4 * und) == UnivariateFunctions.Undefined_Function
typeof(und * f4) == UnivariateFunctions.Undefined_Function

typeof(-1*f4) == UnivariateFunctions.Piecewise_Function

test_result(f1 + first_, 0.0, 505.0, 4)
test_result(first_ + f1, 0.0, 505.0, 4)
test_result(f1 - first_, 0.0, 495.0, 4)
test_result(first_ - f1, 0.0, -495.0, 4)
test_result(f1 * first_, 0.0, 2500.0, 4)
test_result(first_ * f1, 0.0, 2500.0, 4)

test_result(f1 + last_, 0.0, 1005.0, 4)
test_result(last_ + f1, 0.0, 1005.0, 4)
test_result(f1 - last_, 0.0, -5.0, 4)
test_result(last_ - f1, 0.0, 5.0, 4)
test_result(f1 * last_, 0.0, 252500.0, 4)
test_result(last_ * f1, 0.0, 252500.0, 4)

test_result(f1 + f4, 0.0, 1000.0, 8)
test_result(f1 - f4, 0.0, 0.0, 8)
test_result(f1 * f4, 0.0, 250000.0, 8)

# Testing linear rescaling
#test function will be f4
alpha = 0.95
beta = 3.0
X = [-0.3, 0.0, 1.3, 1.5, 4.0, 7.0, 20.0, 40.0, 70.0, 100.0]
X_converted = (alpha .* X) .+ beta
f4_converted = convert_to_linearly_rescale_inputs(f4, alpha, beta)
y = evaluate.(f4,X)
converted_y = evaluate.(f4_converted, X_converted)
sum(abs.(y .- converted_y)) < 1e-10

# Testing with Undefined
und = Undefined_Function()
typeof(sum0 + und) == Undefined_Function
typeof(sum0 - und) == Undefined_Function
typeof(sum0 / und) == Undefined_Function
typeof(sum0 * und) == Undefined_Function
typeof(und + sum0) == Undefined_Function
typeof(und - sum0) == Undefined_Function
typeof(und / sum0) == Undefined_Function
typeof(und * sum0) == Undefined_Function
typeof(und / f5) == Undefined_Function
typeof(f5 / und) == Undefined_Function


text:
mutable struct UnivariateKalmanState{Fl<:AbstractFloat}
    v::Fl
    F::Fl
    att::Vector{Fl}
    a::Vector{Fl}
    Ptt::Matrix{Fl}
    P::Matrix{Fl}
    llk::Fl
    steady_state::Bool
    # Auxiliary matrices
    P_to_check_steady_state::Matrix{Fl}
    ZP::Vector{Fl}
    TPtt::Matrix{Fl}
    function UnivariateKalmanState(a1::Vector{Fl}, P1::Matrix{Fl}) where Fl
        m = length(a1)
        P_to_check_steady_state = zeros(Fl, m, m)
        ZP = zeros(Fl, m)
        TPtt = zeros(Fl, m, m)

        return new{Fl}(
            zero(Fl),
            zero(Fl),
            zeros(Fl, m),
            a1,
            zeros(Fl, m, m),
            P1,
            zero(Fl),
            false,
            P_to_check_steady_state,
            ZP,
            TPtt,
        )
    end
end

function save_a1_P1_in_filter_output!(
    filter_output::FilterOutput{Fl}, kalman_state::UnivariateKalmanState{Fl}
) where Fl
    filter_output.a[1] = deepcopy(kalman_state.a)
    filter_output.P[1] = deepcopy(kalman_state.P)
    return filter_output
end

function save_kalman_state_in_filter_output!(
    filter_output::FilterOutput{Fl}, kalman_state::UnivariateKalmanState{Fl}, t::Int
) where Fl
    filter_output.v[t] = copy(fill(kalman_state.v, 1))
    filter_output.F[t] = copy(fill(kalman_state.F, 1, 1))
    filter_output.a[t + 1] = copy(kalman_state.a)
    filter_output.att[t] = copy(kalman_state.att)
    filter_output.P[t + 1] = copy(kalman_state.P)
    filter_output.Ptt[t] = copy(kalman_state.Ptt)
    # There is no diffuse part in UnivariateKalmanState
    filter_output.Pinf[t] = copy(fill(zero(Fl), size(kalman_state.Ptt)))
    return filter_output
end

# Univariate Kalman filter with the recursions as described
# in Koopman's book TODO
"""
    UnivariateKalmanFilter{Fl <: AbstractFloat}

A Kalman filter that is tailored to univariate systems, exploiting the fact that the
dimension of the observations at any time period is 1.

# TODO equations and descriptions of a1 and P1
"""
mutable struct UnivariateKalmanFilter{Fl<:AbstractFloat} <: KalmanFilter
    steadystate_tol::Fl
    a1::Vector{Fl}
    P1::Matrix{Fl}
    skip_llk_instants::Int
    kalman_state::UnivariateKalmanState

    function UnivariateKalmanFilter(
        a1::Vector{Fl},
        P1::Matrix{Fl},
        skip_llk_instants::Int=length(a1),
        steadystate_tol::Fl=Fl(1e-5),
    ) where Fl
        kalman_state = UnivariateKalmanState(copy(a1), copy(P1))
        return new{Fl}(steadystate_tol, a1, P1, skip_llk_instants, kalman_state)
    end
end

function reset_filter!(kf::UnivariateKalmanFilter{Fl}) where Fl
    copyto!(kf.kalman_state.a, kf.a1)
    copyto!(kf.kalman_state.P, kf.P1)
    set_state_llk_to_zero!(kf.kalman_state)
    fill!(kf.kalman_state.P_to_check_steady_state, zero(Fl))
    fill!(kf.kalman_state.ZP, zero(Fl))
    fill!(kf.kalman_state.TPtt, zero(Fl))
    kf.kalman_state.steady_state = false
    return kf
end

function set_state_llk_to_zero!(kalman_state::UnivariateKalmanState{Fl}) where Fl
    kalman_state.llk = zero(Fl)
    return kalman_state
end

# TODO this should either have ! or return something
function check_steady_state(kalman_state::UnivariateKalmanState{Fl}, tol::Fl) where Fl
    @inbounds for j in axes(kalman_state.P, 2), i in axes(kalman_state.P, 1)
        if abs(
            (kalman_state.P[i, j] - kalman_state.P_to_check_steady_state[i, j]) /
            kalman_state.P[i, j],
        ) > tol
            # Update the P_to_check_steady_state matrix
            copyto!(kalman_state.P_to_check_steady_state, kalman_state.P)
            return nothing
        end
    end
    kalman_state.steady_state = true
    return nothing
end

function update_v!(
    kalman_state::UnivariateKalmanState{Fl}, y::Fl, Z::Vector{Fl}, d::Fl
) where Fl
    kalman_state.v = y - dot(Z, kalman_state.a) - d
    return kalman_state
end

function update_F!(kalman_state::UnivariateKalmanState{Fl}, Z::Vector{Fl}, H::Fl) where Fl
    LinearAlgebra.BLAS.gemv!('N', one(Fl), kalman_state.P, Z, zero(Fl), kalman_state.ZP)
    kalman_state.F = H
    kalman_state.F += dot(kalman_state.ZP, Z)
    return kalman_state
end

function update_att!(kalman_state::UnivariateKalmanState{Fl}, Z::Vector{Fl}) where Fl
    copyto!(kalman_state.att, kalman_state.a)
    @inbounds for i in axes(kalman_state.P, 1), j in axes(kalman_state.P, 2)
        kalman_state.att[i] +=
            (kalman_state.v / kalman_state.F) * kalman_state.P[i, j] * Z[j]
    end
    return kalman_state
end

function repeat_a_in_att!(kalman_state::UnivariateKalmanState{Fl}) where Fl
    for i in eachindex(kalman_state.att)
        kalman_state.att[i] = kalman_state.a[i]
    end
    return kalman_state
end

function update_a!(
    kalman_state::UnivariateKalmanState{Fl}, T::Matrix{Fl}, c::Vector{Fl}
) where Fl
    LinearAlgebra.BLAS.gemv!('N', one(Fl), T, kalman_state.att, zero(Fl), kalman_state.a)
    kalman_state.a .+= c
    return kalman_state
end

function update_Ptt!(kalman_state::UnivariateKalmanState{Fl}) where Fl
    LinearAlgebra.BLAS.gemm!(
        'N', 'T', one(Fl), kalman_state.ZP, kalman_state.ZP, zero(Fl), kalman_state.Ptt
    )
    @. kalman_state.Ptt = kalman_state.P - (kalman_state.Ptt / kalman_state.F)
    return kalman_state
end

function repeat_P_in_Ptt!(kalman_state::UnivariateKalmanState{Fl}) where Fl
    for i in axes(kalman_state.P, 1), j in axes(kalman_state.P, 2)
        kalman_state.Ptt[i, j] = kalman_state.P[i, j]
    end
    return kalman_state
end

function update_P!(
    kalman_state::UnivariateKalmanState{Fl}, T::Matrix{Fl}, RQR::Matrix{Fl}
) where Fl
    LinearAlgebra.BLAS.gemm!(
        'N', 'N', one(Fl), T, kalman_state.Ptt, zero(Fl), kalman_state.TPtt
    )
    LinearAlgebra.BLAS.gemm!(
        'N', 'T', one(Fl), kalman_state.TPtt, T, zero(Fl), kalman_state.P
    )
    kalman_state.P .+= RQR
    return kalman_state
end

function update_P!(
    kalman_state::UnivariateKalmanState{Fl}, T::Matrix{Fl}, R::Matrix{Fl}, Q::Matrix{Fl}
) where Fl
    LinearAlgebra.BLAS.gemm!(
        'N', 'N', one(Fl), T, kalman_state.Ptt, zero(Fl), kalman_state.TPtt
    )
    LinearAlgebra.BLAS.gemm!(
        'N', 'T', one(Fl), kalman_state.TPtt, T, zero(Fl), kalman_state.P
    )
    kalman_state.P .+= R * Q * R'
    return kalman_state
end

function update_llk!(kalman_state::UnivariateKalmanState{Fl}) where Fl
    kalman_state.llk -= (
        HALF_LOG_2_PI + (log(kalman_state.F) + kalman_state.v^2 / kalman_state.F) / 2
    )
    return kalman_state
end

function update_kalman_state!(
    kalman_state::UnivariateKalmanState{Fl},
    y::Fl,
    Z::Vector{Fl},
    T::Matrix{Fl},
    H::Fl,
    RQR::Matrix{Fl},
    d::Fl,
    c::Vector{Fl},
    skip_llk_instants::Int,
    tol::Fl,
    t::Int,
) where Fl
    if isnan(y)
        kalman_state.v = NaN
        update_F!(kalman_state, Z, H)
        repeat_a_in_att!(kalman_state)
        update_a!(kalman_state, T, c)
        repeat_P_in_Ptt!(kalman_state)
        update_P!(kalman_state, T, RQR)
        kalman_state.steady_state = false # Not on steadystate anymore
    elseif kalman_state.steady_state
        update_v!(kalman_state, y, Z, d)
        update_att!(kalman_state, Z)
        update_a!(kalman_state, T, c)
        if t > skip_llk_instants
            update_llk!(kalman_state)
        end
    else
        update_v!(kalman_state, y, Z, d)
        update_F!(kalman_state, Z, H)
        update_att!(kalman_state, Z)
        update_a!(kalman_state, T, c)
        update_Ptt!(kalman_state)
        update_P!(kalman_state, T, RQR)
        check_steady_state(kalman_state, tol)
        if t > skip_llk_instants
            update_llk!(kalman_state)
        end
    end
    return kalman_state
end

function update_kalman_state!(
    kalman_state::UnivariateKalmanState{Fl},
    y::Fl,
    Z::Vector{Fl},
    T::Matrix{Fl},
    H::Fl,
    R::Matrix{Fl},
    Q::Matrix{Fl},
    d::Fl,
    c::Vector{Fl},
    skip_llk_instants::Int,
    t::Int,
) where Fl
    if kalman_state.steady_state
        update_v!(kalman_state, y, Z, d)
        update_att!(kalman_state, Z)
        update_a!(kalman_state, T, c)
    elseif isnan(y)
        kalman_state.v = NaN
        update_F!(kalman_state, Z, H)
        repeat_a_in_att!(kalman_state)
        update_a!(kalman_state, T, c)
        repeat_P_in_Ptt!(kalman_state)
        update_P!(kalman_state, T, R, Q)
    else
        update_v!(kalman_state, y, Z, d)
        update_F!(kalman_state, Z, H)
        update_att!(kalman_state, Z)
        update_a!(kalman_state, T, c)
        update_Ptt!(kalman_state)
        update_P!(kalman_state, T, R, Q)
    end
    if t > skip_llk_instants
        update_llk!(kalman_state)
    end
    return kalman_state
end

function optim_kalman_filter(
    sys::StateSpaceSystem, filter::UnivariateKalmanFilter{Fl}
) where Fl
    return filter_recursions!(
        filter.kalman_state, sys, filter.steadystate_tol, filter.skip_llk_instants
    )
end

function kalman_filter!(
    filter_output::FilterOutput, sys::StateSpaceSystem, filter::UnivariateKalmanFilter{Fl}
) where Fl
    filter_recursions!(
        filter_output,
        filter.kalman_state,
        sys,
        filter.steadystate_tol,
        filter.skip_llk_instants,
    )
    return filter_output
end

function filter_recursions!(
    kalman_state::UnivariateKalmanState{Fl},
    sys::LinearUnivariateTimeInvariant,
    steadystate_tol::Fl,
    skip_llk_instants::Int,
) where Fl
    RQR = sys.R * sys.Q * sys.R'
    @inbounds for t in eachindex(sys.y)
        update_kalman_state!(
            kalman_state,
            sys.y[t],
            sys.Z,
            sys.T,
            sys.H,
            RQR,
            sys.d,
            sys.c,
            skip_llk_instants,
            steadystate_tol,
            t,
        )
    end
    return kalman_state.llk
end

# TODO this doesn't use steadystate_tol
function filter_recursions!(
    kalman_state::UnivariateKalmanState{Fl},
    sys::LinearUnivariateTimeVariant,
    steadystate_tol::Fl,
    skip_llk_instants::Int,
) where Fl
    @inbounds for t in eachindex(sys.y)
        update_kalman_state!(
            kalman_state,
            sys.y[t],
            sys.Z[t],
            sys.T[t],
            sys.H[t],
            sys.R[t],
            sys.Q[t],
            sys.d[t],
            sys.c[t],
            skip_llk_instants,
            t,
        )
    end
    return kalman_state.llk
end

function filter_recursions!(
    filter_output::FilterOutput,
    kalman_state::UnivariateKalmanState{Fl},
    sys::LinearUnivariateTimeInvariant,
    steadystate_tol::Fl,
    skip_llk_instants::Int,
) where Fl
    RQR = sys.R * sys.Q * sys.R'
    save_a1_P1_in_filter_output!(filter_output, kalman_state)
    @inbounds for t in eachindex(sys.y)
        update_kalman_state!(
            kalman_state,
            sys.y[t],
            sys.Z,
            sys.T,
            sys.H,
            RQR,
            sys.d,
            sys.c,
            skip_llk_instants,
            steadystate_tol,
            t,
        )
        save_kalman_state_in_filter_output!(filter_output, kalman_state, t)
    end
    return filter_output
end

# TODO this doesn't use steadystate_tol
function filter_recursions!(
    filter_output::FilterOutput,
    kalman_state::UnivariateKalmanState{Fl},
    sys::LinearUnivariateTimeVariant,
    steadystate_tol::Fl,
    skip_llk_instants::Int,
) where Fl
    save_a1_P1_in_filter_output!(filter_output, kalman_state)
    @inbounds for t in eachindex(sys.y)
        update_kalman_state!(
            kalman_state,
            sys.y[t],
            sys.Z[t],
            sys.T[t],
            sys.H[t],
            sys.R[t],
            sys.Q[t],
            sys.d[t],
            sys.c[t],
            skip_llk_instants,
            t,
        )
        save_kalman_state_in_filter_output!(filter_output, kalman_state, t)
    end
    return filter_output
end


text:
"""
    module test_StaticArrays

Evaluate and test properties of StaticArrays
"""
module test_StaticArrays

import ..ModiaMath

function A_add_B!(C, A, B)
    for i = 1:length(A)
        C[i] = A[i] + B[i]
    end
end

function rot(n::AbstractVector, angleRad::Number)
    s  = sin(angleRad)
    c  = cos(angleRad)
    e  = normalize(n)
    es = e * s
    ee = e * e'
    R  = ee + (eye(3) - ee) * c - [ 0.0   -es[3]  es[2];
                                    es[3]  0.0   -es[1];
                                   -es[2]  es[1]  0.0 ]
end

mutable struct ModelVariables
    R1::AbstractMatrix
    R2::MMatrix{3,3,Float64,9}
    R3::SMatrix{3,3,Float64,9}
    R4::SMatrix{3,3,Float64,9}
    v1a::AbstractVector
    v1b::AbstractVector
    v1c::AbstractVector
    v2a::SVector{3,Float64}
    v2b::SVector{3,Float64}
    v2c::SVector{3,Float64}
    v3a::SVector{3,Float64}
    v3b::SVector{3,Float64}
    v3c::SVector{3,Float64}

    function ModelVariables()
        R1  = rot([1,2,-1], 0.5)
        R2  = MMatrix{3,3,Float64,9}(2 * R1)
        R3  = SMatrix{3,3,Float64,9}(3 * R1)
        R4  = SMatrix{3,3,Float64,9}(4 * R1)

        v1a = [1.1, 2.2, 3.3]
        v1b = v1a
        v1c = 3 * copy(v1a)

        v2a = SVector{3,Float64}(v1a)
        v2b = v2a
        v2c = 3 * copy(v2a)

        v3a = SVector{3,Float64}(v1a)
        v3b = v3a
        v3c = 3 * copy(v3a)

        new(R1, R2, R3, R4, v1a, v1b, v1c, v2a, v2b, v2c, v3a, v3b, v3c)
    end
end

const var = ModelVariables()
const nmax = 100000

# Evaluate overhead of for-loop
function test_for_loop(var)
    j = 0
    for i = 1:nmax
        j += 1
    end
    return nothing
end

(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed test_for_loop(var)
(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed test_for_loop(var)
println("\ntest_for_loop: elapsedTime = ", elapsedTime, ", bytesAllocated = ", bytesAllocated, "\n")

# Evaluate matrix multiplication and assignment
function mult_StandardArray_assign(var::ModelVariables)::Nothing
    for i = 1:nmax
        var.v1a = var.v1c + i * var.R1 * var.v1c
    end
    return nothing
end

function mult_StandardArray_dotAssign(var::ModelVariables)::Nothing
    for i = 1:nmax
        var.v1a .= var.v1c + i * var.R1 * var.v1c
    end
    return nothing
end

function mult_StandardArray_elementAssign(var::ModelVariables)::Nothing
    for i = 1:nmax
        var.v1a[1:3] = var.v1c + i * var.R1 * var.v1c
    end
    return nothing
end

function mult_StandardArray_mult(var::ModelVariables)::Nothing
    for i = 1:nmax
        A_mul_B!(var.v1a, var.R1, var.v1c)
        A_add_B!(var.v1a, var.v1c, i * var.v1a)
    end
    return nothing
end

function mult_MArray_assign(var::ModelVariables)::Nothing
    for i = 1:nmax
        var.v2a = var.v2c + i * var.R2 * var.v2c
    end
    return nothing
end

function mult_MArray_dotAssign(var::ModelVariables)::Nothing
    for i = 1:nmax
        var.v2a .= var.v2c + i * var.R2 * var.v2c
    end
    return nothing
end

function mult_MArray_elementAssign(var::ModelVariables)::Nothing
    for i = 1:nmax
        var.v2a[1:3] = var.v2c + i * var.R2 * var.v2c
    end
    return nothing
end

function mult_MArray_mult(var::ModelVariables)::Nothing
    for i = 1:nmax
        A_mul_B!(var.v2a, var.R2, var.v2c)
        A_add_B!(var.v2a, var.v2c, i * var.v2a)
    end
    return nothing
end

function mult_SArray_assign(var::ModelVariables)::Nothing
    for i = 1:nmax
        var.v3a = var.v3c + i * var.R3 * var.v3c
    end
    return nothing
end

function mult_SArray_assignWithMMatrix(var::ModelVariables)::Nothing
    for i = 1:nmax
        var.v3a = var.v3c + i * var.R2 * var.v3c
    end
    return nothing
end

#function mult_SArray_dotAssign(var::ModelVariables)
#   for i = 1:nmax
#      var.v3a .= var.R1*var.v3c
#   end
#end

#function mult_SArray_elementAssign(var::ModelVariables)
#   for i = 1:nmax
#      var.v3a[1:3] = var.R1*var.v3c
#   end
#end

table = DataFrames.DataFrame(assignment=String[], relTime=Float64[], elapsedTime=Float64[], bytesAllocated=Int[], link=Bool[])

(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_StandardArray_assign(var)
var.v1b = var.v1a
(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_StandardArray_assign(var)
push!(table, ["Vector: va = vc+R*vc", 0.0, elapsedTime, bytesAllocated, var.v1b === var.v1a])

(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_StandardArray_dotAssign(var)
var.v1b = var.v1a
(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_StandardArray_dotAssign(var)
push!(table, ["Vector: va .= vc+R*vc", 0.0, elapsedTime, bytesAllocated, var.v1b === var.v1a])

(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_StandardArray_elementAssign(var)
var.v1b = var.v1a
(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_StandardArray_elementAssign(var)
push!(table, ["Vector: va[1:3] = vc+R*vc", 0.0, elapsedTime, bytesAllocated, var.v1b === var.v1a])

(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_StandardArray_mult(var)
var.v1b = var.v1a
(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_StandardArray_mult(var)
push!(table, ["Vector: A_mult_B!(va,R,vc)", 0.0, elapsedTime, bytesAllocated, var.v1b === var.v1a])

(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_MArray_assign(var)
var.v2b = var.v2a
(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_MArray_assign(var)
push!(table, ["SVector: va = vc+R*vc", 0.0, elapsedTime, bytesAllocated, var.v2b === var.v2a])

(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_MArray_dotAssign(var)
var.v2b = var.v2a
(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_MArray_dotAssign(var)
push!(table, ["SVector: va .= vc+R*vc", 0.0, elapsedTime, bytesAllocated, var.v2b === var.v2a])

(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_MArray_elementAssign(var)
var.v2b = var.v2a
(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_MArray_elementAssign(var)
push!(table, ["SVector: va[1:3] = vc+R*vc", 0.0, elapsedTime, bytesAllocated, var.v2b === var.v2a])

(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_MArray_mult(var)
var.v2b = var.v2a
(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_MArray_mult(var)
push!(table, ["SVector: A_mult_B!(va,R,vc)", 0.0, elapsedTime, bytesAllocated, var.v2b === var.v2a])

var.v3b = var.v3a
(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_SArray_assignWithMMatrix(var)
(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_SArray_assignWithMMatrix(var)
push!(table, ["SVector/MMatrix: va = vc+R*vc", 0.0, elapsedTime, bytesAllocated, var.v3b === var.v3a])

var.v3b = var.v3a
(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_SArray_assign(var)
(val, elapsedTime, bytesAllocated, dummy1, dummy2) = @timed mult_SArray_assign(var)
push!(table, ["SVector: va = vc+R*vc", 0.0, elapsedTime, bytesAllocated, var.v3b === var.v3a])

elapsed = table[:,3]
imin = indmin(elapsed)
vmin = elapsed[imin]
for i in 1:length(elapsed)
    table[i,2] = table[i,3] / vmin
end

println("\nnmax = ", nmax, "\n\ntable = ", table)

end


text:
using DataDecomposition;
using MLKernels;
nc = 3;
p = 3;
k = 4;
y = repeat(Int64[i for i = 1:k],inner=[nc],outer=[1]);
X = vcat([rand(nc,p) .+ 10rand(1,3) for i = 1:k]...);
freq = Float64[1/k for i=1:k];
lda!(copy(X), y, k, freq)


text:
function ProperBigFractionFormat()
    return ProperBigFractionFormat(())
end

function ProperBigFractionFormat(arg0::NumberFormat)
    return ProperBigFractionFormat((NumberFormat,), arg0)
end

function ProperBigFractionFormat(arg0::NumberFormat, arg1::NumberFormat, arg2::NumberFormat)
    return ProperBigFractionFormat((NumberFormat, NumberFormat, NumberFormat), arg0, arg1, arg2)
end

function format(obj::ProperBigFractionFormat, arg0::BigFraction, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (BigFraction, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function get_whole_format(obj::ProperBigFractionFormat)
    return jcall(obj, "getWholeFormat", NumberFormat, ())
end

function parse(obj::ProperBigFractionFormat, arg0::JString, arg1::ParsePosition)
    return jcall(obj, "parse", BigFraction, (JString, ParsePosition), arg0, arg1)
end


text:
export RandomHypergraph

"""
`RandomHypergraph(n::Int,k::Int,p::Real=0.5)` creates a random
`k`-uniform hypergraph with `n` vertices. a `k`-element
subset of `{1,2,...,n}` is an edge (independently) with
probility `p`.
"""
function RandomHypergraph(n::Int, k::Int, p::Real = 0.5)
    H = IntHypergraph(n)
    for e in subsets(1:n, k)
        if rand() < p
            add!(H, Set(e))
        end
    end
    return H
end


text:
# Extract some info from participants' data
# (i.e SARSPEICZVG_all_fMRI.csv)
using DelimitedFiles

function csv_extractallevents()
    # Process participants' data and create the files that are now in
    # /src/mcmc_rl_fit/projects/fmri/someresults/fits/all_events

    data_path = "./mcmc_rl_fit/projects/fmri/data/SARSPEICZVG_all_fMRI.csv"

    base_path = makehomesavepath("all_events")
    mkdir(base_path)

    all_data, all_stats = RL_Fit.load_SARSPEI(data_path)

    all_subject_ids = all_stats.all_subject_ids

    for subj_id in all_subject_ids

        subj_data = RL_Fit.filter_by_person(all_data, subj_id)

        subj_path = joinpath(base_path, "subj_"*string(subj_id))

        mkdir(subj_path)
        # mkpath(subj_path)
        writedlm(subj_path*"/states.csv", subj_data[:, Int(RL_Fit.c_state)])
        writedlm(subj_path*"/actions.csv", subj_data[:, Int(RL_Fit.c_action)])
        writedlm(subj_path*"/rewards.csv", subj_data[:, Int(RL_Fit.c_reward)])
        writedlm(subj_path*"/is_surprise_trial.csv", subj_data[:, Int(RL_Fit.c_surprise_trial)])
        writedlm(subj_path*"/graph_version.csv", subj_data[1, Int(RL_Fit.c_graph_version)])

    end

end

# NOTE: the nr_states are actually the number of actions
#(ie excluding start states)
# And then, the % of surprise trials is the nr_surprise_trials/nr_states
# (ie number of actions = transitions).
function getnrofevents()
    data_path = "./mcmc_rl_fit/projects/fmri/data/SARSPEICZVG_all_fMRI.csv"

    all_data, all_stats = RL_Fit.load_SARSPEI(data_path)

    all_subject_ids = all_stats.all_subject_ids

    nr_states = []
    nr_episodes = []
    nr_surprise_trials = []
    for subj_id in all_subject_ids

        subj_data = RL_Fit.filter_by_person(all_data, subj_id)

        push!(nr_states, length(subj_data[:, Int(RL_Fit.c_state)]))
        push!(nr_episodes, length(findall(subj_data[:, Int(RL_Fit.c_reward)] .== 1)))
        push!(nr_surprise_trials, length(findall(subj_data[:, Int(RL_Fit.c_surprise_trial)] .== 1)))
    end
    nr_states, nr_episodes, nr_surprise_trials
end


text:
module SQLiteGraph

using SQLite: SQLite
import DBInterface: execute
using JSON3: JSON3
using EasyConfig

export DB, Node, Edge

#-----------------------------------------------------------------------------# utils
function single_result_execute(db, stmt, args...)
    ex = execute(db, stmt, args...)
    isempty(ex) ? nothing : first(first(ex))
end

function print_props(io::IO, o::Config)
    for (i,(k,v)) in enumerate(pairs(o))
        if i < 5
            print(io, k, '=', repr(v))
            i == length(o) || print(io, ", ")
        end
    end
    length(o) > 5 && print(io, "…")
end

#-----------------------------------------------------------------------------# Model
struct Node
    id::Int
    labels::Vector{String}
    props::Config
end
Node(id::Int, labels::String...; props...) = Node(id, collect(labels), Config(props))
Node(row::SQLite.Row) = Node(row.id, split(row.labels, ';', keepempty=false), JSON3.read(row.props, Config))
function Base.show(io::IO, o::Node)
    print(io, "Node($(o.id)")
    !isempty(o.labels) && print(io, ", ", join(repr.(o.labels), ", "))
    !isempty(o.props) && print(io, "; "); print_props(io, o.props)
    print(io, ')')
end
args(n::Node) = (n.id, isempty(n.labels) ? "" : join(n.labels, ';'), JSON3.write(n.props))

struct Edge
    source::Int
    target::Int
    type::String
    props::Config
end
Edge(src::Int, tgt::Int, type::String; props...) = Edge(src, tgt, type, Config(props))
Edge(row::SQLite.Row) = Edge(row.source, row.target, row.type, JSON3.read(row.props, Config))
function Base.show(io::IO, o::Edge)
    print(io, "Edge($(o.source), $(o.target), ", repr(o.type))
    !isempty(o.props) && print(io, "; "); print_props(io, o.props)
    print(io, ')')
end
args(e::Edge) = (e.source, e.target, e.type, JSON3.write(e.props))

#-----------------------------------------------------------------------------# Base methods
Base.:(==)(a::Node, b::Node) = all(getfield(a,f) == getfield(b,f) for f in fieldnames(Node))
Base.:(==)(a::Edge, b::Edge) = all(getfield(a,f) == getfield(b,f) for f in fieldnames(Edge))

Base.pairs(o::T) where {T<: Union{Node, Edge}} = (f => getfield(o,f) for f in fieldnames(T))

Base.NamedTuple(o::Union{Node,Edge}) = NamedTuple(pairs(o))

#-----------------------------------------------------------------------------# DB
struct DB
    sqlitedb::SQLite.DB

    function DB(file::String = ":memory:")
        db = SQLite.DB(file)
        foreach(x -> execute(db, x), [
            "PRAGMA foreign_keys = ON;",
            # nodes
            "CREATE TABLE IF NOT EXISTS nodes (
                id INTEGER NOT NULL UNIQUE PRIMARY KEY,
                labels TEXT NOT NULL,
                props TEXT NOT NULL
            );",
            # edges
            "CREATE TABLE IF NOT EXISTS edges (
                source INTEGER NOT NULL REFERENCES nodes(id),
                target INTEGER NOT NULL REFERENCES nodes(id),
                type TEXT NOT NULL,
                props TEXT NOT NULL,
                PRIMARY KEY (source, target, type)
            );",
            "CREATE INDEX IF NOT EXISTS source_idx ON edges(source);",
            "CREATE INDEX IF NOT EXISTS target_idx ON edges(target);",
            "CREATE INDEX IF NOT EXISTS type_idx ON edges(type);",
        ])
        new(db)
    end
end
function Base.show(io::IO, db::DB)
    print(io, "SQLiteGraph.DB(\"$(db.sqlitedb.file)\") ($(n_nodes(db)) nodes, $(n_edges(db)) edges)")
end

execute(db::DB, args...; kw...) = execute(db.sqlitedb, args...; kw...)

n_nodes(db::DB) = single_result_execute(db, "SELECT Count(*) FROM nodes")
n_edges(db::DB) = single_result_execute(db, "SELECT Count(*) FROM edges")

Base.length(db::DB) = n_nodes(db)
Base.size(db::DB) = (nodes=n_nodes(db), edges=n_edges(db))
Base.lastindex(db::DB) = length(db)
Base.axes(db::DB, i) = size(db)[i]

Broadcast.broadcastable(db::DB) = Ref(db)

#-----------------------------------------------------------------------------# insert!
function Base.insert!(db::DB, node::Node)
    execute(db, "INSERT INTO nodes VALUES(?, ?, json(?))", args(node))
    db
end
function Base.insert!(db::DB, edge::Edge)
    execute(db, "INSERT INTO edges VALUES(?, ?, ?, json(?))", args(edge))
    db
end

#-----------------------------------------------------------------------------# replace!
function Base.replace!(db::DB, node::Node)
    execute(db, "INSERT INTO nodes VALUES(?, ?, json(?)) ON CONFLICT(id) DO UPDATE SET labels=excluded.labels, props=excluded.props", args(node))
    db
end
function Base.replace!(db::DB, edge::Edge)
    execute(db, "INSERT INTO edges VALUES(?, ?, ?, json(?)) ON CONFLICT(source,target,type) DO UPDATE SET props=excluded.props", args(edge))
    db
end

#-----------------------------------------------------------------------------# query
function query(db::DB, select::String, from::String, whr::String, args=nothing)
    stmt = "SELECT $select FROM $from WHERE $whr"
    # @info stmt
    res = isnothing(args) ? execute(db, stmt) : execute(db, stmt, args)
    if isempty(res)
        error("No $from found where: $whr")
    else
        return res
    end
end

#-----------------------------------------------------------------------------# getindex (Node)
Base.getindex(db::DB, i::Integer) = Node(first(query(db, "*", "nodes", "id=$i")))
Base.getindex(db::DB, ::Colon) = (Node(row) for row in query(db, "*", "nodes", "TRUE"))

#-----------------------------------------------------------------------------# getindex (Edge)
# all specified
function Base.getindex(db::DB, i::Integer, j::Integer, type::AbstractString)
    Edge(first(query(db, "*", "edges", "source=$i AND target=$j AND type LIKE '$type'")))
end

# one colon
function Base.getindex(db::DB, i::Integer, j::Integer, ::Colon)
    (Edge(row) for row in query(db, "*", "edges", "source=$i AND target=$j"))
end
function Base.getindex(db::DB, i::Integer, ::Colon, type::AbstractString)
    (Edge(row) for row in query(db, "*", "edges", "source=$i AND type LIKE '$type'"))
end
function Base.getindex(db::DB, ::Colon, j::Integer, type::AbstractString)
    (Edge(row) for row in query(db, "*", "edges", "target=$j AND type LIKE '$type'"))
end

# two colons
function Base.getindex(db::DB, i::Integer, ::Colon, ::Colon)
    (Edge(row) for row in query(db, "*", "edges", "source=$i"))
end
function Base.getindex(db::DB, i::Colon, j::Integer, ::Colon)
    (Edge(row) for row in query(db, "*", "edges", "target=$j"))
end
function Base.getindex(db::DB, ::Colon, ::Colon, type::AbstractString)
    (Edge(row) for row in query(db, "*", "edges", "type LIKE '$type'"))
end

# all colons
Base.getindex(db::DB, ::Colon, ::Colon, ::Colon) = (Edge(row) for row in query(db,"*", "edges", "TRUE"))

#-----------------------------------------------------------------------------# adjacency_matrix
"""
    adjacency_matrix(db, type)

Create the adjacency matrix for a given edge `type`.  If `A[i,j] == true`, there exists an
edge from node `i` to node `j` with type `type`.
"""
function adjacency_matrix(db::DB, type)
    n = n_nodes(db)
    out = falses(n, n)
    for row in execute(db, "SELECT DISTINCT source, target FROM edges WHERE type=?;", (type,))
        out[row.source, row.target] = true
    end
    out
end

end


text:
# This file is a part of project JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using Documenter
using Materials

makedocs(
    modules = [Materials],
    checkdocs = :all,
    strict = true)


text:
include("struct.jl")
include("phase1.jl")
include("phase2.jl")

max_stars_repo_name: fbastin/Jasmin

max_issues_repo_name: fbastin/Jasmin

max_forks_repo_name: fbastin/Jasmin


text:
using StaticArrays, PyCall

import Base.convert

export mat, vecs
export SVec, SMat, JVec, JVecs, JVecF, JVecsF
export zerovecs, maxdist, maxnorm

export unsafe_pyarrayref, safe_pyarrayref

"typealias for a static vector type; currently `StaticArrays.SVector`"
typealias SVec SVector
"typealias for a static matrix type; currently `StaticArrays.SMatrix`"
typealias SMat SMatrix

"`JVec{T}` : 3-dimensional immutable vector"
typealias JVec{T} SVec{3,T}
typealias JVecF JVec{Float64}
typealias JVecI JVec{Int}

"`JVecs{T}` : List of 3-dimensional immutable vectors"
typealias JVecs{T} Vector{JVec{T}}
typealias JVecsF JVecs{Float64}
typealias JVecsI JVecs{Int}

"""
`vecs(V::Matrix)` : convert (as reference) a 3 x N matrix representing
N vectors (e.g. forces) in R³ to a list (vector) of fixed-size-array vectors.

`vecs(V::Vector)` : assumes that V is morally a 3 x N matrix, stored in a long
vector

`vecs(V::Array{T,N})` : If `V` has dimensions 3 x n2 x ... x nN then
it gets converted to an n2 x ... x nN array with JVec{T} entries.
"""
vecs{T}(V::Matrix{T}) = reinterpret(JVec{T}, V, (size(V,2),))
vecs{T}(V::Vector{T}) = reinterpret(JVec{T}, V, (length(V) ÷ 3,))
vecs{T,N}(V::Array{T,N}) = reinterpret(JVec{T}, V, tuple(size(V)[2:end]...))

"""
`mat(X::JVecs)`: convert (as reference) a list (Vector) of
fixed-size-array points or vecs to a 3 x N matrix representing those.

### Usage:
```
X = positions(at)          # returns a Vector{JVec}
X = positions(at) |> mat   # returns a Matrix
```
"""
mat{T}(V::JVecs{T}) = reinterpret(T, V, (3, length(V)))
mat{T}(X::AbstractVector{JVec{T}}) = mat(collect(X))

# rewrite all of this in terms of `convert` (TODO: is this needed?)
convert{T}(::Type{Matrix{T}}, V::JVecs{T}) = mat(V)
convert{T}(::Type{JVecs{T}}, V::Matrix{T}) = vec(V)

# initialise a vector of vecs or points
zerovecs(n::Integer) = zerovecs(Float64, n)
zerovecs(T::Type, n::Integer) = zeros(T, 3, n) |> vecs

"""
maximum of distances between two sets of JVec's, usually positions;
`maximum(a - b for (a,b) in zip(x,y))`
"""
function maxdist{T}(x::JVecs{T}, y::JVecs{T})
   @assert length(x) == length(y)
   return maximum( norm(a-b)  for (a,b) in zip(x,y) )
end

"`maximum(norm(y) for y in x);` typically, x is a vector of forces"
maxnorm{T}(x::JVecs{T}) = maximum( norm.(x) )

# The next function is conversion of python arrays  _by reference_
pyarrayref(a::PyObject; own=false) = pyarrayref(PyArray(a); own=own)
pyarrayref(a::PyArray; own=false) = unsafe_wrap(Array, a.data, reverse(a.dims), own)


text:
using Gadfly
using Colors
using MLLibrary

path=joinpath(Pkg.dir("MLLibrary"),"resources","datasets")
seed=11
threshold =.1
delimiter = ","

#resources = ["3percent-miscategorization","10percent-miscatergorization","abalone.data","Seperable.csv"]
resources=["res"=>"SpamDataPruned.csv" ,"proxy"=>[],"path"=>Pkg.dir("MLLibrary","resources","datasets")]

training, test = MLLibrary.DataSampling.ThresholdSampling.getSets(joinpath(resources["path"],resources["res"]),",",threshold,seed)

# println(size(training), size(test))

training = MLLibrary.DataManipulation.DataTransform.convertToFloat(training)
test = MLLibrary.DataManipulation.DataTransform.convertToFloat(test)

training_data = training[:,1:end-1]
test_data = test[:,1:end-1]
training_labels = training[:,end]
test_labels = test[:,end]

dict, dictProb = MLLibrary.LearningAlgorithms.MultinomialBayes.getDictionary(training_data,training_labels)

predicted_labels = MLLibrary.LearningAlgorithms.MultinomialBayes.classifyDocuments(test_data,dictProb)

# println(size(predicted_labels),size(test_labels))
println(MLLibrary.Utils.ErrorAnalysis.generateConfusionMatrix(test_labels,predicted_labels)
)

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
function compress7(inio, outio)
    nextwritebyte = read(inio, UInt8) & 0x7f
    filled = 7
    while !eof(inio)
        inbyte = read(inio, UInt8)
        write(outio, UInt8(nextwritebyte | inbyte << filled))
        nextwritebyte = inbyte >> (8 - filled)
        filled = (filled + 7) % 8
        if filled == 0
            if eof(inio)
                break
            end
            nextwritebyte = read(inio, UInt8) & 0x7f
            filled = 7
        end
    end
    if filled != 0
        write(outio, UInt8(nextwritebyte))
    end
end

function expand7(inio, outio)
    newbyte = read(inio, UInt8)
    write(outio, UInt8(newbyte & 0x7f))
    residualbyte::UInt8 = newbyte >> 7
    filled = 1
    while !eof(inio)
        inbyte = read(inio, UInt8)
        write(outio, UInt8((residualbyte | inbyte << filled) & 0x7f))
        residualbyte = inbyte >> (7 - filled)
        filled = (filled + 1) % 7
        if filled == 0
            write(outio, UInt8(residualbyte & 0x7f))
            residualbyte = 0
        end
    end
end

str = b"These bit oriented I/O functions can be used to implement compressors and decompressors."
ins = IOBuffer(str)
outs = IOBuffer()
newouts = IOBuffer()
compress7(ins, outs)
seek(outs,0)
expand7(outs, newouts)
println("Initial string of length $(length(str)): ", String(ins.data))
println("Compressed to length $(length(outs.data)) on line below:\n", String(outs.data))
println("Decompressed string: ", String(newouts.data))


text:
mutable struct Cached{T}
    value::Nullable{T}

    Cached{T}() where T = new(Nullable{T}())
    Cached{T}(val::T) where T = new(Nullable{T}(val))
end

function set!(cv::Cached{T}, val::T) where T
    cv.value = Nullable(val)
    return val
end
iscached(cv::Cached) = !isnull(cv.value)

Base.get(cv::Cached) = get(cv.value)
function Base.get(f::Function, cv::Cached, args...)
    if iscached(cv)
        return get(cv.value)
    else
        return set!(cv, f(args...))
    end
end