text:
### press-rybicki.jl ---  Plan the Lomb–Scargle periodogram
#
# Copyright (C) 2017 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: periodogram, lomb scargle
#
# This file is a part of LombScargle.jl.
#
# License is BSD 3-clause "New" or "Revised".
#
### Code:

# Switches to select the appropriate algorithm to compute the periodogram.
function _plan_no_fast(times::AbstractVector{R1}, signal::AbstractVector{R2}, sumw::Real,
                       w::AbstractVector{R3}, frequencies::AbstractVector{R4},
                       with_errors::Bool, center_data::Bool,
                       fit_mean::Bool, noise::Real,
                       normalization::Symbol) where {R1<:Real,R2<:Real,R3<:Real,R4<:Real}
    P_type = promote_type(float(R1), float(R2), float(R3), float(R4))
    P = Vector{P_type}(undef, length(frequencies))
    if fit_mean || with_errors
        # If "center_data" or "fit_mean" keywords are true,
        # subtract the weighted mean from each point.
        if center_data || fit_mean
            y = signal .- (w ⋅ signal)
        else
            y = signal
        end
        YY = w ⋅ (y .^ 2)
        if fit_mean
            Y   = w ⋅ y
            YY -= Y * Y
            return GLSPlan_fit_mean(times, signal, frequencies, sumw, w, y, Y,
                                    YY, noise, normalization, P)
        else
            return GLSPlan(times, signal, frequencies, sumw, w, y,
                           YY, noise, normalization, P)
        end
    else
        if center_data
            X = signal .- mean(signal)
        else
            X = signal
        end
        return LSPlan(times, signal, frequencies, sumw, w, X, X ⋅ X, noise, normalization, P)
    end
end

# There are two "_plan" methods, the only different argument being "frequencies".  When it
# is a `Range' object (first method) it could be possible to use the fast method, provided
# that fast == true, otherwise we can only use the non fast methods.
function _plan(times::AbstractVector{<:Real}, signal::AbstractVector{R1}, sumw::Real,
               w::AbstractVector{R2}, frequencies::AbstractRange{<:Real}, fast::Bool,
               with_errors::Bool, center_data::Bool, fit_mean::Bool, flags::Integer,
               timelimit::Real, oversampling::Integer, Mfft::Integer, noise::Real,
               normalization::Symbol) where {R1<:Real,R2<:Real}
    if fast
        @assert Mfft > 0
        @assert step(frequencies) > 0
        if center_data || fit_mean
            y = signal .- w⋅signal
        else
            y = signal
        end
        YY = w ⋅ (y .^ 2)
        N = length(frequencies)
        Nfft = nextpow(2, N * oversampling)
        T = promote_type(float(R1), float(R2))
        bfft_vect = Vector{Complex{T}}(undef, Nfft)
        bfft_grid = Vector{Complex{T}}(undef, Nfft)
        fftgrid   = Vector{Complex{T}}(undef, N)
        bfft_plan = FFTW.plan_bfft(bfft_vect, flags = flags, timelimit = timelimit)
        if fit_mean
            return FastGLSPlan_fit_mean(times, signal, frequencies, sumw, w, y, YY,
                                        fftgrid, bfft_vect, bfft_grid, bfft_plan, Mfft,
                                        noise, normalization, Vector{T}(undef, N))
        else
            return FastGLSPlan(times, signal, frequencies, sumw, w, y, YY,
                               fftgrid, bfft_vect, bfft_grid, bfft_plan, Mfft,
                               noise, normalization, Vector{T}(undef, N))
        end
    else
        return _plan_no_fast(times, signal, sumw, w, frequencies, with_errors,
                             center_data, fit_mean, noise, normalization)
    end
end

_plan(times::AbstractVector{<:Real}, signal::AbstractVector{<:Real}, sumw::Real,
      w::AbstractVector{<:Real}, frequencies::AbstractVector{<:Real},
      fast::Bool, with_errors::Bool, center_data::Bool, fit_mean::Bool, flags::Integer,
      timelimit::Real, oversampling::Integer, Mfft::Integer,
      noise::Real, normalization::Symbol) =
          _plan_no_fast(times, signal, sumw, w, frequencies,
                        with_errors, center_data, fit_mean, noise, normalization)

function _plan(times::AbstractVector{<:Real},
               signal::AbstractVector{<:Real},
               with_errors::Bool,
               sumw::Real=length(signal),
               w::AbstractVector{<:Real}=fill!(similar(signal), one(eltype(signal)))/sumw;
               normalization::Symbol=:standard,
               noise_level::Real=1,
               center_data::Bool=true,
               fit_mean::Bool=true,
               flags::Integer=FFTW.ESTIMATE,
               timelimit::Real=Inf,
               oversampling::Integer=5,
               Mfft::Integer=4,
               samples_per_peak::Integer=5,
               nyquist_factor::Integer=5,
               minimum_frequency::Real=NaN,
               maximum_frequency::Real=NaN,
               frequencies::AbstractVector{<:Real}=
               autofrequency(times,
                             samples_per_peak=samples_per_peak,
                             nyquist_factor=nyquist_factor,
                             minimum_frequency=minimum_frequency,
                             maximum_frequency=maximum_frequency),
               fast::Bool=(length(frequencies) > 200))
    @assert length(times) == length(signal) == length(w)
    return _plan(times, signal, sumw, w, frequencies, fast, with_errors, center_data,
                 fit_mean, flags, timelimit, oversampling, Mfft, noise_level, normalization)
end

### Main interface functions

# No uncertainties
plan(times::AbstractVector{<:Real}, signal::AbstractVector{<:Real}; kwargs...) =
    _plan(times, signal, false; kwargs...)

# Uncertainties provided
function plan(times::AbstractVector{<:Real}, signal::AbstractVector{<:Real},
              errors::AbstractVector{<:Real}; kwargs...)
    # Compute weights vector
    w = 1 ./ errors .^ 2
    sumw = sum(w)
    w ./= sumw
    return _plan(times, signal, true, sumw, w; kwargs...)
end

# Uncertainties provided via Measurement type
plan(times::AbstractVector{<:Real}, signal::AbstractVector{<:Measurement}; kwargs...) =
         plan(times, Measurements.value.(signal), Measurements.uncertainty.(signal); kwargs...)

"""
    LombScargle.plan(times::AbstractVector{<:Real}, signal::AbstractVector{<:Real},
                     [errors::AbstractVector{<:Real}];
                     normalization::Symbol=:standard,
                     noise_level::Real=1,
                     center_data::Bool=true,
                     fit_mean::Bool=true,
                     fast::Bool=true,
                     flags::Integer=FFTW.ESTIMATE,
                     timelimit::Real=Inf,
                     oversampling::Integer=5,
                     Mfft::Integer=4,
                     samples_per_peak::Integer=5,
                     nyquist_factor::Integer=5,
                     minimum_frequency::Real=NaN,
                     maximum_frequency::Real=NaN,
                     frequencies::AbstractVector{Real}=
                     autofrequency(times,
                                   samples_per_peak=samples_per_peak,
                                   nyquist_factor=nyquist_factor,
                                   minimum_frequency=minimum_frequency,
                                   maximum_frequency=maximum_frequency))

Pre-plan the Lomb–Scargle periodogram of the `signal` vector, observed at
`times`.  The periodogram can then be computed by passing the result of this
function to `lombscargle`.

You can also specify the uncertainties for each signal point with `errors`
argument.  All these vectors must have the same length.

Optional keywords arguments are:

* `normalization`: how to normalize the periodogram.  Valid choices are:
  `:standard`, `:model`, `:log`, `:psd`, `:Scargle`, `:HorneBaliunas`,
  `:Cumming`
* `noise_level`: the noise level used to normalize the periodogram when
  `normalization` is set to `:Scargle`
* `fit_mean`: if `true`, fit for the mean of the signal using the Generalised
  Lomb–Scargle algorithm (see Zechmeister & Kürster paper below).  If this is
  `false` and no uncertainty on the signal is provided, the original algorithm
  by Lomb and Scargle will be employed (see Townsend paper below)
* `center_data`: if `true`, subtract the weighted mean of `signal` from `signal`
  itself before performing the periodogram.  This is especially important if
  `fit_mean` is `false`
* `frequencies`: the frequecy grid (not angular frequencies) at which the
  periodogram will be computed, as a vector.  If not provided, it is an evenly
  spaced grid of type `Range`, automatically determined with
  `LombScargle.autofrequency` function, which see.  See below for other
  available keywords that can be used to affect the frequency grid without
  directly setting `frequencies`

You can explicitely require to use or not the fast method by Press & Rybicki,
overriding the default choice, by setting the `fast` keyword.  In any case,
`frequencies` must be a `Range` object (this is the default) in order to
actually use this method.  A few other keywords are available to adjust the
settings of the periodogram when the fast method is used (otherwise they are
ignored):

* `fast`: whether to use the fast method.
* `flags`: this integer keyword is a bitwise-or of FFTW planner flags,
  defaulting to `FFTW.ESTIMATE`.  Passing `FFTW.MEASURE` or `FFTW.PATIENT` will
  instead spend several seconds (or more) benchmarking different possible FFT
  algorithms and picking the fastest one; see the FFTW manual for more
  information on planner flags.
* `timelimit`: specifies a rough upper bound on the allowed planning time, in seconds.
* `oversampling`: oversampling the frequency factor for the approximation;
  roughly the number of time samples across the highest-frequency sinusoid.
  This parameter contains the tradeoff between accuracy and speed.
* `Mfft`: the number of adjacent points to use in the FFT approximation.

In addition, you can use all optional keyword arguments of
[`LombScargle.autofrequency`](@ref) function in order to tune the `frequencies`.

If the signal has uncertainties, the `signal` vector can also be a vector of
`Measurement` objects (from
[`Measurements.jl`](https://github.com/giordano/Measurements.jl) package), in
which case you don’t need to pass a separate `errors` vector for the
uncertainties of the signal.
"""
plan


text:
import GitHub

function all_specified_statuses_passed(registry::GitHub.Repo,
                                       pr::GitHub.PullRequest,
                                       sha::AbstractString,
                                       specified_status_context_list::AbstractVector{<:AbstractString};
                                       auth::GitHub.Authorization)
    specified_status_context_did_fail = Dict{String, Bool}()
    for context in specified_status_context_list
        specified_status_context_did_fail[context] = true
    end
    combined_status = GitHub.status(registry, sha; auth = auth)
    all_statuses = combined_status.statuses
    for status in all_statuses
        context = status.context
        if haskey(specified_status_context_did_fail, context)
            status_was_success = status.state == "success"
            status_was_failure = !status_was_success
            specified_status_context_did_fail[context] = status_was_failure
        end
    end
    if any(values(specified_status_context_did_fail))
        return false
    end
    return true
