text:
module PlutoSliderServer

using FromFile

@from "./MoreAnalysis.jl" import bound_variable_connections_graph
@from "./FileHelpers.jl" import find_notebook_files_recursive, list_files_recursive
@from "./Export.jl" import generate_index_html
@from "./Actions.jl" import process,
    should_shutdown, should_update, should_launch, will_process
@from "./Types.jl" import NotebookSession
@from "./Lock.jl" import withlock
@from "./Configuration.jl" import PlutoDeploySettings,
    ExportSettings, SliderServerSettings, get_configuration
@from "./ConfigurationDocs.jl" import @extract_docs,
    get_kwdocs, list_options_md, list_options_toml
@from "./ReloadFolder.jl" import update_sessions!, select
@from "./HTTPRouter.jl" import make_router
@from "./gitpull.jl" import fetch_pull

import Pluto
import Pluto:
    ServerSession,
    Firebasey,
    Token,
    withtoken,
    pluto_file_extensions,
    without_pluto_file_extension
using HTTP
using Base64
using SHA
using Sockets
import BetterFileWatching: watch_folder
using TerminalLoggers: TerminalLogger
using Logging: global_logger
using GitHubActions: GitHubActionsLogger

export export_directory, run_directory, run_git_directory, github_action
export export_notebook, run_notebook

export show_sample_config_toml_file

function __init__()
    if get(ENV, "GITHUB_ACTIONS", "false") == "true"
        global_logger(GitHubActionsLogger())
    else
        global_logger(try
            TerminalLogger(; margin=1)
        catch
            TerminalLogger()
        end)
    end
end

const sample_config_toml_file = """
# WARNING: this sample configuration file contains settings for **all options**, to demonstrate what is possible. For most users, we recommend keeping the configuration file small, and letting PlutoSliderServer choose the default settings automatically.

# This means: DO NOT use this file in your setup, instead, create an empty toml file and **add only the settings that you want to change**.

[Export]
$(list_options_toml(ExportSettings))

[SliderServer]
$(list_options_toml(SliderServerSettings))

[Pluto]
[Pluto.compiler]
threads = 1

# See documentation for `Pluto.Configuration` for the full list of options. You need specify the categories within `Pluto.Configuration.Options` (`compiler`, `evaluation`, etc).
"""

function show_sample_config_toml_file()
    Text(sample_config_toml_file)
end

merge_recursive(a::AbstractDict, b::AbstractDict) = mergewith(merge_recursive, a, b)
merge_recursive(a, b) = b

merge_recursive(a::PlutoDeploySettings, b::PlutoDeploySettings) = Configurations.from_dict(
    PlutoDeploySettings,
    merge_recursive(Configurations.to_dict(a), Configurations.to_dict(b)),
)

with_kwargs(original::PlutoDeploySettings; kwargs...) =
    merge_recursive(original, Configurations.from_kwargs(PlutoDeploySettings; kwargs...))

showall(xs) = Text(join(string.(xs), "\n"))

default_config_path() = joinpath(Base.active_project() |> dirname, "PlutoDeployment.toml")

macro ignorefailure(x)
    quote
        try
            $(esc(x))
        catch e
            showerror(stderr, e, catch_backtrace())
        end
    end
end

"""
    export_directory(start_dir::String="."; kwargs...)

Search recursively for all Pluto notebooks in the current folder, and for each notebook:
- Run the notebook and wait for all cells to finish
- Export the state object
- Create a .html file with the same name as the notebook, which has:
  - The JS and CSS assets to load the Pluto editor
  - The state object embedded
  - Extra functionality enabled, such as hidden UI, binder button, and a live bind server

# Keyword rguments
$(list_options_md(ExportSettings; prefix="Export"))
- `notebook_paths::Vector{String}=find_notebook_files_recursive(start_dir)`: If you do not want the recursive save behaviour, then you can set this to a vector of absolute paths. In that case, `start_dir` is ignored, and you should set `Export_output_dir`.
"""
function export_directory(args...; kwargs...)
    run_directory(args...; Export_enabled=true, SliderServer_enabled=false, kwargs...)
end
"""
    export_notebook(notebook_filename::String; kwargs...)

A single-file version of [`export_directory`](@ref).
"""
export_notebook(p; kwargs...) = run_notebook(
    p;
    Export_enabled=true,
    SliderServer_enabled=false,
    Export_create_index=false,
    kwargs...,
)
github_action = export_directory

"""
    run_notebook(notebook_filename::String; kwargs...)

A single-file version of [`run_directory`](@ref).
"""
function run_notebook(path::String; kwargs...)
    path = Pluto.tamepath(path)
    @assert isfile(path)
    run_directory(dirname(path); notebook_paths=[basename(path)], kwargs...)
end

"""
    run_directory(start_dir::String="."; export_options...)

Run the Pluto bind server for all Pluto notebooks in the given directory (recursive search).

# Keyword arguments
$(list_options_md(SliderServerSettings; prefix="SliderServer"))
- `notebook_paths::Union{Nothing,Vector{String}}=nothing`: If you do not want the recursive save behaviour, then you can set this to a vector of absolute paths. In that case, `start_dir` is ignored, and you should set `Export_output_dir`.
- `Export_enabled::Bool=true`: Also export HTML files?

---

## Export keyword arguments

If `Export_enabled` is `true`, then additional `Export_` keywords can be given:
$(list_options_md(ExportSettings; prefix="Export"))
"""
function run_directory(
    start_dir::String=".";
    notebook_paths::Union{Nothing,Vector{String}}=nothing,
    on_ready::Function=((args...) -> ()),
    config_toml_path::Union{String,Nothing}=default_config_path(),
    kwargs...,
)

    @assert joinpath("a", "b") == "a/b" "PlutoSliderServer does not work on Windows yet!"

    start_dir = Pluto.tamepath(start_dir)
    @assert isdir(start_dir)

    settings = get_configuration(config_toml_path; kwargs...)
    output_dir = something(
        settings.Export.output_dir,
        settings.SliderServer.enabled ? mktempdir() : start_dir,
    )
    mkpath(output_dir)

    if joinpath("a", "b") != "a/b"
        @error "PlutoSliderServer.jl is only designed to work on unix systems."
        exit()
    end

    function getpaths()
        all_nbs =
            notebook_paths !== nothing ? notebook_paths :
            find_notebook_files_recursive(start_dir)

        if settings.Export.enabled
            if settings.SliderServer.enabled
                setdiff(all_nbs, settings.SliderServer.exclude ∩ settings.Export.exclude)
            else
                setdiff(all_nbs, settings.Export.exclude)
            end
        else
            s = setdiff(all_nbs, settings.SliderServer.exclude)
            filter(s) do f
                occursin("@bind", read(joinpath(start_dir, f), String))
            end
        end
    end

    to_run = getpaths()

    @info "Settings" Text(settings)

    settings.SliderServer.enabled &&
        @warn "Make sure that you run this slider server inside a containerized environment -- it is not intended to be secure. Assume that users can execute arbitrary code inside your notebooks."

    # if to_run != notebook_paths
    #     @info "Excluded notebooks:" showall(setdiff(notebook_paths, to_run))
    # end

    settings.Pluto.server.disable_writing_notebook_files = true
    settings.Pluto.evaluation.lazy_workspace_creation = true

    server_session = Pluto.ServerSession(; options=settings.Pluto)

    notebook_sessions = NotebookSession[]
    # notebook_sessions = NotebookSession[QueuedNotebookSession(;path, current_hash=myhash(read(joinpath(start_dir, path)))) for path in to_run]

    if settings.SliderServer.enabled
        static_dir =
            (settings.Export.enabled && settings.SliderServer.serve_static_export_folder) ?
            output_dir : nothing
        router = make_router(notebook_sessions, server_session; settings, static_dir)

        # This is boilerplate HTTP code, don't read it
        host = settings.SliderServer.host
        port = settings.SliderServer.port

        # This is boilerplate HTTP code, don't read it
        hostIP = parse(Sockets.IPAddr, host)
        if port === nothing
            port, serversocket = Sockets.listenany(hostIP, UInt16(1234))
        else
            serversocket = try
                Sockets.listen(hostIP, UInt16(port))
            catch e
                @error "Port with number $port is already in use."
                return
            end
        end

        address = let
            host_str = string(hostIP)
            host_pretty = if isa(hostIP, Sockets.IPv6)
                if host_str == "::1"
                    "localhost"
                else
                    "[$(host_str)]"
                end
            elseif host_str == "127.0.0.1" # Assuming the other alternative is IPv4
                "localhost"
            else
                host_str
            end
            port_pretty = Int(port)
            "http://$(host_pretty):$(port_pretty)/"
        end

        @info "# Starting server..." address

        # This is boilerplate HTTP code, don't read it
        # We start the HTTP server before launching notebooks so that the server responds to heroku/digitalocean garbage fast enough
        http_server_task = @async HTTP.serve(
            hostIP,
            UInt16(port),
            stream=true,
            server=serversocket,
        ) do http::HTTP.Stream
            request::HTTP.Request = http.message
            request.body = read(http)
            HTTP.closeread(http)

            params = HTTP.queryparams(HTTP.URI(request.target))

            response_body = HTTP.handle(router, request)

            request.response::HTTP.Response = response_body
            request.response.request = request
            try
                HTTP.setheader(http, "Referrer-Policy" => "origin-when-cross-origin")
                HTTP.startwrite(http)
                write(http, request.response.body)
                HTTP.closewrite(http)
            catch e
                if isa(e, Base.IOError) || isa(e, ArgumentError)
                    # @warn "Attempted to write to a closed stream at $(request.target)"
                else
                    rethrow(e)
                end
            end
        end
    else
        http_server_task = @async 1 + 1
        serversocket = nothing
    end

    if (
        settings.Export.enabled &&
        settings.Export.create_index &&
        # If `settings.SliderServer.serve_static_export_folder`, then we serve a dynamic index page (inside HTTPRouter.jl), so we don't want to create a static index page.
        !(settings.SliderServer.enabled && settings.SliderServer.serve_static_export_folder)
    )

        exists = any([
            "index.html",
            "index.md",
            ("index" * e for e in pluto_file_extensions)...,
        ]) do f
            joinpath(output_dir, f) |> isfile
        end
        if !exists
            write(
                joinpath(output_dir, "index.html"),
                generate_index_html((
                    without_pluto_file_extension(path) =>
                        without_pluto_file_extension(path) * ".html" for path in to_run
                )),
            )
        end
    end

    function refresh_until_synced(check_dir_on_every_step::Bool, did_something::Bool=false)
        should_continue = withlock(notebook_sessions) do

            check_dir_on_every_step &&
                update_sessions!(notebook_sessions, getpaths(); start_dir)

            # todo: try catch to release lock?
            to_shutdown = select(should_shutdown, notebook_sessions)
            to_update = select(should_update, notebook_sessions)
            to_launch = select(should_launch, notebook_sessions)

            s = something(to_shutdown, to_update, to_launch, "not found")

            if s != "not found"

                progress = "[$(
                    count(!will_process, notebook_sessions) + 1
                )/$(length(notebook_sessions))]"

                new = process(s; server_session, settings, output_dir, start_dir, progress)
                if new !== s
                    if new isa NotebookSession{Nothing,Nothing,<:Any}
                        # remove it
                        filter!(!isequal(s), notebook_sessions)
                    elseif s !== new
                        replace!(notebook_sessions, s => new)
                    end
                end

                true
            else
                did_something && @info "# ALL NOTEBOOKS READY"
                false
            end
        end

        should_continue && refresh_until_synced(check_dir_on_every_step, true)
    end

    # RUN ALL NOTEBOOKS AND KEEP THEM RUNNING
    update_sessions!(notebook_sessions, getpaths(); start_dir)
    refresh_until_synced(false)

    should_watch = settings.SliderServer.enabled && settings.SliderServer.watch_dir

    watch_dir_task = Pluto.@asynclog if should_watch
        @info "Watching directory for changes..."
        debounced = kind_of_debounced() do _
            @debug "File change detected!"
            sleep(0.5)
            refresh_until_synced(true)
        end
        watch_folder(debounced, start_dir)
    end

    if should_watch
        # todo: skip first watch_folder so that we dont need this sleepo
        sleep(2)
    end

    on_ready((; serversocket, server_session, notebook_sessions))

    try
        wait(http_server_task)
    catch e
        @ignorefailure close(serversocket)
        @ignorefailure schedule(watch_dir_task, e; error=true)
        e isa InterruptException || rethrow(e)
    end
