text:
"""
EventSimulation is an event-based discrete event simulation engine.
"""
module EventSimulation

using Random

export # core.jl
       Action, AbstractState, EmptyState, Scheduler,
       register!, repeat_register!,
       bulk_register!, repeat_bulk_register!,
       interrupt!, terminate!, go!,

       # resource.jl, queue.jl
       AbstractReservoir, # defined in this file
       ResourceRequest, SimResource, SimQueue,
       request!, waive!, provide!, withdraw!,

       # prioritytime.jl
       PriorityTime

include("core.jl")

"""
Abstract class for reservoirs.
`SimQueue` and `SimResource` are concrete types implementing it.
"""
abstract type AbstractReservoir end

"""
    dispatch!(s, r)
    dispatch!(s, q)

Internal function used for dispatching requests in `SimResource` and
`SimQueue`. Puts appropriate `Action`s in `s` immediately.
"""
function dispatch! end

"""
    request!(s, r, quantity, request)
    request!(s, q, request)

Function used to register request for resource in `SimResource`
or object from `SimQueue`.

In `SimResource` requested `quantity` must be provided and
`request` accepts only `Scheduler` argument (it must know what it wanted).
Returns tuple of:
* `true` if successfull and `false` when too many requests were made
* `ResourceRequest` object created

In `SimResource` function `request` must accept one argument `Scheduler`.
In `SimQueue` function `request` must accept two arguments `Scheduler` and object.
"""
function request! end

"""
    waive!(r, res_request)
    waive!(q, request)

Allows to remove first occurence that would be served
of `res_request` from `SimResource` or `request` from `SimQueue`.

Returns `true` on success and `false` if `res_request`
or `request` respectively was not found.
"""
function waive! end

"""
    provide!(s, r, quantity)
    provide!(s, q, object)

Allows to fill `SimResource` with `quantity` or `SimQueue` with `object`.

In `SimResource` changes the balance of `r.quantity`. Given quantity may be
any number, but the balance of `SimResource` will be changed
only in `lo`-`hi` range. Returns the actual change in `SimResource` balance.

In `SimQueue` adds `object` to `q.queue`. Returns `true` on success and
`false` if there were too many objects in queue already.
"""
function provide! end

include("resource.jl")
include("queue.jl")
include("prioritytime.jl")

end # module


text:
"""
Return a vector of stops `v` with positive demand from stop `u`.

It ensures that `np.odmatrix[u,v] > 0`, and that `u` itself does not appear.
"""
nonzerodests(np::TransitNetwork, u::Int) =
    filter!(i -> i != u, findall(np.odmatrix[u,:] .> 0))

"""
Return a vector of commutes with nonzero demand.
"""
commutes(np::TransitNetwork) =
    unique(
        vcat([[(min(u,v),max(u,v)) for v in nonzerodests(np,u)] for u in 1:np.nstations]...)
    )

"""
Return a vector of commutes with demand > lb
"""
function commutes(np::TransitNetwork, lb::Int)
    coms = Vector{Tuple{Int,Int}}()
    for (u,v) in commutes(np)
        if demand(np, (u,v)) > lb
            push!(coms, (u,v))
        end
    end
    coms
end

"""
Return a vector of commutes with nonzero dual variables.
"""
function commutes(rmp::MasterProblem)
    p = JuMP.getdual(rmp.model[:choseline])
    alluvs = commutes(rmp.np)
    smalluvs = Vector{Tuple{Int,Int}}()
    for (u,v) in alluvs
        if p[(u,v)] > 0
            push!(smalluvs, (u,v))
        end
    end
    smalluvs
end

"Total demand between u and v"
demand(np::TransitNetwork, uv::Tuple{Int,Int}) =
    np.odmatrix[uv[1],uv[2]] + np.odmatrix[uv[2],uv[1]]

"""
Return the unique lines of `linelist`.

It ensures that both a line and its reverse will not appear.
"""
function uniquelines(linelist::Vector{Vector{Int}})
    uniquelines = Vector{Int}[]
    for line in linelist
        isunique = true
        for line2 in uniquelines
            if length(line) == length(line2) && all(sort(line) .== sort(line2))
                isunique = false
                break
            end
        end
        isunique && push!(uniquelines, copy(line))
    end

    uniquelines
end

"""
Return whether `line` or its reversal appears in `linelist`.
"""
linein(line::Vector{Int}, linelist::Vector{Vector{Int}}) = in(line, linelist) || in(reverse(line), linelist)

"""
The total distance along a line
"""
linecost(np::TransitNetwork, line::Vector{Int}) =
    sum(edgecost(np, line[i], line[i+1]) for i in 1:(length(line)-1))

"""
Returns distance between stations `u` and `v` based on the `gridtype`.

If `gridtype` is `:latlong`, it returns the distance in miles. If gridtype is
`:euclidean`, it returns the euclidean distance in latlng space.
"""
function edgecost(
        np::TransitNetwork,
        u::Int,
        v::Int
    )
    if np.dists != nothing
        return haskey(np.dists, (u,v)) ? np.dists[u,v] : Inf
    elseif np.gridtype == :latlong
        return TN.haversinedistance(np, u, v)
    elseif np.gridtype == :euclidean
        return norm(np.latlon[u,:] - np.latlon[v,:])
    else
        error("Unrecognized gridtype: $np.gridtype. Use :latlong or :euclidean.")
    end
end

"""
Returns the distance of the shortest path between u and v.
If distance is not specified in TransitNetwork, calculates a default
value based on directly going from u to v via haversine or euclidean distance.
"""
function spcost(
    np::TransitNetwork,
    u::Int,
    v::Int)
    if np.spdists != nothing
        return haskey(np.spdists, (u,v)) ? np.spdists[u,v] : Inf
    else
        return edgecost(np,u,v)
    end
end

"""
Returns the direction of edge (u, v) projected onto cartesian coordinates
"""
function dir(np::TransitNetwork, u::Int, v::Int)
    if np.gridtype == :latlong
        latu, lonu = np.latlon[u,:]
        latv, lonv = np.latlon[v,:]
        x = cosd(latv) * sind(lonv-lonu)
        y = cosd(latu) * sind(latv) - sind(latu) * cosd(latv) * cosd(lonv-lonu)
        return [x,y]
    elseif gridtype == :euclidean
        return np.latlon[v,:] - np.latlon[u,:]
    else
        error("Unrecognized gridtype: $gridtype. Use :latlong or :euclidean.")
    end
end

ext: jl


text:
using OrdinaryDiffEq, ParameterizedFunctions, Test, Random

f = @ode_def LotkaVolterra begin
  dx = 1.5x - x*y
  dy = -3y + x*y
end

prob = ODEProblem(f,big.([1.0;1.0]),(big(0.0),big(10.0)))

using DiffEqDevTools

dts = big(1/2).^(8:-1:4)
test_setup = Dict(:alg=>Vern9(),:reltol=>1e-25,:abstol=>1e-25)
sim1 = analyticless_test_convergence(dts,prob,Tsit5(),test_setup)
sim2 = analyticless_test_convergence(dts,prob,Vern9(),test_setup)

@test sim1.𝒪est[:final]-5 < 0.2
@test sim2.𝒪est[:final]-9 < 0.2

function f2(du,u,p,t)
  du .= 1.01u
end
function g2(du,u,p,t)
  du .= 1.01u
end
prob = SDEProblem(f2,g2,[1.0;1.0],(0.0,10.0))

using StochasticDiffEq

dts = (1/2).^(6:-1:3)
test_dt = 1/2^8
Random.seed!(100)
sim1 = analyticless_test_convergence(dts,prob,SRIW1(),test_dt,trajectories=400)
@test sim1.𝒪est[:final]-1.5 < 0.3


text:
@doc raw"""
    ExponentialKernel([α=1])

The exponential kernel is given by the formula:

```math
\kappa(\mathbf{x},\mathbf{y}) = \exp\left(-\alpha ||\mathbf{x} - \mathbf{y}||\right)
\qquad \alpha > 0
```

where ``\alpha`` is a scaling parameter of the Euclidean distance. The exponential kernel,
also known as the Laplacian kernel, is an isotropic Mercer kernel. The constructor is
aliased by `LaplacianKernel`, so both names may be used:
"""
struct ExponentialKernel{T<:AbstractFloat} <: MercerKernel{T}
    alpha::HyperParameter{T}
    ExponentialKernel{T}(α::Real) where {T<:AbstractFloat} = new{T}(
        HyperParameter(convert(T,α), interval(OpenBound(zero(T)), nothing))
    )
end
ExponentialKernel(α::T=1.0) where {T<:Real} = ExponentialKernel{floattype(T)}(α)
LaplacianKernel = ExponentialKernel

@inline exponentialkernel(z::T, α::T) where {T<:AbstractFloat} = exp(-α*sqrt(z))

@inline basefunction(::ExponentialKernel) = SquaredEuclidean()
@inline function kappa(κ::ExponentialKernel{T}, z::T) where {T<:AbstractFloat}
    exponentialkernel(z, getvalue(κ.alpha))
end


text:
"Particle Filter module."
module BaytesFilters

############################################################################################
#Import modules

import BaytesCore:
    BaytesCore,
    update!,
    infer,
    results,
    init,
    init!,
    propose,
    propose!,
    propagate,
    propagate!,
    result!,
    get_result,
    generate_showvalues,
    ParameterWeighting,
    ResamplingMethod,
    resample!,
    weight!,
    shuffle!,
    shuffle_forward!,
    shuffle_backward!

using BaytesCore:
    BaytesCore,
    AbstractAlgorithm,
    AbstractTune,
    AbstractConfiguration,
    AbstractDiagnostics,
    AbstractKernel,
    AbstractConstructor,
    Updater,
    Iterator,
    Accumulator,
    UpdateBool,
    UpdateTrue,
    UpdateFalse,
    logsumexp,
    logaddexp,
    logmeanexp,
    issmaller,
    grab,
    ArrayConfig,
    to_NamedTuple,
    update,
    ParameterWeights,
    draw!,
    normalize!,
    computeESS,
    randcat,
    ChainsTune,
    ParameterBuffer,
    SampleDefault,
    to_Tuple

import ModelWrappers: ModelWrappers, predict, dynamics
using ModelWrappers:
    ModelWrappers, ModelWrapper, Tagged, Objective, DiffObjective, ℓObjectiveResult, generate

import Base: max, push!, resize!
import Random: Random, rand!

using DocStringExtensions:
    DocStringExtensions, TYPEDEF, TYPEDFIELDS, FIELDS, SIGNATURES, FUNCTIONNAME
using ArgCheck: ArgCheck, @argcheck
using UnPack: UnPack, @unpack, @pack!

using Random: Random, AbstractRNG, GLOBAL_RNG
using Distributions: Distributions, logpdf, Categorical
using ElasticArrays: ElasticArrays, ElasticMatrix
using Statistics: Statistics, mean, var

############################################################################################
# Import sub-folder
include("ParticleFilter/ParticleFilter.jl")
include("KalmanFilter/KalmanFilter.jl")

############################################################################################
# Export
export
    # BaytesCore
    UpdateBool,
    UpdateTrue,
    UpdateFalse,
    propose,
    propose!,
    propagate,
    propagate!,
    resample!,
    ParameterWeighting,
    ResamplingMethod,
    shuffle!,
    shuffle_forward!,
    shuffle_backward!,
    SampleDefault

    # ModelWrappers
    predict,
    dynamics

end


text:
#!/usr/bin/env julia

# Raytracer.jl
# Raytracing for the generation of photorealistic images in Julia
# Copyright (c) 2021 Samuele Colombo, Paolo Galli

# CLI tool for to manage through Raytracer.jl package the generation
# and rendering of photorealistic images

using Pkg
Pkg.activate(normpath(@__DIR__))

using Raytracer

using ArgParse, ProgressMeter
using FileIO:
    File, @format_str, query

###########
# ArgParse

function print_error(args...)
    printstyled("Error: ", bold=true, color=:red)
    print(args..., "\n")
end

function parse_commandline_error_handler(settings::ArgParseSettings, err, err_code::Int = 1)
    if "out of range input for" |> occursin(err.text)
        s = reduce(replace, ["out of range input for " => "", " " => ""], init=err.text)
        parameter, value = split(s, ":")
        if parameter == "input-file"
            print_error("the parameter \"$(parameter)\" (value: $(value)) must be an existing PFM file")
        elseif parameter == "input-script"
            print_error("the parameter \"$(parameter)\" (value: $(value)) must be an existing SceneLang script")
        elseif parameter == "output-file"
            print_error("the parameter \"$(parameter)\" (value: $(value)) must be a file name, not a path")
        elseif parameter ∈ ["--alpha", "-a", "--gamma", "-g", "--luminosity", "-l", "--delta-t", "--n-frames", "--fps", "-f"]
            print_error("the parameter \"$(parameter)\" (value: $(value)) must be > 0")
        elseif parameter == "time-limits"
            print_error("the parameter \"$(parameter)\" (value: $(value)) must be composed of 2 Float32 values divided by a ':' (e.g., \"0:10\")")
        else
            print_error(err.text)
        end
    else
        print_error(err.text)
    end

    println(stderr, "\n", usage_string(settings), "\n")
    exit(err_code)
end