end

function all_specified_check_runs_passed(registry::GitHub.Repo,
                                         pr::GitHub.PullRequest,
                                         sha::AbstractString,
                                         specified_check_run_name_list::AbstractVector{<:AbstractString};
                                         auth::GitHub.Authorization)
    specified_check_run_name_did_fail = Dict{String, Bool}()
    for context in specified_check_run_name_list
        specified_check_run_name_did_fail[context] = true
    end
    endpoint = "/repos/$(registry.full_name)/commits/$(sha)/check-runs"
    check_runs = GitHub.gh_get_json(GitHub.DEFAULT_API,
                                    endpoint;
                                    auth = auth,
                                    headers = Dict("Accept" =>
                                                   "application/vnd.github.antiope-preview+json"))
    for check_run in check_runs["check_runs"]
        name = check_run["name"]
        check_run_was_success = (check_run["status"] == "completed") & (check_run["conclusion"] == "success")
        check_run_was_failure = !check_run_was_success
        if haskey(specified_check_run_name_did_fail, name)
            specified_check_run_name_did_fail[name] = check_run_was_failure
        end
    end
    if any(values(specified_check_run_name_did_fail))
        return false
    end
    return true
end

function pr_comment_is_blocking(c::GitHub.Comment)
    c_body = body(c)
    if occursin("[noblock]", c_body)
        return false
    else
        return true
    end
end

function pr_has_no_blocking_comments(registry::GitHub.Repo,
                                     pr::GitHub.PullRequest;
                                     auth::GitHub.Authorization)
    all_pr_comments = get_all_pull_request_comments(registry, pr; auth = auth)
    if isempty(all_pr_comments)
        return true
    else
        num_comments = length(all_pr_comments)
        comment_is_blocking = BitVector(undef, num_comments)
        for i = 1:num_comments
            comment_is_blocking[i] = pr_comment_is_blocking(all_pr_comments[i])
        end
        if any(comment_is_blocking)
            return false
        else
            return true
        end
    end
end

function pr_is_old_enough(pr_type::Symbol,
                          pr_age::Dates.Period;
                          pkg::AbstractString,
                          new_package_waiting_period::Dates.Period,
                          new_jll_package_waiting_period::Dates.Period,
                          new_version_waiting_period::Dates.Period,
                          new_jll_version_waiting_period::Dates.Period)
    if is_jll_name(pkg)
        if pr_type == :NewPackage
            return pr_age > new_jll_package_waiting_period
        elseif pr_type == :NewVersion
            return pr_age > new_jll_version_waiting_period
        else
            throw(ArgumentError("pr_type must be either :NewPackage or :NewVersion"))
        end
    else
        if pr_type == :NewPackage
            return pr_age > new_package_waiting_period
        elseif pr_type == :NewVersion
            return pr_age > new_version_waiting_period
        else
            throw(ArgumentError("pr_type must be either :NewPackage or :NewVersion"))
        end
    end
end

function _get_all_pr_statuses(repo::GitHub.Repo,
                              pr::GitHub.PullRequest;
                              auth::GitHub.Authorization)
    combined_status = GitHub.status(repo, pr.head.sha)
    all_statuses = combined_status.statuses
    return all_statuses
end

function _get_status_description(status::GitHub.Status)::String
    if hasproperty(status, :description)
        if status.description === nothing
            return ""
        else
            return status.description
        end
    else
        return ""
    end
end

function _postprocess_automerge_decision_status(status::GitHub.Status;
                                                whoami)
    @debug("status: ", status)
    @debug("status.creator: ", status.creator)
    new_package_passed_regex = r"New package. Approved. name=\"(\w*)\". sha=\"(\w*)\""
    new_version_passed_regex = r"New version. Approved. name=\"(\w*)\". sha=\"(\w*)\""
    status_description = _get_status_description(status)
    if status.state == "success" && occursin(new_package_passed_regex,
                                             status_description)
        m = match(new_package_passed_regex,
                  status_description)
        passed_pkg_name = m[1]
        passed_pr_head_sha = m[2]
        return true, passed_pkg_name, passed_pr_head_sha, :NewPackage
    end
    if status.state == "success" && occursin(new_version_passed_regex,
                                             status_description)
        m = match(new_version_passed_regex,
                  status_description)
        passed_pkg_name = m[1]
        passed_pr_head_sha = m[2]
        return true, passed_pkg_name, passed_pr_head_sha, :NewVersion
    end
    return false, "", "", :failing
end

function pr_has_passing_automerge_decision_status(repo::GitHub.Repo,
                                                  pr::GitHub.PullRequest;
                                                  auth::GitHub.Authorization,
                                                  whoami)
    all_statuses = _get_all_pr_statuses(repo, pr; auth = auth)
    for status in all_statuses
        if status.context == "automerge/decision"
            return _postprocess_automerge_decision_status(status;
                                                          whoami = whoami)
        end
    end
    return false, "", "", :failing
end

function cron_or_api_build(registry::GitHub.Repo;
                           auth::GitHub.Authorization,
                           authorized_authors::Vector{String},
                           merge_new_packages::Bool,
                           merge_new_versions::Bool,
                           new_package_waiting_period,
                           new_jll_package_waiting_period,
                           new_version_waiting_period,
                           new_jll_version_waiting_period,
                           whoami::String,
                           all_statuses::AbstractVector{<:AbstractString},
                           all_check_runs::AbstractVector{<:AbstractString})
    # first, get a list of ALL open pull requests on this repository
    # then, loop through each of them.
    all_currently_open_pull_requests = my_retry(() -> get_all_pull_requests(registry, "open"; auth = auth))
    reverse!(all_currently_open_pull_requests)
    at_least_one_exception_was_thrown = false
    num_retries = 0
    if isempty(all_currently_open_pull_requests)
        @info("There are no open pull requests.")
    else
        for pr in all_currently_open_pull_requests
            try
                my_retry(() -> cron_or_api_build(pr,
                                                 registry::GitHub.Repo;
                                                 auth = auth,
                                                 authorized_authors = authorized_authors,
                                                 merge_new_packages = merge_new_packages,
                                                 merge_new_versions = merge_new_versions,
                                                 new_package_waiting_period = new_package_waiting_period,
                                                 new_jll_package_waiting_period = new_jll_package_waiting_period,
                                                 new_version_waiting_period = new_version_waiting_period,
                                                 new_jll_version_waiting_period = new_jll_version_waiting_period,
                                                 whoami = whoami,
                                                 all_statuses = all_statuses,
                                                 all_check_runs = all_check_runs),
                         num_retries)
            catch ex
                at_least_one_exception_was_thrown = true
                showerror(stderr, ex)
                Base.show_backtrace(stderr, catch_backtrace())
                println(stderr)
            end
        end
        if at_least_one_exception_was_thrown
            throw(AutoMergeCronJobError("At least one exception was thrown. Check the logs for details."))
        end
    end
    return nothing
end

function cron_or_api_build(pr::GitHub.PullRequest,
                           registry::GitHub.Repo;
                           auth::GitHub.Authorization,
                           authorized_authors::Vector{String},
                           merge_new_packages::Bool,
                           merge_new_versions::Bool,
                           new_package_waiting_period,
                           new_jll_package_waiting_period,
                           new_version_waiting_period,
                           new_jll_version_waiting_period,
                           whoami::String,
                           all_statuses::AbstractVector{<:AbstractString},
                           all_check_runs::AbstractVector{<:AbstractString})
    #       first, see if the author is an authorized author. if not, then skip.
    #       next, see if the title matches either the "New Version" regex or
    #               the "New Package regex". if it is not either a new
    #               package or a new version, skip.
    #       next, see if it is old enough. if it is not old enough, then skip.
    #       then, get the `automerge/decision` status and make sure it is passing
    #       then, get all of the pull request comments. if there is any comment that is
    #               (1) not by me, and (2) does not contain the text [noblock], then skip
    #       if all of the above criteria were met, then merge the pull request
    pr_number = number(pr)
    @info("Now examining pull request $(pr_number)")
    pr_author = author_login(pr)
    if pr_author in authorized_authors
        if is_new_package(pr) || is_new_version(pr)
            if is_new_package(pr) # it is a new package
                pr_type = :NewPackage
                pkg, version = parse_pull_request_title(NewPackage(), pr)
            else # it is a new version
                pr_type = :NewVersion
                pkg, version = parse_pull_request_title(NewVersion(), pr)
            end
            pr_age = time_since_pr_creation(pr)
            this_pr_is_old_enough = pr_is_old_enough(pr_type,
                                                     pr_age;
                                                     pkg = pkg,
                                                     new_package_waiting_period = new_package_waiting_period,
                                                     new_jll_package_waiting_period = new_jll_package_waiting_period,
                                                     new_version_waiting_period = new_version_waiting_period,
                                                     new_jll_version_waiting_period = new_jll_version_waiting_period)
            if this_pr_is_old_enough
                i_passed_this_pr,
                    passed_pkg_name,
                    passed_pr_head_sha,
                    status_pr_type = pr_has_passing_automerge_decision_status(registry,
                                                                              pr;
                                                                              auth = auth,
                                                                              whoami = whoami)
                if i_passed_this_pr
                    always_assert(pkg == passed_pkg_name)
                    always_assert(pr.head.sha == passed_pr_head_sha)
                    always_assert(all_specified_statuses_passed(registry,
                                                                pr,
                                                                passed_pr_head_sha,
                                                                all_statuses;
                                                                auth = auth))
                    always_assert(all_specified_check_runs_passed(registry,
                                                                  pr,
                                                                  passed_pr_head_sha,
                                                                  all_check_runs;
                                                                  auth = auth))
                    if pr_has_no_blocking_comments(registry, pr; auth = auth)
                        "Pull request: $(pr_number). "
                        "Type: $(pr_type). "
                        "Decision: merge. "
                        if pr_type == :NewPackage # it is a new package
                            always_assert(status_pr_type == :NewPackage)
                            if merge_new_packages
                                my_comment = comment_text_merge_now()
                                @info(string("Pull request: $(pr_number). ",
                                             "Type: $(pr_type). ",
                                             "Decision: merge now."))
                                my_retry(() -> merge!(registry, pr, passed_pr_head_sha; auth = auth))
                            else
                                @info(string("Pull request: $(pr_number). ",
                                             "Type: $(pr_type). ",
                                             "Decision: do not merge. ",
                                             "Reason: ",
                                             "This is a new package pull request. ",
                                             "All of the criteria for automerging ",
                                             "were met. ",
                                             "However, merge_new_packages is false, ",
                                             "so I will not merge. ",
                                             "If merge_new_packages had been set to ",
                                             "true, I would have merged this ",
                                             "pull request right now."))
                            end
                        else # it is a new version
                            always_assert(pr_type == :NewVersion)
                            always_assert(status_pr_type == :NewVersion)
                            if merge_new_versions
                                my_comment = comment_text_merge_now()
                                @info(string("Pull request: $(pr_number). ",
                                             "Type: $(pr_type). ",
                                             "Decision: merge now."))
                                my_retry(() -> merge!(registry, pr, passed_pr_head_sha; auth = auth))
                            else
                                @info(string("Pull request: $(pr_number). ",
                                             "Type: $(pr_type). ",
                                             "Decision: do not merge. ",
                                             "Reason: merge_new_versions is false",
                                             "This is a new version pull request. ",
                                             "All of the criteria for automerging ",
                                             "were met. ",
                                             "However, merge_new_versions is false, ",
                                             "so I will not merge. ",
                                             "If merge_new_versions had been set to ",
                                             "true, I would have merged this ",
                                             "pull request right now."))
                            end
                        end
                    else
                        @info(string("Pull request: $(pr_number). ",
                                     "Type: $(pr_type). ",
                                     "Decision: do not merge. ",
                                     "Reason: pull request has one or more blocking comments."))
                    end
                else
                    @info(string("Pull request: $(pr_number). ",
                                 "Type: $(pr_type). ",
                                 "Decision: do not merge. ",
                                 "Reason: automerge/decision status is not passing"),
                          whoami)
                end
            else
                @info(string("Pull request: $(pr_number). ",
                             "Type: $(pr_type). ",
                             "Decision: do not merge. ",
                             "Reason: mandatory waiting period has not elapsed."),
                      pr_type,
                      pr_age,
                      pkg,
                      is_jll_name(pkg),
                      new_package_waiting_period,
                      new_jll_package_waiting_period,
                      new_version_waiting_period,
                      new_jll_version_waiting_period)
            end
        else
            @info(string("Pull request: $(pr_number). ",
                         "Decision: do not merge. ",
                         "Reason: pull request is neither a new package nor a new version."),
                  title(pr))
        end
    else
        @info(string("Pull request: $(pr_number). ",
                     "Decision: do not merge. ",
                     "Reason: pull request author is not authorized to automerge."),
              pr_author,
              authorized_authors)
    end
    return nothing
