text:
"""
This file contains the code for the work by

Y. Chandak, S. Niekum, B. Castro da Silva, E. Learned-Miller, E. Brunskill, and P. S. Thomas.
Universal Off-Policy Evaluation.
In Advances in Neural Information Processing Systems, 2021

Code repo: https://github.com/yashchandak/UnO
"""

using Statistics
using Random
using Bootstrap
using BlackBoxOptim

include("./UnO_params.jl")
include("./UnO_utils.jl")

#######################################################################
# CDF Estimators
#######################################################################

function IS(data, props, ξG=0.0)
    """
    Extract returns from behavior data
    and their associated IS weights wrt eval policy

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    ξG: Control variate to ensure G is always non-positive or non-negative
    """

    # Initialize arrays to store returns and IS weights
    n, horizon, d = size(data)
    γ = props["gamma"]
    returns = fill(0.0, n)
    ISweights = fill(0.0, n)

    for i in 1:1:n
        ρ = 1.0
        G = - ξG                        # Initialize with control-variate, if any
        for j in 1:1:horizon
            ρ = ρ*data[i,j,1]           # Full trajectory IS
            G += γ^(j-1) * data[i,j,2]
        end
        returns[i] = G
        ISweights[i] = ρ
    end

    return returns, ISweights
end

function IS_CDF(data, props)
    """
    Immportance sampling CDF

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    """

    returns, ISweights = IS(data, props)                # Get IS ratios
    order = sortperm(returns)                           # Create order statistic map

    xs = returns[order]                                 # Location of steps
    ys = cumsum(ISweights[order]) / length(ISweights)   # Step heights

    return xs, ys
end

function WIS_CDF(data, props)
    """
    Weighted Immportance sampling CDF

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    """

    returns, ISweights = IS(data, props)             # Get IS ratios
    order = sortperm(returns)                        # Create order statistics map

    xs = returns[order]                              # Location of steps
    ys = cumsum(ISweights[order]) / sum(ISweights)   # Step heights

    return xs, ys
end

#######################################################################
# Guarnateed Coverage Bounds
#######################################################################

function UnO_CI_lower(estimates, δ=0.05)
    """
    Construct a lower bound for CDF for a given keypoint
    This uses the bound by Thomas et al. (2015) as a sub-routine

    Args
    ----

    estimates: Sequence of samples for a keypoint κ: ρ_i 1{G_i ≤ κ}
    ξG: Control variate to ensure G is always non-positive or non-negative
    δ: Acceptable failure rate
    """

    n = length(estimates)
    if any(x -> x<-0.00001, estimates)
        println("Bound will not work, items in array less than 0 (UnO CI Lower)")

        for item in x
            if x < 0.0
                print(x, ",")
            end
        end
    end;

    function lower_bound(D, c, n_post)
        Y = clamp.(D, 0.0, c)
        term1 = mean(Y)
        term2 = (7*c*log(2/δ)) / (3*(n_post - 1))
        term3 = sqrt( (2 * var(Y) * log(2/δ) )/ n_post) # Sample variance form from Maurer(2009)

        lb = term1 - term2 - term3                      # Thm1 from Thomas et al.  (2015)
        return lb
    end;

    #########################################################
    # Prepare data splits into dpre and dpost
    # Size based on the split-value suggested by Thomas et al.  (2015)
    n_pre = Int(floor(0.05 * n))
    n_post = n - n_pre

    perm_idxs = shuffle(collect(1:n))
    D_pre = estimates[perm_idxs[1:n_pre]]
    D_post = estimates[perm_idxs[n_pre + 1:n]]

    ##########################################################
    # Search for best truncation point using basic random search in 1D
    min_obs, max_obs = minimum(D_pre), maximum(D_pre)
    cs = (rand(100) .* (max_obs - min_obs)) .+ min_obs
    v_max, c_max = -Inf, cs[1]

    # Search for c that maximizes the lower bound
    for c in cs
        # Notice the use of n_post instead of n_pre.
        # See the work by Thomas et al.  (2015) for discussion.
        val = lower_bound(D_pre, c, n_post)
        if val > v_max
            v_max, c_max = val, c
        end;
    end;

    ########################################################
    # Bound by Thomas et al. (2015)
    lb = lower_bound(D_post, c_max, n_post)

    # Lower bound on a positive R.V. has should always be >= 0.
    lb = max(lb, 0)

    return lb
end;

function UnO_CI_upper(estimates, δ=0.05)
    """
    Construct a lower bound for CDF for a given keypoint
    This uses the bound by Thomas et al. (2015) as a sub-routine

    Args
    ----

    estimates: Sequence of samples for a keypoint κ: ρ_i 1{G_i ≤ κ}
    ξG: Control variate to ensure G is always non-positive or non-negative
    δ: Acceptable failure rate
    """

    n = length(estimates)
    if any(x -> x>0.00001, estimates)
        println("Bound will not work, items in array more than 0 (UnO CI Upper)")
        for item in x
            if x > 0.0
                print(x, ",")
            end
        end
    end;

    function upper_bound(D, c, n_post)
        Y = clamp.(D, c, 0.0)
        term1 = mean(Y)
        term2 = (7*abs(c)*log(2/δ)) / (3*(n_post - 1))
        term3 = sqrt( (2 * var(Y) * log(2/δ) )/ n_post) # Sample variance form from Maurer(2009)

        ub = term1 + term2 + term3                      # Thm1 from Thomas et al.  (2015)
        return ub
    end;

    #########################################################
    # Prepare data splits into dpre and dpost
    # Size based on the split-value suggested by Thomas et al.  (2015)
    n_pre = Int(floor(0.05 * n))
    n_post = n - n_pre

    perm_idxs = shuffle(collect(1:n))
    D_pre = estimates[perm_idxs[1:n_pre]]
    D_post = estimates[perm_idxs[n_pre + 1:n]]

    ##########################################################
    # Search for best truncation point using basic random search in 1D
    min_obs, max_obs = minimum(D_pre), maximum(D_pre)
    cs = (rand(100) .* (max_obs - min_obs)) .+ min_obs
    v_min, c_min = +Inf, cs[1]

    # Search for c that minimizes the upper bound
    for c in cs
        # Notice the use of n_post instead of n_pre.
        # See the work by Thomas et al.  (2015) for discussion.
        val = upper_bound(D_pre, c, n_post)
        # println(c, val)
        if val < v_min
            v_min, c_min = val, c
        end;
    end;

    ########################################################
    # Bound by Thomas et al. (2015)
    ub = upper_bound(D_post, c_min, n_post)

    # Upper bound on a negative R.V. should always be <= 0
    ub = min(ub, 0)

    return ub
end;

function CDF_band_with_KP_deltas(xs, ys, props, keypoints, deltas)
    """
    This function creates confidence intervals at the given keypoints
    Then under-completes and overcompletes it to get the entire confidence band

    Inputs:
    xs = x values for a CDF
    ys = y values for a CDF
    keypoints = location of the keypoints
    deltas = failure rate for each of the keypoints
    """

    n = length(xs)
    n_keys = length(keypoints)
    gmin, gmax = props["Gmin"], props["Gmax"]
    ys_pdf = decumsum(ys)

    temp = fill(0.0, n)
    # 2 extra positions for Gmin and Gmax
    bound_xs = fill(0.0, n_keys + 2)
    lb_ys = fill(0.0, n_keys + 2)       # Lower bounds
    ub_ys = fill(0.0, n_keys + 2)       # Upper bounds

    # Pre-specified bounds at the corner-points of CDF box
    bound_xs[1], bound_xs[n_keys+2] = gmin, gmax
    lb_ys[1] = 0
    ub_ys[n_keys+1] = 1
    lb_ys[n_keys+2] = 1
    ub_ys[n_keys+2] = 1

    for k in 1:1:n_keys
        key = keypoints[k]
        bound_xs[k+1] = key

        # Get lower bound at key
        idx = searchsortedlast(xs, key)
        temp[1:idx] = ys_pdf[1:idx] .* n        # Multiply by n to undo CDF normalization
        temp[idx+1:n] .= 0
        lb_ys[k+1] = UnO_CI_lower(temp, deltas[k]/2)

        # Get upper bound at key
        idx = searchsortedlast(xs, key)
        temp[1:idx] .= 0
        temp[idx+1:n] = - ys_pdf[idx+1:n] .* n      # Control variate form: ρ[1{G ≤ ν} - 1] = - ρ[1{G > ν}]

        ub_ys[k] = UnO_CI_upper(temp, deltas[k]/2) + 1  # Add back the expected control value
                                                    # Also, note ub_ys[k] instead of ub_ys[k+1]
                                                    # This left shift by 1 does overcompletion.
    end

    # Enforce monotonic consistency
    # Because CI used above has some randomness CI might not be monotonic
    # CI- for x2 > x1 can be enforced to be no lesser than the CI- for x1.
    # Similarly CI+ for all x1 < x2 can be enforced to be no greater than the CI+ for x2.
    for k in 1:1:n_keys
        lb_ys[k+1] = max(lb_ys[k], lb_ys[k+1])

        rev_i = n_keys - k + 1
        ub_ys[rev_i] = min(ub_ys[rev_i + 1], ub_ys[rev_i])
    end

    return bound_xs, lb_ys, ub_ys
end

function CDF_band(xs, ys, props, δ)
    """
    This function chooses equispaced keypoints where the CIs should be constructed.
    Then it creates CIs at those keypoints.
    It then under-completes and overcompletes it to get the confidence band.

    Inputs:
    xs = x values for a CDF
    ys = y values for a CDF
    keypoints = location of the keypoints
    deltas = failure rate for each of the keypoints
    """

    n_keys = ceil(Int32, log(length(xs)))
    # n_keys = ceil(Int32, sqrt(length(xs)))

    gmin, gmax = props["Gmin"], props["Gmax"]
    step = (gmax - gmin) / (n_keys + 1)

    keypoints = collect(gmin:step:gmax)[2:n_keys+1]     # Equi-spaced keypoints
    deltas = fill(δ/n_keys, n_keys)                 # Equal failure rates at each keypoint

    return CDF_band_with_KP_deltas(xs, ys, props, keypoints, deltas)
end

##################################################################
# Optimized CDF Bands with guaranteed coverage
##################################################################

function Optim_CI_lower(estimates, δ, n_post)
    """
    Construct a lower bound for CDF for a given keypoint
    This uses the bound by Thomas et al. (2015) as a sub-routine

        > This function is similar to CDF_CI_lower but the difference is
          that this function is meant to be used during the trainig phase, and hence
          does not split the data any further.

    Args
    ----

    estimates: Sequence of samples for a keypoint κ: ρ_i 1{G_i ≤ κ}
    δ: Acceptable failure rate
    """

    # Use all the data in the training mode
    # No further training and evaluation split

    if any(x -> x<-0.00001, estimates)
        println("Bound will not work, items in array less than 0 (Optim CI Lower)")
    end;

    function lower_bound(D, c, n_post)
        Y = clamp.(D, 0.0, c)
        term1 = mean(Y)
        term2 = (7*c*log(2/δ)) / (3*(n_post - 1))
        term3 = sqrt( (2 * var(Y) * log(2/δ) )/ n_post) # Sample variance form from Maurer(2009)

        lb = term1 - term2 - term3                      # Thm1 from Thomas et al.  (2015)
        return lb
    end;

    ##########################################################
    # Search for best truncation point using basic random search in 1D
    min_obs, max_obs = minimum(estimates), maximum(estimates)
    cs = (rand(100) .* (max_obs - min_obs)) .+ min_obs
    v_max, c_max = -Inf, cs[1]

    # Search for c that maximizes the lower bound
    for c in cs
        # Notice the use of n_post instead of n_pre.
        # See the work by Thomas et al.  (2015) for discussion.
        val = lower_bound(estimates, c, n_post)
        if val > v_max
            v_max, c_max = val, c
        end;
    end;

    # Lower bound on a positive R.V. should always be >= 0.
    lb = max(v_max, 0)

    return lb
end;