function parse_commandline()
    s = ArgParseSettings()
    s.description = "Raytracing for the generation of photorealistic images in Julia."
    s.exc_handler = parse_commandline_error_handler
    s.version = "Raytracer.jl version: $(@project_version)"
    s.add_version = true

    # main

    @add_arg_table! s begin
        "render"
            action = :command
            help = "render an image from a SceneLang script"
        "tonemapping"
            action = :command
            help = "apply tone mapping to a pfm image and save it to file"
        "docs"
            action = :command
            help = "show documentation links"
    end

    # render

    s["render"].description = "Render an image from a SceneLang script."
    @add_arg_table! s["render"] begin
        "image"
            action = :command
            help = "render an image from a SceneLang script"
        "animation"
            action = :command
            help = "create an animation as mp4 video from a SceneLang script (require ffmpeg)"
    end

    # render image

    s["render"]["image"].description = "Render an image from a SceneLang script."
    @add_arg_table! s["render"]["image"] begin
        "input-script"
            help = "path to the input SceneLang script"
            arg_type = String
            range_tester = x -> x[end-2:end] == ".sl" && isfile(x)
            required = true
        "--output-file", "-O"
            help = "output image name, without extension"
            arg_type = String
            default = "out"
        "--time", "-t"
            help = "time for script"
            arg_type = Float32
            default = 0f0
        "--vars", "-v"
            help = "add variables using the SET SceneLang syntax; multiple variables separated with commas (e.g., \"a_variable 2, a_material Material()\")"
            arg_type = String
            default = ""
        "--force", "-f"
            help = "force overwrite"
            action = :store_true
    end
    add_arg_group!(s["render"]["image"], "tonemapping settings");
    @add_arg_table! s["render"]["image"] begin
        "--with-tonemapping"
            help = "apply the tone mapping process"
            action = :store_true
        "--ldr-extension", "-e"
            help = "only with \"--with-tonemapping\": extension of the generated ldr image (e.g., \"jpg\" or \"png\")"
            arg_type = String
            default = "jpg"
        "--alpha", "-a"
            help = "only with \"--with-tonemapping\": scaling factor for the normalization process"
            arg_type = Float32
            default = 0.75f0
            range_tester = x -> x > 0
        "--gamma", "-g"
            help = "only with \"--with-tonemapping\": gamma value for the tone mapping process"
            arg_type = Float32
            default = 1f0
            range_tester = x -> x > 0
        "--luminosity", "-l"
            help = "only with \"--with-tonemapping\": luminosity for the tone mapping process (-1 = auto)"
            arg_type = Float32
            default = -1f0
            range_tester = x -> x == -1 || x > 0
    end

    # render animation

    s["render"]["animation"].description = "Create an animation as mp4 video from a SceneLang script (require ffmpeg)."
    @add_arg_table! s["render"]["animation"] begin
        "input-script"
            help = "path to the input SceneLang script"
            arg_type = String
            range_tester = x -> x[end-2:end] == ".sl"
            required = true
        "time-limits"
            help = "time limits as \"t_start:t_end\" (e.g., \"0:10\")"
            arg_type = String
            range_tester = x -> (y=Tuple(parse.(Float32, split(x, ":"))); length(y) == 2)
            required = true
        "--output-dir", "-F"
            help = "output directory"
            arg_type = String
            default = "animation"
        "--output-file", "-O"
            help = "name of saved frames and video, without extension"
            arg_type = String
            default = "out"
            range_tester = x -> isempty(dirname(x))
        "--vars", "-v"
            help = "add variables using the SET SceneLang syntax; multiple variables separated with commas (e.g., \"a_variable 2, a_material Material()\")"
            arg_type = String
            default = ""
        "--force", "-f"
            help = "force overwrite"
            action = :store_true
        "--fps", "-r"
            help = "FPS (frame-per-second) of the output video"
            arg_type = Int
            default = 15
            range_tester = x -> x > 0
    end
    add_arg_group!(s["render"]["animation"], "delta animation (mutually exclusive)", exclusive=true, required=true);
    @add_arg_table! s["render"]["animation"] begin
        "--delta-t"
            help = "proceed from \"t_start\" to \"t_end\" in steps of \"delta-t\""
            arg_type = Float32
            range_tester = x -> x > 0
        "--n-frames"
            help = "number of frames between \"t_start\" and \"t_end\""
            arg_type = Int
            range_tester = x -> x > 0
    end
    add_arg_group!(s["render"]["animation"], "tonemapping settings");
    @add_arg_table! s["render"]["animation"] begin
        "--ldr-extension", "-e"
            help = "extension of the generated ldr image (e.g., \"jpg\" or \"png\")"
            arg_type = String
            default = "jpg"
        "--alpha", "-a"
            help = "scaling factor for the normalization process"
            arg_type = Float32
            default = 0.75f0
            range_tester = x -> x > 0
        "--gamma", "-g"
            help = "gamma value for the tone mapping process"
            arg_type = Float32
            default = 1f0
            range_tester = x -> x > 0
        "--luminosity", "-l"
            help = "luminosity for the tone mapping process (-1 = auto)"
            arg_type = Float32
            default = -1f0
            range_tester = x -> x == -1 || x > 0
    end

    # tonemapping

    s["tonemapping"].description = "Apply tone mapping to a pfm image and save it to file."
    @add_arg_table! s["tonemapping"] begin
        "input-file"
            help = "path to input file, it must be a PFM file"
            arg_type = String
            range_tester = x -> typeof(query(x))<:File{format"PFM"} && isfile(x)
            required = true
        "output-file"
            help = "output file name"
            arg_type = String
            required = true
        "--force", "-f"
            help = "force overwrite"
            action = :store_true
    end
    add_arg_group!(s["tonemapping"], "tonemapping settings");
    @add_arg_table! s["tonemapping"] begin
        "--alpha", "-a"
            help = "scaling factor for the normalization process"
            arg_type = Float32
            default = 0.5f0
            range_tester = x -> x > 0
        "--gamma", "-g"
            help = "gamma value for the tone mapping process"
            arg_type = Float32
            default = 1f0
            range_tester = x -> x > 0
        "--luminosity", "-l"
            help = "only with \"--with-tonemapping\": luminosity for the tone mapping process (-1 = auto)"
            arg_type = Float32
            default = -1f0
            range_tester = x -> x == -1 || x > 0
    end

    # docs

    s["docs"].description = "Show documentation links."
    @add_arg_table! s["docs"] begin
        "--dev"
            help = "documentation of the dev version"
            action = :store_true
    end

    # end

    parse_args(s)
end

####################
# Utility functions

function renderimage(options::Dict{String, Any})
    printstyled("Raytracer.jl image rendering\n\n", bold=true)
    println("Number of threads: $(Threads.nthreads())")

    options["input-script"] = normpath(options["input-script"])
    options["output-file"] = normpath(options["output-file"])
    output_hdr_file = join([options["output-file"], "pfm"], ".")
    output_ldr_file = join([options["output-file"], options["ldr-extension"]], ".")

    cond1 = isfile(output_hdr_file)
    cond2 = options["with-tonemapping"] && isfile(output_ldr_file)
    if !options["force"] && (cond1 || cond2)
        img_str = cond1 ? cond2 ? output_hdr_file*"\" and \""*output_ldr_file : output_hdr_file : output_ldr_file
        printstyled("\nWarning!", bold=true, color=:yellow)
        print(" Image \"$(img_str)\" existing: overwrite? [y|n] ")
        if readline() != "y"
            printstyled("Aborting.\n", bold=true, color=:red)
            exit(1)
        end
    end
    mkpath(dirname(options["input-script"]))

    Raytracer.render_from_script(
        options["input-script"],
        output_file = output_hdr_file,
        time = options["time"],
        vars = options["vars"]
    )

    if options["with-tonemapping"]
        if options["luminosity"] == -1
            options["luminosity"] = nothing
        end
        Raytracer.tonemapping(
            output_hdr_file,
            output_ldr_file,
            α = options["alpha"],
            γ = options["gamma"],
            luminosity = options["luminosity"]
        )
    end
end

function renderanimation_loop(elem::Tuple{Int, Float32}, total_elem::Int, options::Dict{String, Any})
    index, t = elem
    # create file name
    filename = joinpath(options["output-dir"], "$(options["output-file"])_$(lpad(repr(index), trunc(Int, log10(total_elem))+1, '0'))")
    output_hdr_file = "$(filename).pfm"
    output_ldr_file = "$(filename).$(options["ldr-extension"])"
    # render
    Raytracer.render_from_script(
        options["input-script"],
        output_file = output_hdr_file,
        time = t,
        vars = options["vars"],
        use_threads = false,
        disable_output = true
    )
    # apply tonemapping
    if options["luminosity"] == -1
        options["luminosity"] = nothing
    end
    Raytracer.tonemapping(
        output_hdr_file,
        output_ldr_file,
        α = options["alpha"],
        γ = options["gamma"],
        luminosity = options["luminosity"],
        disable_output = true
    )
    # remove hdr image
    rm(output_hdr_file)
end

function renderanimation(options::Dict{String, Any})
    printstyled("Raytracer.jl animation rendering\n\n", bold=true)
    println("Number of threads: $(Threads.nthreads())\n")

    if Sys.which("ffmpeg") === nothing
        println()
        print_error("ffmpeg not found in PATH", "\n")
        exit(1)
    end

    options["input-script"] = normpath(options["input-script"])
    options["output-dir"] = normpath(options["output-dir"])
    options["output-file"] = normpath(options["output-file"])

    demodir = options["output-dir"]
    if !options["force"] && isdir(demodir)
        printstyled("Warning!", bold=true, color=:yellow)
        print(" Directory \"$(demodir)\" existing: overwrite content? [y|n] ")
        if readline() != "y"
            printstyled("Aborting.\n", bold=true, color=:red)
            exit(1)
        end
    end

    print("Creating directory \"$(options["output-dir"])\"...")
    rm(demodir, force=true, recursive=true)
    mkpath(demodir)
    println(" done!")

    println("\n-> RENDERING")
    t_start, t_end = Tuple(parse.(Float32, split(options["time-limits"], ":")))
    t_list = (if options["n-frames"] === nothing
        range(t_start, t_end, step=ifelse(t_start<t_end,1,-1)*options["delta-t"])
    else
        range(t_start, t_end, length=options["n-frames"]+1)
    end)[begin:end-1]
    l_t_list = length(t_list)
    println("Generating $(l_t_list) frames...")
    p = Progress(l_t_list, dt=2, color=:white)
    Threads.@threads for elem in collect(enumerate(t_list))
        renderanimation_loop(elem, l_t_list, options)
        next!(p)
    end

    println("\n-> CREATE ANIMATION")
    print("Generating animation...")
    padding = trunc(Int, log10(l_t_list)) + 1
    file_pattern = joinpath(options["output-dir"], "$(options["output-file"])_%0$(padding)d.$(options["ldr-extension"])")
    framerate = ["-r", options["fps"]]
    encoding_options = ["-c:v", "libx264", "-preset", "slow", "-vf", "format=yuv420p", "-movflags", "+faststart"]
    out_file = joinpath(options["output-dir"], "$(options["output-file"]).mp4")
    run(pipeline(
        `ffmpeg -y -i $(file_pattern) $(framerate) $(encoding_options) $(out_file)`,
        stdout=devnull,
        stderr=devnull
    ))
    println(" done!")
end

function tonemapping(options::Dict{String, Any})
    printstyled("Raytracer.jl tone mapping process\n", bold=true)
    # normalize path
    options["input-file"] = normpath(options["input-file"])
    options["output-file"] = normpath(options["output-file"])
    # check if output file exist
    if !options["force"] && isfile(options["output-file"])
        printstyled("\nWarning!", bold=true, color=:yellow)
        print(" Image \"$(options["output-file"])\" existing: overwrite? [y|n] ")
        if readline() != "y"
            printstyled("Aborting.\n", bold=true, color=:red)
            exit(1)
        end
    end
    mkpath(dirname(options["output-file"]))
    # apply tonemapping
    if options["luminosity"] == -1
        options["luminosity"] = nothing
    end
    Raytracer.tonemapping(
        options["input-file"],
        options["output-file"],
        α = options["alpha"],
        γ = options["gamma"],
        luminosity = options["luminosity"]
    )
end

function docs(options::Dict{String, Any})
    doc_type = options["dev"] ? "dev" : "stable"
    printstyled("Raytracer.jl $(doc_type) documentation\n", bold=true)
    printstyled("\nDocumentation\n", bold=true)
    println("  - https://paolo97gll.github.io/Raytracer.jl/$(doc_type)")
    printstyled("\nCLI quickstart\n", bold=true)
    println("  - https://paolo97gll.github.io/Raytracer.jl/$(doc_type)/quickstart/cli")
    printstyled("\nSceneLang quickstart\n", bold=true)
    println("  - https://paolo97gll.github.io/Raytracer.jl/$(doc_type)/quickstart/scenelang")
    println()
end

#######
# main

function main()
    printstyled(raw"""
               _                            _ _    _ _
 _ _ __ _ _  _| |_ _ _ __ _ __ ___ _ _   __| (_)  (_) |
| '_/ _` | || |  _| '_/ _` / _/ -_) '_| / _| | |_ | | |
|_| \__,_|\_, |\__|_| \__,_\__\___|_| __\__|_|_(_)/ |_|
          |__/                       |__|       |__/

""", bold=true, color=:green)

    parsed_args = parse_commandline()

    parsed_command = parsed_args["%COMMAND%"]
    parsed_args = parsed_args[parsed_command]

    if parsed_command == "render"
        parsed_subcommand = parsed_args["%COMMAND%"]
        parsed_command *= parsed_subcommand
        parsed_args = parsed_args[parsed_subcommand]
    end

    printstyled("raytracer_cli.jl : ", bold=true)

    (Symbol(parsed_command) |> eval)(parsed_args)
end

main()


text:
function fit_cumulative(cum,time)
    model(x,p) = @. (1 - exp(-p[1]*x)) * p[2]
    p0 = [0.5,1]
    fitting = curve_fit(model, cum, time, p0)
    # return DataFrame(Coef = fitting.param[1], Scale = fitting.param[2])
    # return (fitting.param...,)
    return fitting.param
end

ParetoModel(x,p) = @. (1 - ((p[1]/x)^p[2])) * p[3]
function fit_pareto(x,y)
    p0 = [0.5,0.5,0.5]
    fitting = curve_fit(ParetoModel, x, y, p0)
    return fitting.param
end

ExpoModel(x,p) = @. ((1 - exp(-p[1]*x)) * p[2]) + p[3]
function fit_ExpoModel(x,y)
    p0 = [0.5,0.5,0.5]
    fitting = curve_fit(ExpoModel, x, y, p0)
    return fitting.param
end