end

ext: jl


text:
using AssigmentProblem
using Test

#test matrix
m = [10  2 10 10 10 10;
	  1 10 10 10 10 10;
	 10 10  3 10 10 10;
	 10 10 10 10  5 10;
	 10 10 10  4 10 10;
	 10 10 10 10 10  6];

#small matrix
sm = rand(Int,2,2);

#optimal combination
oc = [1,2,3,4,5,6];

@testset "AssigmentProblem.jl" begin
    println("Test Matrix");
    display(m);
    println();
    print("Possible combinations: ")
    println(factorial(6));
    println("Optimal Combination");
    println(oc);
    println();
    println("Brute Force Solution (test all the posible combinations)");	#Brute force
    solve_with_brute_force(sm);		#Precompiling
    @time bfs = solve_with_brute_force(m)[1]; 	#Solve
    @test bfs == oc 	#Test with optimal
   	println(bfs);   #Print
    println();
    println("Linear Programming Solution");		#Lineal programming
    solve_with_lineal_programming(sm); 		#precompiling
    @time lps = solve_with_lineal_programming(m)[1];	#solve
    @test lps == oc		#Test with optimal
    println(lps);
    println();

end


text:
# Copyright 2020 Google LLC
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

struct SmoothPlusProximalFunction
  smooth_value::Function
  smooth_gradient::Function
  proximal_value::Function
  proximal_operator::Function
  dim::Int64 # 0 indicates the problem size can change
end

function value(f::SmoothPlusProximalFunction, point::Vector{Float64})
  return f.smooth_value(point) + f.proximal_value(point)
end

"""
Create a differentiable function (with no proximal component).
"""
function differentiable_function(value::Function, gradient::Function, dim::Int64)
  function identity_proximal_operator(x::Vector, step_size::Float64)
    return x
  end
  function zero_proximal_value(x::Vector)
    return 0.0
  end
  return SmoothPlusProximalFunction(
    value,
    gradient,
    zero_proximal_value,
    identity_proximal_operator,
    dim,
  )
end

function create_linear_regression_function(
    data::RegressionData;
    regularizer::Float64=0.0,
  )
  n = size(data.coefficient_matrix, 1)
  function regression_value(x::Vector)
    residual = data.coefficient_matrix * x - data.target_values
    return 0.5 * sum(residual.^2) + 0.5 * regularizer * sum(x.^2)
  end
  function regression_gradient(x::Vector)
    residual = data.coefficient_matrix * x - data.target_values
    return data.coefficient_matrix' * residual + regularizer * x
  end
  return differentiable_function(
    regression_value,
    regression_gradient,
    size(data.coefficient_matrix, 2)
  )
end

mutable struct UnivariateFunction
  value::Function
  gradient::Function
end

mutable struct UnivariateProximalOperator
  value::Function
  prox::Function
end

function create_general_regression_function(;
    data::RegressionData,
    loss::UnivariateFunction,
    regularizer::UnivariateProximalOperator,
  )
  n = size(data.coefficient_matrix, 1)
  function regression_value(x::Vector)
    residual = data.coefficient_matrix * x - data.target_values
    return sum(loss.value.(residual))
  end
  function regression_gradient(x::Vector)
    residual = data.coefficient_matrix * x - data.target_values
    return data.coefficient_matrix' * loss.gradient.(residual)
  end
  function proximal_value(x::Vector)
    return sum(regularizer.value.(x))
  end
  function proximal_operator(x::Vector, step_size::Float64)
    return regularizer.prox.(x, step_size)
  end
  return SmoothPlusProximalFunction(
    regression_value,
    regression_gradient,
    proximal_value,
    proximal_operator,
    size(data.coefficient_matrix, 2),
  )
end

function create_squared_function(lambda::Float64)
  function squared(x::Float64)
    return 0.5 * lambda * x^2
  end
  function squared_gradient(x::Float64)
    return lambda * x
  end
  return UnivariateFunction(squared, squared_gradient)
end

function create_L1_proximal_operator(lambda::Float64)
  function value(x::Float64)
    return lambda * abs(x)
  end
  function prox(x::Float64, step_size::Float64)
    return sign(x) * max( abs(x) - step_size * lambda, 0.0)
  end
  return UnivariateProximalOperator(value, prox)
end

function log_logistic_loss(x::Float64)
  return log(1 + exp(x))
end

function gradient_of_log_logistic_loss(x::Float64)
  return exp(x) / (1 + exp(x))
end

function create_logistic_regression_problem(;
    data::RegressionData,
    regularizer::UnivariateProximalOperator,
  )
  positive_label = maximum(data.target_values)
  negative_label = minimum(data.target_values)

  # identify the subsets of rows that correspond to the positive and negative
  # labels
  positive_labels = positive_label .== data.target_values
  negative_labels = negative_label .== data.target_values

  if sum(positive_labels) + sum(negative_labels) != length(data.target_values)
    error("More than two labels in data. This code can only handle two labels. ")
  end
  coefficient_matrix_positive_labels = data.coefficient_matrix[positive_labels,:]
  coefficient_matrix_negative_labels = data.coefficient_matrix[negative_labels,:]
  target_values_positive_labels = data.target_values[positive_labels,:]
  target_values_negative_labels = data.target_values[negative_labels,:]

  function logistic_value(x::Vector)
    positive_residuals = coefficient_matrix_positive_labels * x
    negative_residuals = -coefficient_matrix_negative_labels * x

    return sum(log_logistic_loss.(positive_residuals)) +
      sum(log_logistic_loss.(negative_residuals))
  end

  function logistic_gradient(x::Vector)
    positive_residuals = coefficient_matrix_positive_labels * x
    positive_gradient = coefficient_matrix_positive_labels' * gradient_of_log_logistic_loss.(positive_residuals)
    negative_residuals = -coefficient_matrix_negative_labels * x
    negative_gradient = -coefficient_matrix_negative_labels' * gradient_of_log_logistic_loss.(negative_residuals)

    return (positive_gradient + negative_gradient)[:]
  end

  function proximal_value(x::Vector)
    return sum(regularizer.value.(x))
  end
  function proximal_operator(x::Vector, step_size::Float64)
    return regularizer.prox.(x, step_size)
  end
  return SmoothPlusProximalFunction(
    logistic_value,
    logistic_gradient,
    proximal_value,
    proximal_operator,
    size(data.coefficient_matrix, 2),
  )
end


text:
using Metaheuristics

if VERSION < v"0.7.0"
    using Base.Test
    srand(31415926534)
else
    using Test
end