function Optim_CI_upper(estimates, δ, n_post)
    """
    Construct an upper bound for CDF for a given keypoint
    This uses the bound by Thomas et al. (2015) as a sub-routine

        > This function is similar to CDF_CI_lower but the difference is
          that this function is meant to be used during the trainig phase, and hence
          does not split the data any further.

    Args
    ----

    estimates: Sequence of samples for a keypoint κ: ρ_i 1{G_i ≤ κ}
    δ: Acceptable failure rate
    """

    if any(x -> x>0.00001, estimates)
        println("Bound will not work, items in array more than 0 (Optim CI Upper)")
    end;

    function upper_bound(D, c, n_post)
        Y = clamp.(D, c, 0.0)
        term1 = mean(Y)
        term2 = (7*abs(c)*log(2/δ)) / (3*(n_post - 1))
        term3 = sqrt( (2 * var(Y) * log(2/δ) )/ n_post) # Sample variance form from Maurer(2009)

        ub = term1 + term2 + term3                      # Thm1 from Thomas et al.  (2015)
        return ub
    end;

    ##########################################################
    # Search for best truncation point using basic random search in 1D
    min_obs, max_obs = minimum(estimates), maximum(estimates)
    cs = (rand(100) .* (max_obs - min_obs)) .+ min_obs
    v_min, c_min = +Inf, cs[1]

    # Search for c that minimizes the upper bound
    for c in cs
        # Notice the use of n_post instead of n_pre.
        # See the work by Thomas et al.  (2015) for discussion.
        val = upper_bound(estimates, c, n_post)
        # println(c, val)
        if val < v_min
            v_min, c_min = val, c
        end;
    end;

    # Upper bound on a negative R.V. should always be <= 0
    ub = min(v_min, 0)

    return ub
end;

function Optim_band_with_KP_deltas(xs, ys, props, keypoints, deltas, n_post)
    """
    This function creates confidence intervals at the given keypoints
    Then under-completes and overcompletes it to get the entire confidence band

        > This function is similar to CDF_band_with_KP_deltas but the difference is
          that this function is meant to be used during the trainig phase, and hence
          does not split the data any further for creating the upper and lower CIs

    Inputs:
    xs = x values for a CDF
    ys = y values for a CDF
    keypoints = location of the keypoints
    deltas = failure rate for each of the keypoints
    """

    # Obtain meta-properties
    n = length(xs)
    n_keys = length(keypoints)
    gmin, gmax = props["Gmin"], props["Gmax"]
    ys_pdf = decumsum(ys)

    # Initialize the bound variables
    temp = fill(0.0, n)
    bound_xs = fill(0.0, n_keys + 2)    # 2 extra positions for Gmin and Gmax
    lb_ys = fill(0.0, n_keys + 2)       # Lower bounds
    ub_ys = fill(0.0, n_keys + 2)       # Upper bounds

    # Pre-specified bounds at the corner-points of CDF box
    bound_xs[1], bound_xs[n_keys+2] = gmin, gmax
    lb_ys[1] = 0
    ub_ys[n_keys+1] = 1
    lb_ys[n_keys+2] = 1
    ub_ys[n_keys+2] = 1

    # FOr each keypoint get the upper and lower bounds
    for k in 1:1:n_keys
        key = keypoints[k]
        bound_xs[k+1] = key

        # Get lower bound at key
        idx = searchsortedlast(xs, key)
        # Estimates using samples less than key = ρ_i G_i, and estimates after key = 0
        temp[1:idx] = ys_pdf[1:idx] .* n        # Multiply by n to undo CDF normalization
        temp[idx+1:n] .= 0
        lb_ys[k+1] = Optim_CI_lower(temp, deltas[k]/2, n_post)

        # Get upper bound at key
        idx = searchsortedlast(xs, key)
        temp[1:idx] .= 0
        temp[idx+1:n] = - ys_pdf[idx+1:n] .* n      # Control variate form: ρ[1{G ≤ ν} - 1] = - ρ[1{G > ν}]
        ub_ys[k] = Optim_CI_upper(temp, deltas[k]/2, n_post) + 1    # Add back the expected control value
                                                                    # Also, note ub_ys[k] instead of ub_ys[k+1]
                                                                    # This left shift by 1 does overcompletion.
    end

    # Enforce monotonic consistency
    # Because CI used above has some randomness CI might not be monotonic
    # CI- for x2 > x1 can be enforced to be no lesser than the CI- for x1.
    # Similarly CI+ for all x1 < x2 can be enforced to be no greater than the CI+ for x2.
    for k in 1:1:n_keys
        lb_ys[k+1] = max(lb_ys[k], lb_ys[k+1])

        rev_i = n_keys - k + 1
        ub_ys[rev_i] = min(ub_ys[rev_i + 1], ub_ys[rev_i])
    end

    return bound_xs, lb_ys, ub_ys
end

function CDF_band_optimized(xs, ys, props, δ)
    """
    This function optmizes for the keypoints where the CIs should be constructed.
    Then it creates CIs at those keypoints.
    It then under-completes and overcompletes it to get the confidence band.

    IMPORTANT: If bounds for only a single parameter is required,
    then the following optimization procedure can be improved. Instead of optimizing for the area enlcosed,
    one can directly optimize for the desired parameter by changing the area_fn() below.
    See Appendix E.3 in the work by Chandak et al. (2021) for more details.

    Inputs:
    xs = x values for a CDF
    ys = y values for a CDF
    keypoints = location of the keypoints
    deltas = failure rate for each of the keypoints
    """

    # Obtain meta-properties
    n_keys = ceil(Int32, log(length(xs)))       # Number of keypoints to be used (hyper-parameter)
    # n_keys = ceil(Int32, sqrt(length(xs)))
    gmin, gmax = props["Gmin"], props["Gmax"]
    step = (gmax - gmin) / (n_keys + 1)

    # Split data into train and test set
    # Train set is used for doing all the optimization
    # Test set is used to get the actual CDF band
    n = length(xs)
    n_train = Int(floor(0.05 * n))
    n_post = n - n_train

    ys_unnorm_pdf = decumsum(ys) .* n       # Basically get back IS ratios

    perm_idxs = shuffle(collect(1:n))
    x_train = xs[perm_idxs[1:n_train]]
    y_train = ys_unnorm_pdf[perm_idxs[1:n_train]]

    x_post = xs[perm_idxs[n_train+1:n]]
    y_post = ys_unnorm_pdf[perm_idxs[n_train+1:n]]

    ##################
    # Training phase
    ##################

    order = sortperm(x_train)
    xs_train = x_train[order]
    ys_train = cumsum(y_train[order]) / n_train

    function constraints(candidate)
        """
        Helper function to project the candidate solution for location and the delta
        at each keypoint within the desired limits.

        (This is one choice. There could be other ways to enforce these constraints
        that make the optimization procedure better in some sense)
        """
        key_locs = candidate[1:length(candidate)÷2]
        key_dels = candidate[(length(candidate)÷2+1):length(candidate)]

        # Rescale individual leypoint locs to ensure they cover the desired range
        temp1 = cumsum(key_locs)
        if sum(key_locs) > (gmax - gmin)
            temp1 = gmin .+ (gmax - gmin) .* temp1 / sum(key_locs)
        else
            temp1 =  gmin .+ temp1
        end

        # Ensure sum of individual deltas is never more than δ
        temp2 = cumsum(key_dels)
        if sum(key_dels) > δ
            temp2 = δ .* temp2 / sum(temp2)
        end

        return temp1, temp2
    end

    function area_fn(candidate)
        """
        Helper function to compute the total area enclosed within the CDF band
        """

        # As this is training phase, use `optim' version of KP delta which does not create further data splits.
        candidate_keys, candidate_deltas = constraints(candidate)
        bound_xs, lb_ys, ub_ys = Optim_band_with_KP_deltas(xs_train, ys_train, props,
                                                           candidate_keys, candidate_deltas, n_post)

        # Compute the area within the CDF band
        area = 0
        for i in 1:1:length(bound_xs) - 1
            height = ub_ys[i] - lb_ys[i]
            width = bound_xs[i+1] - bound_xs[i]
            area = area + height * width
        end

        return area
    end

    # Use blackbox search to find a good set of keypoint locs and deltas
    # For blackbox search, define the range of domain
    # (Notice, through constrain projection, locs get rescaled to: gmin to gmax)
    lb_range = fill(0.0, n_keys*2)
    ub_range = fill(δ, n_keys*2)
    ub_range[1:n_keys] .= abs(gmax-gmin)
    res = bboptimize(area_fn, SearchRange = collect(zip(lb_range, ub_range)), NumDimensions = n_keys*2,
                                            method=:random_search, TraceMode=:silent)

    # Get the optimal key points
    best_keypoints, best_deltas = constraints(best_candidate(res))

    ########################
    # Final evaluation phase
    #########################

    order = sortperm(x_post)
    xs_post = x_post[order]
    ys_post = cumsum(y_post[order]) / n_post

    # For final eval, use non-optim version of KP_delta function such that multi-comparison is avoided
    return CDF_band_with_KP_deltas(xs_post, ys_post, props, best_keypoints, best_deltas)
end

#######################################################################
# Approximate Bounds Using Statistical Bootstrap
#######################################################################

function CDF_boot_band(xs, ys, props, delta)
    """
    This was a preliminary attempt at developing bootstrap based Off-policy CDF bounds.

    We tried a functional bootstrap based approach below, but it does not work well.
    Avoid using it. Instead, one can easily bound any parameter directly when using bootstrap,
    instead of bounding CDF as an intermediate step (which was helpful when using concentration inequalities)
    """

    function get_distance(x, y_pdf)
        order = sortperm(x)                 # Create order statistic map
        xs_boot = x[order]                  # Location of steps
        ys_boot = cumsum(y_pdf[order])      # Step heights
        return KS_distance(xs, ys, xs_boot, ys_boot)
    end

    n = length(xs)
    idxs = collect(1:n)
    ys_pdf = decumsum(ys)

    bs = bootstrap(boot_id -> get_distance(xs[boot_id], ys_pdf[boot_id]), idxs, BasicSampling(100))   # Bootstrap KS distance's distribution
    # bci = confint(bs, BCaConfInt(1 - 2*delta))            # BCa Bootstrap CI on the KS distance
    bci = confint(bs, PercentileConfInt(1 - 2*delta))       # Percentile Bootstrap CI on the KS distance
                                                            # Compute with 2*delta because it internally divides
                                                            # by 2 to get bounds for both the tails
                                                            # However, we only need the bound for the upper tail

    _, _, ub_dist = bci[1]      # return format: estimate, lower, upper

    # 2 extra positions for Gmin and Gmax
    gmin, gmax = props["Gmin"], props["Gmax"]
    bound_xs = fill(0.0, n+2)
    lb_ys = fill(0.0, n+2)       # Lower bounds
    ub_ys = fill(0.0, n+2)       # Upper bounds

    # Pre-specified bounds at the corner-points of CDF box
    bound_xs[1], bound_xs[n+2] = gmin, gmax
    bound_xs[2:n+1] = xs
    lb_ys[1] = 0
    ub_ys[n+1] = 1
    lb_ys[n+2] = 1
    ub_ys[n+2] = 1

    # Create the upper (over-complete) and lower (under-complete) CDF
    lb_ys[2:n+1] = ys .- ub_dist
    ub_ys[1:n] = ys .+ ub_dist      # Note ub_ys[1:n] instead of ub_ys[2:n+1]
                                    # This left shift by 1 does overcompletion.

    # Bounds can be clipped between 0 and 1
    clamp!(lb_ys, 0, 1)
    clamp!(ub_ys, 0, 1)

    return bound_xs, lb_ys, ub_ys
end

#######################################################################
# Approximate Bounds for Non-stationary Setting Using Wild Bootstrap
#######################################################################

function CDF_NS_band_with_KP_deltas(eps, returns, ISweights, props, keypoints, deltas, L, mode)
    """
    This function creates wild-bootstrap based confidence intervals at the given keypoints
    Then under-completes and overcompletes it to get the entire confidence band

    Inputs:
    xs = x values for a CDF
    ys = y values for a CDF
    keypoints = location of the keypoints
    deltas = failure rate for each of the keypoints
    """

    n = length(returns)
    n_keys = length(keypoints)
    gmin, gmax = props["Gmin"], props["Gmax"]

    temp = fill(0.0, n)
    # 2 extra positions for Gmin and Gmax
    bound_xs = fill(0.0, n_keys + 2)
    lb_ys = fill(0.0, n_keys + 2)       # Lower bounds
    ub_ys = fill(0.0, n_keys + 2)       # Upper bounds

    # Pre-specified bounds at the corner-points of CDF box
    bound_xs[1], bound_xs[n_keys+2] = gmin, gmax
    lb_ys[1] = 0
    ub_ys[n_keys+1] = 1
    lb_ys[n_keys+2] = 1
    ub_ys[n_keys+2] = 1

    for k in 1:1:n_keys
        key = keypoints[k]
        bound_xs[k+1] = key

        # Get lower bound at key
        idx = returns .<= key
        temp = ISweights .* idx
        lb_ys[k+1], ub_ys[k] = nswildbs_CI(eps, temp, deltas[k], L, mode)
        # lb_ys[k+1], ub_ys[k] = nswildbst_CI(eps, temp, deltas[k], n)
    end

    # Enforce monotonic consistency
    # Because CI used above has some randomness CI might not be monotonic
    # CI- for x2 > x1 can be enforced to be no lesser than the CI- for x1.
    # Similarly CI+ for all x1 < x2 can be enforced to be no greater than the CI+ for x2.
    for k in 1:1:n_keys
        lb_ys[k+1] = min(1, max(lb_ys[k], lb_ys[k+1]))

        rev_i = n_keys - k + 1
        ub_ys[rev_i] = max(0, min(ub_ys[rev_i + 1], ub_ys[rev_i]))
    end

    return bound_xs, lb_ys, ub_ys