function MVTtangent(coef, P)
    # P is the point where the line pass through
    # cumulative to be tangential to
    c(x) = @. (1 - exp(-coef[1]*x))*coef[2]
    #derivative of the cum, equal to the slope
    d(x) = @. coef[1] * exp(-coef[1]*x) *coef[2]
    #solve using slope-intercept form of the line
    # y-y0 = m(x-x0)
    # x0 has to be a point on the curve so express it in this way
    # y-c(x0) = d(x0)(x-x0)
    # solve passing for the point P
    # yp -c(x0) = d(x0)(xp-x0)
    # yp -c(x0) - d(x0)(xp-x0) = 0
    # f(x) = P[2] - (1 - exp(-coef[1]*x))*coef[2] - (coef[1] * exp(-coef[1]*x))*coef[2]*(P[1]-x)
    f(x) = P[2] - c(x) - d(x)*(P[1]-x)
    x0 = Roots.find_zero(f,0) # search x0 in a reasonable range of trial time
    #find inercept using P and m = d(x0)
    # y = mx +q
    # P[2] = d(x0)*P[1] + q
    # q = -d(x0)*P[1] - P[2]
    q = -d(x0)*P[1] - P[2]
    # sol(x) = @. x * d(x0) + q
    # return DataFrame(OptLeave = x0, Slope = d(x0), Intercept = q)
    return (x0,d(x0),q)
end

GainModel(x,p) = @. (1 - exp(-p[1]*x)) *p[2]

# control = filter(r->r.Treatment == "Control", timedf)
# coef = fit_cumulative(control.Time, control.CumRew)
# control2 = filter(r->r.Treatment == "Control", timedf2)
# coef2 = fit_cumulative(control2.Time, control2.CumRew)
#
# xaxis = 0:0.5:35
# plot(xaxis,GainModel(xaxis,coef), label = "TrialTime", legend = :bottomright)
# plot!(xaxis,GainModel(xaxis,coef2), label = "ForagingTime", linecolor = :purple)


text:
##
cd(@__DIR__)
foreach(eachline("input.txt")) do line
    println(line)
end


text:
"""
Package functions for interacting with Travis.

$(EXPORTS)
"""
module Travis

using DocStringExtensions
using Base64
import LibGit2: GITHUB_REGEX
using ..DocumenterTools: package_devpath

export genkeys

"""
    genkeys(; user="\$USER", repo="\$REPO")

Generates the SSH keys necessary for the automatic deployment of documentation with
Documenter from Travis to GitHub Pages.

By default the links in the instructions need to be modified to correspond to actual URLs.
The optional `user` and `repo` keyword arguments can be specified so that the URLs in the
printed instructions could be copied directly. They should be the name of the GitHub user or
organization where the repository is hosted and the full name of the repository,
respectively.

This method of [`genkeys`](@ref) requires the following command lines programs to be
installed:

- `which`
- `ssh-keygen`

# Examples

```julia-repl
julia> using DocumenterTools

julia> Travis.genkeys()
[ Info: add the public key below to https://github.com/\$USER/\$REPO/settings/keys with read/write access:

ssh-rsa AAAAB3NzaC2yc2EAAAaDAQABAAABAQDrNsUZYBWJtXYUk21wxZbX3KxcH8EqzR3ZdTna0Wgk...jNmUiGEMKrr0aqQMZEL2BG7 username@hostname

[ Info: add a secure environment variable named 'DOCUMENTER_KEY' to https://travis-ci.org/\$USER/\$REPO/settings with value:

LS0tLS1CRUdJTiBSU0EgUFJJVkFURSBLRVktLS0tLQpNSUlFb3dJQkFBS0NBUUVBNnpiRkdXQVZpYlIy...QkVBRWFjY3BxaW9uNjFLaVdOcDU5T2YrUkdmCi0tLS0tRU5EIFJTQSBQUklWQVRFIEtFWS0tLS0tCg==

julia> Travis.genkeys(user="JuliaDocs", repo="DocumenterTools.jl")
[Info: add the public key below to https://github.com/JuliaDocs/DocumenterTools.jl/settings/keys with read/write access:

ssh-rsa AAAAB3NzaC2yc2EAAAaDAQABAAABAQDrNsUZYBWJtXYUk21wxZbX3KxcH8EqzR3ZdTna0Wgk...jNmUiGEMKrr0aqQMZEL2BG7 username@hostname

[ Info: add a secure environment variable named 'DOCUMENTER_KEY' to https://travis-ci.org/JuliaDocs/DocumenterTools.jl/settings with value:

LS0tLS1CRUdJTiBSU0EgUFJJVkFURSBLRVktLS0tLQpNSUlFb3dJQkFBS0NBUUVBNnpiRkdXQVZpYlIy...QkVBRWFjY3BxaW9uNjFLaVdOcDU5T2YrUkdmCi0tLS0tRU5EIFJTQSBQUklWQVRFIEtFWS0tLS0tCg==
```
"""
function genkeys(; user="\$USER", repo="\$REPO")
    # Error checking. Do the required programs exist?
    success(`which which`)      || error("'which' not found.")
    success(`which ssh-keygen`) || error("'ssh-keygen' not found.")

    directory = pwd()
    filename  = "documenter-private-key"

    isfile(filename) && error("temporary file '$(filename)' already exists in working directory")
    isfile("$(filename).pub") && error("temporary file '$(filename).pub' already exists in working directory")

    # Generate the ssh key pair.
    success(`ssh-keygen -N "" -f $filename`) || error("failed to generate a SSH key pair.")

    # Prompt user to add public key to github then remove the public key.
    let url = "https://github.com/$user/$repo/settings/keys"
        @info("add the public key below to $url with read/write access:")
        println("\n", read("$filename.pub", String))
        rm("$filename.pub")
    end

    # Base64 encode the private key and prompt user to add it to travis. The key is
    # *not* encoded for the sake of security, but instead to make it easier to
    # copy/paste it over to travis without having to worry about whitespace.
    let url = "https://travis-ci.org/$user/$repo/settings"
        @info("add a secure environment variable named 'DOCUMENTER_KEY' to $url with value:")
        println("\n", base64encode(read(filename, String)), "\n")
        rm(filename)
    end
end

"""
    genkeys(package::Module; remote="origin")

Like the other method, this generates the SSH keys necessary for the automatic deployment of
documentation with Documenter from Travis to GitHub Pages, but attempts to guess the package
URLs from the Git remote.

`package` needs to be the top level module of the package. The `remote` keyword argument can
be used to specify which Git remote is used for guessing the repository's GitHub URL.

This method requires the following command lines programs to be installed:

- `which`
- `git`
- `ssh-keygen`

!!! note
    The package must be in development mode. Make sure you run
    `pkg> develop pkg` from the Pkg REPL, or `Pkg.develop(\"pkg\")`
    before generating the SSH keys.

# Examples

```julia-repl
julia> using DocumenterTools

julia> Travis.genkeys(DocumenterTools)
[Info: add the public key below to https://github.com/JuliaDocs/DocumenterTools.jl/settings/keys with read/write access:

ssh-rsa AAAAB3NzaC2yc2EAAAaDAQABAAABAQDrNsUZYBWJtXYUk21wxZbX3KxcH8EqzR3ZdTna0Wgk...jNmUiGEMKrr0aqQMZEL2BG7 username@hostname

[ Info: add a secure environment variable named 'DOCUMENTER_KEY' to https://travis-ci.org/JuliaDocs/DocumenterTools.jl/settings with value:

LS0tLS1CRUdJTiBSU0EgUFJJVkFURSBLRVktLS0tLQpNSUlFb3dJQkFBS0NBUUVBNnpiRkdXQVZpYlIy...QkVBRWFjY3BxaW9uNjFLaVdOcDU5T2YrUkdmCi0tLS0tRU5EIFJTQSBQUklWQVRFIEtFWS0tLS0tCg==
```
"""
function genkeys(package::Module; remote="origin")
    # Error checking. Do the required programs exist?
    success(`which which`)      || error("'which' not found.")
    success(`which git`)        || error("'git' not found.")

    path = package_devpath(package)

    # Are we in a git repo?
    user, repo = cd(path) do
        success(`git status`) || error("Failed to run `git status` in $(path). 'Travis.genkeys' only works with Git repositories.")

        let r = readchomp(`git config --get remote.$remote.url`)
            m = match(GITHUB_REGEX, r)
            m === nothing && error("no remote repo named '$remote' found.")
            m[2], m[3]
        end
    end

    # Generate the ssh key pair.
    genkeys(; user=user, repo=repo)
end

end # module


text:
function add_points_in_z_and_interpolate(potential::Array{T, 3}, grid::CylindricalGrid{T}, idcs::Vector{Int})::Tuple{Array{T, 3}, CylindricalGrid{T}} where {T}
    axz::Vector{T} = collect(grid.axes[3])
    n = length(idcs)
    if n == 0 return potential, grid end
    if idcs[end] >= length(axz)
        error("Can not append point in z at the end. Only in between.")
    end

    ranges = [ 1:idcs[1] ]
    for i in (2:n)
        push!(ranges, idcs[i-1]+1:idcs[i])
    end
    push!(ranges, idcs[end]+1:length(axz))

    axz = cat(axz, zeros(T, n), dims=1)

    for (i, range) in enumerate(reverse(ranges))
        axz[range .+ (n+1-i)] = axz[range]
    end
    for (i, range) in enumerate(reverse(ranges[2:end]))
        axz[first(range .+ (n+1-i))-1] = 0.5*(axz[first(range .+ (n+1-i))]+axz[first(range .+ (n+1-i))-2])
    end

    # 1: extend arrays
     potential = cat( potential, zeros(T, size( potential, 1), size( potential, 2), n), dims=3)
    # 2: copy old values to new locations
    for (i, range) in enumerate(reverse(ranges))
        potential[:, :, range .+ (n+1-i)] = potential[:, :, range]
    end
     # 3: interpolate new positions
    for (i, range) in enumerate(reverse(ranges[2:end]))
        potential[:, :, first(range .+ (n+1-i))-1] = 0.5 * ( potential[:, :, first(range .+ (n+1-i))] + potential[:, :, first(range .+ (n+1-i))-2])
    end

    BL::Symbol = typeof(grid.axes[3]).parameters[2]
    BR::Symbol = typeof(grid.axes[3]).parameters[3]
    daz::DiscreteAxis{T, BL, BR} = DiscreteAxis{T, BL, BR}( grid.axes[3].interval, axz )
    grid::CylindricalGrid{T} = CylindricalGrid{T}( (grid.axes[1], grid.axes[2], daz) )

    return potential, grid
end

function add_points_in_r_and_interpolate(potential::Array{T, 3}, grid::CylindricalGrid{T}, idcs::Vector{Int})::Tuple{Array{T, 3}, CylindricalGrid{T}} where {T}
    axr::Vector{T} = collect(grid.axes[1])
    n = length(idcs)
    # idcs = collect(idcs)
    if n == 0 return potential, grid end
    # if !issorted(idcs) idcs = sort(idcs) end
    if idcs[end] >= length(axr)
        error("Can not append point in r at the end. Only in between.")
    end

    ranges = [ 1:idcs[1] ]
    for i in (2:n)
        push!(ranges, idcs[i-1]+1:idcs[i])
    end
    push!(ranges, idcs[end]+1:length(axr))

    axr = cat(axr, zeros(T, n), dims=1)

    for (i, range) in enumerate(reverse(ranges))
        axr[range .+ (n+1-i)] = axr[range]
    end
    for (i, range) in enumerate(reverse(ranges[2:end]))
        axr[first(range .+ (n+1-i))-1] = 0.5*(axr[first(range .+ (n+1-i))]+axr[first(range .+ (n+1-i))-2])
    end

    # 1: extend arrays
    potential = cat( potential, zeros(T, n, size( potential, 2), size( potential, 3)), dims=1)
    # 2: copy old values to new locations
    for (i, range) in enumerate(reverse(ranges))
        potential[range .+ (n+1-i), :, :] = potential[range, :, :]
    end
     # 3: interpolate new positions
    for (i, range) in enumerate(reverse(ranges[2:end]))
        potential[first(range .+ (n+1-i))-1, :, :] = 0.5 * ( potential[first(range .+ (n+1-i)), :, :] + potential[first(range .+ (n+1-i))-2, :, :])
    end

    BL::Symbol = typeof(grid.axes[1]).parameters[2]
    BR::Symbol = typeof(grid.axes[1]).parameters[3]
    da::DiscreteAxis{T, BL, BR} = DiscreteAxis{T, BL, BR}( grid.axes[1].interval, axr )
    grid::CylindricalGrid{T} = CylindricalGrid{T}( (da, grid.axes[2], grid.axes[3]) )

    return potential, grid
end

function add_points_in_φ_and_interpolate(potential::Array{T, 3}, grid::CylindricalGrid{T}, idcs::Vector{Int})::Tuple{Array{T, 3}, CylindricalGrid{T}} where {T}
    axφ::Vector{T} = collect(grid.axes[2])
    n = length(idcs)
    if n == 0 return potential, grid end
    cyclic::T = grid.axes[2].interval.right

    if idcs[end] > length(axφ)
        error("Can not append array in φ at index $(idcs[end]). (Cylcic: $(cyclic)°)")
    end

    append_after_last_index = false
    if idcs[end] == length(axφ)
        append_after_last_index = true
        idcs = idcs[1:end-1]
        n -= 1
    end

    ranges = []
    if length(idcs) > 0
        ranges = [ 1:idcs[1] ]
        for i in (2:n)
            push!(ranges, idcs[i-1]+1:idcs[i])
        end
        if idcs[end] < length(axφ)
            push!(ranges, idcs[end]+1:length(axφ))
        end
    end

    if append_after_last_index
        axφ = cat( axφ, zeros(T, n+1), dims=1)
    else
        axφ = cat( axφ, zeros(T, n), dims=1)
    end

    for (i, range) in enumerate(reverse(ranges))
        axφ[range .+ (n+1-i)] = axφ[range]
    end
    for (i, range) in enumerate(reverse(ranges[2:end]))
        axφ[first(range .+ (n+1-i))-1] = 0.5*(axφ[first(range .+ (n+1-i))]+axφ[first(range .+ (n+1-i))-2])
    end

    if append_after_last_index
        axφ[end] = 0.5 * (cyclic + axφ[end-1])
    end

    # 1: extend arrays
    if append_after_last_index
        potential = cat( potential, zeros(T, size( potential, 1), n + 1, size( potential, 3)), dims=2 )
    else
        potential = cat( potential, zeros(T, size( potential, 1), n, size( potential, 3)), dims=2)
    end
    # 2: copy old values to new locations
    for (i, range) in enumerate(reverse(ranges))
        potential[:, range .+ (n+1-i), :] = potential[:, range, :]
    end
     # 3: interpolate new positions
    for (i, range) in enumerate(reverse(ranges[2:end]))
        potential[: ,first(range .+ (n+1-i))-1, :] = 0.5 * ( potential[:, first(range .+ (n+1-i)), :] + potential[:, first(range .+ (n+1-i))-2, :])
    end
    if append_after_last_index
        potential[: ,end, :] = 0.5 * ( potential[:, end-1, :] + potential[:, 1, :])
    end

    BL::Symbol = typeof(grid.axes[2]).parameters[2]
    BR::Symbol = typeof(grid.axes[2]).parameters[3]
    da::DiscreteAxis{T, BL, BR} = DiscreteAxis{T, BL, BR}( grid.axes[2].interval, axφ )
    grid::CylindricalGrid{T} = CylindricalGrid{T}( (grid.axes[1], da, grid.axes[3]) )

    return potential, grid