# write your own tests here
@testset "Multi objective" begin

    function run_methods(problem)

        ff, bounds, pf = Metaheuristics.TestProblems.get_problem(problem)
        D = size(bounds, 2)

        # number of function evaluations
        f_calls = 0

        f(x) = begin
            f_calls += 1
            ff(x)
        end

        options = Options( seed = 1, iterations=10000, f_calls_limit = 25000)
        nobjectives = length(pf[1].f)
        npartitions = nobjectives == 2 ? 100 : 12

        methods = [
                SPEA2(options=options),
                SMS_EMOA(N = 50, n_samples=500, options=options),
                NSGA2(options=options),
                MOEAD_DE(gen_ref_dirs(nobjectives, npartitions), options=Options( seed = 1, iterations = 500)),
                NSGA3(options=options),
              ]

        for method in methods
            f_calls = 0
            result = ( optimize(f, bounds, method) )
            show(IOBuffer(), "text/html", result)
            show(IOBuffer(), "text/plain", result.population)
            show(IOBuffer(), "text/html", result.population)
            show(IOBuffer(), result.population[1])
            @test Metaheuristics.PerformanceIndicators.igd(result.population, pf) <= 0.2
            @test Metaheuristics.PerformanceIndicators.spacing(result) < 0.5
            @test Metaheuristics.PerformanceIndicators.covering(pf, result.population) <= 1.0
            @test Metaheuristics.PerformanceIndicators.covering(result, result) ≈ 0.0

            # number of function evaluations should be reported correctly
            @test f_calls == result.f_calls

            # test obtaining non-dominated solutions
            pf1 = pareto_front(result)
            pf2 = pareto_front(result.population)

            @test size(pf1, 1) == size(pf2,1) &&
                  Metaheuristics.PerformanceIndicators.igd(pf1, pf2) ≈ 0.0
        end
    end

    function run_methods_parallel_eval(problem)
        ff, bounds, pf = Metaheuristics.TestProblems.get_problem(problem)
        D = size(bounds, 2)

        # number of function evaluations
        f_calls = 0

        f(X) = begin
            f_calls += size(X,1)
            fx, gx, hx = ff(X[1,:])
            F = zeros(size(X,1), length(fx)); F[1,:] = fx;
            G = zeros(size(X,1), length(gx)); G[1,:] = gx;
            H = zeros(size(X,1), length(hx)); H[1,:] = hx;

            for i = 2:size(X,1)
                fx, gx, hx = ff(X[i,:])
                F[i,:] = fx
                G[i,:] = gx
                H[i,:] = hx
            end
            F, G, H
        end

        options = Options( seed = 1, iterations=10000, f_calls_limit = 5000, parallel_evaluation=true)

        methods = [
                NSGA3(options=options),
                SMS_EMOA(N = 50, n_samples=500, options=options),
                NSGA2(options=options),
                SPEA2(options=options),
              ]

        for method in methods
            f_calls = 0
            result = optimize(f, bounds, method)
            # number of function evaluations should be reported correctly
            @test f_calls == result.f_calls
        end
    end

    for problem in [:ZDT3, :DTLZ2, :MTP]
        run_methods(problem)
    end

    for problem in [:ZDT1, :MTP]
        # run_methods(problem)
        run_methods_parallel_eval(problem)
    end

end


text:
using Documenter, JetPackTransforms

makedocs(sitename = "JetPackTransforms", modules=[JetPackTransforms])

deploydocs(
    repo = "github.com/ChevronETC/JetPackTransforms.jl.git",
)


text:
using Surrogates
using AbstractGPs
using LinearAlgebra
using Flux
using Flux: @epochs
using Zygote
using PolyChaos
using Test
#using Zygote: @nograd
#=
#FORWARD
###### 1D ######
lb = 0.0
ub = 10.0
n = 5
x = sample(n,lb,ub,SobolSample())
f = x -> x^2
y = f.(x)

#Radials
my_rad = RadialBasis(x,y,lb,ub,x->norm(x),2)
g = x -> ForwardDiff.derivative(my_rad,x)
g(5.0)

#Kriging
p = 1.5
my_krig = Kriging(x,y,p)
g = x -> ForwardDiff.derivative(my_krig,x)
g(5.0)

#Linear Surrogate
my_linear = LinearSurrogate(x,y,lb,ub)
g = x -> ForwardDiff.derivative(my_linear,x)
g(5.0)

#Inverse distance
p = 1.4
my_inverse = InverseDistanceSurrogate(x,y,p,lb,ub)
g = x -> ForwardDiff.derivative(my_inverse,x)
g(5.0)

#Lobachevsky
n = 4
α = 2.4
my_loba = LobachevskySurrogate(x,y,α,n,lb,ub)
g = x -> ForwardDiff.derivative(my_loba,x)
g(5.0)

#Second order polynomial
my_second = SecondOrderPolynomialSurrogate(x,y,lb,ub)
g = x -> ForwardDiff.derivative(my_second,x)
g(5.0)

###### ND ######
lb = [0.0,0.0]
ub = [10.0,10.0]
n = 5
x = sample(n,lb,ub,SobolSample())
f = x -> x[1]*x[2]
y = f.(x)

#Radials
my_rad = RadialBasis(x,y,[lb,ub],z->norm(z),2)
g = x -> ForwardDiff.gradient(my_rad,x)
g([2.0,5.0])

#Kriging
theta = [2.0,2.0]
p = [1.9,1.9]
my_krig = Kriging(x,y,p,theta)
g = x -> ForwardDiff.gradient(my_krig,x)
g([2.0,5.0])

#Linear Surrogate
my_linear = LinearSurrogate(x,y,lb,ub)
g = x -> ForwardDiff.gradient(my_linear,x)
g([2.0,5.0])

#Inverse Distance
p = 1.4
my_inverse = InverseDistanceSurrogate(x,y,p,lb,ub)
g = x -> ForwardDiff.gradient(my_inverse,x)
g([2.0,5.0])

#Lobachevsky
alpha = [1.4,1.4]
n = 4
my_loba_ND = LobachevskySurrogate(x,y,alpha,n,lb,ub)
g = x -> ForwardDiff.gradient(my_loba_ND,x)
g([2.0,5.0])

#Second order polynomial
my_second = SecondOrderPolynomialSurrogate(x,y,lb,ub)
g = x -> ForwardDiff.gradient(my_second,x)
g([2.0,5.0])
=#

##############
### ZYGOTE ###
##############

############
#### 1D ####
############
lb = 0.0
ub = 3.0
n = 10
x = sample(n,lb,ub,SobolSample())
f = x -> x^2
y = f.(x)

#AbstractGP 1D
@testset "AbstractGP 1D" begin
    agp1D = AbstractGPSurrogate(x,y, gp=GP(SqExponentialKernel()), Σy=0.05)
    g = x -> agp1D'(x)
    g([2.0])
end

#Radials
@testset "Radials 1D" begin
    my_rad = RadialBasis(x,y,lb,ub,rad = linearRadial)
    g = x -> my_rad'(x)
    g(5.0)
end

# #Kriging
@testset "Kriging 1D" begin
    my_p = 1.5
    my_krig = Kriging(x,y,lb,ub,p=my_p)
    g = x -> my_krig'(x)
    g(5.0)
end

# #Linear Surrogate
@testset "Linear Surrogate" begin
    my_linear = LinearSurrogate(x,y,lb,ub)
    g = x -> my_linear'(x)
    g(5.0)
end

#Inverse distance
@testset "Inverse Distance" begin
    my_p = 1.4
    my_inverse = InverseDistanceSurrogate(x,y,lb,ub,p=my_p)
    g = x -> my_inverse'(x)
    g(5.0)
end

#Second order polynomial
@testset "Second Order Polynomial" begin
    my_second = SecondOrderPolynomialSurrogate(x,y,lb,ub)
    g = x -> my_second'(x)
    g(5.0)
end

#Lobachevsky
@testset "Lobachevsky" begin
    n = 4
    α = 2.4
    my_loba = LobachevskySurrogate(x,y,lb,ub, alpha = α, n = 4)
    g = x -> my_loba'(x)
    g(0.0)
end

#NN
@testset "NN" begin
    my_model = Chain(Dense(1,1), first)
    my_loss(x, y) = Flux.mse(my_model(x), y)
    my_opt = Descent(0.01)
    n_echos = 1
    my_neural = NeuralSurrogate(x,y,lb,ub,model=my_model,loss=my_loss,opt=my_opt,n_echos=1)
    g = x->my_neural'(x)
    g(3.4)
end

#Wendland
@testset "Wendland" begin
    my_wend = Wendland(x,y,lb,ub)
    g = x -> my_wend'(x)
    g(3.0)
end

# MOE and VariableFidelity for free because they are Linear combinations
# of differentiable surrogates

# #Polynomialchaos
@testset "Polynomial Chaos" begin
    n = 50
    x = sample(n,lb,ub,SobolSample())
    y = f.(x)
    my_poli = PolynomialChaosSurrogate(x,y,lb,ub)
    g = x -> my_poli'(x)
    g(3.0)
end

# #Gek
@testset "GEK" begin
    n = 10
    lb = 0.0
    ub = 5.0
    x = sample(n,lb,ub,SobolSample())
    f = x-> x^2
    y1 = f.(x)
    der = x->2*x
    y2 = der.(x)
    y = vcat(y1,y2)
    my_gek = GEK(x,y,lb,ub)
    g = x-> my_gek'(x)
    g(3.0)
end

################
###### ND ######
################

lb = [0.0,0.0]
ub = [10.0,10.0]
n = 5
x = sample(n,lb,ub,SobolSample())
f = x -> x[1]*x[2]
y = f.(x)

# AbstractGP ND
@testset "AbstractGPSurrogate ND" begin
    my_agp = AbstractGPSurrogate(x,y, gp=GP(SqExponentialKernel()), Σy=0.05)
    g = x ->Zygote.gradient(my_agp, x)
    #g([(2.0,5.0)])
    g((2.0,5.0))
 end

#Radials
@testset "Radials ND" begin
    my_rad = RadialBasis(x,y,lb,ub,rad = linearRadial, scale_factor = 2.1)
    g = x -> Zygote.gradient(my_rad,x)
    g((2.0,5.0))
end

# Kriging
@testset "Kriging ND" begin
    my_theta = [2.0,2.0]
    my_p = [1.9,1.9]
    my_krig = Kriging(x,y,lb,ub,p=my_p,theta=my_theta)
    g = x -> Zygote.gradient(my_krig,x)
    g((2.0,5.0))
end

# Linear Surrogate
@testset "Linear Surrogate ND" begin
    my_linear = LinearSurrogate(x,y,lb,ub)
    g = x -> Zygote.gradient(my_linear,x)
    g((2.0,5.0))
end

#Inverse Distance
@testset "Inverse Distance ND" begin
    my_p = 1.4
    my_inverse = InverseDistanceSurrogate(x,y,lb,ub,p=my_p)
    g = x -> Zygote.gradient(my_inverse,x)
    g((2.0,5.0))
end

#Lobachevsky not working yet weird issue with Zygote @nograd
#=
Zygote.refresh()
alpha = [1.4,1.4]
n = 4
my_loba_ND = LobachevskySurrogate(x,y,alpha,n,lb,ub)
g = x -> Zygote.gradient(my_loba_ND,x)
g((2.0,5.0))
=#

# Second order polynomial mutating arrays
@testset "SecondOrderPolynomialSurrogate ND" begin
    my_second = SecondOrderPolynomialSurrogate(x,y,lb,ub)
    g = x -> Zygote.gradient(my_second,x)
    g((2.0,5.0))
end

#NN
@testset "NN ND" begin
    my_model = Chain(Dense(2,1), first)
    my_loss(x, y) = Flux.mse(my_model(x), y)
    my_opt = Descent(0.01)
    n_echos = 1
    my_neural = NeuralSurrogate(x,y,lb,ub,model=my_model,loss=my_loss,opt=my_opt,n_echos=1)
    g = x -> Zygote.gradient(my_neural, x)
    g((2.0,5.0))