end

function CDF_NS_band(data, props, δ)
    """
    This function constucts equi-spaced keypoints and then creates wild-bootstrap based bounds at those keypoints.
    It then under-completes and overcompletes it to get the confidence band.

        > This function is similar to CDF_band() but used wild-bootstrap instead of CIs

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    δ: Acceptable failure rate
    """

    # Immportance sampling CDF
    returns, ISweights = IS(data, props)                # Get IS ratios

    n_keys = ceil(Int32, log(length(returns)))
    # n_keys = ceil(Int32, sqrt(length(returns)))

    gmin, gmax = props["Gmin"], props["Gmax"]
    step = (gmax - gmin) / (n_keys + 1)

    keypoints = collect(gmin:step:gmax)[2:n_keys+1]     # Equi-spaced keypoints
    deltas = fill(δ/n_keys, n_keys)                 # Equal failure rates at each keypoint

    eps = collect(Float64, 1:length(returns))
    return CDF_NS_band_with_KP_deltas(eps, returns, ISweights, props, keypoints, deltas, last(eps), :normal)
end

function CDF_NS_band_optimized(data, props, δ)
    """
    This function optmizes for the keypoints where the CIs should be constructed.
    Then it creates wild-bootstrap based bounds at those keypoints.
    It then under-completes and overcompletes it to get the confidence band.

    IMPORTANT: If bounds for only a single parameter is required,
    then the following optimization procedure can be improved. Instead of optimizing for the area enlcosed,
    one can directly optimize for the desired parameter by changing the area_fn() below.
    See Appendix E.3 in the work by Chandak et al. (2021) for more details.

    Args
    ----

    data: Trajectories x Steps x (Importance ratio, Reward)
    prop: Meta properties of the data
    δ: Acceptable failure rate
    """

    # Importance sampling CDF
    returns, ISweights = IS(data, props)                # Get IS ratios
    n_keys = ceil(Int32, log(length(returns)))
    # n_keys = ceil(Int32, sqrt(length(returns)))

    # Obtain meta-properties
    gmin, gmax = props["Gmin"], props["Gmax"]
    step = (gmax - gmin) / (n_keys + 1)

    # Variables for data splitting
    n = length(returns)
    n_train = Int(floor(0.5 * n))
    n_post = n - n_train
    perm_idxs = shuffle(collect(1:n))

    #############################
    # Training phase
    # Similar to CDF_band_optimized,
    # but uses wild bootstrap instead of concentration inequalities
    #############################

    eps_train = perm_idxs[1:n_train]
    return_train = returns[eps_train]
    ISweights_train = ISweights[eps_train]

    function constraints(candidate)
        """
        Projects candidate solutions to the feasible/desired set.
        """

        key_locs = candidate[1:length(candidate)÷2]
        key_dels = candidate[(length(candidate)÷2+1):length(candidate)]

        temp1 = cumsum(key_locs)
        if sum(key_locs) > (gmax - gmin)
            temp1 = gmin .+ (gmax - gmin) .* temp1 / sum(key_locs)
        else
            temp1 =  gmin .+ temp1
        end

        temp2 = cumsum(key_dels)
        if sum(key_dels) > δ
            temp2 = δ .* temp2 / sum(temp2)
        end

        return temp1, temp2
    end

    function area_fn(candidate)
        """
        Area of the CDF band
        """

        candidate_keys, candidate_deltas = constraints(candidate)
        bound_xs, lb_ys, ub_ys = CDF_NS_band_with_KP_deltas(eps_train, return_train, ISweights_train, props, candidate_keys, candidate_deltas, n, :fast)

        area = 0
        for i in 1:1:length(bound_xs) - 1
            height = ub_ys[i] - lb_ys[i]
            width = bound_xs[i+1] - bound_xs[i]
            area = area + height * width
        end

        return area
    end

    lb_range = fill(0.0, n_keys*2)
    ub_range = fill(δ, n_keys*2)
    ub_range[1:n_keys] .= abs(gmax-gmin)
    res = bboptimize(area_fn, SearchRange = collect(zip(lb_range, ub_range)), NumDimensions = n_keys*2,
                                            method=:random_search, TraceMode=:silent)

    # Get the optimal key points
    best_keypoints, best_deltas = constraints(best_candidate(res))

    ##############################
    # Final evaluation phase
    ##############################

    eps_post = perm_idxs[n_train+1:n]
    return_post = returns[eps_post]
    ISweights_post = ISweights[eps_post]

    return CDF_NS_band_with_KP_deltas(eps_post, return_post, ISweights_post, props, best_keypoints, best_deltas, n, :normal)
end

max_stars_repo_name: yashchandak/UnO

max_issues_repo_name: yashchandak/UnO

max_forks_repo_name: yashchandak/UnO


text:
include("common.jl")

example = """
1163751742
1381373672
2136511328
3694931569
7463417111
1319128137
1359912421
3125421639
1293138521
2311944581
"""

struct Problem
    risk::Matrix{Int8}
end

function Base.parse(::Type{Problem}, input::AbstractString)
    rows = Vector{Int8}[]
    for line in split(input, '\n')
        isempty(line) && continue
        push!(rows, map(x -> parse(Int8, x), collect(line)))
    end
    return Problem(reduce(vcat, transpose.(rows)))
end

"""Determine the neighbours we can go to from this point."""
function neighbours(problem::Problem, i::Integer, j::Integer)
    i_max, j_max = size(problem.risk)

    result = Tuple{Int64,Int64}[]
    sizehint!(result, 4)

    if i > 1
        push!(result, (i - 1, j))
    end
    if j > 1
        push!(result, (i, j - 1))
    end
    if i < i_max
        push!(result, (i + 1, j))
    end
    if j < j_max
        push!(result, (i, j + 1))
    end

    return result
end

function neighbours(problem::Problem, point::Tuple{Integer,Integer})
    return neighbours(problem, point[1], point[2])
end

function get_total_risk(problem::Problem)
    total_risk = zeros(Int64, size(problem.risk))

    point_stack = [size(total_risk)]
    while !isempty(point_stack)
        point = pop!(point_stack)
        initial_risks = filter(
            !iszero, [total_risk[point_in...] for point_in in neighbours(problem, point)]
        )
        best_initial_risk = isempty(initial_risks) ? 0 : minimum(initial_risks)
        best_risk = best_initial_risk + problem.risk[point...]

        # If we have already visited this point, determine if this path is useful. It is, then
        # we should add the next points.
        is_better = iszero(total_risk[point...]) || total_risk[point...] > best_risk

        if is_better
            total_risk[point...] = best_risk
            # Add the next points onto the stack
            append!(point_stack, neighbours(problem, point))
        end
    end

    return total_risk
end

function get_ans(problem::Problem)
    total_risk = get_total_risk(problem)
    # The starting position is never entered, so its risk is not counted
    return total_risk[1, 1] - problem.risk[1, 1]
end

get_ans(input::AbstractString) = get_ans(parse(Problem, input))

@assert get_ans(example) == 40
get_ans(load_input(15))

# Part 2

function pluswrap(x::Integer, y::Integer)
    z = x + y
    return z > 9 ? z - 9 : z
end

@assert pluswrap(8, 1) == 9
@assert pluswrap(8, 2) == 1

"""Perform the appropriate tiling operation."""
function make_problem2(problem::Problem)
    size_i, size_j = size(problem.risk)
    new_risk = Matrix{Int8}(undef, 5 .* size(problem.risk))

    for i in 0:4
        for j in 0:4
            risk_increment = Int8(i + j)
            range_i = (1 + i * size_i):((i + 1) * size_i)
            range_j = (1 + j * size_j):((j + 1) * size_j)
            for (i_old, i_new) in enumerate(range_i)
                for (j_old, j_new) in enumerate(range_j)
                    new_risk[i_new, j_new] = pluswrap(
                        problem.risk[i_old, j_old], risk_increment
                    )
                end
            end
        end
    end
    return Problem(new_risk)
end

function get_ans2(input::AbstractString)
    problem = parse(Problem, input)
    problem2 = make_problem2(problem)
    return get_ans(problem2)
end

@assert get_ans2(example) == 315
get_ans2(load_input(15))

ext: jl


text:
using Test

include("resistor-color-trio.jl")

@testset "Orange and orange and black" begin
    @test label(["orange", "orange", "black"]) == "33 ohms"
end

@testset "Blue and grey and brown" begin
    @test label(["blue", "grey", "brown"]) == "680 ohms"
end

@testset "Red and black and red" begin
    @test label(["red", "black", "red"]) == "2 kiloohms"
end

@testset "Green and brown and orange" begin
    @test label(["green", "brown", "orange"]) == "51 kiloohms"
end

@testset "Yellow and violet and yellow" begin
    @test label(["yellow", "violet", "yellow"]) == "470 kiloohms"
end

@testset "Orange and orange and red" begin
    @test label(["orange", "orange", "red"]) == "3300 ohms"
end


text:
module Threadings
using Base.Threads
export tmapreduce, tmap!, tmap

function tmapreduce(f::T, op, itr; init=nothing) where {T<:Function}
  @assert length(itr) > 0
  spinlock = SpinLock()
  popable_itr = Vector(deepcopy(itr)) # vector is poppable
  output = init == nothing ? f(pop!(popable_itr)) : deepcopy(init)
  # array of input arguments to f, to store input for each thread
  inputs = Vector{eltype(popable_itr)}(undef, nthreads()) # deprecated soon?
  @threads for i ∈ eachindex(popable_itr)
    lock(spinlock)
    inputs[threadid()] = pop!(popable_itr)
    unlock(spinlock)
    loop_output = f(inputs[threadid()])
    lock(spinlock)
    output = op(output, loop_output)
    unlock(spinlock)
  end
  return output
end
# changes only output in-place
function tmap!(f::T, itr, output::U) where {T<:Function, U<:AbstractVector}
  @assert length(itr) == length(output)
  @assert !isempty(itr)
  spinlock = SpinLock()
  all_indices = Vector{Int}(1:length(itr))
  thread_index = Vector{Int}(undef, nthreads())
  @threads for i ∈ eachindex(itr)
    lock(spinlock)
    thread_index[threadid()] = pop!(all_indices)
    unlock(spinlock)
    thread_output = f(itr[thread_index[threadid()]])
    lock(spinlock)
    output[thread_index[threadid()]] = thread_output
    unlock(spinlock)
  end
  return nothing
end
function tmap(f::T, itr, output_type::Type=Any) where {T<:Function}
  output = Vector{output_type}(undef, length(itr))
  tmap!(f, itr, output)
  return output
end

end # Threadings

max_stars_repo_name: jwscook/Threadings

max_issues_repo_name: jwscook/Threadings

max_forks_repo_name: jwscook/Threadings


text:
import MathOptInterface
const MOI = MathOptInterface

mutable struct ConstraintInfo{F, S}
    func::F
    set::S
end

const AFF = MOI.ScalarAffineFunction{Float64}
const QUA = MOI.ScalarQuadraticFunction{Float64}
const LE = MOI.LessThan{Float64}
const EQ = MOI.EqualTo{Float64}

mutable struct Optimizer <: MOI.AbstractOptimizer
    inner::Union{Opt,Nothing}

    # Problem data.
    variables::MOI.Utilities.VariablesContainer{Float64}
    starting_values::Vector{Union{Nothing,Float64}}
    nlp_data::MOI.NLPBlockData
    sense::MOI.OptimizationSense
    objective::Union{MOI.VariableIndex,AFF,QUA,Nothing}
    linear_le_constraints::Vector{ConstraintInfo{AFF, LE}}
    linear_eq_constraints::Vector{ConstraintInfo{AFF, EQ}}
    quadratic_le_constraints::Vector{ConstraintInfo{QUA, LE}}
    quadratic_eq_constraints::Vector{ConstraintInfo{QUA, EQ}}

    # Parameters.
    silent::Bool
    options::Dict{String, Any}

    # Solution attributes.
    objective_value::Float64
    solution::Vector{Float64}
    constraint_primal_linear_le::Vector{Float64}
    constraint_primal_linear_eq::Vector{Float64}
    constraint_primal_quadratic_le::Vector{Float64}
    constraint_primal_quadratic_eq::Vector{Float64}
    status::Symbol
    solve_time::Float64
end

struct EmptyNLPEvaluator <: MOI.AbstractNLPEvaluator end
MOI.features_available(::EmptyNLPEvaluator) = [:Grad, :Jac, :Hess]
MOI.initialize(::EmptyNLPEvaluator, features) = nothing
MOI.eval_objective(::EmptyNLPEvaluator, x) = NaN
function MOI.eval_constraint(::EmptyNLPEvaluator, g, x)
    @assert length(g) == 0
    return