end

function add_points_and_interpolate(potential::Array{T, 3}, grid::CylindricalGrid{T}, idcs_r::Vector{Int}, idcs_φ::Vector{Int}, idcs_z::Vector{Int})::Tuple{Array{T, 3}, CylindricalGrid{T}} where {T}
    potential::Array{T, 3}, grid::CylindricalGrid{T} = add_points_in_r_and_interpolate(potential, grid, idcs_r)
    potential, grid = add_points_in_φ_and_interpolate(potential, grid, idcs_φ)
    if isodd(length(grid.φ)) && length(length(grid.φ)) != 1
        potential, grid = add_points_in_φ_and_interpolate(potential, grid, [1])
    end
    potential, grid = add_points_in_z_and_interpolate(potential, grid, idcs_z)
    potential, grid
end

function add_points_and_interpolate(potential::Array{T, 3}, grid::CartesianGrid3D{T}, idcs_x::Vector{Int}, idcs_y::Vector{Int}, idcs_z::Vector{Int})::Tuple{Array{T, 3}, CartesianGrid3D{T}} where {T}
    potential::Array{T, 3}, grid::CartesianGrid3D{T} = add_points_in_x_and_interpolate(potential, grid, idcs_x)
    potential, grid = add_points_in_y_and_interpolate(potential, grid, idcs_y)
    potential, grid = add_points_in_z_and_interpolate(potential, grid, idcs_z)
    potential, grid
end

function add_points_in_z_and_interpolate(potential::Array{T, 3}, grid::CartesianGrid3D{T}, idcs::Vector{Int})::Tuple{Array{T, 3}, CartesianGrid3D{T}} where {T}
    ax::Vector{T} = collect(grid.z)
    n = length(idcs)
    if n == 0 return potential, grid end
    if idcs[end] >= length(ax)
        error("Can not append point in z at the end. Only in between.")
    end

    ranges = [ 1:idcs[1] ]
    for i in (2:n)
        push!(ranges, idcs[i-1]+1:idcs[i])
    end
    push!(ranges, idcs[end]+1:length(ax))

    ax = cat(ax, zeros(T, n), dims = 1)

    for (i, range) in enumerate(reverse(ranges))
        ax[range .+ (n + 1 - i)] = ax[range]
    end
    for (i, range) in enumerate(reverse(ranges[2:end]))
        ax[first(range .+ (n + 1 - i)) - 1] = 0.5 * (ax[first(range .+ (n + 1 - i))] + ax[first(range .+ (n + 1 - i)) - 2])
    end

    # 1: extend arrays
    potential = cat( potential, zeros(T, size( potential, 1), size( potential, 2), n), dims = 3)
    # 2: copy old values to new locations
    for (i, range) in enumerate(reverse(ranges))
        potential[:, :, range .+ (n+1-i)] = potential[:, :, range]
    end
     # 3: interpolate new positions
    for (i, range) in enumerate(reverse(ranges[2:end]))
        potential[:, :, first(range .+ (n+1-i))-1] = 0.5 * ( potential[:, :, first(range .+ (n+1-i))] + potential[:, :, first(range .+ (n+1-i))-2])
    end

    BL::Symbol = typeof(grid.axes[3]).parameters[2]
    BR::Symbol = typeof(grid.axes[3]).parameters[3]
    da::DiscreteAxis{T, BL, BR} = DiscreteAxis{T, BL, BR}( grid.z.interval, ax )
    grid::CartesianGrid3D{T} = CartesianGrid3D{T}( (grid[:x], grid[:y], da) )

    return potential, grid
end

function add_points_in_x_and_interpolate(potential::Array{T, 3}, grid::CartesianGrid3D{T}, idcs::Vector{Int})::Tuple{Array{T, 3}, CartesianGrid3D{T}} where {T}
    ax::Vector{T} = collect(grid[:x])
    n = length(idcs)
    if n == 0 return potential, grid end
    if idcs[end] >= length(ax)
        error("Can not append point in z at the end. Only in between.")
    end

    ranges = [ 1:idcs[1] ]
    for i in (2:n)
        push!(ranges, idcs[i-1]+1:idcs[i])
    end
    push!(ranges, idcs[end]+1:length(ax))

    ax = cat(ax, zeros(T, n), dims = 1)

    for (i, range) in enumerate(reverse(ranges))
        ax[range .+ (n + 1 - i)] = ax[range]
    end
    for (i, range) in enumerate(reverse(ranges[2:end]))
        ax[first(range .+ (n + 1 - i)) - 1] = 0.5 * (ax[first(range .+ (n + 1 - i))] + ax[first(range .+ (n + 1 - i)) - 2])
    end

    # 1: extend arrays
    potential = cat( potential, zeros(T, n, size( potential, 2), size( potential, 3)), dims = 1)
    # 2: copy old values to new locations
    for (i, range) in enumerate(reverse(ranges))
        potential[range .+ (n+1-i), :, :] = potential[range, :, :]
    end
     # 3: interpolate new positions
    for (i, range) in enumerate(reverse(ranges[2:end]))
        potential[first(range .+ (n+1-i))-1, :, :] = 0.5 * ( potential[first(range .+ (n+1-i)), :, :] + potential[first(range .+ (n+1-i))-2, :, :])
    end

    BL::Symbol = typeof(grid[:x]).parameters[2]
    BR::Symbol = typeof(grid[:x]).parameters[3]
    da::DiscreteAxis{T, BL, BR} = DiscreteAxis{T, BL, BR}( grid[:x].interval, ax )
    grid::CartesianGrid3D{T} = CartesianGrid3D{T}( (da, grid[:y], grid.z) )

    return potential, grid
end

function add_points_in_y_and_interpolate(potential::Array{T, 3}, grid::CartesianGrid3D{T}, idcs::Vector{Int})::Tuple{Array{T, 3}, CartesianGrid3D{T}} where {T}
    ax::Vector{T} = collect(grid[:y])
    n = length(idcs)
    if n == 0 return potential, grid end
    if idcs[end] >= length(ax)
        error("Can not append point in z at the end. Only in between.")
    end

    ranges = [ 1:idcs[1] ]
    for i in (2:n)
        push!(ranges, idcs[i-1]+1:idcs[i])
    end
    push!(ranges, idcs[end]+1:length(ax))

    ax = cat(ax, zeros(T, n), dims = 1)

    for (i, range) in enumerate(reverse(ranges))
        ax[range .+ (n + 1 - i)] = ax[range]
    end
    for (i, range) in enumerate(reverse(ranges[2:end]))
        ax[first(range .+ (n + 1 - i)) - 1] = 0.5 * (ax[first(range .+ (n + 1 - i))] + ax[first(range .+ (n + 1 - i)) - 2])
    end

    # 1: extend arrays
    potential = cat( potential, zeros(T, size( potential, 1), n, size( potential, 3)), dims = 2)
    # 2: copy old values to new locations
    for (i, range) in enumerate(reverse(ranges))
        potential[:, range .+ (n+1-i), :] = potential[:, range, :]
    end
     # 3: interpolate new positions
    for (i, range) in enumerate(reverse(ranges[2:end]))
        potential[:, first(range .+ (n+1-i))-1, :] = 0.5 * ( potential[:, first(range .+ (n+1-i)), :] + potential[:, first(range .+ (n+1-i))-2, :])
    end

    BL::Symbol = typeof(grid[:y]).parameters[2]
    BR::Symbol = typeof(grid[:y]).parameters[3]
    da::DiscreteAxis{T, BL, BR} = DiscreteAxis{T, BL, BR}( grid[:y].interval, ax )
    grid::CartesianGrid3D{T} = CartesianGrid3D{T}( (grid[:x], da,  grid.z) )

    return potential, grid
end


text:
using Documenter, CarlemanLinearization
using DynamicPolynomials

DocMeta.setdocmeta!(CarlemanLinearization, :DocTestSetup,
                   :(using CarlemanLinearization); recursive=true)

makedocs(
    sitename = "CarlemanLinearization.jl",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",
                             assets = ["assets/aligned.css"]),
    strict = true,
    pages = [
        "Home" => "index.md",
        "API Reference" => Any["Linearization" => "lib/linearization.md",
                               "Error bounds" => "lib/errors.md"],
        "References" => "references.md",
        "About" => "about.md"
    ]
)

deploydocs(
    repo = "github.com/JuliaReach/CarlemanLinearization.jl.git",
    push_preview = true
)


text:
# https://github.com/grasingerm/lbxflow/blob/master/inc/lattice.jl
# http://mafija.fmf.uni-lj.si/seminar/files/2012_2013/Lattice_Boltzmann_method.pdf

# const original_order = [1, 2 ,3, 4]

# From: E3,1,5
struct D2Q5{
    T <: Real,
    Abscissaes <: AbstractMatrix{<:Integer},
    Weights <: AbstractVector{<:Real},
} <: Quadrature
    abscissae::Abscissaes
    weights::Weights
    speed_of_sound_squared::T
end

D2Q5() = D2Q5(
    [
        0 1 0 -1 0
        0 0 1 0 -1
    ],
    [
        4 / 6,
        1 / 12,
        1 / 12,
        1 / 12,
        1 / 12,
        # 1/3, 1/6, 1/6, 1/6, 1/6
    ],
    6.0,
)

order(q::D2Q5) = 3
function opposite(q::D2Q5, idx::Int64)
    if idx == 1
        return 1
    end
    if idx <= 3
        return idx + 2
    end
    return idx - 2
end

Base.show(io::IO, q::D2Q5) = show(io, "D2Q5")
Base.string(q::D2Q5) = "D2Q5"

# The D2Q5 lattice is unable to include temperature
pressure(q::D2Q5, f::VT, ρ::T, u::VT) where {T <: Real, VT <: AbstractVector{T}} = one(T)


text:
using Multibreak

function part1(input)
    scanners = parse_input(input)
    signatures = compute_signature.(scanners)
    overlaps = [length(intersect(s1, s2))
                for s1 in signatures, s2 in signatures]
    align_coordinates!(scanners, signatures, overlaps)
    return length(reduce(union, scanners))
end

function part2(input)
    scanners = parse_input(input)
    signatures = compute_signature.(scanners)
    overlaps = [length(intersect(s1, s2))
                for s1 in signatures, s2 in signatures]
    scanner_positions = align_coordinates!(scanners, signatures, overlaps)
    return maximum(sum(abs.(s1 .- s2)) for s1 in scanner_positions, s2 in scanner_positions)
end

function parse_input(input)
    scanners = Vector{Vector{Int}}[]
    for line in eachline(input)
        if occursin("scanner", line)
            push!(scanners, Vector{Int}[])
        elseif !isempty(line)
            push!(last(scanners), parse.(Int, split(line, ",")))
        end
    end
    return scanners
end

function compute_signature(scanner)
    return [sort(abs.(scanner[i] - scanner[j]))
            for i=1:length(scanner), j=1:length(scanner) if i > j]
end

function compute_signature2(scanner)
    return [sort(abs.(scanner[i] - scanner[j]))
            for i=1:length(scanner), j=1:length(scanner)]
end

function align_coordinates!(scanners, signatures, overlaps)
    aligned = [1]
    scanner_positions = [[0, 0, 0]]
    @multibreak while length(aligned) < length(scanners)
        for i in aligned
            for j in setdiff(1:length(scanners), aligned)
                if overlaps[i, j] >= 66
                    for signature in intersect(signatures[i], signatures[j])
                        if length(unique(signature)) == 3 && prod(signature) != 0
                            align_coordinates2!(scanners[i], scanners[j],
                                                signature, scanner_positions)
                            push!(aligned, j)
                            break; break; break
                        end
                    end
                end
            end
        end
    end
    return scanner_positions
end

function align_coordinates2!(scanner1, scanner2, signature, scanner_positions)
    s1 = compute_signature2(scanner1)
    s2 = compute_signature2(scanner2)
    a, b = Tuple(findfirst(s1 .== Ref(signature)))
    c, d = Tuple(findfirst(s2 .== Ref(signature)))
    d1 = scanner1[a] - scanner1[b]
    d2 = scanner2[c] - scanner2[d]
    for p in ([1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1])
        if abs.(d1) == abs.(d2[p])
            scanner2 .= (x -> x[p]).(scanner2)
            scanner2 .= (x -> x .* d2[p] ./ d1).(scanner2)
            if prod(d2) / prod(d1)  * (p[3] - p[1]) * (p[3] - p[2]) * (p[2] - p[1]) < 0
                scanner2 .= (x -> -x).(scanner2)
            end
            dd2 = scanner2[c] - scanner2[d]
            if d1 == dd2
                Δ = scanner1[a] .- scanner2[c]
            else
                Δ = scanner1[b] .- scanner2[c]
            end
            scanner2 .= (x -> x .+ Δ).(scanner2)
            push!(scanner_positions, Δ)
            return
        end
    end
end


text:
function nearby(; kwargs...)
    resultjson = requestjson("place/nearbysearch"; kwargs...)
    if resultjson["status"] == "OK"
        return resultjson["results"]
    else
        return []
    end
end