end

"""
    run_git_directory(start_dir::String="."; export_options...)

Like [`run_directory`](@ref), but will automatically keep running `git pull` in `start_dir` and update the slider server to match changes in the directory. See our README for more info.

If you use a `PlutoDeployment.toml` file, we also keep checking whether this file changed. If it changed, we **exit Julia session**, and you will have to restart it. Open an issue if this is a problem.
"""
function run_git_directory(
    start_dir::String=".";
    notebook_paths::Union{Nothing,Vector{String}}=nothing,
    on_ready::Function=((args...) -> ()),
    config_toml_path::Union{String,Nothing}=default_config_path(),
    kwargs...,
)

    start_dir = Pluto.tamepath(start_dir)
    @assert isdir(start_dir)

    get_settings() =
        get_configuration(config_toml_path; SliderServer_watch_dir=true, kwargs...)
    old_settings = get_settings()

    run_dir_task = Pluto.@asynclog begin
        run_directory(
            start_dir;
            notebook_paths,
            on_ready,
            config_toml_path,
            SliderServer_watch_dir=true,
            kwargs...,
        )
    end
    pull_loop_task = Pluto.@asynclog while true
        new_settings = get_settings()

        if old_settings != new_settings
            @error "Configuration changed. Shutting down!"

            println(stderr, "Old settings:")
            println(stderr, repr(old_settings))
            println(stderr, "")
            println(stderr, "New settings:")
            println(stderr, repr(new_settings))

            # @ignorefailure schedule(run_dir_task[], InterruptException(); error=true)
            exit()
        end

        sleep(5)
        fetch_pull(start_dir)
    end

    waitall([run_dir_task, pull_loop_task])
end

function waitall(tasks)
    killing = Ref(false)
    @sync for t in tasks
        try
            wait(t)
        catch e
            if !(e isa InterruptException)
                showerror(stderr, e, catch_backtrace())
            end
            if !killing[]
                killing[] = true
                for t2 in tasks
                    try
                        schedule(t2, e; error=true)
                    catch e
                    end
                end
            end
        end
    end
end

function kind_of_debounced(f)
    update_waiting = Ref(true)
    running = Ref(false)

    function go(args...)
        # @info "trigger" update_waiting[] running[]
        update_waiting[] = true
        if !running[]
            running[] = true
            update_waiting[] = false

            f(args...)
            running[] = false
        end
        update_waiting[] && go(args...)
    end

    return go
end

end


text:
using LapSolve
using Test

@testset "StiffWrapper" begin
     A = [3.0 5.0 6.0 9.0 9.0;
          2.0 4.0 3.0 5.0 6.0;
          4.0 1.0 4.0 9.0 9.0;
          4.0 7.0 8.0 5.0 10.0;
          1.0 1.0 10.0 10.0 5.0]
     M = [2.0 4.0 5.0 8.0 8.0 18.0 Inf Inf Inf Inf;
          1.0 3.0 2.0 4.0 5.0 Inf 18.0 Inf Inf Inf;
          3.0 0.0 3.0 8.0 8.0 Inf Inf 18.0 Inf Inf;
          3.0 6.0 7.0 4.0 9.0 Inf Inf Inf 18.0 Inf;
          0.0 0.0 9.0 9.0 4.0 Inf Inf Inf Inf 18.0;
          18.0 Inf Inf Inf Inf 0.0 0.0 0.0 0.0 0.0;
          Inf 18.0 Inf Inf Inf 0.0 0.0 0.0 0.0 0.0;
          Inf Inf 18.0 Inf Inf 0.0 0.0 0.0 0.0 0.0;
          Inf Inf Inf 18.0 Inf 0.0 0.0 0.0 0.0 0.0;
          Inf Inf Inf Inf 18.0 0.0 0.0 0.0 0.0 0.0]
      W=LapSolve.StiffWrapper(A,2)
      okay=true
      for i in eachindex(M,W)
            M[i]==W[i] || (okay=false; break)
      end
      @test typeof(first(eachindex(M,W))) <: CartesianIndex
      @test okay
end