end
function MOI.eval_objective_gradient(::EmptyNLPEvaluator, g, x)
    fill!(g, 0.0)
    return
end
MOI.jacobian_structure(::EmptyNLPEvaluator) = Tuple{Int64,Int64}[]
function MOI.eval_constraint_jacobian(::EmptyNLPEvaluator, J, x)
    return
end

empty_nlp_data() = MOI.NLPBlockData([], EmptyNLPEvaluator(), false)

const DEFAULT_OPTIONS = Dict{String, Any}(
    "algorithm" => :none,
    "stopval" => NaN,
    "ftol_rel" => 1e-7,
    "ftol_abs" => NaN,
    "xtol_rel" => 1e-7,
    "xtol_abs" => nothing,
    "constrtol_abs" => 1e-7,
    "maxeval" => 0,
    "maxtime" => 0,
    "initial_step" => nothing,
    "population" => 0,
    "seed" => nothing,
    "vector_storage" => 0,
    "local_optimizer" => nothing,
)

function Optimizer()
    return Optimizer(nothing, MOI.Utilities.VariablesContainer{Float64}(), Union{Nothing,Float64}[],
                     empty_nlp_data(), MOI.FEASIBILITY_SENSE,
                     nothing, [], [], [], [],
                     false, copy(DEFAULT_OPTIONS), NaN, Float64[], Float64[], Float64[], Float64[], Float64[], :NOT_CALLED, NaN)
end

MOI.supports(::Optimizer, ::MOI.NLPBlock) = true

MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{MOI.VariableIndex}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{AFF}) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{QUA}) = true

MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true

MOI.supports(::Optimizer, ::MOI.Silent) = true

MOI.supports(::Optimizer, p::MOI.RawOptimizerAttribute) = p.name in DEFAULT_OPTIONS

function MOI.supports(::Optimizer, ::MOI.VariablePrimalStart,
                      ::Type{MOI.VariableIndex})
    return true
end

const BOUNDS = Union{MOI.LessThan{Float64},MOI.GreaterThan{Float64},MOI.EqualTo{Float64},MOI.Interval{Float64}}

MOI.supports_constraint(::Optimizer, ::Type{MOI.VariableIndex}, ::Type{<:BOUNDS}) = true
MOI.supports_constraint(::Optimizer, ::Type{AFF}, ::Type{LE}) = true
MOI.supports_constraint(::Optimizer, ::Type{AFF}, ::Type{EQ}) = true
MOI.supports_constraint(::Optimizer, ::Type{QUA}, ::Type{LE}) = true
MOI.supports_constraint(::Optimizer, ::Type{QUA}, ::Type{EQ}) = true

MOI.supports_incremental_interface(::Optimizer) = true

function MOI.copy_to(model::Optimizer, src::MOI.ModelLike)
    return MOI.Utilities.default_copy_to(model, src)
end

MOI.get(::Optimizer, ::MOI.SolverName) = "NLopt"

MOI.get(model::Optimizer, ::MOI.ObjectiveFunctionType) = typeof(model.objective)

MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{AFF, LE}) = length(model.linear_le_constraints)
MOI.get(model::Optimizer, ::MOI.NumberOfConstraints{AFF, EQ}) = length(model.linear_eq_constraints)
MOI.get(model::Optimizer, attr::Union{MOI.NumberOfVariables,MOI.ListOfVariableIndices,MOI.NumberOfConstraints{MOI.VariableIndex}}) = MOI.get(model.variables, attr)

function MOI.get(model::Optimizer, attr::MOI.ListOfConstraintTypesPresent)
    constraints = Set{Tuple{Type, Type}}(MOI.get(model.variables, attr))

    # handling model constraints separately
    if !isempty(model.linear_le_constraints)
        push!(constraints, (AFF, LE))
    end
    if !isempty(model.linear_eq_constraints)
        push!(constraints, (AFF, EQ))
    end
    if !isempty(model.quadratic_le_constraints)
        push!(constraints, (QUA, LE))
    end
    if !isempty(model.quadratic_eq_constraints)
        push!(constraints, (QUA, EQ))
    end

    return collect(constraints)
end

function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{AFF, LE}
)
    return MOI.ConstraintIndex{AFF, LE}.(eachindex(model.linear_le_constraints))
end

function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{AFF, EQ}
)
    return MOI.ConstraintIndex{AFF, EQ}.(eachindex(model.linear_eq_constraints))
end

function MOI.get(model::Optimizer, attr::MOI.ListOfConstraintIndices{MOI.VariableIndex})
    return MOI.get(model.variables, attr)
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintFunction,
    c::MOI.ConstraintIndex{AFF, LE}
)
    return model.linear_le_constraints[c.value].func
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintFunction,
    c::MOI.ConstraintIndex{AFF, EQ}
)
    return model.linear_eq_constraints[c.value].func
end

function MOI.get(
    model::Optimizer,
    attr::Union{MOI.ConstraintFunction,MOI.ConstraintSet},
    ci::MOI.ConstraintIndex{MOI.VariableIndex}
)
    return MOI.get(model.variables, attr, ci)
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{AFF, LE}
)
    return model.linear_le_constraints[c.value].set
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{AFF, EQ}
)
    return model.linear_eq_constraints[c.value].set
end

function MOI.get(
    model::Optimizer,
    ::MOI.ObjectiveFunction
)
    return model.objective
end

function MOI.set(model::Optimizer, ::MOI.ObjectiveSense,
                 sense::MOI.OptimizationSense)
    model.sense = sense
    return
end

MOI.get(model::Optimizer, ::MOI.ObjectiveSense) = model.sense

function MOI.set(model::Optimizer, ::MOI.Silent, value)
    model.silent = value
    return
end

MOI.get(model::Optimizer, ::MOI.Silent) = model.silent

const TIME_LIMIT = "max_cpu_time"
MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
function MOI.set(model::Optimizer, ::MOI.TimeLimitSec, value::Real)
    MOI.set(model, MOI.RawOptimizerAttribute(TIME_LIMIT), Float64(value))
end
function MOI.set(model::Optimizer, attr::MOI.TimeLimitSec, ::Nothing)
    delete!(model.options, TIME_LIMIT)
end
function MOI.get(model::Optimizer, ::MOI.TimeLimitSec)
    return get(model.options, TIME_LIMIT, nothing)
end

function MOI.set(model::Optimizer, p::MOI.RawOptimizerAttribute, value)
    model.options[p.name] = value
    return
end

function MOI.get(model::Optimizer, p::MOI.RawOptimizerAttribute)
    if haskey(model.options, p.name)
        return model.options[p.name]
    end
    error("RawOptimizerAttribute with name $(p.name) is not set.")
end

MOI.get(model::Optimizer, ::MOI.SolveTimeSec) = model.solve_time

function MOI.empty!(model::Optimizer)
    model.inner = nothing
    MOI.empty!(model.variables)
    empty!(model.starting_values)
    model.nlp_data = empty_nlp_data()
    model.sense = MOI.FEASIBILITY_SENSE
    model.objective = nothing
    empty!(model.linear_le_constraints)
    empty!(model.linear_eq_constraints)
    empty!(model.quadratic_le_constraints)
    empty!(model.quadratic_eq_constraints)
    model.status = :NOT_CALLED
    return
end

function MOI.is_empty(model::Optimizer)
    return MOI.is_empty(model.variables) &&
           isempty(model.starting_values) &&
           model.nlp_data.evaluator isa EmptyNLPEvaluator &&
           model.sense == MOI.FEASIBILITY_SENSE &&
           isempty(model.linear_le_constraints) &&
           isempty(model.linear_eq_constraints) &&
           isempty(model.quadratic_le_constraints) &&
           isempty(model.quadratic_eq_constraints)
end

function MOI.add_variable(model::Optimizer)
    push!(model.starting_values, nothing)
    return MOI.add_variable(model.variables)
end
function MOI.is_valid(model::Optimizer, index::Union{MOI.VariableIndex,MOI.ConstraintIndex{MOI.VariableIndex}})
    return MOI.is_valid(model.variables, index)
end
function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{AFF, LE})
    return ci.value in eachindex(model.linear_le_constraints)
end
function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{AFF, EQ})
    return ci.value in eachindex(model.linear_eq_constraints)
end
function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{QUA, LE})
    return ci.value in eachindex(model.quadratic_le_constraints)
end
function MOI.is_valid(model::Optimizer, ci::MOI.ConstraintIndex{QUA, EQ})
    return ci.value in eachindex(model.quadratic_eq_constraints)
end

function check_inbounds(model::Optimizer, var::MOI.VariableIndex)
    return MOI.throw_if_not_valid(model, var)
end

function check_inbounds(model::Optimizer, aff::AFF)
    for term in aff.terms
        MOI.throw_if_not_valid(model, term.variable)
    end
end

function check_inbounds(model::Optimizer, quad::QUA)
    for term in quad.affine_terms
        MOI.throw_if_not_valid(model, term.variable)
    end
    for term in quad.quadratic_terms
        MOI.throw_if_not_valid(model, term.variable_1)
        MOI.throw_if_not_valid(model, term.variable_2)
    end
end

function MOI.add_constraint(model::Optimizer, vi::MOI.VariableIndex, set::BOUNDS)
    return MOI.add_constraint(model.variables, vi, set)
end

function MOI.set(model::Optimizer, attr::MOI.ConstraintSet, ci::MOI.ConstraintIndex{MOI.VariableIndex}, set)
    return MOI.set(model.variables, attr, ci, set)
end

function MOI.delete(model::Optimizer, ci::MOI.ConstraintIndex{MOI.VariableIndex})
    return MOI.delete(model.variables, ci)
end

function _fill_jacobian_terms(jac::Matrix, x, row, terms::Vector{MOI.ScalarAffineTerm{Float64}})
    for term in terms
        jac[term.variable.value, row] += term.coefficient
    end
end
function _fill_jacobian_terms(jac::Matrix, x, row, terms::Vector{MOI.ScalarQuadraticTerm{Float64}})
    for term in terms
        jac[term.variable_1.value, row] += term.coefficient * x[term.variable_2.value]
        if term.variable_1 != term.variable_2
            jac[term.variable_2.value, row] += term.coefficient * x[term.variable_1.value]
        end
    end
end
function _fill_jacobian(jac::Matrix, x, offset, constraints::Vector{<:ConstraintInfo{AFF}})
    for (k, con) in enumerate(constraints)
        _fill_jacobian_terms(jac, x, offset + k, con.func.terms)
    end
end
function _fill_jacobian(jac::Matrix, x, offset, constraints::Vector{<:ConstraintInfo{QUA}})
    for (k, con) in enumerate(constraints)
        _fill_jacobian_terms(jac, x, offset + k, con.func.affine_terms)
        _fill_jacobian_terms(jac, x, offset + k, con.func.quadratic_terms)
    end
end

function _fill_result(result::Vector, x, offset, constraints::Vector)
    for (k, con) in enumerate(constraints)
        result[offset + k] = MOI.Utilities.eval_variables(vi -> x[vi.value], con.func) - MOI.constant(con.set)
    end
end

function _fill_primal(x, constraints::Vector)
    result = zeros(length(constraints))
    for (k, con) in enumerate(constraints)
        result[k] = MOI.Utilities.eval_variables(vi -> x[vi.value], con.func)
    end
    return result
end

macro define_add_constraint(function_type, set_type, prefix)
    array_name = Symbol(string(prefix) * "_constraints")
    quote
        function MOI.add_constraint(model::Optimizer, func::$function_type, set::$set_type)
            check_inbounds(model, func)
            push!(model.$(array_name), ConstraintInfo(func, set))
            return MOI.ConstraintIndex{$function_type, $set_type}(length(model.$(array_name)))
        end
    end
end

@define_add_constraint(AFF, LE,
                       linear_le)
@define_add_constraint(AFF, EQ,
                       linear_eq)
@define_add_constraint(QUA,
                       LE, quadratic_le)
@define_add_constraint(QUA,
                       EQ, quadratic_eq)

function starting_value(optimizer::Optimizer, i)
    if optimizer.starting_values[i] !== nothing
        return optimizer.starting_values[i]
    else
        v = optimizer.variables
        return min(max(0.0, v.lower[i]), v.upper[i])
    end
end

function MOI.set(model::Optimizer, ::MOI.VariablePrimalStart,
                 vi::MOI.VariableIndex, value::Union{Real, Nothing})
    MOI.throw_if_not_valid(model, vi)
    model.starting_values[vi.value] = value
    return
end