"""
Returns a list of places around a location within the given radius.

## Description

A Nearby Search lets you search for places within a specified area. You can
refine your search request by supplying keywords or specifying the type of
place you are searching for.

https://developers.google.com/places/web-service/search#PlaceSearchRequests

## Parameters

* `key`: your google maps API key
* `lon`: the longitude of the location
* `lat`: the latitude of the location
* `radius`: the distance (in meters) within which to return place results.

## Keyword arguments

* `keyword`: A term to be matched against all content that Google has indexed
    for this place, including but not limited to name, type, and address, as
    well as customer reviews and other third-party content.
* `opennow`: return only places that are open for business at the time the
    query is sent. Places that do not specify opening hours in the Google Places
    database will not be returned if you include this parameter in your query.
* `rankby`: Specifies the order in which results are listed. Note that rankby
    must not be included if radius (described under Required parameters above)
    is specified. Possible values are:

    * `prominence` (default). This option sorts results based on their
        importance. Ranking will favor prominent places within the specified
        area. Prominence can be affected by a place's ranking in Google's index,
        global popularity, and other factors.
    * `distance`. This option biases search results in ascending order by their
        distance from the specified location. When distance is specified, one or
        more of keyword, name, or type is required.
* `type`: Restricts the results to places matching the specified type.

"""
function nearby(key::String, lon::Real, lat::Real, radius::Real; kwargs...)
    radius > 50_000 && error("The maximum allowed radius is 50,000 meters.")
    nearby(key = key, location = string(lat,",",lon), radius=radius; kwargs...)
end

nearby(key::String, lonlat::Tuple{Float64,Float64}, radius::Real; kwargs...) =
    nearby(key, lonlat[1], lonlat[2], radius; kwargs...)


text:
##

using ACE
using Printf, Test, LinearAlgebra, StaticArrays
using ACE: evaluate, evaluate_d, Rn1pBasis, Ylm1pBasis,
      PositionState, Product1pBasis, State, ACEConfig,
      SymmetricBasis, get_spec
using Random: shuffle
using ACEbase.Testing: dirfdtest, fdtest, print_tf, test_fio
using ACE.OrthPolys: transformed_jacobi

##

maxdeg = 10
trans = ACE.Transforms.IdTransform()
P = transformed_jacobi(maxdeg, trans, 1.0, 0.0; pin = 0, pcut = 0)
Bsel = ACE.SimpleSparseBasis(3, maxdeg)
B1p = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg)
symB = SymmetricBasis(ACE.Invariant(), B1p, Bsel)

##
@info("Attempt to sparsify a basis")
symB2 = deepcopy(symB)
spec = get_spec(symB2)
keep = findall(bb -> all(b.n <= 4 && b.l <= 3 for b in bb), spec)
ACE.sparsify!(symB2, keep = keep)

@show length(symB), length(symB2)
@show length(symB.pibasis), length(symB2.pibasis)
@show length(symB.pibasis.basis1p), length(symB2.pibasis.basis1p)

@info("some basic consistency checks")
print_tf(@test length(symB2) == length(keep))
print_tf(@test length(symB.pibasis) > length(symB2.pibasis))
print_tf(@test length(symB.pibasis.basis1p) > length(symB2.pibasis.basis1p))
println()

##

@info("check evaluation is consistent")

for ntest = 1:30
   local cfg , BB1
   cfg = [ State(rr = ACE.rand_vec3(B1p["Rn"])) for _=1:10 ]
   BB1 = evaluate(symB, cfg)
   BB2 = evaluate(symB2, cfg)
   print_tf(@test BB1[keep] ≈ BB2)
end
println()


text:
function _make_kkt_solver!(ws::COSMO.Workspace)

	ws.kkt_solver = ws.settings.kkt_solver(SparseMatrixCSC(ws.p.P),
						SparseMatrixCSC(ws.p.A),
						ws.settings.sigma,
						ws.ρvec)
end

function setup!(ws::COSMO.Workspace)

  	allocate_set_memory!(ws)
	# scale problem data
	if ws.settings.scaling != 0
		if ws.settings.verbose_timing
			ws.times.scaling_time = @elapsed scale_ruiz!(ws)
		else
			scale_ruiz!(ws)
		end
	else
		ws.times.scaling_time = 0.
	end

	# construct a set index --> row range map
	ws.row_ranges = get_set_indices(ws.p.C.sets)
	classify_constraints!(ws)

	set_rho_vec!(ws)

	# create a KKT solver object populated with our data
	if(ws.flags.FACTOR_LHS)
		if ws.settings.verbose_timing
			ws.times.init_factor_time = @elapsed _make_kkt_solver!(ws)
		else
			_make_kkt_solver!(ws)
		end

	end
end

function allocate_set_memory!(ws::COSMO.Workspace)
  for (k, C) in enumerate(ws.p.C.sets)
    allocate_memory!(C)
  end
end

"For inequality and box constraints, check if any constraints are obviously loose or tight."
function classify_constraints!(ws::COSMO.Workspace)
	for (k, C) in enumerate(ws.p.C.sets)
		if C isa Nonnegatives
			b = view(ws.p.b, ws.row_ranges[k])
			classify_constraints!(C.constr_type, b, ws.settings.COSMO_INFTY, ws.settings.MIN_SCALING)
		elseif C isa Box
			classify_box_constraints!(C, ws.settings.COSMO_INFTY, ws.settings.MIN_SCALING, ws.settings.RHO_TOL)
		end
	end
	return nothing
end


text:
module TaglessTrees

struct TraverseWith{F,X}
    call :: F
    term :: X
end

(tw::TraverseWith)(tree) = tree(tw)

export call

function call(f, args...)
    function(tw::TraverseWith)
        tw.call(f, (arg(tw) for arg in args)...)
    end
end

export term

function term(x)
    function(tw::TraverseWith)
        tw.term(x)
    end
end

export evaluate

evaluate = let
    call = function(f, args...)
        return f(args...)
    end

    term = function(x)
        return x
    end

    TraverseWith(call, term)
end

export trace

trace = let
    call = function(f, args...)
        println("call(",nameof(f),",...)")
        return f(args...)
    end

    term = function(x)
        println("term(",x,")")
        return x
    end

    TraverseWith(call, term)
end

end # module


text:
MakieLayout.@Layoutable LTeX

function MakieLayout.default_attributes(::Type{LTeX}, scene)
    Attributes(
        tex = raw"\int \mathbf E \cdot d\mathbf a = \frac{Q_{encl}}{4\pi\epsilon_0}",
        visible = true,
        padding = (0f0, 0f0, 0f0, 0f0),
        height = Auto(),
        width = Auto(),
        alignmode = Inside(),
        valign = :center,
        halign = :center,
        tellwidth = true,
        tellheight = true,
        dpi = 72.0,
        textsize = 12,
        rotation = 0
    )
end

function MakieLayout.layoutable(::Type{LTeX}, fig_or_scene, tex; kwargs...)
    MakieLayout.layoutable(LTeX, fig_or_scene; tex = tex, kwargs...)
end

function MakieLayout.layoutable(::Type{LTeX}, fig_or_scene; bbox = nothing, kwargs...)
    # this is more or less copied from MakieLayout.Label
    topscene = MakieLayout.get_topscene(fig_or_scene)
    default_attrs = MakieLayout.default_attributes(LTeX, topscene).attributes
    theme_attrs = MakieLayout.subtheme(topscene, :LTeX)
    attrs = MakieLayout.merge!(MakieLayout.merge!(Attributes(kwargs), theme_attrs), Attributes(default_attrs))

    # @extract attrs (tex, textsize, font, color, visible, halign, valign,
    #     rotation, padding)
    @extract attrs (tex, dpi, textsize, visible, padding, halign, valign, rotation)

    layoutobservables = LayoutObservables{LTeX}(
        attrs.width, attrs.height, attrs.tellwidth, attrs.tellheight,
        halign, valign, attrs.alignmode; suggestedbbox = bbox
    )

    # This is Point3f0 in Label
    textpos = Node(Point3f0(0))

    cached_tex = @lift CachedTeX($tex, $dpi)

    # this is just a hack until boundingboxes in abstractplotting are perfect
    alignnode = lift(halign, rotation) do h, rot
        # left align the text if it's not rotated and left aligned
        if rot == 0 && (h == :left || h == 0.0)
            (:left, :center)
        else
            (:center, :center)
        end
    end

    t = teximg!(
        topscene, cached_tex; position = textpos, visible = visible, raw = true,
        textsize = textsize, dpi = dpi, align = alignnode, rotation = rotation
    )

    textbb = Ref(BBox(0, 1, 0, 1))

    # Label
    onany(cached_tex, textsize, rotation, padding) do _, textsize, rotation, padding
        # rotation is applied via a model matrix which isn't used in the bbox calculation
        # so we need to deal with it here
        bb = boundingbox(t)
        R = Makie.Mat3f0(
             cos(rotation), sin(rotation), 0,
            -sin(rotation), cos(rotation), 0,
            0, 0, 1
        )
        points = map(p -> R * p, unique(Makie.coordinates(bb)))
        new_bb = Makie.xyz_boundingbox(identity, points)
        textbb[] = FRect2D(new_bb)
        # textbb[] = FRect2D(boundingbox(t))
        autowidth  = Makie.width(textbb[]) + padding[1] + padding[2]
        autoheight = Makie.height(textbb[]) + padding[3] + padding[4]
        layoutobservables.autosize[] = (autowidth, autoheight)
    end

    onany(layoutobservables.computedbbox, padding) do bbox, padding
        tw = Makie.width(textbb[])
        th = Makie.height(textbb[])

        box = bbox.origin[1]
        boy = bbox.origin[2]

        # this is also part of the hack to improve left alignment until
        # boundingboxes are perfect
        tx = if rotation[] == 0 && (halign[] == :left || halign[] == 0.0)
            box + padding[1]
        else
            box + padding[1] + 0.5 * tw
        end
        ty = boy + padding[3] + 0.5 * th

        textpos[] = Point3f0(tx, ty, 0)
    end

    # trigger first update, otherwise bounds are wrong somehow
    padding[] = padding[]
    # trigger bbox
    layoutobservables.suggestedbbox[] = layoutobservables.suggestedbbox[]

    lt = LTeX(fig_or_scene, layoutobservables, attrs, Dict(:tex => t))

    lt
end


text:
#!/usr/bin/env julia
#
# Jerzy Wroczyński (nr. indeksu 250075)
# Obliczenia Naukowe, Laboratorium
# Lista 2, Zadanie 6.
#

# to avoid writing `global` wrap the entire program inside of a function
function program()

    # default case
    k = "1"

    # the first argument correlates to the considered case in the exercise
    if length(ARGS) != 1
        println("usage: ./ex-6.jl <number between 1-7>")
        exit(1)
    end
    k = ARGS[1]

    # if k == "1"
    c = -2
    x₀= 1

    if k == "2"
        c = -2
        x₀= 2
    elseif k == "3"
        c = -2
        x₀= 1.99999999999999
    elseif k == "4"
        c = -1
        x₀= 1
    elseif k == "5"
        c = -1
        x₀= -1
    elseif k == "6"
        c = -1
        x₀= 0.75
    elseif k == "7"
        c = -1
        x₀= 0.25
    end

    # iterate and print out the results
    x = x₀
    for n in 1:40

        x = x^2 + c

        print(
            "(", n, ",", x, ") "
        )

    end
    # put a newline at the end of the program
    println()

end

program()

max_stars_repo_name: jerry-sky/academic-notebook

max_issues_repo_name: jerry-sky/academic-notebook

max_forks_repo_name: jerry-sky/academic-notebook


text:
# make your scripts automatically re-activate your project
cd(@__DIR__)
using Pkg; Pkg.activate("."); Pkg.instantiate()

using Flux, DiffEqFlux, DiffEqSensitivity
using DifferentialEquations
using DiffEqSensitivity
using DiffEqNoiseProcess
using Plots
using Printf, BSON
using QuantumOptics
using Statistics
using Zygote
using StaticArrays
using LinearAlgebra

#################################################
#read parameters from external file
include("parameters_qb.jl")
using Main.parameters

struct Parameters
	dim::Int64
	w::Float32
	force_mag::Float32
	n_steps::Int64
	dt::Float32
	n_substeps::Int64
	gamma::Float32
end

MyParameters=Parameters(
	parameters.parameters_para["N"],
	parameters.parameters_para["w"],
	parameters.parameters_para["force_mag"],
	parameters.parameters_para["max_episode_steps"],
	parameters.parameters_para["dt"],
	parameters.parameters_para["n_substeps"],
	parameters.parameters_para["gamma"])

n_par = 64
# loss hyperparameters
C1 = 1.2f0 # evolution state fid
C2 = 0.0005f0 # action amplitudes
C3 = 0.8f0*MyParameters.n_steps/50 #enhance last 50steps

########################
using Random
Random.seed!(2)