@testset "simple examples" begin
      @testset "simple examples" begin
          A = [ 0.891171  0.0320582   0.564188  0.8999    0.620615;
                0.166402  0.861136    0.201398  0.911772  0.0796335;
                0.77272   0.782759    0.905982  0.800239  0.297333;
                0.561423  0.170607    0.615941  0.960503  0.981906;
                0.748248  0.00799335  0.554215  0.745299  0.42637]

          assign, cost = @inferred solve_lap(A)
          @test assign == [2, 3, 5, 1, 4]

          B = [ 24     1     8;
                 5     7    14;
                 6    13    20;
                12    19    21;
                18    25     2]

          assign, cost = solve_lap(B)
          @test assign == [2, 1, 0, 0, 3]
          @test cost == 8

          assign, cost = solve_lap(B')
          @test assign == [2, 1, 5]
          @test cost == 8

          assign, cost = solve_lap([ifelse(i-j==0,0,i*j) for i in 1:5, j in 1:5])
          @test assign == [1, 2, 3, 4, 5]
          @test cost == 0
      end
end

@testset "UInt16" begin
    M=UInt16[28092 44837 19882 39481 59139;
             26039 46258 38932 51057 9;
             11527 59487 61993 29072 8734;
             10691 16977 12796 16370 14266;
             5199  42319 34194 41332 16472]
    assign,cost=solve_lap(M)
    @test assign == [3, 5, 4, 2, 1]
    @test cost   == 71139
end

@testset "UInt8" begin
    M=UInt8[67  228 135 197 244;
            112 44  84  206 31;
            225 103 231 225 227;
            170 37  135 9   130;
            110 22  133 77  96]
    assign,cost=solve_lap(M)
    @test assign == [1, 5, 2, 4, 3]
    @test cost   == 343
end

@testset "stiff Problem" begin
    M=[0.08886586609282232 0.5001403306848344;
       0.6848622386100172 0.1379247179005112;
       0.9410463734407415 0.30774601623574327]
    A=[Inf Inf 1000.0;
     0.5266269523801999 Inf Inf;
      0.9238847081919377 Inf 0.6017336762941108]
    @test solve_stiff_lap(M)[1:3]==[(1, 1), (2, 2), (3, -1)]
    @test solve_stiff_lap(A)[1:3]== [(1, -1), (2, 1), (3, 3)]
    @test_throws ErrorException solve_lap(A)
end

ext: jl

max_stars_repo_name: francescoalemanno/SolveLAP


text:
using RedBlackGraph
using Test

function test_floyd_warshall(T::Type)
	r = red_one(T)

	M = [
	r 2 3 0 0;
	0 r 0 2 0;
	0 0 1 0 0;
	0 0 0 r 0;
	2 0 0 0 1
	]
    W = floyd_warshall_transitive_closure(M)
    W_expected = [
	r 2 3 4 0;
	0 r 0 2 0;
	0 0 1 0 0;
	0 0 0 r 0;
	2 4 5 8 1
	]
	@test M == [
	r 2 3 0 0;
	0 r 0 2 0;
	0 0 1 0 0;
	0 0 0 r 0;
	2 0 0 0 1
	]
	@test W == W_expected

	W = floyd_warshall_transitive_closure!(M)
	@test M == W_expected
	@test W == W_expected

end

@testset "shortestpaths" begin
	test_floyd_warshall(AInt8)
	test_floyd_warshall(AInt16)
	test_floyd_warshall(AInt32)
	test_floyd_warshall(AInt64)
	test_floyd_warshall(AInt128)
end


text:
# ---
# title: 1491. Average Salary Excluding the Minimum and Maximum Salary
# id: problem1491
# author: zhwang
# date: 2022-03-18
# difficulty: Easy
# categories: Array, Sort
# link: <https://leetcode.com/problems/average-salary-excluding-the-minimum-and-maximum-salary/description/>
# hidden: true
# ---
#
# Given an array of **unique** integers `salary` where `salary[i]` is the salary
# of the employee `i`.
#
# Return the average salary of employees excluding the minimum and maximum
# salary.
#
#
#
# **Example 1:**
#
#
#
#     Input: salary = [4000,3000,1000,2000]
#     Output: 2500.00000
#     Explanation: Minimum salary and maximum salary are 1000 and 4000 respectively.
#     Average salary excluding minimum and maximum salary is (2000+3000)/2= 2500
#
#
# **Example 2:**
#
#
#
#     Input: salary = [1000,2000,3000]
#     Output: 2000.00000
#     Explanation: Minimum salary and maximum salary are 1000 and 3000 respectively.
#     Average salary excluding minimum and maximum salary is (2000)/1= 2000
#
#
# **Example 3:**
#
#
#
#     Input: salary = [6000,5000,4000,3000,2000,1000]
#     Output: 3500.00000
#
#
# **Example 4:**
#
#
#
#     Input: salary = [8000,9000,2000,3000,6000,1000]
#     Output: 4750.00000
#
#
#
#
# **Constraints:**
#
#   * `3 <= salary.length <= 100`
#   * `10^3 <= salary[i] <= 10^6`
#   * `salary[i]` is unique.
#   * Answers within `10^-5` of the actual value will be accepted as correct.
#
#
## @lc code=start
using LeetCode

average_salary_excluding_mami(salary::Vector{Int}) = (sum(salary) - minimum(salary) - maximum(salary)) / (length(salary) - 2)
## @lc code=end

## @lc test=start
@testset "1491.average-salary-excluding-the-minimum-and-maximum-salary.jl" begin
    @test average_salary_excluding_mami([4000,3000,1000,2000]) ≈ 2500.00000
    @test average_salary_excluding_mami([1000,2000,3000]) ≈ 2000.00000
    @test average_salary_excluding_mami([6000,5000,4000,3000,2000,1000]) ≈ 3500.00000
    @test average_salary_excluding_mami([8000,9000,2000,3000,6000,1000]) ≈ 4750.00000
end
## @lc test=end

## @lc info=start
# link: [solution to 1491]()
## @lc info=end


text:
using BenchmarkTools, PyPlot
using Base.Math.@horner, ParPoly

bench(f::Function, X) = mapreduce(x -> reinterpret(Unsigned,x), |, f(x) for x in X)
X = randn(1_000_000)

nn = 5:25
time_h = Float64[]
time_hs = Float64[]
time_hss = Float64[]
for n in nn
    println("running bench for order $n")
    c = randn(Float64,n+1)

    @gensym f1
    @gensym f2
    @gensym f3
    @eval begin
        $f1(x) = @horner x $(c...)
        $f2(x) = @horner_split x $(c...)
        $f3(x) = @horner_split_simd x $(c...)
    end
    @eval begin
        $f1(2.); $f2($2.); $f3($2.) # precompile them
    end
    t_h    = @benchmark bench($f1, $X)
    t_hs   = @benchmark bench($f2, $X)
    t_hss  = @benchmark bench($f3, $X)

    push!(time_h,   mean(t_h).time)
    push!(time_hs,  mean(t_hs).time)
    push!(time_hss, mean(t_hss).time)
end

plot(nn, time_h*0.000001,   label="@horner")
plot(nn, time_hs*0.000001,  label="@horner_split")
plot(nn, time_hss*0.000001, label="@horner_split_simd")

legend(loc="upper left")
xlabel("polynomial order")
ylabel("time (ms)")
savefig("out.jpg")


text:
using Revise

using IndividualRewrite

using Catlab.Present, Catlab.CategoricalAlgebra
using Plots, Random

Random.seed!(5);

# Define structure of a simulation state
@present ThSIR(FreeSchema) begin
    (S,I,R)::Ob
end
@acset_type SIR(ThSIR)

# Helper function for plotting SIR over time.
function plot_states(states::Vector{SIR})
    ss, is, rs = [[nparts(s, x) for s in states] for x in [:S,:I,:R]]
    x = 1:length(states)
    plot(x,ss; label="S")
    plot!(x, is; label="I")
    plot!(x, rs; label="R")
end

# ACSets for rewrite rules
I = @acset SIR begin I=1 end
R = @acset SIR begin R=1 end
I2 = @acset SIR begin I=2 end
SI = @acset SIR begin S=1; I=1 end

# Rewrite rule morphisms
L_infect =  homomorphism(I, SI)
R_infect = homomorphism(I, I2)
L_recov = homomorphism(SIR(), I)
R_recov = homomorphism(SIR(), R)

# Rewrite rules
rules = Dict(pairs((
  infect=Rule(L_infect,R_infect, 2, 1),
  recover=Rule(L_recov, R_recov, 1, 2)
))...)

# Initial state
init = @acset SIR begin S=10; I=3; R=1 end

# Run the simulation and visualize
states = sim(init, rules, 50; verbose=true)
plot_states(states)

using Debugger
Debugger.@enter sim(init, rules, 50; verbose=true)


text:
module CopEnt

using DocStringExtensions
using SpecialFunctions: logabsgamma, digamma
using Distances: pairwise, Euclidean, Chebyshev
using StatsBase: tiedrank

"""
$(SIGNATURES)

Compute the emperical copula.
"""
empirical_copula(x::AbstractVector) = tiedrank(x) / length(x)
function empirical_copula(x::AbstractMatrix)
    d, n = size(x)
    c = similar(x)
    for i in 1:d
        c[i,:] = tiedrank(x[i,:]) / n
    end
    return c
end

function _entropy_knn_info(x::AbstractVector, dist)
    pwd = pairwise(dist, x)
    return 1, length(x), pwd
end
function _entropy_knn_info(x::AbstractMatrix, dist)
    d, n = size(x)
    pwd = pairwise(dist, x, dims=2)
    return d, n, pwd
end

_entropy_knn_constant(::Euclidean, d) =
    d * log(π) / 2 - d * log(2) - first(logabsgamma(1 + d / 2))
_entropy_knn_constant(::Chebyshev, d) = 0

"""
$(SIGNATURES)

Estimate the entropy using the Kraskov method [1]

# References

1. Alexander Kraskov, Harald Stögbauer and Peter Grassberger. "Estimating mutual information." Physical review, 2004.
"""
function entropy_knn(
    x::AbstractVecOrMat; k::Int=3, dist::Union{Euclidean, Chebyshev}=Euclidean()
)
    d, n, pwd = _entropy_knn_info(x, dist)
    logd = map(1:n) do i
        log(2 * sort(pwd[i,:])[k+1])
    end |> sum
    return digamma(n) - digamma(k) + _entropy_knn_constant(dist, d) + logd * d / n
end

"""
$(SIGNATURES)

Estimate the copula entropy by

1. Compute the emperical copula (see also [`empirical_copula`](@ref));
2. Estimate the entropy using the Kraskov method (see [`entropy_knn`](@ref) for keyword parameters) [1].

The returned copula entropy is the negative mutual information [2].

# References

1. Alexander Kraskov, Harald Stögbauer and Peter Grassberger. "Estimating mutual information." Physical review, 2004.
2. Jian Ma and Zengqi Sun. "Mutual information is copula entropy." Tsinghua Science & Technology, 2011.
"""
function copula_entropy(x::AbstractVecOrMat; kwargs...)
    c = empirical_copula(x)
    return entropy_knn(c; kwargs...)
end

export empirical_copula, entropy_knn, copula_entropy
export Euclidean, Chebyshev # reexport two distance structs

end # module


text:
using RvSpectMLBase
using Test

#=
 This only tests functions from src/utils that don't depend on data structures.
 test/data_structure.jl should test functions that depend on the data structure.
 Functions using multiple data structures should be tested in the more complex data structure's file.
=#

@testset "File Utilities" begin
    @test_nowarn read_data_paths(paths_to_search=[pwd(),joinpath(pkgdir(RvSpectMLBase),"test")])
    @test_nowarn code_to_include_param_jl(paths_to_search=[pwd(),joinpath(pkgdir(RvSpectMLBase),"test")])
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
using Distances, StatsBase, ProgressMeter
export pairwise_inbounds, pairwise_inbounds!, progress_pairwise

"""
In-place distance matrix calculation between elements of Vectors.
"""
function Distances.pairwise!(
    A::AbstractArray,
    metric::Metric,
    a::Vector{T} where T,
    b::Vector{T} where T
    )

    for j in 1:length(b)
        for i in 1:length(a)
            A[i,j] = metric(a[i],b[j])
        end
    end
end

function pairwise_inbounds!(
    A::AbstractArray,
    metric::Metric,
    a::Vector{T},
    b::Vector{T}
    ) where {T}
    @inbounds begin
        for j in eachindex(b)
            for i in eachindex(a)
                A[i,j] = metric(a[i],b[j])
            end
        end
    end
end

"""
In-place distance matrix calculation between elements of Vectors. (For SubArray)
"""
function Distances.pairwise!(
    A::SubArray,
    metric::Metric,
    a::Vector{T} where T,
    b::Vector{T} where T
    )

    for j in 1:length(b)
        for i in 1:length(a)
            A[i,j] = metric(a[i],b[j])
        end
    end
end

# function pairwise_inbounds!(
#     A::SubArray,
#     metric::Metric,
#     a::Vector{T},
#     b::Vector{T}
#     ) where {T}
#     @inbounds begin
#         for j in eachindex(b)
#             for i in eachindex(a)
#                 A[i,j] = metric(a[i],b[j])
#             end
#         end
#     end
# end

"""
In-place distance matrix calculation between elements of Vectors. (For SubArray)
"""
function Distances.pairwise!(
    A::AbstractMatrix,
    metric::Metric,
    a::Vector{T} where T,
    b::Vector{T} where T
    )

    for j in 1:length(b)
        for i in 1:length(a)
            A[i,j] = metric(a[i],b[j])
        end
    end
end

"""
Distance matrix calculation between elements of Vectors. This is a custom extension
of the function in the Distances.jl package to allow vectors of general type. The
function in Distances.jl is designed for univariate/multivariate data and so takes
as input either vectors or matrices (data points as rows).
"""
function Distances.pairwise(
    metric::SemiMetric,
    a::Vector{T},
    b::Vector{T}
    ) where {T}
    D = Array{Float64,2}(undef, length(a), length(b))
    for j in 1:length(b)
        for i in 1:length(a)
            D[i,j] = metric(a[i],b[j])
        end
    end
    return D
end

function pairwise_inbounds(
    metric::SemiMetric,
    a::Vector{T},
    b::Vector{T}
    ) where {T}
    D = Array{Float64,2}(undef, length(a), length(b))
    @inbounds begin
        for j in eachindex(b)
            for i in eachindex(a)
                D[i,j] = metric(a[i],b[j])
            end
        end
    end
    return D
end

"""
Distance matrix calculation between elements of Vectors with progress bar.
"""
function progress_pairwise(
    d::SemiMetric,
    a::Vector{T}
    ) where {T}

    D = zeros(length(a), length(a))
    iter = Progress(Int(length(a)*(length(a)-1)/2), 1)
    for j in 1:length(a)
        for i in 1:(j-1)
            D[i,j] = d(a[i],a[j])
            next!(iter)
        end
    end
    D += D'
    return D

end


text:
# ---
# title: 991. Broken Calculator
# id: problem991
# author: Indigo
# date: 2021-05-31
# difficulty: Medium
# categories: Math, Greedy
# link: <https://leetcode.com/problems/broken-calculator/description/>
# hidden: true
# ---
#
# On a broken calculator that has a number showing on its display, we can
# perform two operations:
#
#   * **Double** : Multiply the number on the display by 2, or;
#   * **Decrement** : Subtract 1 from the number on the display.
#
# Initially, the calculator is displaying the number `X`.
#
# Return the minimum number of operations needed to display the number `Y`.
#
#
#
# **Example 1:**
#
#
#
#     Input: X = 2, Y = 3
#     Output: 2
#     Explanation: Use double operation and then decrement operation {2 -> 4 -> 3}.
#
#
# **Example 2:**
#
#
#
#     Input: X = 5, Y = 8
#     Output: 2
#     Explanation: Use decrement and then double {5 -> 4 -> 8}.
#
#
# **Example 3:**
#
#
#
#     Input: X = 3, Y = 10
#     Output: 3
#     Explanation: Use double, decrement and double {3 -> 6 -> 5 -> 10}.
#
#
# **Example 4:**
#
#
#
#     Input: X = 1024, Y = 1
#     Output: 1023
#     Explanation: Use decrement operations 1023 times.
#
#
#
#
# **Note:**
#
#   1. `1 <= X <= 10^9`
#   2. `1 <= Y <= 10^9`
#
#
## @lc code=start
using LeetCode

function broken_calc(x::Integer, y::Integer)
    cnt = 0
    while x < y
        isodd(y) ? (y += 1) : y >>= 1
        cnt += 1
    end
    return cnt + x - y
end
## @lc code=end


text:
# watch the directory cheese saves webcam images in and classify the
# new images as they appear.

using Knet,MAT
isdefined(:VGG) || include(Knet.dir("examples/vgg.jl"))
using VGG: weights, data, predict
const vggpath = Knet.dir("data","imagenet-vgg-verydeep-16.mat")
const imgpath = "/mnt/ai/home/dyuret/.gnome2/cheese/media"
isdefined(:vgg) || (vgg = matread(vggpath))
model = weights(vgg["layers"])
averageImage = convert(Array{Float32},vgg["meta"]["normalization"]["averageImage"])
description = vgg["meta"]["classes"]["description"]

function vg()
    imgdir1 = readdir(imgpath)
    while true
        sleep(2)
        imgdir2 = readdir(imgpath); print(".")
        if !isequal(imgdir1, imgdir2)
            for img in setdiff(imgdir2, imgdir1)
                x1 = data(joinpath(imgpath,img), averageImage)
                y1 = predict(model, x1)
                z1 = vec(Array(y1))
                s1 = sortperm(z1,rev=true)
                p1 = exp(logp(z1))
                println(img)
                display(hcat(p1[s1[1:5]], description[s1[1:5]]))
                println()
            end
            imgdir1 = imgdir2
        end
    end
end

# x1 = data(o[:image], averageImage)
# info("Classifying")
# @time y1 = predict(w,x1)
# z1 = vec(Array(y1))
# s1 = sortperm(z1,rev=true)
# p1 = exp(logp(z1))
# display(hcat(p1[s1[1:o[:top]]], description[s1[1:o[:top]]]))
# println()

max_stars_repo_name: JuliaTagBot/play

max_issues_repo_name: JuliaTagBot/play

max_forks_repo_name: JuliaTagBot/play


text:
"""
    findpeaks(signal, threshold, distance)
Find all the local maxima of `signal` above `threshold` and at least `distance` apart.
The maxima closer than `distance` are removed starting from the ones with
higher neighbour maxima up until the condition is met.

Return two vectors one with indices and one with heights of the peaks.

The implementation is taken from [`scipy.signal`](https://github.com/scipy/scipy/blob/\
8a64c938ddf1ae4c02a08d2c5e38daeb8d061d38/scipy/signal/_peak_finding.py#L723-L1003).

# Examples

```jldoctest
julia> findpeaks(abs.(sin.(0:0.01pi:2pi)), 0.5, 0)
([51, 151], [1.0, 1.0])
```
"""
function findpeaks(signal::AbstractVector, threshold::Number, distance::Number)
    start = firstindex(signal) + 1
    stop = lastindex(signal) - 1
    peaks = Vector{eltype(axes(signal, 1))}(undef, 0)
    heights = Vector{eltype(signal)}(undef, 0)
    for n in start:stop
        if signal[n-1] < signal[n]
            height = signal[n]
            n_ahead = n + 1
            while height == signal[n_ahead] && n_ahead < stop
                n_ahead += 1
            end
            if signal[n_ahead] < height
                midpoint = div(n + n_ahead, 2)
                push!(peaks, midpoint)
                push!(heights, height)
            end
        end
    end

    if threshold > 0 || distance > 0
        tokeep = trues(length(peaks))
        if threshold > 0
            tokeep .&= selectbypeaksheight(heights, threshold)
        end
        if distance > 0
            tokeep .&= selectbypeaksdistance(peaks, heights, distance)
        end
        peaks[tokeep], heights[tokeep]
    else
        peaks, heights
    end
end

function selectbypeaksheight(heights, threshold)
    heights .> threshold
end

function selectbypeaksdistance(peaks, heights, distance)
    npeaks = length(peaks)
    tokeep = fill(true, npeaks)
    for j in sortperm(heights, rev=true)

        if !tokeep[j]
            continue
        end

        k = j - 1
        while k >= 1 && (peaks[j] - peaks[k]) <= distance
            tokeep[k] = false
            k -= 1
        end

        k = j + 1
        while k <= npeaks && (peaks[k] - peaks[j]) <= distance
            tokeep[k] = false
            k += 1
        end
    end

    tokeep
end

"""
    mad_test(xs, [r=3])

Return a vector of booleans whose elements are `true` if the absolute deviation
from the median of `xs` of the corresponding elements (of `xs`) is at least
`r` times the median absolute deviation, and `false` otherwise.

# Examples
```jldoctest
julia> TemplateMatching.mad_test([1, 2, 3, 100])
4-element BitVector:
 0
 0
 0
 1
```
"""
function mad_test(xs::AbstractVector{<:Number}, r=3.0)
    if isempty(xs)
        BitVector()
    else
        deviations = abs.(xs .- median(xs))
        deviations .> r * median(deviations)
    end
end

function series_magnitude(series, waveform, index)
    series_amp = maximum(abs.(view(series, index:index + length(waveform) - 1)))
    template_amp = maximum(abs.(waveform))
    if iszero(template_amp)
        NaN
    elseif iszero(series_amp)
        NaN
    else
        log10(series_amp / template_amp)
    end
end

"""
    magnitude(data, template, indices, [outlier=TemplateMatching.mad], [avg=StatsBase.mean])

Return the average relative magnitude of the view of series in `data` starting at `indices` and
having the same length of the corresponding `template` series. The outliers, i.e. the elements
corresponding to those of `outlier(magnitudes)` evaluates that are true, are excluded. The average
is calculated using `avg(magnitudes)`.
"""
function magnitude(data, template, indices, outlier=mad_test, avg=mean)
    magnitudes = series_magnitude.(data, template, indices)
    avg(magnitudes[.!outlier(magnitudes)])
end


text:
############### Kaggle FER Competition ###################
# See https://www.kaggle.com/c/challenges-in-representation-learning-facial-expression-recognition-challenge/ for details

immutable KaggleFERDataset <: FaceDataset
end

function load_images(::Type{KaggleFERDataset}, datadir::AbstractString)
    datafile = joinpath(datadir, "fer2013.csv")
    if !isfile(datafile)
        throw(ArgumentError("Expected to have CSV file at $datafile"))
    end
    n_images = 35887
    imgs = Vector{Matrix{Float64}}(n_images)
    open(datafile) do f
        readline(f) # skip header
        for (i, line) in zip(1:n_images, eachline(f))
            label, data, _ = split(line, ",")
            arr = [parse(Float64, s) for s in split(data, " ")] / 255
            img = reshape(arr, 48, 48)'
            imgs[i] = img
            if i % 1000 == 0
                info("Loaded $i images")
            end
        end
    end
    return imgs
end

function load_labels(::Type{KaggleFERDataset}, datadir::AbstractString)
    datafile = joinpath(datadir, "fer2013.csv")
    if !isfile(datafile)
        throw(ArgumentError("Expected to have CSV file at $datafile"))
    end
    n_images = 35887
    labels = Vector{Float64}(n_images)
    open(datafile) do f
        readline(f) # skip header
        for (i, line) in zip(1:n_images, eachline(f))
            label, data, _ = split(line, ",")
            labels[i] = parse(Float64, label)
        end
    end
    return labels
end


text:
using Test

using ClimateMachine.SurfaceFluxes
const SF = SurfaceFluxes

using CLIMAParameters: AbstractEarthParameterSet
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

@testset "SurfaceFluxes - FMS Profiles" begin
    FT = Float32

    ## Discretisation altitude z
    z = FT[
        29.432779269303,
        30.0497139076724,
        31.6880000418153,
        34.1873479240475,
    ]
    ## Virtual Dry Static Energy at height z
    pt = FT[
        268.559120403867,
        269.799228886728,
        277.443023238556,
        295.79192777341,
    ]
    ## Surface Pottemp
    pt0 = FT[
        273.42369841804,
        272.551410044203,
        278.638168565727,
        298.133068766049,
    ]
    ## Roughness lengths
    z0 = FT[
        5.86144925739178e-05,
        0.0001,
        0.000641655193293549,
        3.23383768877187e-05,
    ]
    zt = FT[
        3.69403636275411e-05,
        0.0001,
        1.01735489109205e-05,
        7.63933834969505e-05,
    ]
    zq = FT[
        5.72575636226887e-05,
        0.0001,
        5.72575636226887e-05,
        5.72575636226887e-05,
    ]
    ## Speed
    speed =
        FT[2.9693638452068, 2.43308757772094, 5.69418282305367, 9.5608693754561]
    ## Scale velocity and moisture
    u_star = FT[
        0.109462510724615,
        0.0932942802513508,
        0.223232887323184,
        0.290918439028557,
    ]
    q_star = FT[
        0.000110861442197537,
        9.44983279664197e-05,
        4.17643828631936e-05,
        0.000133135421415819,
    ]
    # No explicit buoyancy terms in ClimateMachine
    #b_star =  [0.00690834676781433, 0.00428178089592372, 0.00121229800895103, 0.00262353784027441, -0.000570314880866852]

    for ii in 1:length(u_star)
        # Data at first interior node (x_ave)
        qt_ave = FT(0)
        z_ave = z[ii]
        vdse_ave = pt[ii]
        u_ave = speed[ii]
        x_ave = FT[u_ave, vdse_ave, qt_ave]

        ## Initial guesses for MO parameters
        LMO_init = eps(FT)
        u_star_init = FT(0.1)
        th_star_init = -FT(0.1)
        qt_star_init = -FT(1e-5)
        x_init = FT[LMO_init, u_star_init, th_star_init, qt_star_init]

        # Surface values for variables
        u_sfc = FT(0)
        thv_sfc = pt0[ii]
        qt_sfc = q_star[ii]
        z_sfc = FT(0)
        x_s = FT[u_sfc, thv_sfc, qt_sfc]

        # Roughness
        z_rough = FT[z0[ii], zt[ii], zq[ii]]

        # Constants
        a = FT(4.7)
        Δz = z[ii]

        # F_exchange
        F_exchange = [FT(0.01), -FT(0.01), -FT(0.000001)]

        args = (
            param_set,
            x_init,
            x_ave,
            x_s,
            z_rough,
            F_exchange,
            vdse_ave,
            qt_ave,
            Δz,
            z_ave / 2,
            a,
        )

        ## Assuming surface fluxes are not given
        result = surface_conditions(args..., SF.DGScheme())
        x_star = result.x_star
        @test (abs((x_star[1] - u_star[ii]) / u_star[ii])) <= FT(0.15)

        # Keeping FV-based scheme in case needed
        result = surface_conditions(args..., SF.FVScheme())
    end
end

include("test_universal_functions.jl")

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
export Selector, Class, Object
export @sel_str
export classexists

# import Base: show, convert, unsafe_convert, super, methods
import Base: show, unsafe_convert

int2bool(x::Integer) = x != 0

## Selector ###############################################

struct Selector
  ptr::Ptr{Cvoid}
  # Selector(ptr::Ptr{CVoid}) = new(ptr)
end

unsafe_convert(::Type{Ptr{Cvoid}}, sel::Selector) = sel.ptr

"Turn a selector name into a selector"
function Selector(name::Union{AbstractString, Symbol})
  Selector(ccall(:sel_registerName, Ptr{Cvoid}, (Ptr{Cchar},),
                 pointer(string(name))))
end

macro sel_str(name)
  Selector(name)
end

selname(s::Ptr{Cvoid}) =
  ccall(:sel_getName, Ptr{Cchar}, (Ptr{Cvoid},),
        s) |> unsafe_string

name(sel::Selector) = selname(sel.ptr)

function show(io::IO, sel::Selector)
  print(io, "sel")
  show(io, string(name(sel)))
end

## Class ###############################################

struct Class
  ptr::Ptr{Cvoid}
  # Class(ptr::Ptr{Void}) = new(ptr)
end

"Get class pointer"
classptr(name) = ccall(:objc_getClass, Ptr{Cvoid}, (Ptr{Cchar},),
                       pointer(string(name)))

unsafe_convert(::Type{Ptr{Cvoid}}, class::Class) = class.ptr

function Class(name::Union{AbstractString, Symbol})
  ptr = classptr(name)
  ptr == C_NULL && error("Couldn't find class $name")
  return Class(ptr)
end

"""
    classexists(name::Union{AbstractString, Symbol}) -> Bool

Check if class named `name` is already registered with the system. Useful to
do before registering a new custom class.
"""
function classexists(name::Union{AbstractString, Symbol})
    classptr(name) ≠ C_NULL
end

function name(class::Class)
  cname = ccall(:class_getName, Ptr{Cchar}, (Ptr{Cvoid},), class)
  @assert cname != nil
  Symbol(unsafe_string(cname))
end

function ismeta(class::Class)
  ismeta = ccall(:class_isMetaClass, Cint, (Ptr{Cvoid},), class)
  int2bool(ismeta)
end

function show(io::IO, class::Class)
  ismeta(class) && print(io, "^")
  print(io, name(class))
end

## Object ###############################################
function release end

mutable struct Object
  ptr::Ptr{Cvoid}
  # function Object(ptr)
  #     obj = new(ptr)
  #     function f(t)
  #         # msg = string("About to release ", repr(UInt(t.ptr)))
  #         # TODO: Remove println, just for debugging. Call release directly
  #         @async begin
  #             println("Releasing $t")
  #             release(t)
  #         end
  #     end
  #     finalizer(f, obj)
  # end
end

unsafe_convert(::Type{Ptr{Cvoid}}, obj::Object) = obj.ptr

class(obj) =
  ccall(:object_getClass, Ptr{Cvoid}, (Ptr{Cvoid},),
        obj) |> Class

methods(obj::Object) = methods(class(obj))

# function show(io::IO, obj::Object)
#     addr = repr(UInt(obj.ptr))
#     println(io, class(obj), " Object $addr")
# end


text:
#list of added packages

#routines for adding the packages
add Atom Juno Revise
add Flux CuArrays CUDAapi CUDAnative GPUArrays
add DataFrames Revise Distributions StatsBase GLM Dates Formatting
add DataStructures LinearAlgebra Measures Formatting XLSX
add CSV CategoricalArrays Serialization Distributed Random
add StaticArrays VegaLite CodecZlib CodecLz4 CodecZstd
add BenchmarkTools StatsModels
add LeastSquaresOptim UnPack BlackBoxOptim
add https://github.com/clintonTE/Finometrics

(Temporary)
add Flux#master CuArrays#master CUDAapi#master GPUArrays#master

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#using Convex
#using SCS

using Printf
abstract type Solver end

struct QRSolver <: Solver end
struct DefaultSolver <: Solver end
struct CvxSolver <: Solver
    verbose
    eps
end

mutable struct ProxGradientSolver <: Solver
    verbose
    eps
    max_iters
    min_iters
    gamma_initial
    thetas
    gammas
    gradfs
    fs
    gs
    progress
    stochastic
    storedata
end

name(S::ProxGradientSolver) = "ProxGradientSolver"
name(S::CvxSolver) = "CvxSolver"
name(S::QRSolver) = "QRSolver"

import Base.print
print(io, S::Solver) = print(io, name(S))
print(io::IO, S::Solver) = print(io, name(S))

ProxGradientSolver(;verbose=false, eps=1e-6,
                   maxiters=1000, miniters=1,
                   gamma_initial=0.1,
                   progress=nothing, stochastic=false,
                   storedata = true) = ProxGradientSolver(verbose, eps, maxiters, miniters,
                                                          gamma_initial, nothing, nothing,
                                                          nothing, nothing, nothing, progress,
                                                          stochastic, storedata)

CvxSolver(;verbose=false, eps=1e-5) = CvxSolver(verbose, eps)

# we query for the solver when the user has
# asked for the default. We don't use Julia dispatch on type Solver
# since we want to use one solver for the entire
# regularization path
getsolver(L::SquareLoss, R::L2Reg) = QRSolver()
getsolver(L::LossDiff, R) = ProxGradientSolver()
getsolver(L::LossNonDiff, R) = CvxSolver()

# non convex losses or regularizers
getsolver(L::LogHuberLoss, R) = ProxGradientSolver()
getsolver(L::LossDiff, R::SqrtReg) = ProxGradientSolver()

##############################################################################
# CVX for huber

function solve(S::CvxSolver, L, R, regweights, X, Y, lambda;
                  theta_guess=nothing)
    n,d = size(X)
    m = size(Y,2)
    theta = Convex.Variable(d,m)
    # we shouldn't need to specify Positive here, but
    # Convex.jl is unreliable without it
    losses = Convex.Variable(n, Positive())
    regs = Convex.Variable(d, m, Positive())
    problem = Convex.minimize( sum(losses)/n  + sum(regs))
    for i=1:n
        problem.constraints +=    losses[i] >= cvxloss(L, theta'*X[i,:],  Y[i,:])
    end
    if isa(R, NonnegReg)
        for i=1:d
            if regweights[i] > 0
                problem.constraints +=    theta[i] >= 0
            end
            problem.constraints +=    regs[i] == 0
        end
    else
        for i=1:d
            for j=1:m
                problem.constraints +=    regs[i]   >= lambda * regweights[i] * cvxreg(R, theta[i,j])
            end
        end
    end
    solve!(problem, SCSSolver(verbose=S.verbose, eps=S.eps))
    println("theta.value = ", theta.value)
    return theta.value
end

##############################################################################
# QR for quadratic case

function solve(s::QRSolver, L::SquareLoss, R::L2Reg, regweights,
               X, Y, lambda; theta_guess=nothing)
    return lsreg(X, Y, lambda, regweights)
end

# inputs:
#    X        n by d
#    y        n by 1
#
# returns:
#    theta:   d by 1 vector
#
function ls(X,Y)
    F = qr(X)
    Q = Matrix(F.Q)
    R = Matrix(F.R)
    theta = R\(Q'*Y)
    return theta
end

function lsreg(X, Y, lmb, regweights)
    d = size(X,2)
    if d>5000
        return svdlsreg(X,Y,lmb)
    end
    m = size(Y,2)
    A = sqrt(lmb) * LinearAlgebra.diagm(0=>sqrt.(regweights))
    Xh = [X; A]
    Yh = [Y; zeros(d,m)]
    theta = ls(Xh,Yh)
end

# always regularizes all of theta
function svdlsreg(X, Y, lmb)
    U,S,V=svd(X)
    d = S./(S.*S + lmb)
    th = V*(d.*(U'*Y))
    return th
end

##############################################################################

#
# minimize (1/n) * sum loss(theta^T*x^i, y^i) + lambda*\sum_i regweights[i]*reg(theta_i)
#

function solve(S::ProxGradientSolver, L::Loss, R::Regularizer,
               regweights, X, Y, regparam; theta_guess = nothing)

    # assume Y is n by 1
    d = size(X,2)
    n = size(X,1)
    m = size(Y,2)
    Xt = X'

    # useful stuff
    # f = loss(theta)
    # g = sum(regparam*regweights*reg(theta))
    #
    # so we minimize f + g
    #
    # reshape to deal with matrix case
    #
    R1 = repmat(regweights, 1,m)
    rw = regparam*R1[:]
    f(theta) = loss(L, X*reshape(theta,d,m) , Y)
    gradf(theta) = derivlosstheta(L, X*reshape(theta,d,m), Y, Xt, S.stochastic)[:]
    g(theta) = reg(R, rw, theta[:])
    proxg(gamma, v) = prox(R, gamma, rw, v)

    tg = theta_guess
    if tg != nothing
        tg = tg[:]
    end

    theta =  proxgradient(d*m, f, gradf, g, proxg, S; theta_guess = tg)
    return reshape(theta,d,m)
end

# generic solver will work for any f and g
# gets parameters from S and stores results in S
function proxgradient(d, f, gradf, g, proxg, S; theta_guess=nothing)

    max_iters = S.max_iters
    min_iters = S.min_iters
    verbose = S.verbose
    eps = S.eps
    gamma_initial = S.gamma_initial

    # buffers
    if S.storedata
        thetas = zeros(d, max_iters)
        gammas = zeros(max_iters)
        gradfs = zeros(d, max_iters)
        fs = zeros(max_iters)
        gs = zeros(max_iters)
    end

    # initial conditions
    if theta_guess != nothing
        theta_initial = theta_guess
    else
        theta_initial = zeros(d)
    end

    # initialize
    theta = theta_initial
    gamma = gamma_initial
    gradfvalue = zeros(d)
    fvalue = f(theta)
    gvalue = g(theta)

    # store first step
    if S.storedata
        thetas[:,1] = theta
        gammas[1] = gamma
        fs[1] = fvalue
        gs[1] = gvalue
    end

    # make inner loop variables accessible
    k=1

    # loop
    for k=1:max_iters-1

        # printing
        if verbose
            @printf("%d   gamma: %f   loss: %f   reg: %f \n", k, gamma, fvalue, gvalue)
        end

        gradfvalue = gradf(theta)

        if S.storedata
            gradfs[:,k] = gradfvalue
        end

        if S.stochastic
            v = theta - gradfs[:,k]/(2*gamma)
            theta_next = proxg(gamma, v)
            f_next = f(theta_next)
            g_next = g(theta_next)
            gamma_next = gamma_initial*sqrt(k+1)

        else
            # line search
            while true
                v = theta - gradfvalue/(2*gamma)
                theta_next = proxg(gamma, v)
                f_next = f(theta_next)
                g_next = g(theta_next)

                # should be <= not < else can get stuck if g is an indicator function
                if f_next + g_next <= fvalue + gvalue
                    # increase the step size and move to next step
                    gamma = gamma/1.2
                    break
                else
                    # decrease the step size and try again
                    gamma = gamma*2
                end
            end
        end

        if k > min_iters && ! S.stochastic && fvalue + gvalue - (f_next + g_next) < eps
            break
        end

        fvalue = f_next
        gvalue = g_next
        theta = theta_next

        # save the variables
        if S.storedata
            thetas[:,k+1] = theta
            gammas[k+1] = gamma
            fs[k+1] = fvalue
            gs[k+1] = gvalue
        end

        if S.progress != nothing
            S.progress(k+1,theta)
        end

    end
    if S.storedata
        S.thetas = thetas[:,1:k+1]
        S.gammas = gammas[1:k+1]
        S.gradfs = gradfs[:,1:k+1]
        S.fs = fs[1:k+1]
        S.gs = gs[1:k+1]
    end

    return theta

end

##############################################################################


text:
"""
based on  https://github.com/Visceral-Project/EvaluateSegmentation/blob/master/source/MutualInformationMetric.h
"""
module InformationTheorhetic
using ..BasicStructs, Parameters,  CUDA

"""
calculate mutual information  based on precalulated constants
"""
function mutualInformationMetr(tn,tp,fp, fn) ::Float64
    H1=0
    H2=0
    H12=0
    if(tn>0)
        # print("tn $(tn) tp $(tp) fp $(fp) fn $(fn)\n ")
        n = tn+tp+fp+fn
        # print("tn $(tn) tp $(tp) fp $(fp) fn $(fn) n $(n) \n ")
        row1 = tn + fn
        row2 = fp + tp
        # print("row1 $(row1) row2 $(row2) row1/n $(row1/n) row2/n $(row2/n)")

        H1 = - ( (row1/n)*log2(row1/n) + (row2/n)*log2(row2/n))

        col1 = tn + fp
        col2 = fn + tp
        H2 = - ( (col1/n)*log2(col1/n) + (col2/n)*log2(col2/n))

        p00::Float32 = tn==0 ? 1 : (tn/n)
        p01::Float32 = fn==0 ? 1 : (fn/n)
        p10::Float32 = fp==0 ? 1 : (fp/n)
        p11::Float32 = tp==0 ? 1 : (tp/n)
        H12= - ( (tn/n)*log2(p00) + (fn/n)*log2(p01) +  (fp/n)*log2(p10) + (tp/n)*log2(p11) )
    end
return H1+H2-H12

end #calculateVolumeMetric

"""
The Variation of Information (VOI) measures the
amount of information lost (or gained) when changing
from one variable to the other
"""
function variationOfInformation(tn,tp,fp, fn)
    H1=0
    H2=0
    Mi=0
    H12=0
    n = tn+tp+fp+fn
    trueVoxels = fn + tp;
    retrievedVoxels = fp + tp;
    if(tn>0)

    H1 = - ( (trueVoxels/n)*log2(trueVoxels/n) + (1- trueVoxels/n)*log2(1- trueVoxels/n));
    H2 = - ( (retrievedVoxels/n)*log2(retrievedVoxels/n) + (1- retrievedVoxels/n)*log2(1- retrievedVoxels/n) );
    p00::Float32 = tn==0 ? 1.0 : (tn/n);
     p01::Float32 = fn==0 ? 1.0 : (fn/n);
     p10::Float32 = fp==0 ? 1.0 : (fp/n);
     p11::Float32 = tp==0 ? 1.0 : (tp/n);
   H12= - ( (tn/n)*log2(p00)
                + (fn/n)*log2(p01) +  (fp/n)*log2(p10) + (tp/n)*log2(p11) )
    end
    MI=H1+H2-H12;

    return H1+H2-2*MI;

end#variationOfInformation

end #InformationTheorhetic


text:
mutable struct Inputs
    word::DubinsPathType
    α::Float64
    β::Float64
    d::Float64
end

mutable struct Outputs
    errcode::Int
    params::Vector{Float64}
    len::Float64
end

mutable struct Param
    inputs::Inputs
    outputs::Outputs
end

function Param(w::DubinsPathType, α, β, d, code, p1, p2, p3, len)
    p = Param(Inputs(w, float(α), float(β), float(d)), Outputs(code, [float(p1), float(p2), float(p3)], len))
    return p
end

params = [
    Param(LSL, -1.57079632679490, -1.57079632679490, 0.200000000000000, 0, 1.57079632679490, 0.200000000000000, 4.71238898038469, 6.48318530717959),
    Param(LSR, -1.57079632679490, -1.57079632679490, 0.200000000000000, 1, 0, 0, 0, 0),
    Param(RSL, -1.57079632679490, -1.57079632679490, 0.200000000000000, 0, 5.85348564102816, 0.916515138991168, 5.85348564102816, 12.6234864210475),
    Param(RSR, -1.57079632679490, -1.57079632679490, 0.200000000000000, 0, 4.71238898038469, 0.200000000000000, 1.57079632679490, 6.48318530717959),
    Param(RLR, -1.57079632679490, -1.57079632679490, 0.200000000000000, 0, 1.52077546998913, 6.18314359356805, 4.66236812357892, 12.3662871871361),
    Param(LRL, -1.57079632679490, -1.57079632679490, 0.200000000000000, 0, 4.66236812357892, 6.18314359356805, 1.52077546998913, 12.3662871871361),
    Param(LSL, -1.57079632679490, -1.57079632679490, 0.600000000000000, 0, 1.57079632679490, 0.600000000000000, 4.71238898038469, 6.88318530717959),
    Param(LSR, -1.57079632679490, -1.57079632679490, 0.600000000000000, 1, 0, 0, 0, 0),
    Param(RSL, -1.57079632679490, -1.57079632679490, 0.600000000000000, 0, 5.59002539960283, 1.66132477258361, 5.59002539960283, 12.8413755717893),
    Param(RSR, -1.57079632679490, -1.57079632679490, 0.600000000000000, 0, 4.71238898038469, 0.600000000000000, 1.57079632679490, 6.88318530717959),
    Param(RLR, -1.57079632679490, -1.57079632679490, 0.600000000000000, 0, 1.42022805401821, 5.98204876162621, 4.56182070760800, 11.9640975232524),
    Param(LRL, -1.57079632679490, -1.57079632679490, 0.600000000000000, 0, 4.56182070760800, 5.98204876162621, 1.42022805401821, 11.9640975232524),
    Param(LSL, -1.57079632679490, -1.57079632679490, 1.20000000000000, 0, 1.57079632679490, 1.20000000000000, 4.71238898038469, 7.48318530717959),
    Param(LSR, -1.57079632679490, -1.57079632679490, 1.20000000000000, 1, 0, 0, 0, 0),
    Param(RSL, -1.57079632679490, -1.57079632679490, 1.20000000000000, 0, 5.38752051332172, 2.49799919935936, 5.38752051332172, 13.2730402260028),
    Param(RSR, -1.57079632679490, -1.57079632679490, 1.20000000000000, 0, 4.71238898038469, 1.20000000000000, 1.57079632679490, 7.48318530717959),
    Param(RLR, -1.57079632679490, -1.57079632679490, 1.20000000000000, 0, 1.26610367277950, 5.67379999914879, 4.40769632636929, 11.3475999982976),
    Param(LRL, -1.57079632679490, -1.57079632679490, 1.20000000000000, 0, 4.40769632636929, 5.67379999914879, 1.26610367277950, 11.3475999982976),
    Param(LSL, -1.57079632679490, -1.57079632679490, 10.0000000000000, 0, 1.57079632679490, 10.0000000000000, 4.71238898038469, 16.2831853071796),
    Param(LSR, -1.57079632679490, -1.57079632679490, 10.0000000000000, 0, 1.82347658193698, 7.74596669241483, 1.82347658193698, 11.3929198562888),
    Param(RSL, -1.57079632679490, -1.57079632679490, 10.0000000000000, 0, 4.87983705960438, 11.8321595661992, 4.87983705960438, 21.5918336854080),
    Param(RSR, -1.57079632679490, -1.57079632679490, 10.0000000000000, 0, 4.71238898038469, 10.0000000000000, 1.57079632679490, 16.2831853071796),
    Param(RLR, -1.57079632679490, -1.57079632679490, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LRL, -1.57079632679490, -1.57079632679490, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LSL, -1.57079632679490, 0, 0.200000000000000, 0, 3.81633359581335, 1.28062484748657, 4.03764803816114, 9.13460648146105),
    Param(LSR, -1.57079632679490, 0, 0.200000000000000, 1, 0, 0, 0, 0),
    Param(RSL, -1.57079632679490, 0, 0.200000000000000, 1, 0, 0, 0, 0),
    Param(RSR, -1.57079632679490, 0, 0.200000000000000, 0, 5.40712725658139, 1.56204993518133, 5.58844703098288, 12.5576242227456),
    Param(RLR, -1.57079632679490, 0, 0.200000000000000, 0, 1.86434638999434, 5.48080888118507, 2.04566616439583, 9.39082143557524),
    Param(LRL, -1.57079632679490, 0, 0.200000000000000, 0, 0.348846568576992, 5.63139655988647, 0.570161010924784, 6.55040413938824),
    Param(LSL, -1.57079632679490, 0, 0.600000000000000, 0, 3.52209903070216, 1.07703296142690, 4.33188260327232, 8.93101459540138),
    Param(LSR, -1.57079632679490, 0, 0.600000000000000, 1, 0, 0, 0, 0),
    Param(RSL, -1.57079632679490, 0, 0.600000000000000, 1, 0, 0, 0, 0),
    Param(RSR, -1.57079632679490, 0, 0.600000000000000, 0, 5.27098829572825, 1.88679622641132, 5.72458599183602, 12.8823705139756),
    Param(RLR, -1.57079632679490, 0, 0.600000000000000, 0, 1.63817896180931, 5.30075194652129, 2.09177665791708, 9.03070756624769),
    Param(LRL, -1.57079632679490, 0, 0.600000000000000, 0, 0.107883633503051, 5.73793981996096, 0.917667206073218, 6.76349065953723),
    Param(LSL, -1.57079632679490, 0, 1.20000000000000, 0, 2.94419709373991, 1.01980390271856, 4.90978454023457, 8.87378553669304),
    Param(LSR, -1.57079632679490, 0, 1.20000000000000, 1, 0, 0, 0, 0),
    Param(RSL, -1.57079632679490, 0, 1.20000000000000, 0, 5.26059043786122, 1.35646599662505, 0.548201457476526, 7.16525789196279),
    Param(RSR, -1.57079632679490, 0, 1.20000000000000, 0, 5.13901647351157, 2.41660919471891, 5.85655781405271, 13.4121834822832),
    Param(RLR, -1.57079632679490, 0, 1.20000000000000, 0, 1.34872217225806, 4.98578201185216, 2.06626351279921, 8.40076769690943),
    Param(LRL, -1.57079632679490, 0, 1.20000000000000, 0, 5.82799274420666, 5.76759130093350, 1.51039488352173, 13.1059789286619),
    Param(LSL, -1.57079632679490, 0, 10.0000000000000, 0, 1.68145354796879, 9.05538513813742, 6.17252808600569, 16.9093667721119),
    Param(LSR, -1.57079632679490, 0, 10.0000000000000, 0, 1.68283839261674, 8.83176086632785, 0.112042065821841, 10.6266413247664),
    Param(RSL, -1.57079632679490, 0, 10.0000000000000, 0, 4.80380493701101, 10.8627804912002, 0.0914159566263198, 15.7580013848375),
    Param(RSR, -1.57079632679490, 0, 10.0000000000000, 0, 4.80304886758544, 11.0453610171873, 6.19252541997884, 22.0409353047515),
    Param(RLR, -1.57079632679490, 0, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LRL, -1.57079632679490, 0, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LSL, -1.57079632679490, 1.57079632679490, 0.200000000000000, 0, 4.71238898038469, 1.80000000000000, 4.71238898038469, 11.2247779607694),
    Param(LSR, -1.57079632679490, 1.57079632679490, 0.200000000000000, 1, 0, 0, 0, 0),
    Param(RSL, -1.57079632679490, 1.57079632679490, 0.200000000000000, 1, 0, 0, 0, 0),
    Param(RSR, -1.57079632679490, 1.57079632679490, 0.200000000000000, 0, 4.71238898038469, 2.20000000000000, 4.71238898038469, 11.6247779607694),
    Param(RLR, -1.57079632679490, 1.57079632679490, 0.200000000000000, 0, 0.988432088926153, 5.11845683144210, 0.988432088926153, 7.09532100929441),
    Param(LRL, -1.57079632679490, 1.57079632679490, 0.200000000000000, 0, 1.10403098774760, 5.34965462908499, 1.10403098774760, 7.55771660458019),
    Param(LSL, -1.57079632679490, 1.57079632679490, 0.600000000000000, 0, 4.71238898038469, 1.40000000000000, 4.71238898038469, 10.8247779607694),
    Param(LSR, -1.57079632679490, 1.57079632679490, 0.600000000000000, 1, 0, 0, 0, 0),
    Param(RSL, -1.57079632679490, 1.57079632679490, 0.600000000000000, 1, 0, 0, 0, 0),
    Param(RSR, -1.57079632679490, 1.57079632679490, 0.600000000000000, 0, 4.71238898038469, 2.60000000000000, 4.71238898038469, 12.0247779607694),
    Param(RLR, -1.57079632679490, 1.57079632679490, 0.600000000000000, 0, 0.863211890069541, 4.86801643372887, 0.863211890069541, 6.59444021386796),
    Param(LRL, -1.57079632679490, 1.57079632679490, 0.600000000000000, 0, 1.21322522314939, 5.56804309988857, 1.21322522314939, 7.99449354618734),
    Param(LSL, -1.57079632679490, 1.57079632679490, 1.20000000000000, 0, 4.71238898038469, 0.800000000000000, 4.71238898038469, 10.2247779607694),
    Param(LSR, -1.57079632679490, 1.57079632679490, 1.20000000000000, 1, 0, 0, 0, 0),
    Param(RSL, -1.57079632679490, 1.57079632679490, 1.20000000000000, 1, 0, 0, 0, 0),
    Param(RSR, -1.57079632679490, 1.57079632679490, 1.20000000000000, 0, 4.71238898038469, 3.20000000000000, 4.71238898038469, 12.6247779607694),
    Param(RLR, -1.57079632679490, 1.57079632679490, 1.20000000000000, 0, 0.643501108793284, 4.42859487117636, 0.643501108793284, 5.71559708876293),
    Param(LRL, -1.57079632679490, 1.57079632679490, 1.20000000000000, 0, 1.36943840600457, 5.88046946559892, 1.36943840600457, 8.61934627760806),
    Param(LSL, -1.57079632679490, 1.57079632679490, 10.0000000000000, 0, 1.57079632679490, 8.00000000000000, 1.57079632679490, 11.1415926535898),
    Param(LSR, -1.57079632679490, 1.57079632679490, 10.0000000000000, 0, 1.77215424758523, 9.79795897113271, 4.91374690117502, 16.4838601198930),
    Param(RSL, -1.57079632679490, 1.57079632679490, 10.0000000000000, 0, 4.91374690117502, 9.79795897113271, 1.77215424758523, 16.4838601198930),
    Param(RSR, -1.57079632679490, 1.57079632679490, 10.0000000000000, 0, 4.71238898038469, 12.0000000000000, 4.71238898038469, 21.4247779607694),
    Param(RLR, -1.57079632679490, 1.57079632679490, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LRL, -1.57079632679490, 1.57079632679490, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LSL, 0, -1.57079632679490, 0.200000000000000, 0, 5.58844703098288, 1.56204993518133, 5.40712725658139, 12.5576242227456),
    Param(LSR, 0, -1.57079632679490, 0.200000000000000, 1, 0, 0, 0, 0),
    Param(RSL, 0, -1.57079632679490, 0.200000000000000, 1, 0, 0, 0, 0),
    Param(RSR, 0, -1.57079632679490, 0.200000000000000, 0, 4.03764803816114, 1.28062484748657, 3.81633359581335, 9.13460648146105),
    Param(RLR, 0, -1.57079632679490, 0.200000000000000, 0, 0.570161010924784, 5.63139655988647, 0.348846568576992, 6.55040413938824),
    Param(LRL, 0, -1.57079632679490, 0.200000000000000, 0, 2.04566616439583, 5.48080888118507, 1.86434638999434, 9.39082143557524),
    Param(LSL, 0, -1.57079632679490, 0.600000000000000, 0, 5.72458599183602, 1.88679622641132, 5.27098829572825, 12.8823705139756),
    Param(LSR, 0, -1.57079632679490, 0.600000000000000, 1, 0, 0, 0, 0),
    Param(RSL, 0, -1.57079632679490, 0.600000000000000, 1, 0, 0, 0, 0),
    Param(RSR, 0, -1.57079632679490, 0.600000000000000, 0, 4.33188260327232, 1.07703296142690, 3.52209903070216, 8.93101459540138),
    Param(RLR, 0, -1.57079632679490, 0.600000000000000, 0, 0.917667206073218, 5.73793981996096, 0.107883633503051, 6.76349065953723),
    Param(LRL, 0, -1.57079632679490, 0.600000000000000, 0, 2.09177665791708, 5.30075194652129, 1.63817896180931, 9.03070756624769),
    Param(LSL, 0, -1.57079632679490, 1.20000000000000, 0, 5.85655781405271, 2.41660919471891, 5.13901647351157, 13.4121834822832),
    Param(LSR, 0, -1.57079632679490, 1.20000000000000, 1, 0, 0, 0, 0),
    Param(RSL, 0, -1.57079632679490, 1.20000000000000, 0, 0.548201457476526, 1.35646599662505, 5.26059043786122, 7.16525789196279),
    Param(RSR, 0, -1.57079632679490, 1.20000000000000, 0, 4.90978454023457, 1.01980390271856, 2.94419709373991, 8.87378553669304),
    Param(RLR, 0, -1.57079632679490, 1.20000000000000, 0, 1.51039488352173, 5.76759130093350, 5.82799274420666, 13.1059789286619),
    Param(LRL, 0, -1.57079632679490, 1.20000000000000, 0, 2.06626351279921, 4.98578201185216, 1.34872217225806, 8.40076769690943),
    Param(LSL, 0, -1.57079632679490, 10.0000000000000, 0, 6.19252541997884, 11.0453610171873, 4.80304886758544, 22.0409353047515),
    Param(LSR, 0, -1.57079632679490, 10.0000000000000, 0, 0.112042065821841, 8.83176086632785, 1.68283839261674, 10.6266413247664),
    Param(RSL, 0, -1.57079632679490, 10.0000000000000, 0, 0.0914159566263198, 10.8627804912002, 4.80380493701101, 15.7580013848375),
    Param(RSR, 0, -1.57079632679490, 10.0000000000000, 0, 6.17252808600569, 9.05538513813742, 1.68145354796879, 16.9093667721119),
    Param(RLR, 0, -1.57079632679490, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LRL, 0, -1.57079632679490, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LSL, 0, 0, 0.200000000000000, 0, 0, 0.200000000000000, 0, 0.200000000000000),
    Param(LSR, 0, 0, 0.200000000000000, 0, 0, 0.200000000000000, 0, 0.200000000000000),
    Param(RSL, 0, 0, 0.200000000000000, 0, 0, 0.200000000000000, 0, 0.200000000000000),
    Param(RSR, 0, 0, 0.200000000000000, 0, 0, 0.200000000000000, 0, 0.200000000000000),
    Param(RLR, 0, 0, 0.200000000000000, 0, 3.09157179678402, 6.18314359356805, 3.09157179678402, 12.3662871871361),
    Param(LRL, 0, 0, 0.200000000000000, 0, 3.09157179678402, 6.18314359356805, 3.09157179678402, 12.3662871871361),
    Param(LSL, 0, 0, 0.600000000000000, 0, 0, 0.600000000000000, 0, 0.600000000000000),
    Param(LSR, 0, 0, 0.600000000000000, 0, 0, 0.600000000000000, 0, 0.600000000000000),
    Param(RSL, 0, 0, 0.600000000000000, 0, 0, 0.600000000000000, 0, 0.600000000000000),
    Param(RSR, 0, 0, 0.600000000000000, 0, 0, 0.600000000000000, 0, 0.600000000000000),
    Param(RLR, 0, 0, 0.600000000000000, 0, 2.99102438081311, 5.98204876162621, 2.99102438081311, 11.9640975232524),
    Param(LRL, 0, 0, 0.600000000000000, 0, 2.99102438081311, 5.98204876162621, 2.99102438081311, 11.9640975232524),
    Param(LSL, 0, 0, 1.20000000000000, 0, 0, 1.20000000000000, 0, 1.20000000000000),
    Param(LSR, 0, 0, 1.20000000000000, 0, 0, 1.20000000000000, 0, 1.20000000000000),
    Param(RSL, 0, 0, 1.20000000000000, 0, 0, 1.20000000000000, 0, 1.20000000000000),
    Param(RSR, 0, 0, 1.20000000000000, 0, 0, 1.20000000000000, 0, 1.20000000000000),
    Param(RLR, 0, 0, 1.20000000000000, 0, 2.83689999957440, 5.67379999914879, 2.83689999957440, 11.3475999982976),
    Param(LRL, 0, 0, 1.20000000000000, 0, 2.83689999957440, 5.67379999914879, 2.83689999957440, 11.3475999982976),
    Param(LSL, 0, 0, 10.0000000000000, 0, 0, 10.0000000000000, 0, 10.0000000000000),
    Param(LSR, 0, 0, 10.0000000000000, 0, 0, 10.0000000000000, 0, 10.0000000000000),
    Param(RSL, 0, 0, 10.0000000000000, 0, 0, 10.0000000000000, 0, 10.0000000000000),
    Param(RSR, 0, 0, 10.0000000000000, 0, 0, 10.0000000000000, 0, 10.0000000000000),
    Param(RLR, 0, 0, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LRL, 0, 0, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LSL, 0, 1.57079632679490, 0.200000000000000, 0, 4.03764803816114, 1.28062484748657, 3.81633359581335, 9.13460648146105),
    Param(LSR, 0, 1.57079632679490, 0.200000000000000, 1, 0, 0, 0, 0),
    Param(RSL, 0, 1.57079632679490, 0.200000000000000, 1, 0, 0, 0, 0),
    Param(RSR, 0, 1.57079632679490, 0.200000000000000, 0, 5.58844703098288, 1.56204993518133, 5.40712725658139, 12.5576242227456),
    Param(RLR, 0, 1.57079632679490, 0.200000000000000, 0, 2.04566616439583, 5.48080888118507, 1.86434638999434, 9.39082143557524),
    Param(LRL, 0, 1.57079632679490, 0.200000000000000, 0, 0.570161010924784, 5.63139655988647, 0.348846568576992, 6.55040413938824),
    Param(LSL, 0, 1.57079632679490, 0.600000000000000, 0, 4.33188260327232, 1.07703296142690, 3.52209903070216, 8.93101459540138),
    Param(LSR, 0, 1.57079632679490, 0.600000000000000, 1, 0, 0, 0, 0),
    Param(RSL, 0, 1.57079632679490, 0.600000000000000, 1, 0, 0, 0, 0),
    Param(RSR, 0, 1.57079632679490, 0.600000000000000, 0, 5.72458599183602, 1.88679622641132, 5.27098829572825, 12.8823705139756),
    Param(RLR, 0, 1.57079632679490, 0.600000000000000, 0, 2.09177665791708, 5.30075194652129, 1.63817896180931, 9.03070756624769),
    Param(LRL, 0, 1.57079632679490, 0.600000000000000, 0, 0.917667206073218, 5.73793981996096, 0.107883633503051, 6.76349065953723),
    Param(LSL, 0, 1.57079632679490, 1.20000000000000, 0, 4.90978454023457, 1.01980390271856, 2.94419709373991, 8.87378553669304),
    Param(LSR, 0, 1.57079632679490, 1.20000000000000, 0, 0.548201457476526, 1.35646599662505, 5.26059043786122, 7.16525789196279),
    Param(RSL, 0, 1.57079632679490, 1.20000000000000, 1, 0, 0, 0, 0),
    Param(RSR, 0, 1.57079632679490, 1.20000000000000, 0, 5.85655781405271, 2.41660919471891, 5.13901647351157, 13.4121834822832),
    Param(RLR, 0, 1.57079632679490, 1.20000000000000, 0, 2.06626351279921, 4.98578201185216, 1.34872217225806, 8.40076769690943),
    Param(LRL, 0, 1.57079632679490, 1.20000000000000, 0, 1.51039488352173, 5.76759130093350, 5.82799274420666, 13.1059789286619),
    Param(LSL, 0, 1.57079632679490, 10.0000000000000, 0, 6.17252808600569, 9.05538513813742, 1.68145354796879, 16.9093667721119),
    Param(LSR, 0, 1.57079632679490, 10.0000000000000, 0, 0.0914159566263198, 10.8627804912002, 4.80380493701101, 15.7580013848375),
    Param(RSL, 0, 1.57079632679490, 10.0000000000000, 0, 0.112042065821841, 8.83176086632785, 1.68283839261674, 10.6266413247664),
    Param(RSR, 0, 1.57079632679490, 10.0000000000000, 0, 6.19252541997884, 11.0453610171873, 4.80304886758544, 22.0409353047515),
    Param(RLR, 0, 1.57079632679490, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LRL, 0, 1.57079632679490, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LSL, 1.57079632679490, -1.57079632679490, 0.200000000000000, 0, 4.71238898038469, 2.20000000000000, 4.71238898038469, 11.6247779607694),
    Param(LSR, 1.57079632679490, -1.57079632679490, 0.200000000000000, 1, 0, 0, 0, 0),
    Param(RSL, 1.57079632679490, -1.57079632679490, 0.200000000000000, 1, 0, 0, 0, 0),
    Param(RSR, 1.57079632679490, -1.57079632679490, 0.200000000000000, 0, 4.71238898038469, 1.80000000000000, 4.71238898038469, 11.2247779607694),
    Param(RLR, 1.57079632679490, -1.57079632679490, 0.200000000000000, 0, 1.10403098774760, 5.34965462908499, 1.10403098774760, 7.55771660458019),
    Param(LRL, 1.57079632679490, -1.57079632679490, 0.200000000000000, 0, 0.988432088926153, 5.11845683144210, 0.988432088926153, 7.09532100929441),
    Param(LSL, 1.57079632679490, -1.57079632679490, 0.600000000000000, 0, 4.71238898038469, 2.60000000000000, 4.71238898038469, 12.0247779607694),
    Param(LSR, 1.57079632679490, -1.57079632679490, 0.600000000000000, 1, 0, 0, 0, 0),
    Param(RSL, 1.57079632679490, -1.57079632679490, 0.600000000000000, 1, 0, 0, 0, 0),
    Param(RSR, 1.57079632679490, -1.57079632679490, 0.600000000000000, 0, 4.71238898038469, 1.40000000000000, 4.71238898038469, 10.8247779607694),
    Param(RLR, 1.57079632679490, -1.57079632679490, 0.600000000000000, 0, 1.21322522314939, 5.56804309988857, 1.21322522314939, 7.99449354618734),
    Param(LRL, 1.57079632679490, -1.57079632679490, 0.600000000000000, 0, 0.863211890069541, 4.86801643372887, 0.863211890069541, 6.59444021386796),
    Param(LSL, 1.57079632679490, -1.57079632679490, 1.20000000000000, 0, 4.71238898038469, 3.20000000000000, 4.71238898038469, 12.6247779607694),
    Param(LSR, 1.57079632679490, -1.57079632679490, 1.20000000000000, 1, 0, 0, 0, 0),
    Param(RSL, 1.57079632679490, -1.57079632679490, 1.20000000000000, 1, 0, 0, 0, 0),
    Param(RSR, 1.57079632679490, -1.57079632679490, 1.20000000000000, 0, 4.71238898038469, 0.800000000000000, 4.71238898038469, 10.2247779607694),
    Param(RLR, 1.57079632679490, -1.57079632679490, 1.20000000000000, 0, 1.36943840600457, 5.88046946559892, 1.36943840600457, 8.61934627760806),
    Param(LRL, 1.57079632679490, -1.57079632679490, 1.20000000000000, 0, 0.643501108793284, 4.42859487117636, 0.643501108793284, 5.71559708876293),
    Param(LSL, 1.57079632679490, -1.57079632679490, 10.0000000000000, 0, 4.71238898038469, 12.0000000000000, 4.71238898038469, 21.4247779607694),
    Param(LSR, 1.57079632679490, -1.57079632679490, 10.0000000000000, 0, 4.91374690117502, 9.79795897113271, 1.77215424758523, 16.4838601198930),
    Param(RSL, 1.57079632679490, -1.57079632679490, 10.0000000000000, 0, 1.77215424758523, 9.79795897113271, 4.91374690117502, 16.4838601198930),
    Param(RSR, 1.57079632679490, -1.57079632679490, 10.0000000000000, 0, 1.57079632679490, 8.00000000000000, 1.57079632679490, 11.1415926535898),
    Param(RLR, 1.57079632679490, -1.57079632679490, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LRL, 1.57079632679490, -1.57079632679490, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LSL, 1.57079632679490, 0, 0.200000000000000, 0, 5.40712725658139, 1.56204993518133, 5.58844703098288, 12.5576242227456),
    Param(LSR, 1.57079632679490, 0, 0.200000000000000, 1, 0, 0, 0, 0),
    Param(RSL, 1.57079632679490, 0, 0.200000000000000, 1, 0, 0, 0, 0),
    Param(RSR, 1.57079632679490, 0, 0.200000000000000, 0, 3.81633359581335, 1.28062484748657, 4.03764803816114, 9.13460648146105),
    Param(RLR, 1.57079632679490, 0, 0.200000000000000, 0, 0.348846568576992, 5.63139655988647, 0.570161010924784, 6.55040413938824),
    Param(LRL, 1.57079632679490, 0, 0.200000000000000, 0, 1.86434638999434, 5.48080888118507, 2.04566616439583, 9.39082143557524),
    Param(LSL, 1.57079632679490, 0, 0.600000000000000, 0, 5.27098829572825, 1.88679622641132, 5.72458599183602, 12.8823705139756),
    Param(LSR, 1.57079632679490, 0, 0.600000000000000, 1, 0, 0, 0, 0),
    Param(RSL, 1.57079632679490, 0, 0.600000000000000, 1, 0, 0, 0, 0),
    Param(RSR, 1.57079632679490, 0, 0.600000000000000, 0, 3.52209903070216, 1.07703296142690, 4.33188260327232, 8.93101459540138),
    Param(RLR, 1.57079632679490, 0, 0.600000000000000, 0, 0.107883633503051, 5.73793981996096, 0.917667206073218, 6.76349065953723),
    Param(LRL, 1.57079632679490, 0, 0.600000000000000, 0, 1.63817896180931, 5.30075194652129, 2.09177665791708, 9.03070756624769),
    Param(LSL, 1.57079632679490, 0, 1.20000000000000, 0, 5.13901647351157, 2.41660919471891, 5.85655781405271, 13.4121834822832),
    Param(LSR, 1.57079632679490, 0, 1.20000000000000, 0, 5.26059043786122, 1.35646599662505, 0.548201457476526, 7.16525789196279),
    Param(RSL, 1.57079632679490, 0, 1.20000000000000, 1, 0, 0, 0, 0),
    Param(RSR, 1.57079632679490, 0, 1.20000000000000, 0, 2.94419709373991, 1.01980390271856, 4.90978454023457, 8.87378553669304),
    Param(RLR, 1.57079632679490, 0, 1.20000000000000, 0, 5.82799274420666, 5.76759130093350, 1.51039488352173, 13.1059789286619),
    Param(LRL, 1.57079632679490, 0, 1.20000000000000, 0, 1.34872217225806, 4.98578201185216, 2.06626351279921, 8.40076769690943),
    Param(LSL, 1.57079632679490, 0, 10.0000000000000, 0, 4.80304886758544, 11.0453610171873, 6.19252541997884, 22.0409353047515),
    Param(LSR, 1.57079632679490, 0, 10.0000000000000, 0, 4.80380493701101, 10.8627804912002, 0.0914159566263198, 15.7580013848375),
    Param(RSL, 1.57079632679490, 0, 10.0000000000000, 0, 1.68283839261674, 8.83176086632785, 0.112042065821841, 10.6266413247664),
    Param(RSR, 1.57079632679490, 0, 10.0000000000000, 0, 1.68145354796879, 9.05538513813742, 6.17252808600569, 16.9093667721119),
    Param(RLR, 1.57079632679490, 0, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LRL, 1.57079632679490, 0, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LSL, 1.57079632679490, 1.57079632679490, 0.200000000000000, 0, 4.71238898038469, 0.200000000000000, 1.57079632679490, 6.48318530717959),
    Param(LSR, 1.57079632679490, 1.57079632679490, 0.200000000000000, 0, 5.85348564102816, 0.916515138991168, 5.85348564102816, 12.6234864210475),
    Param(RSL, 1.57079632679490, 1.57079632679490, 0.200000000000000, 1, 0, 0, 0, 0),
    Param(RSR, 1.57079632679490, 1.57079632679490, 0.200000000000000, 0, 1.57079632679490, 0.200000000000000, 4.71238898038469, 6.48318530717959),
    Param(RLR, 1.57079632679490, 1.57079632679490, 0.200000000000000, 0, 4.66236812357892, 6.18314359356805, 1.52077546998913, 12.3662871871361),
    Param(LRL, 1.57079632679490, 1.57079632679490, 0.200000000000000, 0, 1.52077546998913, 6.18314359356805, 4.66236812357892, 12.3662871871361),
    Param(LSL, 1.57079632679490, 1.57079632679490, 0.600000000000000, 0, 4.71238898038469, 0.600000000000000, 1.57079632679490, 6.88318530717959),
    Param(LSR, 1.57079632679490, 1.57079632679490, 0.600000000000000, 0, 5.59002539960283, 1.66132477258361, 5.59002539960283, 12.8413755717893),
    Param(RSL, 1.57079632679490, 1.57079632679490, 0.600000000000000, 1, 0, 0, 0, 0),
    Param(RSR, 1.57079632679490, 1.57079632679490, 0.600000000000000, 0, 1.57079632679490, 0.600000000000000, 4.71238898038469, 6.88318530717959),
    Param(RLR, 1.57079632679490, 1.57079632679490, 0.600000000000000, 0, 4.56182070760800, 5.98204876162621, 1.42022805401821, 11.9640975232524),
    Param(LRL, 1.57079632679490, 1.57079632679490, 0.600000000000000, 0, 1.42022805401821, 5.98204876162621, 4.56182070760800, 11.9640975232524),
    Param(LSL, 1.57079632679490, 1.57079632679490, 1.20000000000000, 0, 4.71238898038469, 1.20000000000000, 1.57079632679490, 7.48318530717959),
    Param(LSR, 1.57079632679490, 1.57079632679490, 1.20000000000000, 0, 5.38752051332172, 2.49799919935936, 5.38752051332172, 13.2730402260028),
    Param(RSL, 1.57079632679490, 1.57079632679490, 1.20000000000000, 1, 0, 0, 0, 0),
    Param(RSR, 1.57079632679490, 1.57079632679490, 1.20000000000000, 0, 1.57079632679490, 1.20000000000000, 4.71238898038469, 7.48318530717959),
    Param(RLR, 1.57079632679490, 1.57079632679490, 1.20000000000000, 0, 4.40769632636929, 5.67379999914879, 1.26610367277950, 11.3475999982976),
    Param(LRL, 1.57079632679490, 1.57079632679490, 1.20000000000000, 0, 1.26610367277950, 5.67379999914879, 4.40769632636929, 11.3475999982976),
    Param(LSL, 1.57079632679490, 1.57079632679490, 10.0000000000000, 0, 4.71238898038469, 10.0000000000000, 1.57079632679490, 16.2831853071796),
    Param(LSR, 1.57079632679490, 1.57079632679490, 10.0000000000000, 0, 4.87983705960438, 11.8321595661992, 4.87983705960438, 21.5918336854080),
    Param(RSL, 1.57079632679490, 1.57079632679490, 10.0000000000000, 0, 1.82347658193698, 7.74596669241483, 1.82347658193698, 11.3929198562888),
    Param(RSR, 1.57079632679490, 1.57079632679490, 10.0000000000000, 0, 1.57079632679490, 10.0000000000000, 4.71238898038469, 16.2831853071796),
    Param(RLR, 1.57079632679490, 1.57079632679490, 10.0000000000000, 1, 0, 0, 0, 0),
    Param(LRL, 1.57079632679490, 1.57079632679490, 10.0000000000000, 1, 0, 0, 0, 0)
]

@testset "test for the dubins paths" begin
    ρ = 1.0
    for i in 1:length(params)
        param = params[i]
        errcode, path = dubins_path([0.0, 0.0, param.inputs.α], [param.inputs.d, 0.0, param.inputs.β], ρ, param.inputs.word)
        (errcode != 0) && (errcode = 1)
        @test errcode == param.outputs.errcode
        (errcode == 0) && (@test isapprox(path.params, param.outputs.params, atol=1e-8))
        (errcode == 0) && (@test isapprox(dubins_path_length(path), param.outputs.len, atol=1e-8))
    end
end


text:
"""
"""
function fit!(
        transformer::ImmutablePackageMultiLabelPredictionTransformer,
        varargs...;
        kwargs...
        )
    if length(varargs) == 1
        return varargs[1]
    else
        return varargs
    end
end

"""
"""
function predict(
        transformer::ImmutablePackageMultiLabelPredictionTransformer,
        single_labelpredictions::AbstractMatrix,
        varargs...
        )
    result = DataFrames.DataFrame()
    for i = 1:length(transformer.label_names)
        result[transformer.label_names[i]] = single_labelpredictions[:, i]
    end
    return result
end

"""
"""
function predict_proba(
        transformer::ImmutablePackageMultiLabelPredictionTransformer,
        varargs...;
        kwargs...
        )
    if length(varargs) == 1
        return varargs[1]
    else
        return varargs
    end
end


text:
using Test, Random, LinearAlgebra

using Yao, YaoExtensions

@testset "RotBasis" begin
    rt = RotBasis(0.5, 0.4)
    crt = chain(rt)

    dispatch!(crt, [2., 3.])
    @test nparameters(crt) == 2

    for (t1, t2, t3) in zip(parameters(rt), (2, 3), parameters(crt))
        @test t1 == t2 == t3
    end

    # check consistency
    rb = put(1, 1=>RotBasis(0.1, 0.3))#rot_basis(1)
    angles = randpolar(1)
    # prepair a state in the angles direction.
    psi = angles |> polar2u |> ArrayReg

    # rotate to the same direction for measurements.
    dispatch!(rb, vec(angles))
    @test state(apply!(psi, rb)) ≈ [1, 0]

    @test nparameters(rot_basis(3)) == 6
    dispatch!(rb, :zero)
    @test parameters(rb)[1] == 0
    dispatch!(rb, :random)
    @test parameters(rb)[1] != 0
end

@testset "polar and u" begin
    polar = randpolar(10)
    @test size(polar) == (2, 10)
    @test polar |> polar2u |> u2polar ≈ polar
end

@testset "basis rotor" begin
    mblock = Measure(1)

    reg = repeat(ArrayReg([0im,1]/sqrt(2)), 10)
    reg |> basis_rotor(Z) |> mblock
    @test mblock.results ≈ fill(1,10)

    reg = repeat(ArrayReg([1+0im,1]/sqrt(2)), 10)
    reg |> basis_rotor(X) |> mblock
    @test mblock.results ≈ fill(0,10)

    reg = repeat(ArrayReg([im,1.0]/sqrt(2)), 10)
    reg |> basis_rotor(Y) |> mblock
    @test mblock.results ≈ fill(1,10)
end