function MOI.set(model::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
    model.nlp_data = nlp_data
    return
end

function MOI.set(model::Optimizer, ::MOI.ObjectiveFunction,
                 func::Union{MOI.VariableIndex, AFF,
                             QUA})
    check_inbounds(model, func)
    model.objective = func
    return
end

function eval_objective(model::Optimizer, x)
    # The order of the conditions is important. NLP objectives override regular
    # objectives.
    if model.nlp_data.has_objective
        return MOI.eval_objective(model.nlp_data.evaluator, x)
    elseif model.objective !== nothing
        return MOI.Utilities.eval_variables(vi -> x[vi.value], model.objective)
    else
        # No objective function set. This could happen with FEASIBILITY_SENSE.
        return 0.0
    end
end

function fill_gradient!(grad, x, var::MOI.VariableIndex)
    fill!(grad, 0.0)
    grad[var.value] = 1.0
end

function fill_gradient!(grad, x, aff::AFF)
    fill!(grad, 0.0)
    for term in aff.terms
        grad[term.variable.value] += term.coefficient
    end
end

function fill_gradient!(grad, x, quad::QUA)
    fill!(grad, 0.0)
    for term in quad.affine_terms
        grad[term.variable.value] += term.coefficient
    end
    for term in quad.quadratic_terms
        row_idx = term.variable_1
        col_idx = term.variable_2
        coefficient = term.coefficient
        if row_idx == col_idx
            grad[row_idx.value] += coefficient*x[row_idx.value]
        else
            grad[row_idx.value] += coefficient*x[col_idx.value]
            grad[col_idx.value] += coefficient*x[row_idx.value]
        end
    end
end

function eval_objective_gradient(model::Optimizer, grad, x)
    if model.nlp_data.has_objective
        MOI.eval_objective_gradient(model.nlp_data.evaluator, grad, x)
    elseif model.objective !== nothing
        fill_gradient!(grad, x, model.objective)
    else
        fill!(grad, 0.0)
    end
    return
end

# Refers to local variables in eval_constraint() below.
macro eval_function(array_name)
    escrow = esc(:row)
    quote
        for info in $(esc(array_name))
            $(esc(:g))[$escrow] = eval_function(info.func, $(esc(:x)))
            $escrow += 1
        end
    end
end

function eval_constraint(model::Optimizer, g, x)
    row = 1
    @eval_function model.linear_le_constraints
    @eval_function model.linear_eq_constraints
    @eval_function model.quadratic_le_constraints
    @eval_function model.quadratic_eq_constraints
    nlp_g = view(g, row:length(g))
    MOI.eval_constraint(model.nlp_data.evaluator, nlp_g, x)
    return
end

function fill_constraint_jacobian!(values, start_offset, x, aff::AFF)
    num_coefficients = length(aff.terms)
    for i in 1:num_coefficients
        values[start_offset+i] = aff.terms[i].coefficient
    end
    return num_coefficients
end

function fill_constraint_jacobian!(values, start_offset, x, quad::QUA)
    num_affine_coefficients = length(quad.affine_terms)
    for i in 1:num_affine_coefficients
        values[start_offset+i] = quad.affine_terms[i].coefficient
    end
    num_quadratic_coefficients = 0
    for term in quad.quadratic_terms
        row_idx = term.variable_1
        col_idx = term.variable_2
        coefficient = term.coefficient
        if row_idx == col_idx
            values[start_offset+num_affine_coefficients+num_quadratic_coefficients+1] = coefficient*x[col_idx.value]
            num_quadratic_coefficients += 1
        else
            # Note that the order matches the Jacobian sparsity pattern.
            values[start_offset+num_affine_coefficients+num_quadratic_coefficients+1] = coefficient*x[col_idx.value]
            values[start_offset+num_affine_coefficients+num_quadratic_coefficients+2] = coefficient*x[row_idx.value]
            num_quadratic_coefficients += 2
        end
    end
    return num_affine_coefficients + num_quadratic_coefficients
end

# Refers to local variables in eval_constraint_jacobian() below.
macro fill_constraint_jacobian(array_name)
    esc_offset = esc(:offset)
    quote
        for info in $(esc(array_name))
            $esc_offset += fill_constraint_jacobian!($(esc(:values)),
                                                     $esc_offset, $(esc(:x)),
                                                     info.func)
        end
    end
end

function eval_constraint_jacobian(model::Optimizer, values, x)
    offset = 0
    @fill_constraint_jacobian model.linear_le_constraints
    @fill_constraint_jacobian model.linear_eq_constraints
    @fill_constraint_jacobian model.quadratic_le_constraints
    @fill_constraint_jacobian model.quadratic_eq_constraints

    nlp_values = view(values, 1+offset:length(values))
    MOI.eval_constraint_jacobian(model.nlp_data.evaluator, nlp_values, x)
    return
end

function constraint_bounds(model::Optimizer)
    constraint_lb = Float64[]
    constraint_ub = Float64[]
    for info in model.linear_le_constraints
        push!(constraint_lb, -Inf)
        push!(constraint_ub, info.set.upper)
    end
    for info in model.linear_eq_constraints
        push!(constraint_lb, info.set.value)
        push!(constraint_ub, info.set.value)
    end
    for info in model.quadratic_le_constraints
        push!(constraint_lb, -Inf)
        push!(constraint_ub, info.set.upper)
    end
    for info in model.quadratic_eq_constraints
        push!(constraint_lb, info.set.value)
        push!(constraint_ub, info.set.value)
    end
    for bound in model.nlp_data.constraint_bounds
        push!(constraint_lb, bound.lower)
        push!(constraint_ub, bound.upper)
    end
    return constraint_lb, constraint_ub
end

function MOI.optimize!(model::Optimizer)
    # TODO: Reuse model.inner for incremental solves if possible.
    num_variables = length(model.starting_values)
    model.inner = Opt(model.options["algorithm"], num_variables)
    local_optimizer = model.options["local_optimizer"]
    if local_optimizer !== nothing
        if local_optimizer isa Symbol
            local_optimizer = Opt(local_optimizer, num_variables)
        else
            local_optimizer = Opt(local_optimizer.algorithm, num_variables)
        end
        local_optimizer!(model.inner, local_optimizer)
    end

    # load parameters
    stopval!(model.inner, model.options["stopval"])
    if !isnan(model.options["ftol_rel"])
        ftol_rel!(model.inner, model.options["ftol_rel"])
    end
    if !isnan(model.options["ftol_abs"])
        ftol_abs!(model.inner, model.options["ftol_abs"])
    end
    if !isnan(model.options["xtol_rel"])
        xtol_rel!(model.inner, model.options["xtol_rel"])
    end
    if model.options["xtol_abs"] != nothing
        xtol_abs!(model.inner, model.options["xtol_abs"])
    end
    maxeval!(model.inner, model.options["maxeval"])
    maxtime!(model.inner, model.options["maxtime"])
    if model.options["initial_step"] != nothing
        initial_step!(model.inner, model.options["initial_step"])
    end
    population!(model.inner, model.options["population"])
    if isa(model.options["seed"], Integer)
        NLopt.srand(model.options["seed"])
    end
    vector_storage!(model.inner, model.options["vector_storage"])

    lower_bounds!(model.inner, model.variables.lower)
    upper_bounds!(model.inner, model.variables.upper)

    nleqidx = findall(bound -> bound.lower == bound.upper, model.nlp_data.constraint_bounds) # indices of equalities
    nlineqidx = findall(bound -> bound.lower != bound.upper, model.nlp_data.constraint_bounds)

    num_nl_constraints = length(model.nlp_data.constraint_bounds)

    # map from eqidx/ineqidx to index in equalities/inequalities
    constrmap = zeros(Int, num_nl_constraints)
    for i in eachindex(nleqidx)
        constrmap[nleqidx[i]] = i
    end
    ineqcounter = 1
    for i in eachindex(nlineqidx)
        k = nlineqidx[i]
        constrmap[k] = ineqcounter
        bounds = model.nlp_data.constraint_bounds[k]
        if isinf(bounds.lower) || isinf(bounds.upper)
            ineqcounter += 1
        else # constraint has bounds on both sides, keep room for it
            ineqcounter += 2
        end
    end
    num_nl_ineq = ineqcounter - 1
    num_nl_eq = length(nleqidx)

    isderivativefree = string(model.options["algorithm"])[2] == 'N'
    if isderivativefree
        requested_features = Symbol[]
    else
        requested_features = num_nl_constraints > 0 ? [:Grad, :Jac] : [:Grad]
    end

    MOI.initialize(model.nlp_data.evaluator, requested_features)

    if model.sense == MOI.FEASIBILITY_SENSE
        # If we don't give any objective to NLopt, it throws the error:
        # invalid NLopt arguments: NULL args to nlopt_optimize
        function z(x::Vector, grad::Vector)
            fill!(grad, 0.0)
            return 0.0
        end
        min_objective!(model.inner, z)
    else
        function f(x::Vector, grad::Vector)
            if length(grad) > 0
                eval_objective_gradient(model, grad, x)
            end
            return eval_objective(model, x)
        end
        if model.sense == MOI.MIN_SENSE
            min_objective!(model.inner, f)
        else
            max_objective!(model.inner, f)
        end
    end

    Jac_IJ = num_nl_constraints > 0 ? MOI.jacobian_structure(model.nlp_data.evaluator) : Tuple{Int,Int}[]
    Jac_val = zeros(length(Jac_IJ))
    g_vec = zeros(num_nl_constraints)

    num_eq = num_nl_eq + length(model.linear_eq_constraints) + length(model.quadratic_eq_constraints)

    if num_eq > 0
        function g_eq(result::Vector, x::Vector, jac::Matrix)
            if length(jac) > 0
                fill!(jac, 0.0)
                MOI.eval_constraint_jacobian(model.nlp_data.evaluator, Jac_val, x)
                for k in 1:length(Jac_val)
                    row, col = Jac_IJ[k]
                    bounds = model.nlp_data.constraint_bounds[row]
                    if bounds.lower == bounds.upper
                        jac[col, constrmap[row]] += Jac_val[k]
                    end
                end
                _fill_jacobian(jac, x, num_nl_eq, model.linear_eq_constraints)
                _fill_jacobian(jac, x, num_nl_eq + length(model.linear_eq_constraints), model.quadratic_eq_constraints)
            end
            MOI.eval_constraint(model.nlp_data.evaluator, g_vec, x)
            for (ctr, idx) in enumerate(nleqidx)
                bounds = model.nlp_data.constraint_bounds[idx]
                result[ctr] = g_vec[idx] - bounds.upper
            end
            _fill_result(result, x, num_nl_eq, model.linear_eq_constraints)
            _fill_result(result, x, num_nl_eq + length(model.linear_eq_constraints), model.quadratic_eq_constraints)
        end

        g_eq(zeros(num_eq), zeros(num_variables), zeros(num_variables, num_eq))
        equality_constraint!(model.inner, g_eq, fill(model.options["constrtol_abs"], num_eq))
    end

    # inequalities need to be massaged a bit
    # f(x) <= u   =>  f(x) - u <= 0
    # f(x) >= l   =>  l - f(x) <= 0

    num_ineq = num_nl_ineq + length(model.linear_le_constraints) + length(model.quadratic_le_constraints)

    if num_ineq > 0
        function g_ineq(result::Vector, x::Vector, jac::Matrix)
            if length(jac) > 0
                fill!(jac, 0.0)
                MOI.eval_constraint_jacobian(model.nlp_data.evaluator, Jac_val, x)
                for k in 1:length(Jac_val)
                    row, col = Jac_IJ[k]
                    bounds = model.nlp_data.constraint_bounds[row]
                    bounds.lower == bounds.upper && continue
                    if isinf(bounds.lower) # upper bound
                        jac[col, constrmap[row]] += Jac_val[k]
                    elseif isinf(bounds.upper) # lower bound
                        jac[col, constrmap[row]] -= Jac_val[k]
                    else
                        # boxed
                        jac[col, constrmap[row]] += Jac_val[k]
                        jac[col, constrmap[row] + 1] -= Jac_val[k]
                    end
                end
                _fill_jacobian(jac, x, num_nl_ineq, model.linear_le_constraints)
                _fill_jacobian(jac, x, num_nl_ineq + length(model.linear_le_constraints), model.quadratic_le_constraints)
            end
            MOI.eval_constraint(model.nlp_data.evaluator, g_vec, x)
            for row in 1:num_nl_constraints
                bounds = model.nlp_data.constraint_bounds[row]
                bounds.lower == bounds.upper && continue
                if isinf(bounds.lower)
                    result[constrmap[row]] = g_vec[row] - bounds.upper
                elseif isinf(bounds.upper)
                    result[constrmap[row]] = bounds.lower - g_vec[row]
                else
                    result[constrmap[row]] = g_vec[row] - bounds.upper
                    result[constrmap[row] + 1] = bounds.lower - g_vec[row]
                end
            end
            _fill_result(result, x, num_nl_ineq, model.linear_le_constraints)
            _fill_result(result, x, num_nl_ineq + length(model.linear_le_constraints), model.quadratic_le_constraints)
        end

        g_ineq(zeros(num_ineq), zeros(num_variables), zeros(num_variables, num_ineq))
        inequality_constraint!(model.inner, g_ineq, fill(model.options["constrtol_abs"], num_ineq))
    end

    # If nothing is provided, the default starting value is 0.0.
    model.solution = zeros(num_variables)
    for i in eachindex(model.starting_values)
        model.solution[i] = starting_value(model, i)
    end

    start_time = time()

    model.objective_value, _, model.status = optimize!(model.inner, model.solution)

    model.constraint_primal_linear_le = _fill_primal(model.solution, model.linear_le_constraints)
    model.constraint_primal_linear_eq = _fill_primal(model.solution, model.linear_eq_constraints)
    model.constraint_primal_quadratic_le = _fill_primal(model.solution, model.quadratic_le_constraints)
    model.constraint_primal_quadratic_eq = _fill_primal(model.solution, model.quadratic_eq_constraints)

    model.solve_time = time() - start_time
    return
end

function MOI.get(model::Optimizer, ::MOI.TerminationStatus)
    if model.status == :NOT_CALLED
        return MOI.OPTIMIZE_NOT_CALLED
    elseif model.status == :SUCCESS || model.status == :FTOL_REACHED || model.status == :XTOL_REACHED
        return MOI.LOCALLY_SOLVED
    elseif model.status == :ROUNDOFF_LIMITED
        return MOI.ALMOST_LOCALLY_SOLVED
    elseif model.status == :MAXEVAL_REACHED
        return MOI.ITERATION_LIMIT
    elseif model.status == :MAXTIME_REACHED
        return MOI.TIME_LIMIT
    elseif model.status == :STOPVAL_REACHED
        return MOI.OTHER_LIMIT
    elseif model.status == :FORCED_STOP
        # May be due to an error in the callbacks `f`, `g_eq` and `g_ineq` defined in `MOI.optimize!`
        return MOI.OTHER_ERROR
    elseif model.status == :OUT_OF_MEMORY
        return MOI.MEMORY_LIMIT
    elseif model.status == :INVALID_ARGS
        return MOI.INVALID_OPTION
    elseif model.status == :FAILURE
        return MOI.OTHER_ERROR
    else
        error("Unknown status $(model.status)")
    end
end

function MOI.get(model::Optimizer, ::MOI.RawStatusString)
    return string(model.status)
end

# Ipopt always has an iterate available.
function MOI.get(model::Optimizer, ::MOI.ResultCount)
    return model.status == :NOT_CALLED ? 0 : 1
end

function MOI.get(model::Optimizer, attr::MOI.PrimalStatus)
    if !(1 <= attr.result_index <= MOI.get(model, MOI.ResultCount()))
        return MOI.NO_SOLUTION
    end
    if model.status == :SUCCESS || model.status == :FTOL_REACHED || model.status == :XTOL_REACHED
        return MOI.FEASIBLE_POINT
    elseif model.status == :ROUNDOFF_LIMITED
        return MOI.NEARLY_FEASIBLE_POINT
    elseif model.status in (:STOPVAL_REACHED, :MAXEVAL_REACHED, :MAXTIME_REACHED, :FORCED_STOP, :OUT_OF_MEMORY, :INVALID_ARGS, :FAILURE)
        return MOI.UNKNOWN_RESULT_STATUS
    else
        error("Unknown status $(model.status)")
    end
end

function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(model, attr)
    return model.objective_value
end

function MOI.get(model::Optimizer, attr::MOI.VariablePrimal, vi::MOI.VariableIndex)
    MOI.check_result_index_bounds(model, attr)
    MOI.throw_if_not_valid(model, vi)
    return model.solution[vi.value]
end

function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{MOI.VariableIndex,
                                         <:Union{LE, MOI.GreaterThan{Float64}, EQ}})
    MOI.check_result_index_bounds(model, attr)
    MOI.throw_if_not_valid(model, ci)
    return model.solution[ci.value]