################################
# Hamiltonian and evolution
basis = SpinBasis(1//2) #max n=dim-1
σ_z = sigmaz(basis)
σ_x = sigmax(basis)
σ_p = sigmap(basis)
σ_m = sigmam(basis)
σ_pm = σ_p*σ_m

H_0pre = MyParameters.w/2.0f0*σ_z #drift Hamiltonian
H_1pre = σ_x #control field

#to the dense regime, all real, otherwise Flux.train! complains
H_0 = real(Array{Float32}(H_0pre.data))
H_1 = real(Array{Float32}(H_1pre.data))

Arrayσ_p = real(Array{Float32}(σ_p.data))
Arrayσ_m = real(Array{Float32}(σ_m.data))
Arrayσ_pm = real(Array{Float32}(σ_pm.data))

###################
# DRIFT TERM FOR SDE
function qb_dynamics_dt!(du,u, α, t)  #du/u have a standard dimension
	ψRe = u[1:MyParameters.dim]   #(dim,)
	ψIm = u[MyParameters.dim+1:end]

	HRe = H_0 .+ α[1]*H_1 #size (dim,dim) ,alpha is just a single number

	Reρ = ψRe.*transpose(ψRe) + ψIm.*transpose(ψIm)
	ex_x = real(sum(diag((Arrayσ_p+Arrayσ_m)*Reρ)))

	HIm = -Arrayσ_pm/2 + ex_x*Arrayσ_m
	du[1:MyParameters.dim] = dψRe = HIm*ψRe + HRe*ψIm; #size dim)
	du[MyParameters.dim+1:end] = dψIm = HIm*ψIm - HRe*ψRe;
end

###################
# DIFFUSION TERM FOR SDE
function qb_dynamics_dW!(du,u, α, t)  #last action IS NOT STORED IN u
	ψRe = u[1:MyParameters.dim]   #(dim)
	ψIm = u[MyParameters.dim+1:end]

	HRe = Arrayσ_m

	du[1:MyParameters.dim] = dψRe = HRe*ψRe; #size (dim)
	du[MyParameters.dim+1:end] = dψIm = HRe*ψIm;
end

################################################
#model
@info("Constructing model...")
#state-aware, input = homodyne current values as integrated over the last n_substeps
state_1 = Dense(MyParameters.n_substeps, 256,relu ,initb = Flux.glorot_uniform)
state_2 = Dense(256, 256,relu, initb = Flux.glorot_uniform)
state_3 = Dense(256, 128,relu, initb = Flux.glorot_uniform)
#action-aware, input = last actions provided by the NN
nαin = 8  #number of last actions to be remembered
action_1 = Dense(nαin, 128,relu,initb = Flux.glorot_uniform)
action_2 = Dense(128, 128,relu,initb = Flux.glorot_uniform)
#combined
combine_1 = Dense(256, 64,relu,initb = Flux.glorot_uniform)
combine_2 = Dense(64, 32,relu,initb = Flux.glorot_uniform)
combine_3 = Dense(32, 1,softsign,initb = Flux.glorot_uniform)

state_aware=Chain(state_1, state_2, state_3)
action_aware=Chain(action_1,action_2)
combined=Chain(combine_1, combine_2, combine_3)

#this concatenates the full NN
struct Concat{T}
  catted::T
end
Concat(xs...) = Concat(xs)
Flux.@functor Concat

function (C::Concat)(x)
    mapreduce((f, x) -> f(x), vcat, C.catted, x)
end

model = Chain(Concat(state_aware, action_aware), combined);
p1, re = Flux.destructure(model)

###############################################
# initial state anywhere on the Bloch sphere
u0 = Array{Float32,2}(undef, 2*MyParameters.dim,n_par)
fill!(u0,0.0f0)
u0[2,:] .= 1.0f0 #down state

function prepare_initial!(u0) #random position on the Bloch sphere
	fill!(u0,0.0f0)
	theta = acos.(2*rand(n_par).-[1])  #uniform sampling for cos(theta) between -1 and 1
	phi = rand(n_par)*2*pi
	#real parts
	u0[1,:] += cos.(theta/2)
	u0[2,:] += sin.(theta/2).*cos.(phi)
	#imag parts
	#u0[3,:].+=0
	u0[4,:] += sin.(theta/2).*sin.(phi)
	#normalize initial state -already normalized by definition
	norm_factor = sqrt.(sum(u0[1:2*MyParameters.dim,:].^2,dims=1))
	u0[1:2*MyParameters.dim,:] = u0[1:2*MyParameters.dim,:]./norm_factor
	return u0
end

# target state
# ψtar = |up>
ut_complex = Array{Float32,1}(undef, MyParameters.dim)
fill!(ut_complex,0.0f0)
ut_complex[1,:] .= 1.0f0
Re_ut = real(Array{Float32}(ut_complex))
Im_ut = imag(Array{Float32}(ut_complex))
#To visualize
Fock_target = SVector{MyParameters.dim}(Re_ut.^2+Im_ut.^2)

############################################
#Mutate vector outside gradients!
# Necessary for collectiong results during the training
mut(A,index,x) = (A[index] = x)
mut_row(A,index,x) = (A[index,:] = x)
mut_vec(A,x) = (A .=x )
#tell Zygote not to take gradients of these
Zygote.@nograd mut
Zygote.@nograd mut_row
Zygote.@nograd mut_vec

##################################
# time range for the solver
t_interval = round(MyParameters.n_substeps*MyParameters.dt,digits=5)
tspan = (0.0f0,t_interval)

#########################################
#Static Arrays to collect results
mean_fid_store = zeros(MVector{MyParameters.n_steps+1,Float32})
std_fid_store = zeros(MVector{MyParameters.n_steps+1,Float32})
single_traj_fid_store = zeros(MMatrix{MyParameters.n_steps+1,n_par,Float32})
mean_action_store = zeros(MVector{MyParameters.n_steps,Float32})
std_action_store = zeros(MVector{MyParameters.n_steps,Float32})
single_traj_action_store = zeros(MMatrix{MyParameters.n_steps,n_par,Float32})
Fock_end_example = zeros(MVector{MyParameters.dim,Float32})

#########################################
# compute loss
###########################

#normalize the state using callback function
condition(u,t,integrator) = true
function affect!(integrator)
	integrator.u=integrator.u/norm(integrator.u)
end
cb = DiscreteCallback(condition,affect!,save_positions=(false,true))

CreateGrid(W1) =  (NoiseGrid(Array{Float32}((0.0:MyParameters.dt:(t_interval+MyParameters.dt))),W1))
Zygote.@nograd CreateGrid #avoid taking grads of this function

# set scalar random process
W = sqrt(MyParameters.dt)*randn(Float32,MyParameters.n_substeps+1) #for 1 trajectory
W1 = cumsum([zero(MyParameters.dt); W[1:end-1]], dims=1)
NG = CreateGrid(W1)
# define SDE problem
prepare_initial!(u0)
α = zeros(Float32,1,n_par)
prob = SDEProblem{true}(qb_dynamics_dt!,qb_dynamics_dW!, u0[:,1], tspan,α[1,1], noise=NG )

function compute_ex_x(sol)
	#size sol (2*dim, n_substeps+1, n_par)
	# Tr{σx*Reρ_i}...
	# Trick,  σx =(0  1)
	#             (1  0)... ..the trace effectively sums anti-diagonal terms
	# Reρ =  (     st. unimportant          ψRe[1]ψRe[2]+ψIm[1]ψIm[2])
	#	     (ψRe[1]ψRe[2]+ψIm[1]ψIm[2]        st. unimportant       )
	ψRe = @view sol[1:MyParameters.dim,:,:]   #(dim,n_substeps+1, n_par)
	ψIm = @view sol[MyParameters.dim+1:end,:,:]
	ex_x = 2*(ψRe[1,:,:].*ψRe[2,:,:]+ψIm[1,:,:].*ψIm[2,:,:]) #size(n_substeps+1, n_par)
end

function loss_along_trajectory(p1)
	#initial values
	loss = 0.0f0
	u = u0 #initial state
	#prepare ini Jhom, suppose there is no noise in the input
	ex_x = transpose(compute_ex_x(u0)) # (1,n_par)
	Jhom = reshape(repeat(ex_x,MyParameters.n_substeps),(MyParameters.n_substeps,n_par))
	αin = zeros(Float32,nαin,n_par)  #initially all past actions set to zero

	function prob_func(prob, i, repeat)
		#prepare tge vector of Wiener Process
		NG = CreateGrid(W1j[:,i]) # EM and RKMil
		remake(prob,p=α[i:i],u0=u[:,i],noise=NG)
	end

	#initial state, save the features
	Re_uj = @view u[1:MyParameters.dim,:]
	Im_uj= @view u[MyParameters.dim+1:end,:]
	fid = zeros(Float32,1,n_par)
	fid += abs2.(sum(Re_ut.*Re_uj,dims=1))
#	fid += abs2.(sum(Im_ut.*Im_uj,dims=1)) #commented because Im_ut=0, for general target uncomment
#	fid += abs2.(sum(Im_ut.*Re_uj,dims=1))
	fid += abs2.(sum(Re_ut.*Im_uj,dims=1))
#	fid += 2*(sum(Re_ut.*Re_uj,dims=1)).*sum(Im_ut.*Im_uj,dims=1)
#	fid -= 2*(sum(Re_ut.*Im_uj,dims=1)).*sum(Im_ut.*Re_uj,dims=1)
	mut_row(single_traj_fid_store,1,fid)
	mut(mean_fid_store,1,mean(fid))
	mut(std_fid_store,1,std(fid))

	W = sqrt(MyParameters.dt)*randn(Float32,MyParameters.n_substeps*MyParameters.n_steps+1,n_par) #for 1 trajectory
	W1 = cumsum([zero(MyParameters.dt); W[1:end-1,:]], dims=1)

	for j in 1:MyParameters.n_steps
		global W1j #must be global to be seen in prob_func
		global α
		#update action from NN
		α = re(p1)([Jhom,αin]).*MyParameters.force_mag #dim (1,n_par)
		αin = vcat(αin[2:end,:],α)
		#define ensemble problem
		W1j = @view W1[1+(j-1)*MyParameters.n_substeps:j*MyParameters.n_substeps+1,:]
		ensembleprob = EnsembleProblem(prob,prob_func = prob_func)

		sol = Array(solve(ensembleprob, RKMil(),
		ensemblealg=EnsembleThreads(), # EnsembleCPUArray(), EnsembleDistributed()
		sensealg = ForwardDiffSensitivity(),
		dt=MyParameters.dt, save_everystep=false,
		saveat=[t_interval],
		callback=cb, adaptive=false, trajectories=n_par, batch_size=n_par,
		save_noise=false,
		timeseries_errors=false,weak_timeseries_errors=false,weak_dense_errors=false
		))  #[2*dim,n_substeps+1, n_par]

		#collect <x> for all n_substeps and n_par
		ex_x = compute_ex_x(sol)  # (n_substeps, n_par), leave out the ini state
		dWj = W1j[2:end,:]-W1j[1:end-1,:]
		#NOISY MEASUREMENT at inidividual substeps
		Jhom = ex_x.*MyParameters.dt+dWj
		#current state
		u = @view sol[:,end,:]
		#compute loss function
		Re_uj = @view u[1:MyParameters.dim,:]
		Im_uj = @view u[MyParameters.dim+1:end,:]
		fid = zeros(Float32,1,n_par)
		fid += abs2.(sum(Re_ut.*Re_uj,dims=1))
	#	fid += abs2.(sum(Im_ut.*Im_uj,dims=1))
	#	fid += abs2.(sum(Im_ut.*Re_uj,dims=1))
		fid += abs2.(sum(Re_ut.*Im_uj,dims=1))
	#	fid += 2*(sum(Re_ut.*Re_uj,dims=1)).*sum(Im_ut.*Im_uj,dims=1)
	#	fid -= 2*(sum(Re_ut.*Im_uj,dims=1)).*sum(Im_ut.*Re_uj,dims=1)
		loss += C1*MyParameters.gamma^j*(1-mean(fid))
		mut(mean_fid_store,j+1,mean(fid))
		mut(std_fid_store,j+1,std(fid))
		mut_row(single_traj_fid_store,j+1,fid)
		mut(mean_action_store,j,mean(α))
		mut(std_action_store,j,std(α))
		mut_row(single_traj_action_store,j,α)
		#punish large actions--note, that for j we are pointing to the j-1st action!
		loss += C2*MyParameters.gamma^(j)*(mean(abs2.(α))) #mimic...sum max valus
   		#emphasize the main interval
   		if j>(MyParameters.n_steps-50)
  			loss += C3*MyParameters.gamma^j*(1-mean(fid))
   			mut_vec(Fock_end_example,Re_uj[:,1].^2+Im_uj[:,1].^2);
   		end
	end
	return loss
end

###################################
# training loop parameters
epochs=7000
println("total epochs: ",epochs)
training=zeros(epochs+1) #add extra dimension because of grad rejection condition in the training
maxgrads=zeros(epochs+1)
some_nans=zeros(epochs+1)

data = Iterators.repeated((), epochs)
opt = ADAM(0.0001)

#grad clipping
function clip(x)
    min.(max.(x,-40), 40)
end

using DelimitedFiles  #save txt files

## Save model
using BSON: @save
function qb_train!(loss, p1, data, opt,u0)
	ps = Flux.params(p1)
	iter = 1  #start with iter=1 because of the grad rejection lines below
	for d in data
		iter += 1
		prepare_initial!(u0)  #different initial states!
		#ini_fid=abs2.(sum(Re_ut.*u0[1:2],dims=1))+abs2.(sum(Re_ut.*u0[3:4],dims=1)) #initial mean fidelity
		@show iter
		#@show mean(ini_fid)
		@time gs = gradient(ps) do
			training_loss = loss(Zygote.hook(clip,p1)) #grad clipping does not work now???
			mut(training,iter,training_loss)
			println("loss: ",training_loss)
			return training_loss
		end
		maxgrads[iter]=maximum(abs.(gs[p1]))
		some_nans[iter]=sum(isnan.(gs[p1]))
		println("is nan: ",sum(isnan.(gs[p1])))
		println("max grad: ",maximum(abs.(gs[p1])))
		if iter%1 == 0
			fig1 = plot( [1:MyParameters.n_steps+1,1:MyParameters.n_steps+1],
			  [mean_fid_store mean_fid_store],
			    fillrange=[mean_fid_store.-std_fid_store mean_fid_store.+std_fid_store], fillalpha=0.3, c=:blue,
				xlabel = "steps", ylim=(0,1),xlim=(0,150), lw = 1.5, title="Fidelity" ,legend=false)
			fig1=plot!(1:MyParameters.n_steps+1,single_traj_fid_store[:,1], c=:black  )
			fig2 = plot( [1:MyParameters.n_steps,1:MyParameters.n_steps],
			  [mean_action_store mean_action_store],
			    fillrange=[mean_action_store.-std_action_store mean_action_store.+std_action_store], fillalpha=0.3, c=:orange,
				xlabel = "steps", ylim=(-MyParameters.force_mag,MyParameters.force_mag), xlim=(0,150), lw = 1.5, title="Action" ,legend=false)
			fig2=plot!(1:MyParameters.n_steps,single_traj_action_store[:,1], c=:black  )
			#display(plot(fig1, fig2, layout = (1, 2), legend = false,size=(800,300)));
			if iter == 1 || iter%100 == 0
				png(string("Figures/Figure_",iter)) #for saving figs
				epoch_output=hcat(1:MyParameters.n_steps,mean_fid_store[1:end-1],mean_action_store) #save results for epoch
				writedlm(string("Data/Epoch",iter,".txt"), epoch_output)
			end
			println("iter: ", iter)
		end
		println("++++++++++++++++++++++")
		# Only gradients that don't lead to a rapid increase of the loss function are accepted
		# Empirically it prevents instabilities in learning
		if training[iter]<1.1* training[iter-1]
			Flux.Optimise.update!(opt, ps, gs)
		else
			println("GRADS REJECTED!")
		end
		if iter%1000 == 0 #save the model every 1000 epochs
			@save string("model_Jhom-in_parallel_",iter,".bson")  p1 re opt
		end
	end
	@save string("Data/model_Jhom-in_parallel_END.bson")  p1 re opt
end

#training
Random.seed!(50)
qb_train!(loss_along_trajectory, p1, data, opt,u0)

#plotting of the loss function
plot(training[2:end], xlabel = "epochs", ylabel="loss", lw = 1.5, title = "Training")
savefig("Figures/Loss.pdf")
#save training a txt
writedlm("Data/Loss.txt", training[2:end])


text:
using DashBootstrapComponents

progress = dbc_progress(value=50);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


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

mutable struct ManagedServiceIdentityUserAssignedIdentities <: SwaggerModel
    principalId::Any # spec type: Union{ Nothing, String } # spec name: principalId
    clientId::Any # spec type: Union{ Nothing, String } # spec name: clientId

    function ManagedServiceIdentityUserAssignedIdentities(;principalId=nothing, clientId=nothing)
        o = new()
        validate_property(ManagedServiceIdentityUserAssignedIdentities, Symbol("principalId"), principalId)
        setfield!(o, Symbol("principalId"), principalId)
        validate_property(ManagedServiceIdentityUserAssignedIdentities, Symbol("clientId"), clientId)
        setfield!(o, Symbol("clientId"), clientId)
        o
    end
end # type ManagedServiceIdentityUserAssignedIdentities

const _property_map_ManagedServiceIdentityUserAssignedIdentities = Dict{Symbol,Symbol}(Symbol("principalId")=>Symbol("principalId"), Symbol("clientId")=>Symbol("clientId"))
const _property_types_ManagedServiceIdentityUserAssignedIdentities = Dict{Symbol,String}(Symbol("principalId")=>"String", Symbol("clientId")=>"String")
Base.propertynames(::Type{ ManagedServiceIdentityUserAssignedIdentities }) = collect(keys(_property_map_ManagedServiceIdentityUserAssignedIdentities))
Swagger.property_type(::Type{ ManagedServiceIdentityUserAssignedIdentities }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ManagedServiceIdentityUserAssignedIdentities[name]))}
Swagger.field_name(::Type{ ManagedServiceIdentityUserAssignedIdentities }, property_name::Symbol) =  _property_map_ManagedServiceIdentityUserAssignedIdentities[property_name]