end

#wendland
@testset "Wendland ND" begin
    my_wend_ND = Wendland(x,y,lb,ub)
    g = x -> Zygote.gradient(my_wend_ND,x)
    g((2.0,5.0))
end
# #MOE and VariableFidelity for free because they are Linear combinations
# #of differentiable surrogates

# #PolynomialChaos
@testset "Polynomial Chaos ND" begin
    n = 50
    lb = [0.0,0.0]
    ub = [10.0,10.0]
    x = sample(n,lb,ub,SobolSample())
    f = x -> x[1]*x[2]
    y = f.(x)
    my_poli_ND = PolynomialChaosSurrogate(x,y,lb,ub)
    g = x -> Zygote.gradient(my_poli_ND,x)
    g((1.0,1.0))

    n = 10
    d = 2
    lb = [0.0,0.0]
    ub = [5.0,5.0]
    x = sample(n,lb,ub,SobolSample())
    f = x -> x[1]^2 + x[2]^2
    y1 = f.(x)
    grad1 = x -> 2*x[1]
    grad2 = x -> 2*x[2]
    function create_grads(n,d,grad1,grad2,y)
        c = 0
        y2 = zeros(eltype(y[1]),n*d)
        for i = 1:n
            y2[i+c] = grad1(x[i])
            y2[i+c+1] = grad2(x[i])
            c = c+1
        end
        return y2
    end
    y2 = create_grads(n,d,grad1,grad2,y)
    y = vcat(y1,y2)
    my_gek_ND = GEK(x,y,lb,ub)
    g = x -> Zygote.gradient(my_gek_ND,x)
    g((2.0,5.0))
end

# ###### ND -> ND ######

lb = [0.0, 0.0]
ub = [10.0, 2.0]
n = 5
x = sample(n,lb,ub,SobolSample())
f = x -> [x[1]^2, x[2]]
y = f.(x)

#NN
@testset "NN ND -> ND" begin
    my_model = Chain(Dense(2,2))
    my_loss(x, y) = Flux.mse(my_model(x), y)
    my_opt = Descent(0.01)
    n_echos = 1
    my_neural = NeuralSurrogate(x,y,lb,ub,model=my_model,loss=my_loss,opt=my_opt,n_echos=1)
    Zygote.gradient(x -> sum(my_neural(x)), (2.0, 5.0))

    my_rad = RadialBasis(x,y,lb,ub,rad = linearRadial)
    Zygote.gradient(x -> sum(my_rad(x)), (2.0, 5.0))

    my_p = 1.4
    my_inverse = InverseDistanceSurrogate(x,y,lb,ub,p=my_p)
    my_inverse((2.0, 5.0))
    Zygote.gradient(x -> sum(my_inverse(x)), (2.0, 5.0))

    my_second = SecondOrderPolynomialSurrogate(x,y,lb,ub)
    Zygote.gradient(x -> sum(my_second(x)), (2.0, 5.0))
end


text:
function project_relative_path(xs...; pkg = "StatisticsWithJulia")
    pkg_dir = Base.find_package(pkg) # pathof(StatisticsWithJulia)
    normpath(joinpath(dirname(dirname(pkg_dir)), xs...))
end


text:
projected_pos = PVTSolution(
    ECEF([4.0186911757808267e6, 427026.7592380331, 4.91826139491441e6]),
    -2.1332096090475585e7,
    DOP(2.346836665386327, 0, 0, 0, 0), 1:length(satellite_states), [], [])

prev_PVT = PVTSolution()
min_elevation =  0

@testset "Testing calc_PVT" begin
    user = calc_PVT(satellite_states, prev_PVT, min_elevation)
    @test user.pos ≈ projected_pos.pos
    @test user.receiver_time_correction ≈ projected_pos.receiver_time_correction
    @test get_gdop(user) ≈ get_gdop(projected_pos)

    out = false
    try
        calc_PVT(satellite_states[1:3], prev_PVT, min_elevation)
    catch e
        if typeof(e) == PositionVelocityTime.SmallData
            out = true
        end
    end
    @test out == true

    svs_bad = deepcopy(satellite_states)
    svs_bad[1].decoder_state.subframes_decoded = [0,0,0,0,0]
    svs_bad[2].decoder_state.subframes_decoded = [0,0,0,0,0]

    out = false
    try
        calc_PVT(svs_bad, prev_PVT, min_elevation)
    catch e
        if typeof(e) == PositionVelocityTime.SmallData
            out = true
        end
    end
    @test out == true
end


text:
using DataDrivenDiffEq
using OrdinaryDiffEq
using LinearAlgebra
using Plots
gr()

# Create a linear , time discrete system
function linear_discrete(du, u, p, t)
    du[1] = 0.9u[1]
    du[2] = 0.9u[2] + 0.1u[1]
end

u0 = [10.0; -2.0]
tspan = (0.0, 10.0)
prob = DiscreteProblem(linear_discrete, u0, tspan)
sol = solve(prob, FunctionMap())

# Create Approximation
approx = DMD(sol[:,:])

# Create the associated problem
prob_approx = DiscreteProblem(approx, u0, tspan)
approx_sol = solve(prob_approx, FunctionMap())

# Show solutions
plot(sol)
plot!(approx_sol)