end
function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{AFF, LE})
    MOI.check_result_index_bounds(model, attr)
    MOI.throw_if_not_valid(model, ci)
    return model.constraint_primal_linear_le[ci.value]
end
function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{AFF, EQ})
    MOI.check_result_index_bounds(model, attr)
    MOI.throw_if_not_valid(model, ci)
    return model.constraint_primal_linear_eq[ci.value]
end
function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{QUA, LE})
    MOI.check_result_index_bounds(model, attr)
    MOI.throw_if_not_valid(model, ci)
    return model.constraint_primal_quadratic_le[ci.value]
end
function MOI.get(model::Optimizer, attr::MOI.ConstraintPrimal,
                 ci::MOI.ConstraintIndex{QUA, EQ})
    MOI.check_result_index_bounds(model, attr)
    MOI.throw_if_not_valid(model, ci)
    return model.constraint_primal_quadratic_eq[ci.value]
end


text:
"""
    module DielectricFunctions

for defining dielectric and pump functions either as piecewise constant
or with user-supplied functions
"""
module DielectricFunctions

export DielectricFunction,
PumpFunction,
piecewise_constant_ε,
piecewise_constant_F

"""
    struct DielectricFunction

---

    DielectricFunction(n1=1,n2=0) -> dielectric
    DielectricFunction(n::Complex) -> dielectric

piecewise constant dielectric function with index given by `n1+1im*n2`

---

    DielectricFunction(ε::Function, parameters::Dict) -> dielectric

ε is any function with the signature ε(x,y,parameters::Dict{Symbol,Float64})
which evaluates to the dielectric at position `x,y`. For example:
    ε(x,y,params) = 2+sin(2πx/params[:period])

---

    (::DielectricFunction)(x,y) -> ε
"""
struct DielectricFunction
    ε::Function
    parameters::Dict{Symbol,Float64}

    DielectricFunction(ε::Function=piecewise_constant_ε,
        parameters::Dict{Symbol,T}=Dict(:n₁=>1,:n₂=>0)
        ) where T<:Number = new(ε,parameters)

    DielectricFunction(n::Complex) = DielectricFunction(real(n),imag(n))

    function DielectricFunction(n1::Real,n2::Real=0.0)
        parameters = Dict(:n₁=>n1,:n₂=>n2)
        return DielectricFunction(piecewise_constant_ε,parameters)
    end

    (df::DielectricFunction)(x::Number,y::Number) = df.ε(x,y,df.parameters)

    function Base.show(io::IO,d::DielectricFunction)
        print(io, "DielectricFunction: ")
        if d.ε==piecewise_constant_ε
            print(io, "piecewise_constant_ε: ",d.parameters[:n₁]+1im*d.parameters[:n₂])
        else
            print(io, d.ε,": ", d.parameters)
        end
    end
end

"""
    struct PumpFunction

---

    PumpFunction(F=0) -> pump

piecewise constant pump function with pump parameter given by `F`

---

    PumpFunction(F::Function, parameters::Dict) -> pump

F is any real function with the signature F(x,y,parameters::Dict{Symbol,Float64})
which evaluates to the pump profile at position `x,y`. For example:
    F(x,y,params) = sin(2πx/params[:period])

---

    (::PumpFunction)(x,y) -> F
"""
struct PumpFunction
    F::Function
    parameters::Dict{Symbol,Float64}

    PumpFunction(F::Function=piecewise_constant_F,
        parameters::Dict{Symbol,T}=Dict(:F=>0)
        ) where T<:Number = new(F,parameters)

    function PumpFunction(F::Real)
        parameters = Dict(:F=>F)
        return PumpFunction(piecewise_constant_F,parameters)
    end

    (pf::PumpFunction)(x::Number,y::Number) = pf.F(x,y,pf.parameters)

    function Base.show(io::IO,d::PumpFunction)
        print(io, "PumpFunction: ")
        if d.F==piecewise_constant_F
            print(io, "piecewise_constant_F: ",d.parameters[:F])
        else
            print(io, d.F,": ", d.parameters)
        end
    end
end

"""
    piecewise_constant_ε(x,y,parameters) -> ε
"""
function piecewise_constant_ε(x::Number,y::Number,parameters::Dict)
    n1 = get!(parameters,:n₁,nothing)
    n2 = get!(parameters,:n₂,nothing)
    return complex(n1,n2)^2
end

"""
    piecewise_constant_F(x,y,parameters) -> F
"""
piecewise_constant_F(x::Number,y::Number,parameters::Dict) = get!(parameters,:F,nothing)

end # module


text:
@graphs begin
	################################################################################
	@node select_action_mili = Start() -> threat_check
	@node threat_check = JumpToGraph("threat_exposed") -> a1

	########################################
	@node a1 = Dummy() -> a1_1
	@node a1_1 = SetActiveDie('C', may_use_ring = true) -> [next = a1_cond, no_die = a2]
	@node a1_cond = BinaryCondition("""
									The Which King is in play and not in a *mobile* army but is able to create or join one.
									""") -> [n_true = a1_jump, n_false = a2]
	@node a1_jump = JumpToGraph("character_which_king") -> a2

	########################################
	@node a2 = Dummy() -> a2_1
	@node a2_1 = SetActiveDie('M', may_use_ring = true) -> [next = a2_jump, no_die = a3]
	@node a2_jump = JumpToGraph("muster_minion") -> a3

	########################################
	@node a3 = Dummy() -> a3_1
	@node a3_1 = SetActiveDie('M') -> [next = a3_jump, no_die = a4]
	@node a3_jump = JumpToGraph("muster_politics") -> a4

	########################################
	@node a4 = Dummy() -> a4_1
	@node a4_1 = SetActiveDie('C', may_use_ring = true) -> [next = a4_2, no_die = a5]
	@node a4_2 = SetActiveDie('P', may_use_ring = true) -> [next = a4_cond, no_die = a5]
	@node a4_cond = BinaryCondition("""
									The Fellowship is on the Mordor track or is revealed.
									And, a "Fellowship revealed" character card is held.
									""") -> [n_true = a4_die, n_false = a5]
	@node a4_die = UseActiveDie() -> a4_action
	@node a4_action = PerformAction("""
									Play a "Fellowship revealed" character card.

									Priority:
									1. Ascending order of initiative
									2. Random
									""") -> a4_end
	@node a4_end = End() -> []

	########################################
	@node a5 = Dummy() -> a5_1
	@node a5_1 = SetActiveDie('C', may_use_ring = true) -> [next = a5_2, no_die = a6]
	@node a5_2 = SetActiveDie('A', may_use_ring = true) -> [next = a5_cond, no_die = a6]
	@node a5_cond = BinaryCondition("""
									A *mobile* army is adjacent to its *target* or a Free Peoples' army on the shortest route to its *target*.

									And, any of the following conditions hold:
									- The *mobile* army's *target* gives enough points to win.
									- The *mobile* army's *target* is in a nation at war and is not under siege.
									- The Fellowship is on the Mordor track.
									""") -> [n_true = a5_action, n_false = a6]
	@node a5_action = JumpToGraph("movement_attack_basic") -> a6

	########################################
	@node a6 = Dummy() -> a6_1
	@node a6_1 = SetActiveDie('M') -> [next = a6_cond, no_die = a7]
	@node a6_cond = BinaryCondition("""
									A strategy card that musters is *playable*.
									""") -> [n_true = a6_die, n_false = a6_2]
	@node a6_2 = SetActiveDie('A') -> [next = a6_cond_2, no_die = a7]
	@node a6_cond_2 = BinaryCondition("""
									  A strategy card that musters is *playable*.
									  """) -> [n_true = a6_die, n_false = a7]
	@node a6_die = UseActiveDie() -> a6_action
	@node a6_action = PerformAction("""
									Play a strategy card that musters.

									Priority:
									3. Ascending order of initiative
									4. Random
									""") -> a6_end
	@node a6_end = End() -> []

	########################################
	@node a7 = Dummy() -> a7_1
	@node a7_1 = BinaryCondition("The Shadow player is allowed to pass.") -> [n_true = a7_action, n_false = a8]
	@node a7_action = PerformAction("Pass") -> a7_end
	@node a7_end = End() -> []

	########################################
	@node a8 = Dummy() -> a8_1
	@node a8_1 = SetActiveDie('P') -> [next = a8_jump_1, no_die = a9]
	@node a8_jump_1 = JumpToGraph("event_cards_preferred") -> a8_jump_2
	@node a8_jump_2 = JumpToGraph("event_cards_general") -> a9

	########################################
	@node a9 = Dummy() -> a9_1
	@node a9_1 = SetActiveDie('C') -> [next = a9_2, no_die = a10]
	@node a9_2 = SetActiveDie('A') -> [next = a9_cond, no_die = a10]
	@node a9_cond = BinaryCondition("""
									A *mobile* army is adjacent to its *target*.
									""") -> [n_true = a9_action, n_false = a10]
	@node a9_action = JumpToGraph("movement_attack_basic") -> a10

	########################################
	@node a10 = Dummy() -> a10_1
	@node a10_1 = SetActiveDie('A') -> [next = a10_action, no_die = a11]
	@node a10_action = JumpToGraph("movement_attack_corr") -> a11

	########################################
	@node a11 = Dummy() -> a11_1
	@node a11_1 = SetActiveDie('C') -> [next = a11_action, no_die = a12]
	@node a11_action = JumpToGraph("character_army") -> a12

	########################################
	@node a12 = Dummy() -> a12_start
	@node a12_start = SetActiveDie('M') -> [next = a12_1, no_die = a13]
	@node a12_1 = JumpToGraph("muster_minion") -> a12_2
	@node a12_2 = JumpToGraph("muster_politics") -> a12_3
	@node a12_3 = JumpToGraph("muster_muster") -> a13

	########################################
	@node a13 = Dummy() -> a13_1
	@node a13_1 = CheckStrategy("military") -> [n_true = a13_mili, n_false = a13_corr]
	@node a13_mili = PerformAction("""
								   Queller failed to find an action. Discard a random Character or Event die if possible, otherwise discard a random die (do not discard a die set aside for later use).
								   """) -> a13_dice
	@node a13_corr = PerformAction("""
								   Queller failed to find an action. Discard a random Army, Muster, Muster/Army or Event die if possible, otherwise discard a random die (do not discard a die set aside for later use).
								   """) -> a13_dice
	@node a13_dice = GetAvailableDice("""
									  After discarding the random die, input the remaining available dice here (not counting dice set aside for later use).
									  """) -> a_end
	@node a_end = End() -> []