function check_required(o::ManagedServiceIdentityUserAssignedIdentities)
    true
end

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


text:
using ReachabilityAnalysis
Dict([
    "opts" => Dict(:T=>5.0, :alg=>BOX(δ=0.01), :fixpoint_check=>false),
    "X0"   => [(1, Hyperrectangle(low=[68, 0.], high=[69, 0.]))],
    "info" => Dict(
        "name" => "clocked_thermostat",
        "dim" => 2,
        "linear" => true,
        "hybrid" => true
    )
])

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using WebSockets
const BAREHTML = "<meta http-equiv=\"Content-Type\" content=\"text/html; charset=UTF-8\">
 <title>Empty.html</title>"
import Sockets
const LOCALIP = string(Sockets.getipaddr())
const PORT = 8080
const BODY =  "
Press F12. In console:

ws = new WebSocket(\"ws://$LOCALIP:$PORT\")

ws.onmessage = function(e){console.log(e.data)}

ws.send(\"Browser console says hello!\")
                "

function coroutine(ws)
    @info "Started coroutine for " ws
    while isopen(ws)
        data, = readguarded(ws)
        s = String(data)
        if s == ""
            writeguarded(ws, "Goodbye!")
            break
        end
        @info "Received: $s"
        writeguarded(ws, "Hello! Send empty message to exit, or just leave.")
    end
    @info "Will now close " ws
end

function gatekeeper(req, ws)
    orig = WebSockets.origin(req)
    @info "\nOrigin: $orig   Target: $(req.target)   subprotocol: $(subprotocol(req))"
    if occursin(LOCALIP, orig)
        coroutine(ws)
    elseif orig == ""
        @info "Non-browser clients don't send Origin. We liberally accept the update request in this case:" ws
        coroutine(ws)
    else
        @warn "Inacceptable request"
    end
end

handle(req) = replace(BAREHTML, "" => BODY) |> WebSockets.Response

const server = WebSockets.ServerWS(handle,
                                    gatekeeper)

@info "In browser > $LOCALIP:$PORT , F12> console > ws = new WebSocket(\"ws://$LOCALIP:$PORT\") "
@async WebSockets.with_logger(WebSocketLogger()) do
    WebSockets.serve(server, LOCALIP, PORT)
end


text:
# This file was generated, do not modify it. # hide
mach = machine(Deterministic(), X, y; predict=yhat)


text:
CONFIG = Dict(
    "target.project_dir" => (@__DIR__) |> dirname |> abspath,
    "reporter.use_dataframe" => true,
)

ON_TRAVIS = get(ENV, "TRAVIS", "false") == "true"

if ON_TRAVIS
    BENCHMARK_FILES = [
        "dummy.jl",
        "gdemo.jl",
        "mvnormal.jl",
    ]
else
    BENCHMARK_FILES = [
        "dummy.jl",
        "gdemo.jl",
        "mvnormal.jl",
    ]
end

CONFIG["benchmark_files"] = map(BENCHMARK_FILES) do x joinpath(@__DIR__, x) end
CONFIG["comment_templates.ci_bm_complete_for_commit"] = """

The benchmark job for this commit has completed.
Here is the benchmark [report]( {{{ :report_url }}}) and commit ({{{ :report_repo }}}@{{{ :commit_id }}}).
The benchmark log (printed via `log_report`) can be found below.

```
{{#:report_logs}}
{{.}}
{{/:report_logs}}
```
"""


text:
## Implement the Newton's method update to estimate sqrt(x)
""" the function we are using here is y^2 - x, for which
the Newton update step is y - (y^2 - x)/(2y) """
function mysqrt(x; err = max(1.0,x)*eps(typeof(x/x)), maxiter=100)
  y = one(typeof(x/x))
  iter = 0
  while abs(y^2 - x) >= err && iter <= maxiter
    y = y - (y^2 - x)/(2*y)
    iter += 1
  end
  return y
end
@show mysqrt(500)
##
@show mysqrt(0.25)


text:
using Test, FranklinParser, Pkg
FP = FranklinParser
FPE = FP.FranklinParserException

include("testutils.jl")

@testset "utils" begin
    include("utils/strings.jl")
    include("utils/types.jl")
end

@testset "tokens" begin
    include("tokens/utils.jl")
    include("tokens/md_tokens.jl")
    include("tokens/html_tokens.jl")
end

@testset "blocks" begin
    include("blocks/md_blocks.jl")
    include("blocks/html_blocks.jl")
    include("blocks/utils.jl")
end

@testset "partition" begin
    include("partition/md_partition.jl") # depr
    include("partition/md_specs.jl")
    include("partition/html_partition.jl")
    include("partition/math_partition.jl")
    include("partition/split_args.jl")
end

@testset "cmark" begin
    # all the tests here are directly taken from
    # https://github.com/MichaelHatherly/CommonMark.jl/tree/master/test/samples/cmark
    include("cmark/inlines.jl")
    include("cmark/blocks.jl")
    include("cmark/misc.jl")
end

@testset "misc" begin
    @test FP.content("abc") == "abc"
    @test FP.content(FP.subs("abc")) == "abc"
end

@testset "bugfixes" begin
    include("misc_fixes.jl")
end


text:
function sort_gui(wave_points)

    grid = Grid()

    c_sort = Canvas(100,100)

    @guarded draw(c_sort) do widget
        ctx = Gtk.getgc(c_sort)
        set_source_rgb(ctx,0.0,0.0,0.0)
        paint(ctx)
    end
    show(c_sort)
    grid[1,2]=c_sort
    setproperty!(c_sort,:hexpand,true)
    setproperty!(c_sort,:vexpand,true)

    panel_grid = Grid()
    grid[2,2] = panel_grid

    b1 = Button("Plot")

    panel_grid[1,1]=b1

    clusteropts = MenuItem("_Cluster")
    clustermenu = Menu(clusteropts)
    cluster_km = MenuItem("K means")
    push!(clustermenu,cluster_km)
    cluster_dbscan = MenuItem("DBSCAN")
    push!(clustermenu,cluster_dbscan)

    mb = MenuBar()
    push!(mb, clusteropts)
    grid[1,1]=mb

    col_sb=SpinButton(1:3)
    panel_grid[1,3]=col_sb

    #Event
    popup_axis = Menu()

    popup_x = MenuItem("X Axis")
    push!(popup_axis,popup_x)
    popup_x_menu=Menu(popup_x)
    popup_y = MenuItem("Y Axis")
    push!(popup_axis,popup_y)
    popup_y_menu=Menu(popup_y)

    popup_pca1_x=MenuItem("PCA1")
    push!(popup_x_menu,popup_pca1_x)
    popup_pca2_x=MenuItem("PCA2")
    push!(popup_x_menu,popup_pca2_x)
    popup_pca3_x=MenuItem("PCA3")
    push!(popup_x_menu,popup_pca3_x)
    popup_peak_x=MenuItem("Peak")
    push!(popup_x_menu,popup_peak_x)
    popup_valley_x=MenuItem("Valley")
    push!(popup_x_menu,popup_valley_x)
    popup_pv_x=MenuItem("Peak-Valley")
    push!(popup_x_menu,popup_pv_x)
    popup_area_x=MenuItem("Area")
    push!(popup_x_menu,popup_area_x)

    popup_pca1_y=MenuItem("PCA1")
    push!(popup_y_menu,popup_pca1_y)
    popup_pca2_y=MenuItem("PCA2")
    push!(popup_y_menu,popup_pca2_y)
    popup_pca3_y=MenuItem("PCA3")
    push!(popup_y_menu,popup_pca3_y)
    popup_peak_y=MenuItem("Peak")
    push!(popup_y_menu,popup_peak_y)
    popup_valley_y=MenuItem("Valley")
    push!(popup_y_menu,popup_valley_y)
    popup_pv_y=MenuItem("Peak-Valley")
    push!(popup_y_menu,popup_pv_y)
    popup_area_y=MenuItem("Area")
    push!(popup_y_menu,popup_area_y)

    Gtk.showall(popup_axis)

    win = Window(grid,"Sort View") |> Gtk.showall

    myfeatures=Dict{String,Array{Float64,1}}("PCA-1"=>zeros(Float64,0),"PCA-2"=>zeros(Float64,0),"PCA-3"=>zeros(Float64,0),"Peak"=>zeros(Float64,0),"Valley"=>zeros(Float64,0),"Peak-Valley"=>zeros(Float64,0),"Area"=>zeros(Float64,0))

    handles = SortView(win,c_sort,b1,myfeatures,fit(PCA,rand(Float64,10,10)),false,1,1,popup_axis,1,falses(10,2),["Non" for i=1:20,j=1:2],col_sb,[FeaturePlot() for i=1:10],100.0,100.0,Buffer(wave_points))

    signal_connect(b1_cb,b1,"clicked",Void,(),false,(handles,))
    signal_connect(col_sb_cb,col_sb,"value-changed",Void,(),false,(handles,))

    signal_connect(canvas_press,c_sort,"button-press-event",Void,(Ptr{Gtk.GdkEventButton},),false,(handles,))

    signal_connect(popup_pca1_cb_x,popup_pca1_x,"activate",Void,(),false,(handles,))
    signal_connect(popup_pca2_cb_x,popup_pca2_x,"activate",Void,(),false,(handles,))
    signal_connect(popup_pca1_cb_y,popup_pca1_y,"activate",Void,(),false,(handles,))
    signal_connect(popup_pca2_cb_y,popup_pca2_y,"activate",Void,(),false,(handles,))
    signal_connect(popup_pca3_cb_x,popup_pca3_x,"activate",Void,(),false,(handles,))
    signal_connect(popup_pca3_cb_y,popup_pca3_y,"activate",Void,(),false,(handles,))
    signal_connect(popup_peak_cb_x,popup_peak_x,"activate",Void,(),false,(handles,))
    signal_connect(popup_peak_cb_y,popup_peak_y,"activate",Void,(),false,(handles,))
    signal_connect(popup_valley_cb_x,popup_valley_x,"activate",Void,(),false,(handles,))
    signal_connect(popup_valley_cb_y,popup_valley_y,"activate",Void,(),false,(handles,))
    signal_connect(popup_pv_cb_x,popup_pv_x,"activate",Void,(),false,(handles,))
    signal_connect(popup_pv_cb_y,popup_pv_y,"activate",Void,(),false,(handles,))

    signal_connect(popup_area_cb_x,popup_area_x,"activate",Void,(),false,(handles,))
    signal_connect(popup_area_cb_y,popup_area_y,"activate",Void,(),false,(handles,))

    id = signal_connect(win_resize_cb, win, "size-allocate",Void,(Ptr{Gtk.GdkRectangle},),false,(handles,))

    handles
end

function b1_cb(widget::Ptr,user_data::Tuple{SortView})

    han, = user_data
    replot_sort(han)

    nothing
end

function b2_cb(widget::Ptr,user_data::Tuple{SortView})

    han, = user_data
    recalc_features(han)
    replot_sort(han)

    nothing
end

function recalc_features(han::SortView)

    han.pca_calced=false
    reset_pca(han,1)
    reset_pca(han,2)
    reset_pca(han,3)
    reset_peak(han)
    reset_valley(han)
    reset_pv(han)
    reset_area(han)

    for i=han.n_col*han.n_row

        if han.axes_name[i,1]=="Non"
            han.axes_name[i,1]="Peak"
        end
        if han.axes_name[i,2]=="Non"
            han.axes_name[i,2]="Peak"
        end

        scale_xaxis(han,i)
        scale_yaxis(han,i)
    end

    nothing
end

function col_sb_cb(widget::Ptr,user_data::Tuple{SortView})

    han, = user_data
    han.n_col=getproperty(han.col_sb,:value,Int)

    replot_sort(han)

    nothing
end

popup_pca1_cb_x(widget::Ptr,han::Tuple{SortView})=pca_calc(han[1],1,1)
popup_pca2_cb_x(widget::Ptr,han::Tuple{SortView})=pca_calc(han[1],2,1)
popup_pca3_cb_x(widget::Ptr,han::Tuple{SortView})=pca_calc(han[1],3,1)
popup_pca1_cb_y(widget::Ptr,han::Tuple{SortView})=pca_calc(han[1],1,2)
popup_pca2_cb_y(widget::Ptr,han::Tuple{SortView})=pca_calc(han[1],2,2)
popup_pca3_cb_y(widget::Ptr,han::Tuple{SortView})=pca_calc(han[1],3,2)

popup_peak_cb_x(widget::Ptr,han::Tuple{SortView})=peak_calc(han[1],1)
popup_peak_cb_y(widget::Ptr,han::Tuple{SortView})=peak_calc(han[1],2)

popup_valley_cb_x(widget::Ptr,han::Tuple{SortView})=valley_calc(han[1],1)
popup_valley_cb_y(widget::Ptr,han::Tuple{SortView})=valley_calc(han[1],2)

popup_pv_cb_x(widget::Ptr,han::Tuple{SortView})=pv_calc(han[1],1)
popup_pv_cb_y(widget::Ptr,han::Tuple{SortView})=pv_calc(han[1],2)

popup_area_cb_x(widget::Ptr,han::Tuple{SortView})=area_calc(han[1],1)
popup_area_cb_y(widget::Ptr,han::Tuple{SortView})=area_calc(han[1],2)

function reset_pca(han::SortView,num::Int64)
    if !han.pca_calced
        han.pca = fit(PCA,convert(Array{Float64,2},han.buf.spikes[:,1:han.buf.ind]))
        han.pca_calced=true
    end

    if num>size(han.pca.proj,2)
        han.features[string("PCA-",num)]=zeros(Float64,han.buf.ind)
    else
        han.features[string("PCA-",num)] = squeeze(han.pca.proj[:,num]' * han.buf.spikes[:,1:han.buf.ind],1)

    end
    nothing
end

function pca_calc(han::SortView,num::Int64,myaxis::Int64)

    reset_pca(han,num)
    han.axes_name[han.selected_plot,myaxis]=string("PCA-",num)
    scale_axis(han,myaxis)

    replot_sort(han)
    nothing
end

function reset_peak(han::SortView)
    han.features["Peak"]=zeros(Float64,han.buf.ind)
    for i=1:han.buf.ind
        han.features["Peak"][i]=maximum(han.buf.spikes[:,i])
    end
    nothing
end

function peak_calc(han::SortView,myaxis::Int64)

    reset_peak(han)
    han.axes_name[han.selected_plot,myaxis]=string("Peak")
    scale_axis(han,myaxis)

    replot_sort(han)
end

function reset_valley(han::SortView)
    han.features["Valley"]=zeros(Float64,han.buf.ind)
    for i=1:han.buf.ind
        han.features["Valley"][i]=minimum(han.buf.spikes[:,i])
    end
    nothing
end

function valley_calc(han::SortView,myaxis::Int64)

    reset_valley(han)
    han.axes_name[han.selected_plot,myaxis]=string("Valley")
    scale_axis(han,myaxis)

    replot_sort(han)
end

function reset_pv(han::SortView)
    han.features["Peak-Valley"]=zeros(Float64,han.buf.ind)
    for i=1:han.buf.ind
        han.features["Peak-Valley"][i]=maximum(han.buf.spikes[:,i])-minimum(han.buf.spikes[:,i])
    end
    nothing
end

function pv_calc(han::SortView,myaxis::Int64)

    reset_pv(han)
    han.axes_name[han.selected_plot,myaxis]=string("Peak-Valley")
    scale_axis(han,myaxis)

    replot_sort(han)
    nothing
end

function reset_area(han::SortView)
    han.features["Area"]=zeros(Float64,han.buf.ind)
    for i=1:han.buf.ind
        for j=1:size(han.buf.spikes,1)
            han.features["Area"][i]+=maximum(han.buf.spikes[j,i])
        end
    end
    nothing
end

function area_calc(han::SortView,myaxis::Int64)

    reset_area(han)
    han.axes_name[han.selected_plot,myaxis]=string("Area")
    scale_axis(han,myaxis)

    replot_sort(han)
end

function scale_xaxis(han::SortView,myplot::Int64)
    han.plots[myplot].xmin=minimum(han.features[han.axes_name[myplot,1]])
    han.plots[myplot].xscale=maximum(han.features[han.axes_name[myplot,1]])-han.plots[myplot].xmin
end

function scale_yaxis(han::SortView,myplot::Int64)
    han.plots[myplot].ymin=minimum(han.features[han.axes_name[myplot,2]])
    han.plots[myplot].yscale=maximum(han.features[han.axes_name[myplot,2]])-han.plots[myplot].ymin
end

function scale_axis(han::SortView,myaxis::Int64)

    if myaxis==1
        scale_xaxis(han,han.selected_plot)
    else
        scale_yaxis(han,han.selected_plot)
    end

    han.axes[han.selected_plot,myaxis]=true
    nothing
end

function name_axis(han::SortView,myname)

end

function canvas_press(widget::Ptr,param_tuple,user_data::Tuple{SortView})

    han, = user_data

    event = unsafe_load(param_tuple)

    inaxis = get_axis_bounds(han,event.x,event.y)

    if event.button==1
        rubberband_start(han,event.x,event.y)
    elseif event.button==3
        popup(han.popup_axis,event)
    end

    nothing
end

function get_axis_bounds(han::SortView,x,y)

    xbounds=linspace(0.0,han.w,han.n_col+1)
    ybounds=linspace(0.0,han.h,han.n_row+1)

    count=1
    for yy=2:length(ybounds), xx=2:length(xbounds)
        if (x<xbounds[xx])&(y<ybounds[yy])
            han.selected_plot=count
            break
        end
        count+=1
    end

    nothing
end

function replot_sort(han::SortView)

    ctx=Gtk.getgc(han.c)
    set_source_rgb(ctx,0.0,0.0,0.0)
    paint(ctx)

    #Check if something has changed in underlying data
    if han.axes[1,1]&han.axes[1,2]
        if length(han.features[han.axes_name[1,1]]) != han.buf.ind
            recalc_features(han)
        end
    end

    prepare_plots(han)

    for jj=1:(han.n_col*han.n_row)

        if han.axes[jj,1]&han.axes[jj,2]

            xmin=han.plots[jj].xmin
            ymin=han.plots[jj].ymin
            xscale=han.plots[jj].xscale
            yscale=han.plots[jj].yscale

            xdata = han.features[han.axes_name[jj,1]]
            ydata = han.features[han.axes_name[jj,2]]

            Cairo.translate(ctx,50+han.w/(han.n_col)*(jj-1),1)
            Cairo.scale(ctx,(han.w/(han.n_col)-70)/xscale,(han.h/(han.n_row)-50)/yscale)

            for ii=1:(maximum(han.buf.clus)+1)
                for i=1:han.buf.ind
                    if (han.buf.clus[i]+1 == ii)&(han.buf.mask[i])

                        move_to(ctx,xdata[i]-xmin,ydata[i]-ymin)

                        line_to(ctx,xdata[i]-xmin+10.0,ydata[i]-ymin+10.0)
                    end
                end
                select_color(ctx,ii)
                stroke(ctx)
            end

            identity_matrix(ctx)

            set_source_rgb(ctx,1.0,1.0,1.0)
            move_to(ctx,han.w/(han.n_col*2)+han.w/(han.n_col)*(jj-1),han.h-10.0)
            show_text(ctx,han.axes_name[jj,1])

            move_to(ctx,10.0+han.w/han.n_col*(jj-1),han.h/2)
            rotate(ctx,-pi/2)
            show_text(ctx,han.axes_name[jj,2])

            identity_matrix(ctx)
        end
    end
    reveal(han.c)

    nothing
end

function prepare_plots(han::SortView)

    ctx=Gtk.getgc(han.c)

    xbounds=linspace(0.0,han.w,han.n_col+1)
    ybounds=linspace(0.0,han.h,han.n_row+1)

    for i=2:length(ybounds)-1
        move_to(ctx,0.0,ybounds[i])
        line_to(ctx,han.w,ybounds[i])
    end

    for i=2:length(xbounds)-1
        move_to(ctx,xbounds[i],0.0)
        line_to(ctx,xbounds[i],han.h)
    end

    nothing
end

identity_matrix(ctx)=ccall((:cairo_identity_matrix,Cairo._jl_libcairo),Void, (Ptr{Void},), ctx.ptr)

function select_color(ctx,clus,alpha=1.0)

    if clus==1
        set_source_rgba(ctx,1.0,1.0,1.0,alpha) # white
    elseif clus==2
        set_source_rgba(ctx,1.0,1.0,0.0,alpha) #Yellow
    elseif clus==3
        set_source_rgba(ctx,0.0,1.0,0.0,alpha) #Green
    elseif clus==4
        set_source_rgba(ctx,0.0,0.0,1.0,alpha) #Blue
    elseif clus==5
        set_source_rgba(ctx,1.0,0.0,0.0,alpha) #Red
    else
        set_source_rgba(ctx,1.0,1.0,0.0,alpha)
    end

    nothing
end

#=
Rubber Band functions adopted from GtkUtilities.jl package by Tim Holy 2015
=#

function rb_draw(r::Cairo.CairoContext, rb::RubberBand)
    rb_set(r, rb)
    set_line_width(r, 1)

    set_source_rgb(r, 1, 1, 1)
    stroke_preserve(r)
end

function rb_set(r::Cairo.CairoContext, rb::RubberBand)
    move_to(r, rb.pos1.x, rb.pos1.y)
    rel_line_to(r,rb.pos2.x-rb.pos1.x, rb.pos2.y-rb.pos1.y)
end

function rubberband_start(han::SortView, x, y; minpixels::Int=2)
    r = Gtk.getgc(han.c)
    Cairo.save(r)
    ctxcopy = copy(r)
    rb = RubberBand(Vec2(x,y),Vec2(x,y), Vec2(x,y), [Vec2(x,y)],false, minpixels)
    push!((han.c.mouse, :button1motion),  (c, event) -> rubberband_move(han.c, rb, event.x, event.y, ctxcopy))
    push!((han.c.mouse, :motion), Gtk.default_mouse_cb)
    push!((han.c.mouse, :button1release), (c, event) -> rubberband_stop(han, rb, event.x, event.y, ctxcopy))
    nothing
end

function rubberband_move(c::Canvas, rb::RubberBand, x, y, ctxcopy)
    r = Gtk.getgc(c)
    if rb.moved
        #rb_erase(r, ctxcopy)
    end
    rb.moved = true

    # Draw the new rubberband
    rb.pos2 = Vec2(x, y)
    push!(rb.polygon,rb.pos2)
    rb_draw(r, rb)
    rb.pos1=rb.pos2
    reveal(c, false)
end

function rubberband_stop(han::SortView, rb::RubberBand, x, y, ctxcopy)
    pop!((han.c.mouse, :button1motion))
    pop!((han.c.mouse, :motion))
    pop!((han.c.mouse, :button1release))
    if !rb.moved
        return
    end
    r = Gtk.getgc(han.c)
    rb_set(r, rb)
    restore(r)
    set_line_width(r,3.0)
    move_to(r,rb.polygon[1].x,rb.polygon[1].y)
    for i=2:length(rb.polygon)
        line_to(r,rb.polygon[i].x,rb.polygon[i].y)
    end
    stroke(r)

    if han.buf.selected_clus>0
        inside_polygon(rb.polygon,han)
    end
    replot_sort(han)
    reveal(han.c, false)
    nothing
end

function canvas_to_feature(han,x1,y1,myplot)

    xtrans=50+han.w/(han.n_col)*(myplot-1)
    ytrans=0

    xscale=(han.w/han.n_col-70)/han.plots[myplot].xscale
    yscale=(han.h/han.n_row-50)/han.plots[myplot].yscale

    x1=(x1 - xtrans)/(xscale)+han.plots[myplot].xmin
    y1=(y1 - ytrans)/(yscale)+han.plots[myplot].ymin

    (x1,y1)
end

function inside_polygon(xy::Array{Vec2,1},han::SortView)

    xmin=xy[1].x
    ymin=xy[1].y
    xmax=xy[1].x
    ymax=xy[1].y

    #Find bounds of lasso
    for i=2:length(xy)
        if xy[i].x<xmin
            xmin=xy[i].x
        elseif xy[i].x>xmax
            xmax=xy[i].x
        end

        if xy[i].y<ymin
            ymin=xy[i].y
        elseif xy[i].y>ymax
            ymax=xy[i].y
        end
    end

    #Convert canvas coordinates to PCA space coordinates
    (xmin,ymin)=canvas_to_feature(han,xmin,ymin,han.selected_plot)
    (xmax,ymax)=canvas_to_feature(han,xmax,ymax,han.selected_plot)

    xdata = han.features[han.axes_name[han.selected_plot,1]]
    ydata = han.features[han.axes_name[han.selected_plot,2]]

    han.buf.selected=trues(han.buf.ind)

    for i=1:han.buf.ind

        px=xdata[i]
        py=ydata[i]
        if ((px>xmin)&(px<xmax))&((py>ymin)&(py<ymax))
            han.buf.selected[i]=false
        end
    end

    han.buf.c_changed=true
    han.buf.replot=true
    nothing
end

function win_resize_cb(widget::Ptr,param_tuple,user_data::Tuple{SortView})

    han, = user_data

    ctx=Gtk.getgc(han.c)
    han.h=height(ctx)
    han.w=width(ctx)

    replot_sort(han)

    nothing
end


text:
function epsilon_dominates_clear(u, v, epsilon)
    veps = v .+ epsilon
    all(u .<= veps) && any(u .< veps)
end

epsilon_dominates_fast(u, v, epsilon) = pareto_dominates_fast(u, v .+ epsilon)


text:
using Conda
#Conda.add_channel("conda-forge")
#Conda.add("tensorflow")

cmd = `pip3 install tensorflow --user`
run(cmd)


text:
struct CSVSchema{T₁ <: AbstractVendor,  T₂ <: AbstractProduct, T₃ <: AbstractData} <: AbstractSchema
    vendor::T₁
    product::T₂
    data::T₃
end