# Show error
plot((sol .- approx_sol)')
norm(sol .- approx_sol) # ≈ 2.23e-14

# Eigen Decomposition
eigen(approx)
# Eigenvalues
scatter(eigvals(approx), xlim = (-1, 1), ylim = (-1,1))

# Adapt system with new measurements
function linear_discrete_2(du, u, p, t)
    du[1] = 0.9u[1]
    du[2] = 0.9u[2] + 0.2u[1]
end

# Solve the new system
prob2 = DiscreteProblem(linear_discrete_2, u0, (0.0, 20.0))
sol2 = solve(prob2, FunctionMap())

# Split the data
x = sol2[:,1:20]
y = sol2[:,2:21]

# Update our approximation with new measurements
update!(approx, x, y)

# Lets have a look at the operator, which moves near the true value
operator(approx)

# Add time continouos system
function linear(du, u, p, t)
    du[1] = -0.9*u[1] + 0.1*u[2]
    du[2] = -0.8*u[2]
end

prob_cont = ODEProblem(linear, u0, tspan)
sol_cont = solve(prob_cont, Tsit5())

plot(sol_cont)

X = sol_cont[:,:]
DX = sol_cont(sol_cont.t, Val{1})[:,:]

# Giving the method a time step (which should be sequentially sampled)
# Enables us to get the continouos representation
approx_cont = gDMD(X, DX)
generator(approx_cont)
approx_sys = ODEProblem(approx_cont, u0, tspan)
approx_sol = solve(approx_sys, Tsit5(),  saveat = sol_cont.t)
# Lets have a look at the solution
plot(sol_cont)
plot!(approx_sol)
# And the error
plot(abs.(sol_cont .- approx_sol)')
norm(sol_cont .- approx_sol) # ≈ 1e-13


text:
using Optim
using PyPlot

f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2

optimize(f, [0.0, 0.0], NelderMead())

function g!(storage, x)
    storage[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
    storage[2] = 200.0 * (x[2] - x[1]^2)
end

res = optimize(f, g!, [0.0, 0.0], LBFGS())

figure(1)
x = collect(Float16, range(-1,length=100,stop=1));
y = collect(Float16, range(-1,length=100, stop=1));
z = hcat(x, y); z = mapslices(f, z, dims = 2)
surf(x,y,z);

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
using Distributions

"""
Example of a computation of beta credible interval
"""
function betaCredibleInt()
    S = 47
    N = 100
    a = S+1
    b = (N-S)+1
    α = .05
    β = Beta(a, b)
    CI = quantile(β, [α/2, 1-α/2])
end


text:
function equals(obj::Object, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function event_occurred(obj::AbstractDetector, arg0::SpacecraftState, arg1::jboolean)
    return jcall(obj, "eventOccurred", Action, (SpacecraftState, jboolean), arg0, arg1)
end

function g(obj::AbstractDetector, arg0::SpacecraftState)
    return jcall(obj, "g", jdouble, (SpacecraftState,), arg0)
end

function get_class(obj::Object)
    return jcall(obj, "getClass", Class, ())
end

function get_handler(obj::AbstractDetector)
    return jcall(obj, "getHandler", EventHandler, ())
end

function get_max_check_interval(obj::AbstractDetector)
    return jcall(obj, "getMaxCheckInterval", jdouble, ())
end

function get_max_iteration_count(obj::AbstractDetector)
    return jcall(obj, "getMaxIterationCount", jint, ())
end

function get_threshold(obj::AbstractDetector)
    return jcall(obj, "getThreshold", jdouble, ())
end

function hash_code(obj::Object)
    return jcall(obj, "hashCode", jint, ())
end

function init(obj::AbstractDetector, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end

function is_forward(obj::AbstractDetector)
    return jcall(obj, "isForward", jboolean, ())
end

function notify(obj::Object)
    return jcall(obj, "notify", void, ())
end

function notify_all(obj::Object)
    return jcall(obj, "notifyAll", void, ())
end

function reset_state(obj::AbstractDetector, arg0::SpacecraftState)
    return jcall(obj, "resetState", SpacecraftState, (SpacecraftState,), arg0)
end

function to_string(obj::Object)
    return jcall(obj, "toString", JString, ())
end

function wait(obj::Object)
    return jcall(obj, "wait", void, ())
end

function wait(obj::Object, arg0::jlong)
    return jcall(obj, "wait", void, (jlong,), arg0)
end

function wait(obj::Object, arg0::jlong, arg1::jint)
    return jcall(obj, "wait", void, (jlong, jint), arg0, arg1)
end

function with_handler(obj::AbstractDetector, arg0::EventHandler)
    return jcall(obj, "withHandler", AbstractDetector, (EventHandler,), arg0)
end

function with_max_check(obj::AbstractDetector, arg0::jdouble)
    return jcall(obj, "withMaxCheck", AbstractDetector, (jdouble,), arg0)
end

function with_max_iter(obj::AbstractDetector, arg0::jint)
    return jcall(obj, "withMaxIter", AbstractDetector, (jint,), arg0)
end

function with_threshold(obj::AbstractDetector, arg0::jdouble)
    return jcall(obj, "withThreshold", AbstractDetector, (jdouble,), arg0)
end


text:
using LearnBase
import Base: cat, vcat, hcat

"""

"""
struct ArrayNode{A<:AbstractArray,C} <: AbstractNode
    data::A
    metadata::C
end

ArrayNode(data::AbstractMatrix) = ArrayNode(data, nothing)
# ArrayNode(data::AbstractNode, a...) = data

Flux.@functor ArrayNode

mapdata(f, x::ArrayNode) = ArrayNode(mapdata(f, x.data), x.metadata)

Base.ndims(x::ArrayNode) = Colon()
LearnBase.nobs(a::ArrayNode) = size(a.data, 2)
LearnBase.nobs(a::ArrayNode, ::Type{ObsDim.Last}) = nobs(a)

function reduce(::typeof(catobs), as::Vector{T}) where {T<:ArrayNode}
    data = reduce(catobs, [x.data for x in as])
    metadata = reduce(catobs, [a.metadata for a in as])
    ArrayNode(data, metadata)
end

function reduce(::typeof(vcat), as::Vector{T}) where {T<:ArrayNode}
    data = reduce(vcat, [a.data for a in as])
    metadata = as[1].metadata == nothing ? nothing : as[1].metadata
    ArrayNode(data, metadata)
end

# hcat and vcat only for ArrayNode
function Base.vcat(as::ArrayNode...)
    data = vcat([a.data for a in as]...)
    metadata = as[1].metadata == nothing ? nothing : as[1].metadata
    ArrayNode(data, metadata)
end

Base.hcat(as::ArrayNode...) = reduce(catobs, collect(as))

Base.getindex(x::ArrayNode, i::VecOrRange) = ArrayNode(subset(x.data, i), subset(x.metadata, i))

Base.hash(e::ArrayNode{A,C}, h::UInt) where {A,C} = hash((A, C, e.data, e.metadata), h)
Base.:(==)(e1::ArrayNode{A,C}, e2::ArrayNode{A,C}) where {A,C} = e1.data == e2.data && e1.metadata == e2.metadata


text:
#first add workers
using Distributed
addprocs(5)
nprocs()

#load packages into all worker
using MendelIHT
using SnpArrays
using DataFrames
using Distributions
using BenchmarkTools
using Random
using LinearAlgebra
using BenchmarkTools
using GLM
using DelimitedFiles

function run_cv(n :: Int64, p :: Int64, k :: Int64, debias :: Bool, d::UnionAll, l::Link)
    #construct snpmatrix, covariate files, and true model b
    x, = simulate_random_snparray(n, p, undef)
    xbm = SnpBitMatrix{Float64}(x, model=ADDITIVE_MODEL, center=true, scale=true);
    z = ones(n, 1) # the intercept

    # simulate response, true model b, and the correct non-0 positions of b
    y, true_b, correct_position = simulate_random_response(x, xbm, k, d, l)

    #specify path and folds
    path = collect(1:20)
    num_folds = 5
    folds = rand(1:num_folds, size(x, 1))

    #compute cross validation
    result = @timed cv_iht_distributed(d(), l, x, z, y, 1, path, folds, num_folds, use_maf=false, debias=debias, parallel=true);
    dr = result[1]
    runtime = result[2] #seconds
    memory = result[3] / 1e6  #MB

    return dr, runtime, memory
end

function run()
    # set simulation parameters and seed
    repeats = 30
    n = 5000
    p = 50000
    k = 10 # number of true predictors
    d = Normal
    l = canonicallink(d())
    debias = false

    Random.seed!(2019)

    # initialize variables then run
    drs = zeros(20, repeats)
    run_time = zeros(repeats)
    memory = zeros(repeats)
    for i in 1:repeats
        dr, t, m = run_cv(n, p, k, debias, d, l)
        drs[:, i] .= dr
        run_time[i] = t
        memory[i] = m
    end

    writedlm("normal_nodebias_cv_drs", drs)
    writedlm("normal_nodebias_cv_run_times", run_time)
    writedlm("normal_nodebias_cv_memory", memory)
end

run()


text:
using JointMoments
using StatsBase
using Base.Test

include("data.jl")
include("precomputed.jl")
include("tolerance.jl")

ε = 1e-4


text:
module PartitionedGraphs

using SparseArrays
using LinearAlgebra

const Partition = Vector{Vector{Int}}
const WeightMatrix = SparseMatrixCSC{Float64,Int}

struct PartitionedGraph
    # graph (immutable)
    edge_weight::WeightMatrix
    cardinality::Vector{Int}

    # partition (mutable)
    partition::Partition
    size::Vector{Int}
    membership::Vector{Int}
end

function PartitionedGraph(
        adjmat::AbstractMatrix{<:Real};
        cardinality::Vector{Int} = create_unit_cardinality(nv(adjmat)),
        partition::Partition = create_singleton_partition(nv(adjmat)),)
    check_adjacent_matrix(adjmat)
    n = nv(adjmat)
    check_cardinality(n, cardinality)
    check_partition(n, partition)
    size = zeros(Int, length(partition))
    membership = zeros(Int, n)
    for (i, community) in enumerate(partition)
        size[i] = weighted_sum(community, cardinality)
        membership[community] .= i
    end
    return PartitionedGraph(adjmat, cardinality, partition, size, membership)
end

# number of vertices
nv(A::AbstractMatrix) = size(A, 1)
nv(graph::PartitionedGraph) = nv(graph.edge_weight)

# number of communities
nc(graph::PartitionedGraph) = length(graph.partition)

# neighbors of node `u`
function neighbors(graph::PartitionedGraph, u::Int)
    A = graph.edge_weight
    return @view A.rowval[A.colptr[u]:A.colptr[u+1]-1]
end

# communities connected with `u`
function connected_communities(graph::PartitionedGraph, u::Int)
    A = graph.total_weight
    return @view A.rowval[A.colptr[u]:A.colptr[u+1]-1]
end

# generators of the default parameter
create_unit_cardinality(n::Integer) = fill(1, n)
create_singleton_partition(n::Integer) = [[i] for i in 1:n]

function check_adjacent_matrix(adjmat::AbstractMatrix)
    m, n = size(adjmat)
    if m != n
        throw(ArgumentError("invalid adjacent matrix: not a square matrix"))
    end
    if n == 0
        throw(ArgumentError("invalid adjacent matrix: empty matrix"))
    end
    if !issymmetric(adjmat)
        throw(ArgumentError("invalid adjacent matrix: not an symmetric matrix"))
    end
    if adjmat isa SparseMatrixCSC
        if any(v < 0 for v in adjmat.nzval)
            throw(ArgumentError("invalid adjacent matrix: found negative weight(s)"))
        end
    else
        if any(v < 0 for v in adjmat)
            throw(ArgumentError("invalid adjacent matrix: found negative weight(s)"))
        end
    end
    return nothing
end

function check_cardinality(n::Integer, cardinality::Vector{Int})
    if n != length(cardinality)
        throw(ArgumentError("invalid cardinality: mismatching length"))
    end
    if any(x ≤ 0 for x in cardinality)
        throw(ArgumentError("invalid cardinality: found non-positive value"))
    end
    return nothing
end

function check_partition(n::Integer, partition::Partition)
    found = BitSet()
    for community in partition, u in community
        if u ∈ found
            throw(ArgumentError("invalid partition: found duplicated node"))
        end
        if !(1 ≤ u ≤ n)
            throw(ArgumentError("invalid partition: found out-of-bounds node"))
        end
        push!(found, u)
    end
    if length(found) != n
        throw(ArgumentError("invalid partition: found missing node"))
    end
    return nothing
end

# Move node `u` to `dst`.
function move_node!(graph::PartitionedGraph, (u, dst)::Pair{Int,Int})
    @assert 1 ≤ dst ≤ nc(graph) + 1
    src = graph.membership[u]
    if src == dst
        # no movement
        return graph
    end
    cardinality = graph.cardinality[u]
    community_src = graph.partition[src]
    pos = findfirst(isequal(u), community_src)::Int
    if dst > nc(graph)
        community_dst = Int[]
        push!(graph.partition, community_dst)
        push!(graph.size, 0)
    else
        community_dst = graph.partition[dst]
    end
    deleteat!(community_src, pos)
    push!(community_dst, u)
    graph.size[src] -= cardinality
    graph.size[dst] += cardinality
    graph.membership[u] = dst
    return graph
end

# Drop empty communities from the graph.
function drop_empty_communities!(graph::PartitionedGraph)
    empty = Int[]
    for (i, community) in enumerate(graph.partition)
        @assert isempty(community) == (graph.size[i] == 0)
        if isempty(community)
            push!(empty, i)
        end
    end
    deleteat!(graph.partition, empty)
    deleteat!(graph.size, empty)
    for (i, community) in enumerate(graph.partition)
        graph.membership[community] .= i
    end
    return graph
end

function reset_partition!(graph::PartitionedGraph, partition::Partition)
    check_partition(nv(graph), partition)
    empty!(graph.partition)
    empty!(graph.size)
    for (i, community) in enumerate(partition)
        push!(graph.partition, copy(community))
        push!(graph.size, weighted_sum(community, graph.cardinality))
        graph.membership[community] .= i
    end
    return graph
end

function weighted_sum(xs::Vector{Int}, weights::Vector{Int})
    return isempty(xs) ? 0 : sum(weights[x] for x in xs)
end

end


text:
# Objects here are defined as scratchspace for recurring computations when
# fitting models. The naming indicates the dimensions, e.g. _N, _N2 etc are
# vectors of size N (several as sometimes we need non-overlapping vectors)
# _NC are N*C matrices used in multiclass settings
#
# For now all of these are Float64 scratch spaces.

const SCRATCH_N   = Ref(zeros(0))
const SCRATCH_N2  = Ref(zeros(0))
const SCRATCH_N3  = Ref(zeros(0))
const SCRATCH_P   = Ref(zeros(0))
const SCRATCH_NC  = Ref(zeros(0,0))
const SCRATCH_NC2 = Ref(zeros(0,0))
const SCRATCH_NC3 = Ref(zeros(0,0))
const SCRATCH_NC4 = Ref(zeros(0,0))
const SCRATCH_PC  = Ref(zeros(0,0))

allocate(n, p, c=0) = begin
    SCRATCH_N[]  = zeros(n)
    SCRATCH_N2[] = zeros(n)
    SCRATCH_N3[] = zeros(n)
    SCRATCH_P[]  = zeros(p)
    if !iszero(c)
        SCRATCH_NC[]  = zeros(n, c)
        SCRATCH_NC2[] = zeros(n, c)
        SCRATCH_NC3[] = zeros(n, c)
        SCRATCH_NC4[] = zeros(n, c)
        SCRATCH_PC[]  = zeros(p, c)
    end
end

deallocate() = begin
    SCRATCH_N[]   = zeros(0)
    SCRATCH_N2[]  = zeros(0)
    SCRATCH_N3[]  = zeros(0)
    SCRATCH_P[]   = zeros(0)
    SCRATCH_NC[]  = zeros(0,0)
    SCRATCH_NC2[] = zeros(0,0)
    SCRATCH_NC3[] = zeros(0,0)
    SCRATCH_NC4[] = zeros(0,0)
    SCRATCH_PC[]  = zeros(0,0)
end


text:
T1 = [LinearInterpolator,
      CubicInterpolator,
      CubicSplineInterpolator]
T2 = [BilinearInterpolator,
      BicubicInterpolator]

for T in T1
    ϕ = T(x->1.0, -1, 1, 5)
    @test_throws ErrorException ϕ(-2)
    @test_throws ErrorException ϕ(2)
    ϕ = T(x->1.0, -1, 1, 5, NoBoundaries())
    @test isreal(ϕ(-2))
    @test isreal(ϕ(2))
    ϕ = T(x->1.0, -1, 1, 5, WeakBoundaries())
    @test isreal(ϕ(-1 - 1e-10))
    @test isreal(ϕ(1 + 1e-10))
end

for T in T2
    Φ = T((x,y)->1.0, -1, 1, 5, -1, 1, 5)
    @test_throws ErrorException Φ(-2, 0)
    @test_throws ErrorException Φ(0, 2)
    @test_throws ErrorException Φ(2, 0)
    @test_throws ErrorException Φ(0, -2)
    Φ = T((x,y)->1.0, -1, 1, 5, -1, 1, 5, NoBoundaries())
    @test isreal(Φ(-2, 0))
    @test isreal(Φ(0, 2))
    @test isreal(Φ(2, 0))
    @test isreal(Φ(0, -2))
    Φ = T((x,y)->1.0, -1, 1, 5, -1, 1, 5, WeakBoundaries())
    @test isreal(Φ(-1 - 1e-10, 0))
    @test isreal(Φ(0, 1 + 1e-10))
    @test isreal(Φ(1 + 1e-10, 0))
    @test isreal(Φ(0, -1 - 1e-10))
end


text:
module NGG
export to_type, from_type, show_repr, TypeLevel, TVal, TApp, TCons, TNil
export RuntimeFn, Unset, Argument, mkngg, rmlines
using MLStyle
using CanonicalTraits
using DataStructures
List = LinkedList

rmlines(ex::Expr) = begin
    hd = ex.head
    tl = map(rmlines, filter(!islinenumbernode, ex.args))
    Expr(hd, tl...)
end
rmlines(@nospecialize(a)) = a
islinenumbernode(@nospecialize(x)) = x isa LineNumberNode

include("typeable.jl")
include("runtime_fns.jl")

function vectolist(x::Vector{T}) where T
    foldr(x, init=nil(T)) do e, last
        Cons{T}(e, last)
    end
end

"""
julia> using .NGG
julia> mkngg(
           :f, #fname
           [
               Argument(:a, nothing, Unset()),
               Argument(:b, nothing, Unset())
           ],  # args
           Argument[], # kwargs
           :(a + b) #expression
       )
f = (a, b;) -> a + b

julia> ans(1, 2)
3
"""
function mkngg(
    name::Symbol,
    args::Vector{Argument},
    kwargs::Vector{Argument},
    @nospecialize(ex)
)
    arglist = vectolist(args)
    Args = to_type(arglist)

    kwarglist = vectolist(kwargs)
    Kwargs = to_type(kwarglist)
    Ex = to_type(ex)
    RuntimeFn{Args, Kwargs, Ex, name}()
end

end


text:
@testset "Norms" begin
    x = [2.0im, 3.0 - 1im, 5.0 + 2.0im]
    y = [-2.0im, 3.0 - 1im, 5.0 + 2.0im]
    norm = HC2.InfNorm()
    @test norm isa HC2.AbstractNorm
    @test norm(x) ≈ abs(5.0 + 2.0im)
    @test norm(x, x) ≈ 0.0
    @test norm(x, y) ≈ abs(4im)

    # test overflow check
    huge_x = exp2(700) .* x
    huge_y = exp2(700) .* y
    @test norm(huge_x) ≈ exp2(700) * abs(5.0 + 2.0im)
    @test norm(huge_x, huge_x) ≈ 0.0
    @test norm(huge_x, huge_y) ≈ exp2(700) * abs(4im)

    weighted_norm = HC2.WeightedNorm(norm, x)
    @test HC2.weights(weighted_norm) == ones(3)
    weighted_norm .= 2.0
    @test HC2.weights(weighted_norm) == [2.0, 2, 2]
    weighted_norm[1] = 4.0
    @test weighted_norm isa HC2.AbstractNorm
    @test weighted_norm(x) ≈ 0.5 * abs(5 + 2.0im)
    @test weighted_norm(x, x) ≈ 0.0
    @test weighted_norm(x, y) ≈ 0.25 * abs(4im)
    @test !isempty(sprint(show, weighted_norm))

    weighted_norm = HC2.WeightedNorm(norm, huge_x)
    weighted_norm .= 2.0
    weighted_norm[1] = 4.0
    @test weighted_norm isa HC2.AbstractNorm
    @test weighted_norm(huge_x) ≈ exp2(700) * 0.5 * abs(5 + 2.0im)
    @test weighted_norm(huge_x, huge_x) ≈ exp2(700) * 0.0
    @test weighted_norm(huge_x, huge_y) ≈ exp2(700) * 0.25 * abs(4im)
end


text:
# boxplots w/ outliers

# Makie

using CairoMakie

xs = rand(1:3, 1000)
ys = randn(1000)
dodge = rand(1:2, 1000)

CairoMakie.boxplot(xs, ys, dodge = dodge, show_notch = true, color = dodgeb)

# Jl Plots

using StatsPlots
gr(leg = false, bg = :honeydew)
Plots.plot()
Plots.plot(Plots.boxplot!(rand(1:4, 1000), randn(1000), color=:green))

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
#=
    Contents: Functions common to all backends.
=#

function create_user!(username::AbstractString, password::AbstractString, roles = Set{AbstractString}())
    create_user!(config[:acdata], username, password, roles)
end

function delete_user!(username::AbstractString)
    delete_user!(config[:acdata], username)
end

function add_sessionid_to_user!(username::AbstractString, session_id::AbstractString)
    acdata = config[:acdata]
    c1     = acdata != nothing
    c2     = length(session_id) == config[:session][:id_length]
    c1 && c2 && add_sessionid_to_user!(acdata, username, session_id)
end

function remove_sessionid_from_user!(username::AbstractString, session_id::AbstractString)
    acdata = config[:acdata]
    c1     = acdata != nothing
    c2     = length(session_id) == config[:session][:id_length]
    c1 && c2 && remove_sessionid_from_user!(config[:acdata], username, session_id)
end

function get_n_sessions(username::AbstractString)
    acdata = config[:acdata]
    acdata == nothing ? 0 : get_n_sessions(acdata, username)
end

function set_password!(username::AbstractString, password::AbstractString)
    acdata = config[:acdata]
    acdata == nothing && error("You haven't specified a backend for access control data.")
    set_password!(config[:acdata], username, password)
end

function get_salt_hashedpwd(username::AbstractString)
    get_salt_hashedpwd(config[:acdata], username)
end

function add_role!(username::AbstractString, role::AbstractString)
    add_roles!(username, role)
end

function remove_role!(username::AbstractString, role::AbstractString)
    remove_roles!(username, role)
end

function add_roles!(username::AbstractString, roles...)
    add_roles!(config[:acdata], username, roles...)
end

function remove_roles!(username::AbstractString, roles...)
    remove_roles!(config[:acdata], username, roles...)
end

"Returns: True if username has the role, false otherwise."
function has_role(username::AbstractString, role::AbstractString)
    has_role(config[:acdata], username, role)
end

"""
Returns: role for username.

ASSUMES: There is only one role.
"""
function get_role(username::AbstractString)
    get_role(config[:acdata], username)
end

# EOF


text:
using CUDA
using MOCNeutronTransport
using BenchmarkTools
using Test

# Number of points to use in vectors
N = 2^20
println("Using Point arrays of length $N")

# Check num threads and give warning
nthreads = Threads.nthreads()
if nthreads === 1
    @warn "Only using single-thread for cpu. Try restarting julia with 'julia --threads n'"
else
    println("Using $nthreads threads for CPU multi-threading")
end

# Single threads CPU add
function sequential_add!(x, y, z)
    for i in eachindex(x, y)
        @inbounds z[i] = x[i] + y[i]
    end
    return nothing
end

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

# Single threaded GPU add
function gpu_1_thread_add!(x, y, z)
    for i = 1:length(x)
        @inbounds z[i] = x[i] + y[i]
    end
    return nothing
end

function bench_gpu_1_thread_add!(x, y, z)
    CUDA.@sync begin
        @cuda gpu_1_thread_add!(x, y, z)
    end
end

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

function bench_gpu_1_block_add!(x, y, z)
    CUDA.@sync begin
        @cuda threads=512 gpu_1_block_add!(x, y, z)
    end
end

# Multiple blocks GPU add
function gpu_multiblock_add!(x, y, z)
    index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = gridDim().x * blockDim().x
    for i = index:stride:length(y)
        @inbounds z[index] = x[index] + y[index]
    end
    return nothing
end

function bench_gpu_multiblock_add!(x, y, z)
    numblocks = ceil(Int, length(x)/512)
    CUDA.@sync begin
        @cuda threads=512 blocks=numblocks gpu_multiblock_add!(x, y, z)
    end
end

# Use the occupancy API to determine threads and blocks to saturate the GPU
function bench_gpu_multiblock_autooccupancy!(x, y, z)
    kernel = @cuda launch=false gpu_multiblock_add!(x, y, z)
    config = launch_configuration(kernel.fun)
    threads = min(length(y), config.threads)
    blocks = cld(length(y), threads)

    CUDA.@sync begin
        kernel(x, y, z; threads, blocks)
    end
end

cpu_time = 0.0
for T = [Float64, Float32]
    println("Using Point of type $T")
    x = fill(Point_2D{T}(1, 1), N)
    y = fill(Point_2D{T}(2, 2), N)
    z = fill(Point_2D{T}(0, 0), N)
    time = @belapsed sequential_add!($x, $y, $z)
    μs = 1e6*time
    if T == Float64
        global cpu_time = μs
    end
    speedup = cpu_time/μs
    println("    CPU: single-thread = $μs μs")
    println("         Speed up compared to single-thread CPU & Float64 Points = $speedup")
    @test all(z .== Point_2D{T}(3, 3))

    fill!(z, Point_2D{T}(0, 0))
    time = @belapsed parallel_add!($x, $y, $z)
    μs = 1e6*time
    speedup = cpu_time/μs
    println("    CPU: $nthreads threads = $μs μs.")
    println("         Speed up compared to single-thread CPU & Float64 Points = $speedup")
    @test all(z .== Point_2D{T}(3, 3))

    x_d = CUDA.fill(Point_2D{T}(1, 1), N)
    y_d = CUDA.fill(Point_2D{T}(2, 2), N)
    z_d = CUDA.fill(Point_2D{T}(0, 0), N)

    time = @belapsed bench_gpu_1_thread_add!($x_d, $y_d, $z_d)
    μs = 1e6*time
    speedup = cpu_time/μs
    println("    GPU: single-thread/block, 1 blocks = $μs μs.")
    println("         Speed up compared to single-thread CPU & Float64 Points = $speedup")
    @test all(Array(z_d) .== Point_2D{T}(3, 3))

    fill!(z_d, Point_2D{T}(0, 0))
    time = @belapsed bench_gpu_1_block_add!($x_d, $y_d, $z_d)
    μs = 1e6*time
    speedup = cpu_time/μs
    println("    GPU: 512 threads/block, 1 blocks = $μs μs.")
    println("         Speed up compared to single-thread CPU & Float64 Points = $speedup")
    @test all(Array(z_d) .== Point_2D{T}(3, 3))

    fill!(z_d, Point_2D{T}(0, 0))
    time = @belapsed bench_gpu_multiblock_add!($x_d, $y_d, $z_d)
    μs = 1e6*time
    speedup = cpu_time/μs
    numblocks = ceil(Int64, N/512)
    println("    GPU: 512 threads/block, $numblocks blocks = $μs μs.")
    println("         Speed up compared to single-thread CPU & Float64 Points = $speedup")
    @test all(Array(z_d) .== Point_2D{T}(3, 3))

    fill!(z_d, Point_2D{T}(0, 0))
    time = @belapsed bench_gpu_multiblock_autooccupancy!($x_d, $y_d, $z_d)
    μs = 1e6*time
    speedup = cpu_time/μs
    kernel = @cuda launch=false gpu_multiblock_add!(x_d, y_d, z_d)
    config = launch_configuration(kernel.fun)
    threads = min(N, config.threads)
    blocks = cld(N, threads)
    println("    GPU: $threads threads/block, $blocks blocks = $μs μs.")
    println("         Speed up compared to single-thread CPU & Float64 Points = $speedup")
    @test all(Array(z_d) .== Point_2D{T}(3, 3))
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
# functions on `AbstractFittable`s
function fit! end
function get_history end
function get_underlying end
function parse_functions! end
function predict end
function predict_proba end
function set_feature_contrasts! end
function set_max_epochs! end
function transform end

# functions for calculating metrics
function accuracy end
function auprc end
function aurocc end
function binary_brier_score end
function cohen_kappa end
function f1score end
function false_negative_rate end
function false_positive_rate end
function fbetascore end
function mean_squared_error end
function negative_predictive_value end
function positive_predictive_value end
function prcurve end
function precision end
function r2score end
function recall end
function roccurve end
function root_mean_squared_error end
function sensitivity end
function specificity end
function true_negative_rate end
function true_positive_rate end

# functions for generating plots
function plotlearningcurve end
function plotprcurve end
function plotroccurve end

# functions for saving to file and loading from file
function load_model end
function save_model end
function save_plot end


text:
using MIRTjim
using Documenter
using Literate

# based on:
# https://github.com/jw3126/UnitfulRecipes.jl/blob/master/docs/make.jl

# https://juliadocs.github.io/Documenter.jl/stable/man/syntax/#@example-block
ENV["GKSwstype"] = "100"

# generate tutorials and how-to guides using Literate
lit = joinpath(@__DIR__, "lit")
src = joinpath(@__DIR__, "src")
notebooks = joinpath(src, "notebooks")

ENV["GKS_ENCODING"] = "utf-8"

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

execute = true # Set to true for executing notebooks and documenter!
nb = false # Set to true to generate the notebooks
for (root, _, files) in walkdir(lit), file in files
    splitext(file)[2] == ".jl" || continue
    ipath = joinpath(root, file)
    opath = splitdir(replace(ipath, lit=>src))[1]
    Literate.markdown(ipath, opath, documenter = execute)
    nb && Literate.notebook(ipath, notebooks, execute = execute)
end

# Documentation structure
ismd(f) = splitext(f)[2] == ".md"
pages(folder) =
    [joinpath(folder, f) for f in readdir(joinpath(src, folder)) if ismd(f)]

isci = get(ENV, "CI", nothing) == "true"

format = Documenter.HTML(;
    prettyurls = isci,
    edit_link = "main",
#   canonical = "https://JeffFessler.github.io/MIRTjim.jl/stable",
#   assets = String[],
)

makedocs(;
    modules = [MIRTjim],
    authors = "Jeff Fessler and contributors",
    sitename = "MIRTjim.jl",
    format,
    pages = [
        "Home" => "index.md",
        "Methods" => "methods.md",
        "Examples" => pages("examples")
    ],
)

if isci
    deploydocs(;
        repo = "github.com/JeffFessler/MIRTjim.jl.git",
        devbranch = "main",
        devurl = "dev",
        versions = ["stable" => "v^", "dev" => "dev"],
        forcepush = true,
    #   push_preview = true,
        # see https://JeffFessler.github.io/MIRTjim.jl/previews/PR##
    )
else
    @warn "may need to: rm -r src/examples"
end


text:
import Base
using DataStructures

abstract type ShibaObject end
abstract type Func <: ShibaObject end

const NullableRealValueType = Union{Real,AbstractArray{<:Real},Nothing}
const NullableFuncType = Union{Func,Nothing}

mutable struct Variable <: ShibaObject
    data::NullableRealValueType
    creator::NullableFuncType
    grad::Union{Variable,AbstractArray{Variable},Nothing}
    generation::Int
    name::String
end

const NullableVariableArray = Union{AbstractArray{Variable},Nothing}

macro variable(var)
    if var isa Symbol
        name = string(var)
        esc(:($var = Variable($var, nothing, nothing, 0, $name)))
    elseif (var isa Expr) && (var.head == :(=))
        res = eval(var)
        name = string(var.args[1])
        esc(:($(var.args[1]) = Variable($res, nothing, nothing, 0, $name)))
    else
        error("ignore format.")
    end
end

macro DeepShiba_Func(ex)
    pushfirst!(ex.args[3].args, :(generation::Int))
    pushfirst!(ex.args[3].args, :(outputs::NullableVariableArray))
    pushfirst!(ex.args[3].args, :(inputs::NullableVariableArray))
    return esc(ex)
end

@DeepShiba_Func mutable struct Add <: Func end

@DeepShiba_Func mutable struct Sub <: Func end

@DeepShiba_Func mutable struct Neg <: Func end

@DeepShiba_Func mutable struct Mul <: Func end

@DeepShiba_Func mutable struct Div <: Func end

@DeepShiba_Func mutable struct Pow <: Func
    c::NullableRealValueType
end

_Add(x1, x2) =  Add([x1, x2], nothing, 0)(x1, x2)

_Sub(x1, x2) = Sub([x1, x2], nothing, 0)(x1, x2)

_Neg(x) = Neg([x, ], nothing, 0)(x)

_Mul(x1, x2) = Mul([x1, x2], nothing, 0)(x1, x2)

_Div(x1, x2) = Div([x1, x2], nothing, 0)(x1, x2)

_Pow(x, c) = Pow([x, ], nothing, 0, c)(x)

@inline forward(f::Add, x1, x2)  = x1 + x2
@inline forward(f::Sub, x1, x2) = x1 - x2
@inline forward(f::Neg, x) = -x
@inline forward(f::Mul, x1, x2) = x1 * x2
@inline forward(f::Div, x1, x2) = x1 / x2
@inline forward(f::Pow, x) = x^(f.c)

function backward(f::Add)::AbstractArray{Variable}
    gys = get_gys(f.outputs)
    return [1, 1] .* gys
end

function backward(f::Sub)::AbstractArray{Variable}
    gys = get_gys(f.outputs)
    return [1, -1] .* gys
end

function backward(f::Neg)::AbstractArray{Variable}
    gys = get_gys(f.outputs)
    return -gys
end

function backward(f::Mul)::AbstractArray{Variable}
    gys = get_gys(f.outputs)
    x1, x2 = f.inputs
    return [x2, x1] .* gys
end

function backward(f::Div)::AbstractArray{Variable}
    gys = get_gys(f.outputs)
    x1, x2 = f.inputs
    return [1 / x2, -x1 / (x2^2)] .* gys
end

function backward(f::Pow)::AbstractArray{Variable}
    gys = get_gys(f.outputs)
    x, c = f.inputs[1], f.c
    return [c * (x^(c - 1)),] .* gys
end

Base.:+(x1::Variable, x2::Variable) = _Add(x1, x2)
Base.:+(x1::Variable, x2) = _Add(x1, variable(x2))
Base.:+(x1, x2::Variable) = _Add(variable(x1), x2)

Base.:-(x1::Variable, x2::Variable) = _Sub(x1, x2)
Base.:-(x1::Variable, x2) = _Sub(x1, variable(x2))
Base.:-(x1, x2::Variable) = _Sub(variable(x1), x2)

Base.:*(x1::Variable, x2::Variable) = _Mul(x1, x2)
Base.:*(x1::Variable, x2) = _Mul(x1, variable(x2))
Base.:*(x1, x2::Variable) = _Mul(variable(x1), x2)

Base.:/(x1::Variable, x2::Variable) = _Div(x1, x2)
Base.:/(x1::Variable, x2) = _Div(x1, variable(x2))
Base.:/(x1, x2::Variable) = _Div(variable(x1), x2)

Base.:-(x::Variable) = _Neg(x)

Base.:^(x1::Variable, x2::Variable) = _Pow(x1, x2.data)
Base.:^(x1::Variable, x2) = _Pow(x1, x2)
Base.:^(x1, x2::Variable) = _Pow(variable(x1), x2.data)

Base.log(x::Variable) = _Log(x)
Base.size(x::Variable) = size(x.data)
Base.iterate(x::Variable) = iterate(x.data)