end

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
using ScoreDrivenModels

Model(1, 2, LogNormal, 0.5)

Model(1, 2, LogNormal, 0.5; time_varying_params = [1])

Model([1, 12], [1, 12], LogNormal, 0.5)


text:
function zad39(A,el1)
    a,b = size(A)
    s = 0
    V = []
    for i in 1:a
        for j in 1:b
         if i == j && A[i,j] > el1
            s +=A[i,j]
         end
         if rem(i,2)==1 && rem(A[i,j],el1) == 0
            push!(V,A[i,j])
         end
        end
    end

    return V, s
end


text:
## Mimi UI
using CSVFiles
using DataFrames

function get_sim_results(sim_inst::SimulationInstance, comp_name::Symbol, datum_name::Symbol; model_index::Int = 1, scen_name::Union{Nothing, String} = nothing)

    multiple_results = (length(sim_inst.results) > 1)

    key = (comp_name, datum_name)
    df = (sim_inst.results[model_index])[key]
    if scen_name !== nothing
        if in(:scen, names(df))
            error("The results for this simulation contain a scenario dimension, you must specify the scen_name keyword argument, which is currently set to $(scen_name)")
        end
        filter!(row -> row.scen_name === scen_name, df)
    end
    return df
end

function get_sim_results(sim_inst::SimulationInstance, comp_name::Symbol, datum_name::Symbol, results_output_dir::String; model_index::Int = 1, scen_name::Union{Nothing, String} = nothing)

    # scenario sub-folder
    scen_name !== nothing ? scen_dir = joinpath(results_output_dir,scen_name) : scen_dir = results_output_dir

    # model sub-foder
    if (length(sim_inst.results) > 1)
        model_dir = joinpath(scen_dir, "model_$(model_index)")
    else
        model_dir = scen_dir
    end

    # load results
    filepath = joinpath(model_dir, string(comp_name, "_", datum_name, ".csv"));
    if isfile(filepath)
        results = CSVFiles.load(filepath) |> DataFrame
    else
        error("$filepath does not exist, check arguments which currently set comp_name to $(comp_name), datum_name to $(datum_name), results_output_dir to $outputdir, model_index to $(model_index), and scen_name to $(scen_name)")
    end

    return results
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct VSCDCLine <: DCBranch
        name::String
        available::Bool
        active_power_flow::Float64
        arc::Arc
        rectifier_tap_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        rectifier_xrc::Float64
        rectifier_firing_angle::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        inverter_tap_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        inverter_xrc::Float64
        inverter_firing_angle::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        services::Vector{Service}
        ext::Dict{String, Any}
        time_series_container::InfrastructureSystems.TimeSeriesContainer
        internal::InfrastructureSystemsInternal
    end

As implemented in Milano's Book, Page 397.

# Arguments
- `name::String`
- `available::Bool`
- `active_power_flow::Float64`
- `arc::Arc`
- `rectifier_tap_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`
- `rectifier_xrc::Float64`
- `rectifier_firing_angle::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`
- `inverter_tap_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`
- `inverter_xrc::Float64`
- `inverter_firing_angle::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`
- `services::Vector{Service}`: Services that this device contributes to
- `ext::Dict{String, Any}`
- `time_series_container::InfrastructureSystems.TimeSeriesContainer`: internal time_series storage
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct VSCDCLine <: DCBranch
    name::String
    available::Bool
    active_power_flow::Float64
    arc::Arc
    rectifier_tap_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    rectifier_xrc::Float64
    rectifier_firing_angle::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    inverter_tap_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    inverter_xrc::Float64
    inverter_firing_angle::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Services that this device contributes to"
    services::Vector{Service}
    ext::Dict{String, Any}
    "internal time_series storage"
    time_series_container::InfrastructureSystems.TimeSeriesContainer
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function VSCDCLine(name, available, active_power_flow, arc, rectifier_tap_limits, rectifier_xrc, rectifier_firing_angle, inverter_tap_limits, inverter_xrc, inverter_firing_angle, services=Device[], ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), )
    VSCDCLine(name, available, active_power_flow, arc, rectifier_tap_limits, rectifier_xrc, rectifier_firing_angle, inverter_tap_limits, inverter_xrc, inverter_firing_angle, services, ext, time_series_container, InfrastructureSystemsInternal(), )
end

function VSCDCLine(; name, available, active_power_flow, arc, rectifier_tap_limits, rectifier_xrc, rectifier_firing_angle, inverter_tap_limits, inverter_xrc, inverter_firing_angle, services=Device[], ext=Dict{String, Any}(), time_series_container=InfrastructureSystems.TimeSeriesContainer(), internal=InfrastructureSystemsInternal(), )
    VSCDCLine(name, available, active_power_flow, arc, rectifier_tap_limits, rectifier_xrc, rectifier_firing_angle, inverter_tap_limits, inverter_xrc, inverter_firing_angle, services, ext, time_series_container, internal, )
end

# Constructor for demo purposes; non-functional.
function VSCDCLine(::Nothing)
    VSCDCLine(;
        name="init",
        available=false,
        active_power_flow=0.0,
        arc=Arc(Bus(nothing), Bus(nothing)),
        rectifier_tap_limits=(min=0.0, max=0.0),
        rectifier_xrc=0.0,
        rectifier_firing_angle=(min=0.0, max=0.0),
        inverter_tap_limits=(min=0.0, max=0.0),
        inverter_xrc=0.0,
        inverter_firing_angle=(min=0.0, max=0.0),
        services=Device[],
        ext=Dict{String, Any}(),
        time_series_container=InfrastructureSystems.TimeSeriesContainer(),
    )
end

"""Get [`VSCDCLine`](@ref) `name`."""
get_name(value::VSCDCLine) = value.name
"""Get [`VSCDCLine`](@ref) `available`."""
get_available(value::VSCDCLine) = value.available
"""Get [`VSCDCLine`](@ref) `active_power_flow`."""
get_active_power_flow(value::VSCDCLine) = get_value(value, value.active_power_flow)
"""Get [`VSCDCLine`](@ref) `arc`."""
get_arc(value::VSCDCLine) = value.arc
"""Get [`VSCDCLine`](@ref) `rectifier_tap_limits`."""
get_rectifier_tap_limits(value::VSCDCLine) = value.rectifier_tap_limits
"""Get [`VSCDCLine`](@ref) `rectifier_xrc`."""
get_rectifier_xrc(value::VSCDCLine) = value.rectifier_xrc
"""Get [`VSCDCLine`](@ref) `rectifier_firing_angle`."""
get_rectifier_firing_angle(value::VSCDCLine) = value.rectifier_firing_angle
"""Get [`VSCDCLine`](@ref) `inverter_tap_limits`."""
get_inverter_tap_limits(value::VSCDCLine) = value.inverter_tap_limits
"""Get [`VSCDCLine`](@ref) `inverter_xrc`."""
get_inverter_xrc(value::VSCDCLine) = value.inverter_xrc
"""Get [`VSCDCLine`](@ref) `inverter_firing_angle`."""
get_inverter_firing_angle(value::VSCDCLine) = value.inverter_firing_angle
"""Get [`VSCDCLine`](@ref) `services`."""
get_services(value::VSCDCLine) = value.services
"""Get [`VSCDCLine`](@ref) `ext`."""
get_ext(value::VSCDCLine) = value.ext
"""Get [`VSCDCLine`](@ref) `time_series_container`."""
get_time_series_container(value::VSCDCLine) = value.time_series_container
"""Get [`VSCDCLine`](@ref) `internal`."""
get_internal(value::VSCDCLine) = value.internal

"""Set [`VSCDCLine`](@ref) `name`."""
set_name!(value::VSCDCLine, val) = value.name = val
"""Set [`VSCDCLine`](@ref) `available`."""
set_available!(value::VSCDCLine, val) = value.available = val
"""Set [`VSCDCLine`](@ref) `active_power_flow`."""
set_active_power_flow!(value::VSCDCLine, val) = value.active_power_flow = set_value(value, val)
"""Set [`VSCDCLine`](@ref) `arc`."""
set_arc!(value::VSCDCLine, val) = value.arc = val
"""Set [`VSCDCLine`](@ref) `rectifier_tap_limits`."""
set_rectifier_tap_limits!(value::VSCDCLine, val) = value.rectifier_tap_limits = val
"""Set [`VSCDCLine`](@ref) `rectifier_xrc`."""
set_rectifier_xrc!(value::VSCDCLine, val) = value.rectifier_xrc = val
"""Set [`VSCDCLine`](@ref) `rectifier_firing_angle`."""
set_rectifier_firing_angle!(value::VSCDCLine, val) = value.rectifier_firing_angle = val
"""Set [`VSCDCLine`](@ref) `inverter_tap_limits`."""
set_inverter_tap_limits!(value::VSCDCLine, val) = value.inverter_tap_limits = val
"""Set [`VSCDCLine`](@ref) `inverter_xrc`."""
set_inverter_xrc!(value::VSCDCLine, val) = value.inverter_xrc = val
"""Set [`VSCDCLine`](@ref) `inverter_firing_angle`."""
set_inverter_firing_angle!(value::VSCDCLine, val) = value.inverter_firing_angle = val
"""Set [`VSCDCLine`](@ref) `services`."""
set_services!(value::VSCDCLine, val) = value.services = val
"""Set [`VSCDCLine`](@ref) `ext`."""
set_ext!(value::VSCDCLine, val) = value.ext = val
"""Set [`VSCDCLine`](@ref) `time_series_container`."""
set_time_series_container!(value::VSCDCLine, val) = value.time_series_container = val


text:
module Cusp

using DiffEqBase: ODEProblem
using Parameters: @with_kw, @unpack
using StaticArrays: SVector
using Setfield: @lens

using ...Bifurcations: BifurcationProblem

@with_kw struct CuspParam{B1, B2, S}
    # parameters are chosen so that u=1 is a fixed point.
    β₁::B1 = 1 / 2
    β₂::B2 = 1 / 2
    s::S = -1
end

function f(u::Real, p, t)
    @unpack β₁, β₂, s = p
    return β₁ + β₂ * u + s * u^3
end
# http://www.scholarpedia.org/article/Cusp_bifurcation#One-dimensional_Case

f(u::SVector, p, t) = SVector(f(u[1], p, t))

make_prob(
        p = CuspParam();
        u0 = SVector(1.0),
        tspan = (0.0, 30.0),
        ode = ODEProblem(f, u0, tspan, p),
        param_axis = (@lens _.β₁),
        t_domain = (-2.0, 2.0),
        kwargs...) =
    BifurcationProblem(ode, param_axis, t_domain;
                       kwargs...)

prob = make_prob()
ode = prob.p.de_prob
u0 = ode.u0
param_axis = prob.p.param_axis

end  # module


text:
"""
Selectors indicate that index values are not indices, but points to
be selected from the dimension values, such as DateTime objects on a Time dimension.
"""
abstract type Selector{T} end

val(m::Selector) = m.val

(::Type{T})(args...) where T <: Selector = T{typeof(args)}(args)

"""
    At(x)

Selector that exactly matches the value on the passed-in dimensions, or throws an error.
For ranges and arrays, every value must match an existing value - not just the end points.
"""
struct At{T} <: Selector{T}
    val::T
end

"""
    Near(x)

Selector that selects the nearest index to its contained value(s)
"""
struct Near{T} <: Selector{T}
    val::T
end

"""
    Between(a, b)

Selector that retreive all indices located between 2 values.
"""
struct Between{T<:Union{Tuple{Any,Any},Nothing}} <: Selector{T}
    val::T
end
Between(x::Tuple) = Between{typeof(x)}(x)

@inline dims2indices(dim::AbDim, lookup::AbDim{<:Selector}, emptyval) = sel2indices(dim, val(lookup))

sel2indices(a, lookup) = sel2indices(dims(a), lookup)
sel2indices(dims::Tuple, lookup) = sel2indices(dims, (lookup,))
sel2indices(dims::Tuple, lookup::Tuple) =
    (sel2indices(dims[1], lookup[1]), sel2indices(tail(dims), tail(lookup))...)
sel2indices(dims::Tuple{}, lookup::Tuple{}) = ()
sel2indices(dim::AbDim, sel::At) = at(dim, val(sel))
sel2indices(dim::AbDim, sel::At{<:Tuple}) = [at.(Ref(dim), val(sel))...]
sel2indices(dim::AbDim, sel::At{<:AbstractVector}) = at.(Ref(dim), val(sel))
sel2indices(dim::AbDim, sel::Near) = near(dim, val(sel))
sel2indices(dim::AbDim, sel::Near{<:Tuple}) = [near.(Ref(dim), val(sel))...]
sel2indices(dim::AbDim, sel::Near{<:AbstractVector}) = near.(Ref(dim), val(sel))
sel2indices(dim::AbDim, sel::Between{<:Tuple}) = between(dim, val(sel))

at(dim::AbDim, selval) = at(val(dim), selval)
at(dim, selval) = begin
    ind = findfirst(x -> x == selval, val(dim))
    ind == nothing ? throw(ArgumentError("$selval not found in $dim")) : ind
end

near(list, selval) = begin
    list = val(list)
    ind = searchsortedfirst(list, selval)
    if ind <= firstindex(list)
        ind
    elseif abs(list[ind] - selval) < abs(list[ind-1] - selval)
        ind
    else
        ind - 1
    end
end

between(dim::AbDim, sel) = between(order(dim), val(dim), sel)
between(::Forward, dim, sel) =
    searchsortedfirst(dim, first(sel)):searchsortedlast(dim, last(sel))
between(::Reverse, dim, sel) =
    searchsortedfirst(dim, first(sel); rev=true):-1:searchsortedlast(dim, last(sel); rev=true)

Base.@propagate_inbounds Base.getindex(a::AbstractArray, I::Vararg{Selector}) =
    getindex(a, sel2indices(a, I)...)
Base.@propagate_inbounds Base.setindex!(a::AbstractArray, x, I::Vararg{Selector}) = begin
    println(sel2indices(a, I))
    setindex!(a, x, sel2indices(a, I)...)
end
Base.view(a::AbstractArray, I::Vararg{Selector}) =
    view(a, sel2indices(a, I)...)


text:
import InteractiveUtils: subtypes

"""Returns an array of all concrete subtypes of T."""
function get_all_concrete_subtypes(::Type{T}) where T
    return _get_all_concrete_subtypes(T)
end

function _get_all_concrete_subtypes(::Type{T}, sub_types=[]) where T
    for sub_type in subtypes(T)
        push!(sub_types, sub_type)
        if isabstracttype(sub_type)
            _get_all_concrete_subtypes(sub_type, sub_types)
        end
    end

    return sub_types
end

"""Returns an array of concrete types that are direct subtypes of T."""
function get_concrete_subtypes(::Type{T}) where T
    return [x for x in subtypes(T) if isconcretetype(x)]
end

"""Returns an array of abstract types that are direct subtypes of T."""
function get_abstract_subtypes(::Type{T}) where T
    return [x for x in subtypes(T) if isabstracttype(x)]
end

"""Returns an array of all super types of T."""
function supertypes(::Type{T}, types=[]) where T
    super = supertype(T)
    push!(types, super)
    if super == Any
        return types
    end

    supertypes(super, types)
end

"""Converts a DataType to a Symbol, stripping off the module name(s)."""
function type_to_symbol(data_type::DataType)
    text = string(data_type)
    index = findlast(".", text)
    if !isnothing(index)
        text = text[index.start + 1:end]
    end

    return Symbol(text)
end


text:
function gcd{T<:Integer}(a::T, b::T)
    neg = a < 0
    while b != 0
        t = b
        b = rem(a, b)
        a = t
    end
    g = abs(a)
    neg ? -g : g
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using CameraCalibrations, JLD2
import Pkg.TOML

constfile = joinpath(dirname(Base.find_package("SkyRoom")), "..", "constants.toml")
constants = TOML.parsefile(constfile)

picsfolder = constants["calibration"]["camera"]["picsfolder"]

mount_point = "/home/yakir/mnt"
if !ismount(mount_point)
    run(`sshfs skyroom:/ $mount_point`)
end

calibration_images = joinpath(mount_point, "tmp", picsfolder)
check = constants["calibration"]["camera"]["board"]["check"]
intrinsic = readdir(calibration_images, join = true)
extrinsic = intrinsic[1]
_tform, _itform, errors = CameraCalibrations.calibrationmatrices(check, extrinsic, intrinsic)

x = first.(_tform.tform)
y = last.(_tform.tform)
cols, rows = (_tform.cols, _tform.rows)
tform = Dict(pairs((; x, y, cols, rows)))
x = first.(_itform.itform)
y = last.(_itform.itform)
xs, ys = (_itform.xs, _itform.ys)
itform = Dict(pairs((; x, y, xs, ys)))

calibrationfolder = joinpath(mount_point, "home/yakir/", constants["calibration"]["camera"]["datafolder"])
isdir(calibrationfolder) || mkdir(calibrationfolder)
file = joinpath(calibrationfolder, "calibration_matrices.jld2")
@save file tform itform

run(`fusermount -u $mount_point`)

run(`ssh skyroom 'julia --project=skyroom/Project.toml -e "using SkyRoom; assemble_camera_calibration()"'`)


text:
using TuringModels

ProjDir = @__DIR__

# Can't really set a Uniform[-Inf,Inf] on σ

# Turing model
@model m8_4(y) = begin
    α₁ ~ Uniform(-3000, 1000)
    α₂ ~ Uniform(-1000, 3000)
    σ ~ truncated(Cauchy(0,1), 0, Inf)

    y .~ Normal(α₁ + α₂, σ)
end

# Observations

y = rand(Normal(0,1), 100);

# Sample

chns = sample(m8_4(y), Turing.NUTS(0.65), 2000);

# Results rethinking

m84rethinking = "
         mean      sd     5.5%   94.5% n_eff Rhat
 a1    -861.15 558.17 -1841.89  -31.04     7 1.43
 a2     861.26 558.17    31.31 1842.00     7 1.43
 sigma    0.97   0.07     0.89    1.09     9 1.17
";

# Describe the posterior samples

chns |> display

# End of `08/m8.4t.jl`


text:
@safetestset "iscreditcard" begin
    using Validator

    @test Validator.iscreditcard("375556917985515") === true
    @test Validator.iscreditcard("36050234196908") === true
    @test Validator.iscreditcard("4716461583322103") === true
    @test Validator.iscreditcard("4716-2210-5188-5662") === true
    @test Validator.iscreditcard("4929 7226 5379 7141") === true
    @test Validator.iscreditcard("5398228707871527") === true
    @test Validator.iscreditcard("6283875070985593") === true
    @test Validator.iscreditcard("6263892624162870") === true
    @test Validator.iscreditcard("6234917882863855") === true
    @test Validator.iscreditcard("6234698580215388") === true
    @test Validator.iscreditcard("6226050967750613") === true
    @test Validator.iscreditcard("6246281879460688") === true
    @test Validator.iscreditcard("2222155765072228") === true
    @test Validator.iscreditcard("2225855203075256") === true
    @test Validator.iscreditcard("2720428011723762") === true
    @test Validator.iscreditcard("2718760626256570") === true
    @test Validator.iscreditcard("6765780016990268") === true
    @test Validator.iscreditcard("4716989580001715211") === true
    @test Validator.iscreditcard("8171999927660000") === true
    @test Validator.iscreditcard("8171999900000000021") === true

    @test Validator.iscreditcard("foo") === false
    @test Validator.iscreditcard("5398228707871528") === false
    @test Validator.iscreditcard("2718760626256571") === false
    @test Validator.iscreditcard("2721465526338453") === false
    @test Validator.iscreditcard("2220175103860763") === false
    @test Validator.iscreditcard("375556917985515999999993") === false
    @test Validator.iscreditcard("899999996234917882863855") === false
    @test Validator.iscreditcard("prefix6234917882863855") === false
    @test Validator.iscreditcard("623491788middle2863855") === false
    @test Validator.iscreditcard("6234917882863855suffix") === false
    @test Validator.iscreditcard("4716989580001715213") === false
end


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

Common infrastructure shared by PNML objects and labels.
Some optional incidental bits are shared by most PNML objects are also collected here.
"""
struct ObjectCommon
    graphics::Maybe{Graphics}
    tools::Maybe{Vector{ToolInfo}}
    labels::Maybe{Vector{PnmlLabel}}
end

ObjectCommon(pdict::PnmlDict) = ObjectCommon(
    get(pdict, :graphics, nothing),
    get(pdict, :tools, nothing),
    get(pdict, :labels, nothing)
)
ObjectCommon() = ObjectCommon(nothing, nothing, nothing)

"Return `true` if has a `graphics` element."
has_graphics(oc::ObjectCommon) = oc.graphics !== nothing

"Return `true` if has a `tools` element."
has_tools(oc::ObjectCommon) = oc.tools !== nothing #!&& !isempty(oc.tools)

"Return `true` if there is a `labels` element."
has_labels(oc::ObjectCommon) = oc.labels !== nothing #!&& !isempty(oc.labels)

graphics(oc::ObjectCommon) = oc.graphics
tools(oc::ObjectCommon) = oc.tools
labels(oc::ObjectCommon) = oc.labels

Base.isempty(oc::ObjectCommon) = !(has_graphics(oc) ||
                                   (has_tools(oc) && !isempty(oc.tools)) ||
                                   (has_labels(oc) && !isempty(oc.labels)))

function Base.empty!(oc::ObjectCommon)
    #! isnothing(oc.graphics) || replace with Graphics()
    #has_tools(oc)
    if oc.tools !== nothing
        # JET needs help avoiding the Nothing union split.
        t::Vector{ToolInfo} = oc.tools
        empty!(t)
    end
    if has_labels(oc)
        l::Vector{PnmlLabel} = oc.labels
        empty!(l)
    end
end

function Base.append!(l::ObjectCommon, r::ObjectCommon)
    # In the flatten use-case do not overwrite scalars.
    # How useful is propagating scalars?
    # Note that ObjectCommon is immutable so these error.
    #if !has_name(l);     l.name     = r.name; end
    #if !has_graphics(l); l.graphics = r.graphics; end

    update_maybe!(l, r, :tools)
    update_maybe!(l, r, :labels)
end


text:
#======================================#
#             Random Utils             #
#======================================#

function thinking(_message ::String, _end_message ::String)
  for x in 1:rand(2:6) #amount of times we want the full string dots to show
    for i in 1:4 #the max amount of dots we actually want to show on the screen
      print("$_message $("." ^ (i-1)) \r") # Stack overflow answer SHOULD work
      sleep((rand(0.2:0.3)))
    end
  end
  println("$_end_message")
end

function determine_first_player()
  if rand(1:2) == 1
    print("The first player is you!")
    return 1
  else
    print("The first player is the computer!")
    return 2
  end
end

function display_info(_userCard, _compCard, _userChances, _compChances)
  print(
    """\n\t
    \t#==============#\t#==============#
    \t      You    \t\t     Computer
    \t              \t\t
    \t Health: $(_userCard.Health)   \t\t Health: $(_compCard.Health)
    \t Attack: $(_userCard.Attack)   \t\t Attack: $(_compCard.Attack)
    \t Defence: $(_userCard.Defence)  \t\t Defence: $(_compCard.Defence)
    \t Speed: $(_userCard.Speed)    \t\t Speed: $(_compCard.Speed)
    \t Crit: $(_userCard.Crit)     \t\t Crit: $(_compCard.Crit)
    \t#==============#\t#==============#\n\n
    """)
end

max_stars_repo_name: NarutoXY/austroporean-cards

max_issues_repo_name: NarutoXY/austroporean-cards

max_forks_repo_name: NarutoXY/austroporean-cards


text:
nums = parse.(Int, readlines("day9/input.txt"))

function satisfies_preamble(line_num)
    for i in line_num-25:line_num-1, j in i+1:line_num-1
        if nums[i]+nums[j] == nums[line_num]
            return true
        end
    end
    return false
end

println(nums[25+findfirst(!satisfies_preamble, 26:length(nums))])