text:
# power_to_socp
#
# Functions that recursively simplify inequalities of the form
#
# x^n <= t^{n - m} s^m, t >= 0, s >= 0
#
# into a series of equations of the form x^2 <= st, which is
# second-order cone representable as
#
# norm([2x, t - s]) <= t + s, t >= 0, s >= 0.
#
# This reduction is documented in the pdf available at
# https://github.com/jump-dev/Convex.jl/raw/master/docs/supplementary/rational_to_socp.pdf

module psocp

mutable struct InequalityExpression
    # Represents an inequality of the form
    #
    # x^n <= t^{p_1} s^{p_2} u^{p_3},
    #
    # where we always have p_1 + p_2 + p_3 == n, and we must have p_3
    # belonging to {0, 1, 2}. This allows a particular form of recursive
    # enumeration of inequalities.

    # The powers p_1, p_2, p_3
    power1::Int
    power2::Int
    power3::Int

    # Index for variable represented by x in the above equation
    left_var_ind::Int
    # Index for variable represented by t
    t_var_ind::Int
    # Index for variable represented by s
    s_var_ind::Int
    # Index for variable represented by u (may be -1 if power3 == 0).
    u_var_ind::Int
end

struct SimpleInequalityExpression
    # Represents an inequality of the form
    #
    # x^2 <= t s,
    #
    # which is representable as the SOCP norm([2 x, t - s]) <= t + s. A
    # list of these is the first of the tuple returned by the method
    # ProductToSimpleInequalities.

    # Index for variable represented by x in the above equation
    left_var_ind::Int
    # Index for variable represented by t
    t_var_ind::Int
    # Index for variable represented by s
    s_var_ind::Int
end

# ProductToSimpleInequalities(first_power, second_power)
#
# Computes a sequence of inequalities of the form x^2 <= st, beginning
# from an initial inequality of the form
#
# x^n <= t^{p1} s^{p2},
#
# where n = p1 + p2 and p1 = first_power, p2 = second_power.
#
# Returns the tuple (inequalities, vars). The first of the tuple (the
# list inequalities) of which is an array of
# SimpleInequalityExpressions that consists of all the variable
# inequalities, the second of the tuple (the array vars) is an array
# of all the variable indices in the inequalities.
function ProductToSimpleInequalities(first_power::Int, second_power::Int)
    # Construct the first InequalityExpression, which is an inequality
    # of the form x^n <= t^p1 s^p2.
    @assert first_power > 0 && second_power > 0;
    n = first_power + second_power;
    var_list = [1, 2, 3];
    init_inequality = InequalityExpression(first_power, second_power, 0,
                                                                                1, 2, 3, -1);
    return ReducePowers(init_inequality, Array{SimpleInequalityExpression}(undef, 0),
                                            var_list);
end

# ReducePowers(curr_inequality::InequalityExpression,
#              ineq_array::Array{SimpleInequalityExpression, 1},
#              variable_array::Array{Int64, 1})
#
# Takes an initial inequality (curr_inequality) representing an
# expression of the form
#
# x^n <= t^{p1} s^{p2} u^{p3},
#
# where p3 is one of 0, 1, or 2, and n = p1 + p2 + p3, and recursively
# reduces it (see the file method-notes.tex) until the inequality
# array is fully populated by inequalities of the form x^2 <= t s,
# which are SOCP representable.
#
# The inequality array is populated (should be initially called with
# an empty array of SimpleInequalityExpressions), as is the
# variable_array, and returns the pair (ineq_array, variable_array)
# once both are fully populated.
function ReducePowers(curr_inequality::InequalityExpression,
                      ineq_array::Array{SimpleInequalityExpression, 1},
                      variable_array::Array{Int, 1})
    ReArrangePowers!(curr_inequality);
    p1 = curr_inequality.power1;
    p2 = curr_inequality.power2;
    p3 = curr_inequality.power3;
    n = (p1 + p2 + p3);
    @assert p1 >= 1;  # , "Must have at least 1 on first power");
    @assert p2 >= 1;  # , "Must have at least 1 on second power");
    # Evaluate cases for variables
    if (p3 == 0)
        # Double check if we have p1 == 1, p2 == 1
        if (p1 == 1 && p2 == 1)
            new_ineq = SimpleInequalityExpression(curr_inequality.left_var_ind,
                                                  curr_inequality.t_var_ind,
                                                  curr_inequality.s_var_ind);
            return (vcat(ineq_array, new_ineq), variable_array);
        end
        return ReduceThirdZero(curr_inequality,
                               ineq_array, variable_array);
    elseif (p3 == 1)
        return ReduceThirdOne(curr_inequality, ineq_array, variable_array);
    elseif (p3 == 2)
        return ReduceThirdTwo(curr_inequality, ineq_array, variable_array);
    else
        error("Third power in inequality (", curr_inequality.power3,
              ") must be 0, 1, or 2");
    end
end

# ReArrangePowers!(current_inequality::InequalityExpression)
#
# Function that takes in the current inequality, then re-orders
# variables so that the variable with smallest power is stored as
# u_var_ind. Also divides each power by their gcd.
function ReArrangePowers!(inequality::InequalityExpression)
    p1 = inequality.power1;
    p2 = inequality.power2;
    p3 = inequality.power3;
    divisor = gcd(p1, p2, p3);
    if (divisor > 1)
        inequality.power1 /= divisor;
        inequality.power2 /= divisor;
        inequality.power3 /= divisor;
        p1 = inequality.power1;
        p2 = inequality.power2;
        p3 = inequality.power3;
    end
    if (p1 <= min(p2, p3))
        # Swap t and u
        inequality.power3 = p1;
        inequality.power1 = p3;
        ind1 = inequality.t_var_ind;
        ind3 = inequality.u_var_ind;
        inequality.t_var_ind = ind3;
        inequality.u_var_ind = ind1;
    elseif (p2 <= min(p1, p3))
        # Swap s and u
        inequality.power3 = p2;
        inequality.power2 = p3;
        ind2 = inequality.s_var_ind;
        ind3 = inequality.u_var_ind;
        inequality.s_var_ind = ind3;
        inequality.u_var_ind = ind2;
    end  # In this case, p3 <= min(p1, p2), so no swap necessary
end

# ReduceThirdZero(curr_inequality::InequalityExpression,
#                 ineq_array::Array{SimpleInequalityExpression, 1},
#                 var_array::Array{Int64, 1})
#
# Recursively reduces powers when the third power (u) in
# curr_inequality is equal to 1.
function ReduceThirdZero(curr_inequality::InequalityExpression,
                         ineq_array::Array{SimpleInequalityExpression, 1},
                         var_array::Array{Int, 1})
    p1 = curr_inequality.power1;
    p2 = curr_inequality.power2;
    p3 = curr_inequality.power3;
    n = (p1 + p2);
    @assert p3 == 0;
    @assert n >= 2;
    if (mod(n, 2) == 0)
        # n is even, so check even-ness of power1, power2
        if (p1 == p2)
            # We are seriously done
            new_ineq = SimpleInequalityExpression(curr_inequality.left_var_ind,
                                                  curr_inequality.t_var_ind,
                                                  curr_inequality.s_var_ind);
            return (vcat(ineq_array, new_ineq), var_array);
        end
        if (mod(p1, 2) == 0)
            # Both are even, reduce and repeat
            curr_inequality.power1 /= 2;
            curr_inequality.power2 /= 2;
            return ReducePowers(curr_inequality, ineq_array, var_array);
        else
            # p1 is odd, p2 is odd, and they are unequal. Go to two
            # equations, which are
            #
            # x^(n / 2) <= t^((p1 - 1)/2) s^((p2 - 1)/2) u,
            # u^2 <= ts
            #
            # If n/2 == 2, then do not recurse, but simply catenate the
            # inequality on and return
            new_ind = var_array[end] + 1;
            new_ineq = SimpleInequalityExpression(new_ind,
                                                  curr_inequality.t_var_ind,
                                                  curr_inequality.s_var_ind);
            curr_inequality.power1 = (p1 - 1) / 2;
            curr_inequality.power2 = (p2 - 1) / 2;
            curr_inequality.power3 = 1;
            curr_inequality.u_var_ind = new_ind;
            # Recurse on the reduced inequality
            return ReducePowers(curr_inequality,
                                vcat(ineq_array, new_ineq),
                                vcat(var_array, new_ind));
        end  # if (mod(p1, 2) == 0)
    else    # n is odd, and we must recurse differently.
        if (mod(p1, 2) == 1)
            new_ind = var_array[end] + 1;
            new_ineq = SimpleInequalityExpression(new_ind,
                                                  curr_inequality.t_var_ind,
                                                  curr_inequality.left_var_ind);
            curr_inequality.power1 = (p1 - 1) / 2;
            curr_inequality.power2 = p2 / 2;
            curr_inequality.power3 = 1;
            curr_inequality.u_var_ind = new_ind;
            return ReducePowers(curr_inequality,
                                vcat(ineq_array, new_ineq),
                                vcat(var_array, new_ind));
        else    # mod(p2, 2) == 1
            new_ind = var_array[end] + 1;
            new_ineq = SimpleInequalityExpression(new_ind,
                                                  curr_inequality.s_var_ind,
                                                  curr_inequality.left_var_ind);
            curr_inequality.power1 = p1 / 2;
            curr_inequality.power2 = (p2 - 1) / 2;
            curr_inequality.power3 = 1;
            curr_inequality.u_var_ind = new_ind;
            return ReducePowers(curr_inequality,
                                vcat(ineq_array, new_ineq),
                                vcat(var_array, new_ind));
        end
    end  # if mod(n, 2) == 0
end

# ReduceThirdOne(curr_inequality::InequalityExpression,
#                ineq_array::Array{SimpleInequalityExpression, 1},
#                var_array::Array{Int64, 1})
#
# Recursively reduces powers when the third power (u) in
# curr_inequality is equal to 1.
function ReduceThirdOne(curr_inequality::InequalityExpression,
                        ineq_array::Array{SimpleInequalityExpression, 1},
                        var_array::Array{Int, 1})
    p1 = curr_inequality.power1;
    p2 = curr_inequality.power2;
    p3 = curr_inequality.power3;
    n = (p1 + p2 + p3);
    @assert p3 == 1; # , "Must have third power 1");
    @assert n >= 3; #  "Must have power n >= 3");
    if (mod(n, 2) == 0)
        # Exactly one of p1, p2 is odd, find it and reduce
        if (mod(p1, 2) == 1)
            new_ind = var_array[end] + 1;
            new_ineq = SimpleInequalityExpression(new_ind,
                                                  curr_inequality.t_var_ind,
                                                  curr_inequality.u_var_ind);
            curr_inequality.power1 = (p1 - 1)/2;
            curr_inequality.power2 = (p2 / 2);
            curr_inequality.power3 = 1;
            curr_inequality.u_var_ind = new_ind;
            return ReducePowers(curr_inequality,
                                vcat(ineq_array, new_ineq),
                                vcat(var_array, new_ind));
        else    # mod(p2, 2) == 1
            new_ind = var_array[end] + 1;
            new_ineq = SimpleInequalityExpression(new_ind,
                                                  curr_inequality.s_var_ind,
                                                  curr_inequality.u_var_ind);
            curr_inequality.power1 = p1 / 2;
            curr_inequality.power2 = (p2 - 1) / 2;
            curr_inequality.power3 = 1;
            curr_inequality.u_var_ind = new_ind;
            return ReducePowers(curr_inequality,
                                vcat(ineq_array, new_ineq),
                                vcat(var_array, new_ind));
        end
    else    # mod(n, 2) == 1
        # Either both of p1, p2 are odd or both are even
        if (mod(p1, 2) == 0)
            new_ind = var_array[end] + 1;
            new_ineq = SimpleInequalityExpression(new_ind,
                                                  curr_inequality.left_var_ind,
                                                  curr_inequality.u_var_ind);
            curr_inequality.power1 = p1 / 2;
            curr_inequality.power2 = p2 / 2;
            curr_inequality.power3 = 1;
            curr_inequality.u_var_ind = new_ind;
            return ReducePowers(curr_inequality,
                                vcat(ineq_array, new_ineq),
                                vcat(var_array, new_ind));
        else    # mod(p1, 2) == 1, so mod(p2, 2) == 1
            if (p1 == 1 && p2 == 1)
                # There is a simpler reduction to return: we replace
                #
                # x^3 <= tsu
                #
                # with the three inequalities
                #
                # x^4 <= w^2 v^2,  w^2 <= ts,  v^2 <= ux,
                #
                # which is equivalent to x^2 <= wv, w^2 <= ts, v^2 <= ux
                new_w_ind = var_array[end] + 1;
                new_v_ind = var_array[end] + 2;
                new_x_ineq = SimpleInequalityExpression(curr_inequality.left_var_ind,
                                                        new_w_ind, new_v_ind);
                new_w_ineq = SimpleInequalityExpression(new_w_ind,
                                                        curr_inequality.t_var_ind,
                                                        curr_inequality.s_var_ind);
                new_v_ineq = SimpleInequalityExpression(new_v_ind,
                                                        curr_inequality.u_var_ind,
                                                        curr_inequality.left_var_ind);
                return (vcat(ineq_array, [new_x_ineq, new_w_ineq, new_v_ineq]),
                        vcat(var_array, [new_w_ind, new_v_ind]));
            else
                # Do the full reduction
                new_w_ind = var_array[end] + 1;
                new_v_ind = var_array[end] + 2;
                new_y_ind = var_array[end] + 3;
                new_y_ineq = SimpleInequalityExpression(new_y_ind,
                                                        new_w_ind, new_v_ind);
                new_w_ineq = SimpleInequalityExpression(new_w_ind,
                                                        curr_inequality.t_var_ind,
                                                        curr_inequality.s_var_ind);
                new_v_ineq = SimpleInequalityExpression(new_v_ind,
                                                        curr_inequality.u_var_ind,
                                                        curr_inequality.left_var_ind);
                curr_inequality.power1 = (p1 - 1) / 2;
                curr_inequality.power2 = (p2 - 1) / 2;
                curr_inequality.power3 = 2;
                curr_inequality.u_var_ind = new_y_ind;
                return ReducePowers(curr_inequality,
                                    vcat(ineq_array, [new_y_ineq, new_w_ineq, new_v_ineq]),
                                    vcat(var_array, [new_w_ind, new_v_ind, new_y_ind]));
            end
        end
    end
end

# ReduceThirdTwo(curr_inequality::InequalityExpression,
#                ineq_array::Array{SimpleInequalityExpression, 1},
#                var_array::Array{Int64, 1})
#
# Recursively reduces powers when the third power (u) in
# curr_inequality is equal to 2.
function ReduceThirdTwo(curr_inequality::InequalityExpression,
                        ineq_array::Array{SimpleInequalityExpression, 1},
                        var_array::Array{Int, 1})
    p1 = curr_inequality.power1;
    p2 = curr_inequality.power2;
    p3 = curr_inequality.power3;
    n = (p1 + p2 + p3);
    @assert p3 == 2;  # , "Must have third power 2");
    @assert n >= 6;  # , "Must have power n >= 6");
    @assert p1 >= 2 && p2 >= 2;  # , "Did not rearrange powers properly");
    if (mod(n, 2) == 0)
        # Either both p1, p2 are even or both are odd
        if (mod(p1, 2) == 0)
            # Simply reduce powers on everything by a factor of two
            curr_inequality.power1 /= 2;
            curr_inequality.power2 /= 2;
            curr_inequality.power3 /= 2;
            return ReducePowers(curr_inequality, ineq_array, var_array);
        else
            # Both are odd, so introduce two new variables
            new_y_ind = var_array[end] + 1;
            new_w_ind = var_array[end] + 2;
            new_y_ineq = SimpleInequalityExpression(new_y_ind,
                                                    curr_inequality.u_var_ind,
                                                    new_w_ind);
            new_w_ineq = SimpleInequalityExpression(new_w_ind,
                                                    curr_inequality.s_var_ind,
                                                    curr_inequality.t_var_ind);
            curr_inequality.power1 = (p1 - 1) / 2;
            curr_inequality.power2 = (p2 - 1) / 2;
            curr_inequality.power3 = 2;
            curr_inequality.u_var_ind = new_y_ind;
            return ReducePowers(curr_inequality,
                                vcat(ineq_array, [new_y_ineq, new_w_ineq]),
                                vcat(var_array, [new_y_ind, new_w_ind]));
        end
    else    # mod(n, 2) == 1, so exactly one of p1 and p2 is odd
        if (mod(p1, 2) == 1)
            new_y_ind = var_array[end] + 1;
            new_w_ind = var_array[end] + 2;
            new_y_ineq = SimpleInequalityExpression(new_y_ind,
                                                    curr_inequality.u_var_ind,
                                                    new_w_ind);
            new_w_ineq = SimpleInequalityExpression(new_w_ind,
                                                    curr_inequality.left_var_ind,
                                                    curr_inequality.t_var_ind);
            curr_inequality.power1 = (p1 - 1) / 2;
            curr_inequality.power2 = p2 / 2;
            curr_inequality.power3 = 2;
            curr_inequality.u_var_ind = new_y_ind;
            return ReducePowers(curr_inequality,
                                vcat(ineq_array, [new_y_ineq, new_w_ineq]),
                                vcat(var_array, [new_y_ind, new_w_ind]));
        else    # mod(p2, 2) == 1
            new_y_ind = var_array[end] + 1;
            new_w_ind = var_array[end] + 2;
            new_y_ineq = SimpleInequalityExpression(new_y_ind,
                                                    curr_inequality.u_var_ind,
                                                    new_w_ind);
            new_w_ineq = SimpleInequalityExpression(new_w_ind,
                                                    curr_inequality.left_var_ind,
                                                    curr_inequality.s_var_ind);
            curr_inequality.power1 = p1 / 2;
            curr_inequality.power2 = (p2 - 1) / 2;
            curr_inequality.power3 = 2;
            curr_inequality.u_var_ind = new_y_ind;
            return ReducePowers(curr_inequality,
                                vcat(ineq_array, [new_y_ineq, new_w_ineq]),
                                vcat(var_array, [new_y_ind, new_w_ind]));
        end
    end
end

end  # module PSOCP


text:
using CSV, DataFrames, Statistics, Random
using Gurobi, JuMP
using Test

include("../src/OCTHaGOn.jl")
global PROJECT_ROOT = @__DIR__

function train_transonic_tree()
    @info "Trains tree over transonic airfoil data, testing IAI ORT fitting,
           basic tree data structures, and ORT constraint generation."
    # Getting data from aerospace database.
    vks = ["Re", "thick", "M", "C_L"];
    X = CSV.read("../data/airfoil/airfoil_X.csv",
                 copycols=true, header=vks, delim=",");
    Y = CSV.read("../data/airfoil/airfoil_y.csv",
                 copycols=true, header=["C_D"], delim=",");
    (train_X, train_y), (test_X, test_y) = IAI.split_data(:regression, Matrix(X), Matrix(Y), train_proportion=.999);

    # Training tree
    lnr = OCTHaGOn.base_otr()
    IAI.set_params!(lnr, hyperplane_config=(sparsity=1,))
    IAI.fit!(lnr, train_X, train_y)
    IAI.write_json(string(PROJECT_ROOT, "/data/transonic_tree.json"), lnr)
    return true
end

function transonic_mio_model()
    @info "Generates and solves mio model over transonic airfoil data."
    lnr = IAI.read_json(string(PROJECT_ROOT, "/data/transonic_tree.json"));
    Re = Array(range(10000, stop=35000, step=5000));
    thick = [0.100,0.110,0.120,0.130,0.140,0.145];
    M = [0.4, 0.5, 0.6, 0.7, 0.8, 0.9];
    cl = Array(range(0.35, stop=0.70, step=0.05));
    vks = [Symbol("x",i) for i=1:4];
    pwlDict = OCTHaGOn.pwl_constraint_data(lnr, vks);
    upperDict, lowerDict = OCTHaGOn.trust_region_data(lnr, vks);
    # Generate MIO constraints from aerodynamics data
    m = Model()
    set_optimizer(m, SOLVER_SILENT);
    @variable(m, x[1:4])
    @variable(m, y)
    @constraint(m, log(minimum(Re)) <= x[1])
    @constraint(m, x[1] <= log(maximum(Re)))
    @constraint(m, log(minimum(thick)) <= x[2])
    @constraint(m, x[2] <= log(maximum(thick)))
    @constraint(m, log(minimum(M)) <= x[3])
    @constraint(m, x[3] <= log(maximum(M)))
    @constraint(m, log(minimum(cl)) <= x[4])
    @constraint(m, x[4] <= log(maximum(cl)))
    OCTHaGOn.add_regr_constraints!(m, x, y, lnr, vks, M=100., eq=true);
    # Setting optimization constraints (as a demonstration)
    # Re, thickness, M, cl
    @constraint(m, x[1] >= log(20000))
    @constraint(m, x[2] >= log(0.11))
    @constraint(m, x[3] >= log(0.7))
    @constraint(m, x[4] >= log(0.4))
    @objective(m, Min, y)
    status = optimize!(m)
    println("Actual Cd: ", exp(getvalue.(y)))
    println("Min Cd in data: ", 0.00371)
    return true
end

@test train_transonic_tree()

@test transonic_mio_model()

ext: jl


text:
#### load chemical shift and J-coupling information from JSOn.
function loadcouplinginfojson(save_path::String)::Tuple{Vector{Int}, Vector{Float64}, Vector{Tuple{Int,Int}}, Vector{Float64}}

    J = JSON.parsefile(save_path)

    H_IDs = Vector{Int}(undef, 0)
    H_css = Vector{Float64}(undef, 0)
    for dict in J["chemical shift"]
        id = convert(Int, dict["ID"])
        value = convert(Float64, dict["value"])
        push!(H_IDs, id)
        push!(H_css, value)
    end

    J_IDs = Vector{Tuple{Int,Int}}(undef, 0)
    J_vals = Vector{Float64}(undef, 0)
    for dict in J["J-coupling"]
        id1 = convert(Int, dict["ID1"])
        id2 = convert(Int, dict["ID2"])
        value = convert(Float64, dict["value"])
        push!(J_IDs, (id1, id2))
        push!(J_vals, value)
    end

    return H_IDs, H_css, J_IDs, J_vals
end

function constructspinsystemsgraphforcompound(N_vertices::Int,
    J_inds::Vector{Tuple{Int,Int}})

    g = Graphs.SimpleGraph(N_vertices)

    for k = 1:length(J_inds)
        i, j = J_inds[k]
        Graphs.add_edge!(g, i, j)
    end

    return g
end

function convertJIDstoJinds(J_IDs::Vector{Tuple{Int,Int}},
    H_IDs::Vector{Int})::Vector{Tuple{Int,Int}}

    J_inds = Vector{Tuple{Int,Int}}(undef, 0)
    for k = 1:length(J_IDs)

        i = findfirst(xx->xx==J_IDs[k][1], H_IDs)
        j = findfirst(xx->xx==J_IDs[k][2], H_IDs)

        push!(J_inds, (i,j))
    end

    return J_inds
end

"""
matchJlabels((label_pair_list::Vector{Tuple{Int,Int}}, search_list::Vector{Int})::Vector{Int}

Returns the indices of entries of `label_pair_list` that has both labels in the `search_list`.
"""
function matchJlabels(label_pair_list::Vector{Tuple{Int,Int}}, search_list::Vector{Int})::Vector{Int}

    inds = Vector{Int}(undef, 0)
    for k = 1:length(label_pair_list)
        i, j = label_pair_list[k]

        if any(i .== search_list) && any(j .== search_list)
            push!(inds, k)
        end
    end

    return inds
end

function setupcsJ(H_IDs, H_css::Vector{T}, J_IDs, J_vals;
    unique_cs_tol = 1e-6,
    cs_round_digits = 5) where T

    H_inds = collect(1:length(H_IDs))
    J_inds = convertJIDstoJinds(J_IDs, H_IDs)

    g = constructspinsystemsgraphforcompound(length(H_IDs), J_inds)
    systems_g = Graphs.connected_components(g)

    dict_H_inds_to_css = Dict(H_inds .=> H_css)

    J_inds_sys = Vector{Vector{Tuple{Int,Int}}}(undef, 0)
    J_IDs_sys = Vector{Vector{Tuple{Int,Int}}}(undef, 0)
    J_vals_sys = Vector{Vector{T}}(undef, 0)
    H_inds_sys = Vector{Vector{Int}}(undef, 0)
    cs_sys = Vector{Vector{T}}(undef, 0)

    H_inds_singlets = Vector{Vector{Int}}(undef, 0)
    cs_singlets = Vector{T}(undef, 0)

    for i = 1:length(systems_g)

        H_inds_i = systems_g[i]

        if length(H_inds_i) > 1

            cs = collect( dict_H_inds_to_css[H_inds_i[l]] for l = 1:length(H_inds_i) )

            # check if cs is a vector of the same value, up to an absolute tolerance of `unique_cs_tol`.
            if all(isapprox.(cs, cs[1], atol = unique_cs_tol))

                # add singlet.
                push!(H_inds_singlets, H_inds_i)
                push!(cs_singlets, dict_H_inds_to_css[H_inds_i[1]])
            else

                # add spin system.
                push!(H_inds_sys, H_inds_i)
                push!(cs_sys, cs)

                inds = matchJlabels(J_inds, H_inds_i)
                push!(J_inds_sys, J_inds[inds])
                push!(J_IDs_sys, J_IDs[inds])
                push!(J_vals_sys, J_vals[inds])
            end

        elseif length(H_inds_i) == 1

            # add singlet.
            push!(H_inds_singlets, H_inds_i)
            push!(cs_singlets, dict_H_inds_to_css[H_inds_i[1]])
        end
    end

    # remove singlets that have duplicate chemical shifts.
    H_inds_singlets, cs_singlets = removeredundantsinglets(H_inds_singlets, cs_singlets)

    #
    J_inds_sys_local = convertJindsglobaltolocal(H_inds_sys, J_inds_sys)

    # check.
    for i = 1:length(J_inds_sys)
        @assert length(J_inds_sys[i]) == length(J_inds_sys_local[i]) == length(J_IDs_sys[i]) == length(J_vals_sys[i])
        @assert length(H_inds_sys[i]) == length(cs_sys[i])
    end

    @assert length(H_inds_singlets) == length(cs_singlets)

    return J_inds_sys, J_inds_sys_local, J_IDs_sys, J_vals_sys, H_inds_sys,
        cs_sys, H_inds_singlets, cs_singlets, H_inds, J_inds, g
end

function removeredundantsinglets(H_inds::Vector{Vector{Int}}, cs_singlets::Vector{T}) where T <: Real

    cs_singlets_unique, inds_unique = uniqueinds(cs_singlets)

    H_inds_unique = Vector{Vector{Int}}(undef, 0)
    for i = 1:length(inds_unique)

        inds = inds_unique[i]
        push!(H_inds_unique, combinevectors(H_inds[inds]))
    end

    return H_inds_unique, cs_singlets_unique
end

"""
convertJindsglobaltolocal(H_inds_sys::Vector{Vector{Int}},
    J_inds_sys::Vector{Vector{Tuple{Int,Int}}})

The return type is Vector{Vector{Tuple{Int,Int}}}.

Local: Each spin system will have its own spin nucleui numbering that starts at 1.
Global: The spins systems will together have one single spin nucleui numbering that starts at 1.
"""
function convertJindsglobaltolocal(H_inds_sys::Vector{Vector{Int}},
    J_inds_sys::Vector{Vector{Tuple{Int,Int}}})

    @assert length(H_inds_sys) == length(J_inds_sys)
    N_systems = length(H_inds_sys)

    J_inds_sys_local = Vector{Vector{Tuple{Int,Int}}}(undef, N_systems)

    for i = 1:N_systems

        x = H_inds_sys[i]
        J = J_inds_sys[i]

        conversion_dict = Dict(x .=> collect(1:length(x)))

        J_inds_sys_local[i] = collect( (conversion_dict[J[k][1]], conversion_dict[J[k][2]]) for k = 1:length(J) )
    end

    return J_inds_sys_local
end


text:
"""
    greet()

Greets the user the first time they run JuliaTutor.
"""
greet() = println("Hi my name is JuliaTutor, \n
                Here you can interactively learn Julia through pre-made lessons.")

"""
    inform(lesson::Lesson)

A convenience function to display the prompt of a lesson to an end user.
"""
function inform(lesson::Lesson)
    print( Crayon(foreground = :blue, bold = true), "\n> ",
             Crayon(foreground = :white, bold = false), lesson.prompt)
end

"""
    request(lesson::Lesson)

A convenience function to display what action is requested for an end user to accomplish.
"""
function request(lesson::Lesson)
    print( Crayon( foreground = :red, bold = true ), "\n> " )
    print( Crayon( foreground = :white, bold = false ), lesson.request )
end

"""
    hint(lesson::Lesson)

Displays a hint for the user, if the lesson provides it. Otherwise - it directs the user to Discourse.
"""
function hint(lesson::Lesson)
    if length(lesson.hint) == 0
        print( Crayon( foreground = :red, italics = true ), """
        Sorry... No hints were given for this lesson.
        You can always check in with the Humans Of Julia Discord for assistance:
        https://discord.com/channels/762167454973296644/822525930579689552
        """ )
    else
        print( Crayon( foreground = :blue, italics = true ), "\n " * lesson.hint )
    end
    print( Crayon( foreground = :white, italics = false ), "\n" )
end

"""
    help()

Displays generic tips for how to navigate JuliaTutor.
"""
function help()
    println( "Enter `)` to exit or enter the tutor REPL." )
    println( "Ctrl + `c` will also exit the tutor REPL." )
    println( "Ctrl + `d` will end the REPL instance." )
    println( "Type `hint` to see if the author provided hints for a lesson." )
end


text:
function SparseArrays.sparse(model::MetNet, fields::Vector{Symbol} = [:S])
    net = Dict()
    for f in fields
        net[f] = getfield(model, f) |> sparse
    end
    return MetNet(model; reshape = false, net...)
end


text:
function HelmertTransformation(arg0::AbsoluteDate, arg1::jdouble, arg10::jdouble, arg11::jdouble, arg12::jdouble, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return HelmertTransformation((AbsoluteDate, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

function get_epoch(obj::HelmertTransformation)
    return jcall(obj, "getEpoch", AbsoluteDate, ())
end

function get_transform(obj::HelmertTransformation, arg0::AbsoluteDate)
    return jcall(obj, "getTransform", Transform, (AbsoluteDate,), arg0)
end

function get_transform(obj::HelmertTransformation, arg0::FieldAbsoluteDate)
    return jcall(obj, "getTransform", FieldTransform, (FieldAbsoluteDate,), arg0)
end


text:
function table(cur_matrix::Matrix; header_row=[], title=nothing)
    cur_table = ""
    if title != nothing
        cur_table *= "
$title
"
    end
    cur_table *= "
"
    if !isempty(header_row)
        cur_table *= "

"
        for cur_header in header_row
            cur_table *= "$cur_header	"
        end
        cur_table *= "

"
    end
    cur_table *= "
"
    for ii in 1:size(cur_matrix, 1)
        cur_table *= "
"
        for jj in 1:size(cur_matrix, 2)
            cur_table *= ""
            cur_table *= string(cur_matrix[ii, jj])
            cur_table *= "	"
        end
        cur_table *= "
"
    end
    cur_table *= "
"
    cur_table *= "
"
    return cur_table
end
function show(x::Second)::String
    ps = canonicalize(Dates.CompoundPeriod(Second(x)))
    a = Dict{DataType, Int}(k => 0 for k in [Hour, Minute, Second])
    ts = [Day, Week, Month, Year]
    for p in ps.periods
        if typeof(p) in ts
            a[Hour] += Hour(p).value
        else
            a[typeof(p)] += p.value
        end
    end
    return @sprintf "%i:%02i:%02i" a[Hour] a[Minute] a[Second]
end

function show(x::Time)::String
    j = Hour(x).value
    i, ampm = j > 12 ? (j - 12, :pm) : (j, :am)
    s = Minute(x).value >= 30 ? 0x1F550+(i-1)+12 : 0x1F550+(i-1)
    t = uppercase(num2hex(s)[end-4:end])
    return string("&#x$t; $ampm")
end

function report(folder::String)
    a = loadAssociation(folder)
    delete_empty_metadata!(a)

    rep = counter(Vector{String})
    for r in a.runs
        k = collect(values(r.run.metadata))
        push!(rep, k)
    end

    durs = Dict(poi.name => counter(Vector{String}) for poi in a.pois)
    coun = Dict(poi.name => counter(Vector{String}) for poi in a.pois)
    cosa = Dict(poi.name => Accumulator(Vector{String}, Float64) for poi in a.pois)
    sina = Dict(poi.name => Accumulator(Vector{String}, Float64) for poi in a.pois)
    for (p, r) in a.associations
        k = collect(values(r.run.metadata))
        push!(durs[p.name], k, duration(p, folder))
        push!(coun[p.name], k)
        t = timeofday(p, folder)
        # println(t)
        α = convert(u"rad", t)
        push!(cosa[p.name], k, cos(α))
        push!(sina[p.name], k, sin(α))
    end
    dur = Dict(k1 => Dict(k2 => Second(round(Int, v2/coun[k1][k2])) for (k2, v2) in v1) for (k1, v1) in durs)
    tim = Dict(k1 => Dict(k2 => convert(Time, atan2(u"rad", sina[k1][k2], v2)) for (k2, v2) in v1) for (k1, v1) in cosa)

    m = Matrix{String}[]
    for (k, v) in rep
        run = [k; string(v)]
        poi = [haskey(v1, k) ? show(v1[k])*" "*show(tim[k1][k]) : "" for (k1, v1) in dur]
        push!(m, reshape([run; poi], (1,:)))
    end
    m = vcat(m...)

    txt = table(m, header_row = [String.(collect(keys(a.runs[1].run.metadata))); "Repetition"; collect(keys(dur))], title = "Summary for $folder")
    return """<!DOCTYPE html>      $txt  """
end

function report_gui(folder::String)
    filename = joinpath(folder, "report.html")
    open(filename, "w") do o
        print(o, report(folder))
    end
    openit(filename)
end


text:
"""
Wraps DynamicBranch devices from PowerSystems to handle changes in controls and connection
status, and allocate the required indexes of the state space.
"""
struct BranchWrapper
    branch::PSY.DynamicBranch
    system_base_power::Float64
    system_base_frequency::Float64
    connection_status::Base.RefValue{Float64}
    bus_ix_from::Int
    bus_ix_to::Int
    ix_range::Vector{Int}
    ode_range::Vector{Int}
    global_index::Base.ImmutableDict{Symbol, Int}
    function BranchWrapper(
        branch::PSY.DynamicBranch,
        bus_ix_from::Int,
        bus_ix_to::Int,
        ix_range,
        ode_range,
        sys_base_power,
        sys_base_freq,
    )
        branch_states = PSY.get_states(branch)
        new(
            branch,
            sys_base_power,
            sys_base_freq,
            Base.Ref(1.0),
            bus_ix_from,
            bus_ix_to,
            ix_range,
            ode_range,
            Base.ImmutableDict(Dict(branch_states .=> ix_range)...),
        )
    end
end

get_connection_status(wrapper::BranchWrapper) = wrapper.connection_status[]
get_bus_ix_from(wrapper::BranchWrapper) = wrapper.bus_ix_from
get_bus_ix_to(wrapper::BranchWrapper) = wrapper.bus_ix_to
get_ix_range(wrapper::BranchWrapper) = wrapper.ix_range
get_ode_ouput_range(wrapper::BranchWrapper) = wrapper.ode_range
get_global_index(wrapper::BranchWrapper) = wrapper.global_index

get_system_base_power(wrapper::BranchWrapper) = wrapper.system_base_power
get_system_base_frequency(wrapper::BranchWrapper) = wrapper.system_base_frequency

PSY.get_name(wrapper::BranchWrapper) = PSY.get_name(wrapper.branch)
PSY.get_available(wrapper::BranchWrapper) = PSY.get_available(wrapper.branch)
PSY.get_active_power_flow(wrapper::BranchWrapper) = PSY.get_active_power(wrapper.branch)
PSY.get_reactive_power_flow(wrapper::BranchWrapper) = PSY.get_reactive_power(wrapper.branch)
PSY.get_arc(wrapper::BranchWrapper) = PSY.get_arc(wrapper.branch)
PSY.get_r(wrapper::BranchWrapper) = PSY.get_r(wrapper.branch)
PSY.get_x(wrapper::BranchWrapper) = PSY.get_x(wrapper.branch)
PSY.get_b(wrapper::BranchWrapper) = PSY.get_b(wrapper.branch)
PSY.get_rate(wrapper::BranchWrapper) = PSY.get_rate(wrapper.branch)
PSY.get_angle_limits(wrapper::BranchWrapper) = PSY.get_angle_limits(wrapper.branch)
PSY.get_ext(wrapper::BranchWrapper) = PSY.get_ext(wrapper.branch)
PSY.get_states(wrapper::BranchWrapper) = PSY.get_states(wrapper.branch)
PSY.get_n_states(wrapper::BranchWrapper) = PSY.get_n_states(wrapper.branch)

function _add_dynamic_bus_states!(
    DAE_vector::Vector{Bool},
    voltage_buses_ix::Vector{Int},
    bus_ix::Int,
    n_buses::Int,
)
    push!(voltage_buses_ix, bus_ix)
    DAE_vector[bus_ix] = DAE_vector[bus_ix + n_buses] = true
    return
end


text:
#
# # inimize  trace(S*X) - log det X   subject to support(X) ⊆ non_zero_set
# function covsel_refit!{T<:AbstractFloat}(
#     X::StridedMatrix{T},
#     Z::StridedMatrix{T},
#     U::StridedMatrix{T},
#     S::StridedMatrix{T},
#     non_zero_set::Vector{Int64};
#     options::ADMMOptions = ADMMOptions()
#     )
#
#   maxiter = options.maxiter
#   ρ = options.ρ
#   α = options.α
#   abstol = options.abstol
#   reltol = options.reltol
#
#   p = size(S, 1)
#   tmpStorage = zeros(T, (p, p))
#   Zold = copy(Z)
#
#   for iter=1:maxiter
#     # x-update
#     @simd for i=1:length(tmpStorage)
#       @inbounds tmpStorage[i] = ρ * (Z[i] - U[i]) - S[i]
#     end
#     ef = eigfact(Symmetric(tmpStorage))
#     efVectors = ef[:vectors]::Array{T, 2}
#     efValues = ef[:values]::Array{T, 1}
#     @simd for i=1:p
#       @inbounds t = efValues[i]
#       @inbounds efValues[i] = (t + sqrt(t^2. + 4.*ρ)) / (2.*ρ)
#     end
#     @inbounds for c=1:p, r=1:p
#       X[r, c] = zero(T)
#       for i=1:p
#         X[r, c] = X[r, c] + efVectors[r,i] * efValues[i] * efVectors[c, i]
#       end
#     end
#
#     # z-update with relaxation
#     copy!(Zold, Z)
#     fill!(Z, zero(T))
#     for i in non_zero_set
#       @inbounds Z[i] = α*X[i] + (one(T)-α)*Z[i] + U[i]
#     end
#
#     # u-update
#     @simd for i in eachindex(X)
#       @inbounds U[i] = U[i] + X[i] - Z[i]
#     end
#
#     # check convergence
#     r_norm = _normdiff(X, Z)
#     s_norm = _normdiff(Z, Zold) * sqrt(ρ)
#     eps_pri = p*abstol + reltol * max( vecnorm(X), vecnorm(Z) )
#     eps_dual = p*abstol + reltol * ρ * vecnorm(U)
#     if r_norm < eps_pri && s_norm < eps_dual
#       break
#     end
#   end
#   Z
# end

# # inimize  sum_i trace(S_i*X_i) - log det X_i + lambda * sum_ab (sum_i X_i,ab ^2)
# function solve!{T<:FloatingPoint}(
#     solver::ADMMSolver,
#     X::Array{T, 3},
#     Z::Array{T, 3},
#     U::Array{T, 3},
#     S::Array{T, 3},
#     λ::T;
#     penalize_diag::Bool=true
#     )
#   maxiter = solver.maxiter
#   ρ = solver.ρ
#   α = solver.α
#   abstol = solver.abstol
#   reltol = solver.reltol
#   λρ = λ/ρ

#   K = size(S, 3)
#   p = size(S, 1)

#   tmpStorage = zeros(T, (p,p,K))
#   Zold = copy(Z)

#   for iter=1:maxiter

#     # x-update
#     @simd for i in eachindex(X)
#       tmpStorage[i] = ρ * (Z[i] - U[i]) - S[i]
#     end
#     for k=1:K
#       ef = eigfact(Symmetric(sub(tmpStorage, :, :, k)))
#       efVectors = ef[:vectors]::Array{T, 2}
#       efValues = ef[:values]::Array{T, 1}
#       for i=1:p
#         t = efValues[i]
#         efValues[i] = (t + sqrt(t^2. + 4.*ρ)) / (2.*ρ)
#       end
#       for c=1:p, r=1:p
#         X[r, c, k] = zero(T)
#         @simd for i=1:p
#           X[r, c, k] = X[r, c, k] + efVectors[r,i] * efValues[i] * efVectors[c, i]
#         end
#       end
#     end

#     # z-update with relaxation
#     copy!(Zold, Z)
#     @simd @inbounds for i in eachindex(tmpStorage)
#       tmpStorage[i] = α*X[i] + (one(T)-α)*Z[i] + U[i]
#     end
#     for c=1:p
#       for r=1:c-1
#         tmpSlice = slice(tmpStorage, r, c, :)
#         t = max(one(T) - λρ / vecnorm(tmpSlice), zero(T))
#         if t > zero(T)
#           @simd @inbounds for k=1:K
#             Z[r,c,k] = t * tmpSlice[k]
#             Z[c,r,k] = t * tmpSlice[k]
#           end
#         else
#           @simd @inbounds for k=1:K
#             Z[r,c,k] = zero(T)
#             Z[c,r,k] = zero(T)
#           end
#         end
#       end
#       if penalize_diag
#         tmpSlice = slice(tmpStorage, c, c, :)
#         t = max(one(T) - λρ / vecnorm(tmpSlice), zero(T))
#         if t > zero(T)
#           @simd @inbounds for k=1:K
#             Z[c,c,k] = t * tmpSlice[k]
#           end
#         else
#           @simd @inbounds for k=1:K
#             Z[c,r,k] = zero(T)
#           end
#         end
#       else
#         @simd @inbounds for k=1:K
#           Z[c, c, k] = tmpStorage[c, c, k]
#         end
#       end
#     end

#     # u-update
#     @simd @inbounds for i in eachindex
#       U[i] = tmpStorage[i] - Z[i]
#     end

#     # check convergence
#     r_norm = _normdiff(X, Z)
#     s_norm = _normdiff(Z, Zold) * sqrt(ρ)
#     eps_pri = p*abstol + reltol * max( vecnorm(X), vecnorm(Z) )
#     eps_dual = p*abstol + reltol * ρ * vecnorm(U)
#     if r_norm < eps_pri && s_norm < eps_dual
#       break
#     end
#   end
#   Z
# end


text:
function Convmtx(w,nr)
	w  = vec(w)
	nw = length(w)
	nd = nw + nr - 1
	cm = zeros(nd, nr)
	for i=1:nr
		cm[i:i+nw-1, i] = w
	end
	cm = sparse(cm)
	return cm
end

max_stars_repo_name: fercarozzi/myseismicjulia

max_issues_repo_name: fercarozzi/myseismicjulia

max_forks_repo_name: fercarozzi/myseismicjulia


text:
"""
    ms = ModelingStruct{T}(nt::NamedTuple) where {T}

A mutable struct-like data type with homogeneous data type `T` that can be accessed like a
    vector or a struct. Useful for modeling tasks such as differential equations or
    mathematical optimization that operate internally on data as a flat vector while calling
    user functions that could often be better expressed if the data was a named tuple.
    For a more convenient constructor, see `mstruct`.

# Examples

```julia-repl
julia> c = (a=2, b=[1, 2]);

julia> p = (a=1, b=[2, 1, 4], c=c)

julia> ms = ModelingStruct{Float32}(p)
ModelingStruct(a = 1.0f0, b = Float32[2.0, 1.0, 4.0], c = (a = 2.0f0, b = Float32[1.0, 2.0]))

julia> ms.c.a = 4; ms
ModelingStruct(a = 1.0f0, b = Float32[2.0, 1.0, 4.0], c = (a = 4.0f0, b = Float32[1.0, 2.0]))

julia> ms[5]
4.0f0
```
"""
struct ModelingStruct{T,V<:AbstractVector{T},K,P} <: DenseVector{T} where N
    data::V
    structure::NamedTuple{K,P}
    # The ::Val{false} means don't reattach the views, I'll change this later, probably
    function ModelingStruct(vect::V, nt::NamedTuple{K,P}, make_refs::Val{false}) where {T,V<:AbstractVector{T},K,P}
        return new{T,V,K,P}(vect, nt)
    end
end
function ModelingStruct(vect::AbstractVector{T}, nt::NamedTuple{K}) where {T,K}
    vals = values(nt)
    N = vals .|> recursive_length |> indices
    tup = map((val, n) -> attach(view(vect, n), val), vals, N)
    nt = NamedTuple{K}(tup)
    v = maybe_parent(vect)
    return ModelingStruct(v, nt, Val(false))
end
function ModelingStruct{T}(nt::NamedTuple{K}) where {T,K}
    return ModelingStruct(Vector{T}(undef,recursive_length(nt)), nt)
end
ModelingStruct(nt::NamedTuple) = ModelingStruct{Float64}(nt)
ModelingStruct{T}(; kwargs...) where T = ModelingStruct{T}((;kwargs...))
ModelingStruct(; kwargs...) = ModelingStruct{Float64}((;kwargs...))

# MStruct is a convenience alias for ModelingStruct
MStruct = ModelingStruct

function attach(v::AbstractVector, nt::NamedTuple)
    return ModelingStruct(v, nt)
end
function attach(v::AbstractVector{T}, a::AbstractVector{TT}) where {T,TT}
    # if !isconcretetype(eltype(a))
    #     types = typeof.(a)
    #     error("Vectors in ModelingStructs must be homogeneous, this one has types: \n", join(types, "\n"))
    # end
    N = recursive_length(a[1])
    views = [attach(view(v, 1:N), a[1])]
    idx_start = N+1
    for idx in 2:length(a)
        idx_end = idx_start + N - 1
        push!(views, attach(view(v, idx_start:idx_end), a[idx]))
        idx_start = idx_end+1
    end
    return views
end
function attach(v::AbstractVector{T}, a::AbstractVector{<:N}) where {T,N<:Number}
    v .= a
    return @view v[:]
end
function attach(v::AbstractArray{T}, a::Number) where T
    v[1] = a
    return view(v, 1)
end

max_stars_repo_name: jonniedie/NamedViewVectors


text:
## Step3 is the iteration to find the price q

#module Step3

using LinearAlgebra, Random, Distributions, Statistics, Plots
using BenchmarkTools, Interpolations

include("Step2.jl")
include("Step1Grid.jl")

"""
    Relies on step2 and finds the optimal price.
"""

# Define function that computes the market for savings
function mktclearing(A, E, pol_func, F_update)

    """
    Function that computes market clearing conditions.
    Inputs:
        - A: set of assets
        - E: set of shocks
        - pol_func: optimal policy function
        - F_update: stationary distribution function Ψ
    Output:
        - market: demand for assets
    """

    # retrieve dimensions
    a_size = length(A)
    e_size = size(E)[1]

    # initialize market for update
    market = 0.0

    for a = 2:a_size
        for e = 1:e_size

            market += pol_func[a, e] * (F_update[a, e] - F_update[a-1, e])

        end
    end

    return market

end

# Define function for Step3
function step3(A, E, Π, q0, β, σ; maxT = 100, tol_iterations = 0.01,
                tol_market = 0.0025, tol_q = 0.001, weight = 0.5)

    """
    Function that finds the optimal price q.
    Relies on market in the updating procedure.
    Inputs:
        - A: set of assets
        - E: set of shocks
        - q0: guess price
        - β: discount factor
        - σ: preference parameter
    Output:
        - q: optimal price
        - policy: optimal policy
        - Ψ: stationary distribution Ψ
    """

    # check q > β
    @assert q0 > β

    q = q0
    solStep1 = step1(A, E, Π, q, β, σ; tol = tol_iterations)
    solStep2 = step2(A, E, Π, solStep1; tol = tol_iterations)

    # Start price iteration
    normdiff = Inf
    t = 1

    while normdiff > tol_q && t <= maxT

        # solve step1
        solStep1 = step1(A, E, Π, q, β, σ; tol = tol_iterations)

        # solve step2
        solStep2 = step2(A, E, Π, solStep1; tol = tol_iterations)

        # compute market
        market = mktclearing(A, E, solStep1, solStep2)

        # update q based on market clearing:
        # increase guess if market > 0 (but just by a little bit)
        if market > 0 && abs(market) > tol_market
            q = q + (1 - weight) * tol_market

        # decrease guess if market < 0 (but just by a little bit)
        elseif market < 0 && abs(market) > tol_market
            q = q - (1 - weight) * tol_market

        # if abs(market) < tol_market, then no need to update q
        else
            q = q
        end

        t = t + 1
        normdiff = abs(q - q0)

        q0 = q

    end

    return (q = q, policy = solStep1, Ψ = solStep2)

end

#export mktclearing, step3

#end


text:
function calc_rho( pw::PWGrid, psi::Array{ComplexF64,2} )
    Ω = pw.Ω
    Ns = pw.Ns
    Npoints = pw.Npoints
    Nstates = size(psi)[2]
    #
    ρ = zeros(ComplexF64,Npoints)
    # Transform to real space
    psiR = G_to_R(Ns,psi)
    # orthonormalization in real space
    ortho_gram_schmidt!(Nstates,psiR)
    psiR = psiR*sqrt(Npoints/Ω)
    #scale!(sqrt(Npoints/Ω),psiR)
    for is = 1:Nstates
        for ip = 1:Npoints
            ρ[ip] = ρ[ip] + conj(psiR[ip,is])*psiR[ip,is]
        end
    end
    return real(ρ)
end


text:
function mysum_in(a::AbstractArray)
s = zero(eltype(a))
for i in a
s = s + i
end
end

max_stars_repo_path: Module 2/Chapter06/mysum_in.jl

max_issues_repo_path: Module 2/Chapter06/mysum_in.jl

max_forks_repo_path: Module 2/Chapter06/mysum_in.jl


text:
# Brownian Motion Plots

# Geometric Brownian Motion Plots


text:
using Behavior
using Behavior.Gherkin
using Niancat.Languages
using Niancat.Games.NiancatGames
using Niancat.Games.NiancatGames: Incorrect, Correct, Rejected, PuzzleIs, GetPuzzle, NoPuzzleSet
using Niancat.Games

@when("the puzzle is set to {String}") do context, puzzle
    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    response = gamecommand(context[:game], user, SetPuzzle(puzzle))
    context[:response] = response
end

@when("a user asks for the puzzle") do context
    team = Team(1, "defaultteam", "")
    user = User(1, "name", "displayname", team)
    response = gamecommand(context[:game], user, GetPuzzle())
    context[:response] = response
end

@then("the puzzle response is {String}") do context, puzzle
    @expect context[:response] == PuzzleIs(puzzle)
end

@given("that no puzzle has been set") do context
    # No action needed, the game does not have a puzzle
    # set by default after construction.
end

@then("the game responds that no puzzle has been set") do context
    @expect context[:response] == NoPuzzleSet()
end

@then("it is rejected") do context
    @expect context[:response] isa Rejected
end


text:
# Read and properly scale phase
fn_phase = "data/small/Phase.nii"
phase_nii = readphase(fn_phase)
@test maximum(phase_nii) ≈ π atol=2e-3
@test minimum(phase_nii) ≈ -π atol=2e-3

# Read and normalize mag
fn_mag = "data/small/Mag.nii"
mag_nii = readmag(fn_mag; rescale=true)
@test 1 ≤ maximum(mag_nii) ≤ 2
@test 0 ≤ minimum(mag_nii) ≤ 1

# sample
sample = MriResearchTools.sample
@test length(sample(1:10)) >= 10
@test 10 >= length(sample(1:10; n=3)) >= 3
@test isempty(sample([NaN]))
@test all(isfinite.(sample([1:10..., NaN])))
@test length(sample([1])) == 1
@test isempty(sample([]))

# estimatenoise
@test estimatenoise(mag_nii)[2] ≈ 0.03 atol=1e-2
R = rand(500, 500, 500)
R[:, 251:500, :] .= 10
μ, σ = estimatenoise(R)
@test μ ≈ 0.5 atol=1e-1
@test σ ≈ sqrt(1/12) atol=2e-2
R[1:10,:,:] .= NaN; R[:,1:10,:] .= NaN; R[:,:,1:10] .= NaN;
R[end-9:end,:,:] .= NaN; R[:,end-9:end,:] .= NaN; R[:,:,end-9:end] .= NaN
μ, σ = estimatenoise(R)
#@test μ ≈ 0.5 atol=1e-1
@test σ ≈ sqrt(1/12) atol=1e-2

# robust mask
mag = Float32.(readmag(fn_mag; rescale=true))
@test robustmask(mag) |> m -> count(.!m) / count(m) < 0.01
for i in 1:10
    mag[(end÷2):end,:,:,:] .= i .* 0.025 .* rand.()
    m = robustmask(mag)
    @test 1.1 < count(.!m) / count(m) < 1.2
end

# savenii
fn_temp = tempname()
savenii(mag, fn_temp)
mag2 = niread(fn_temp)
@test mag == mag2

# close mmapped files
GC.gc()

@test estimatequantile(1:1000, 0.8) ≈ 800 atol=1

function header_test(hdr, hdr2)
    @test hdr.scl_inter == 0
    @test hdr.scl_slope == 1
    @test hdr.dim == hdr2.dim
end
# similar
header_test(similar(mag_nii.header), mag_nii.header)
# header
header_test(header(mag_nii), mag_nii.header)
header_test(header(phase_nii), phase_nii.header)

# to_dim
@test [1 2] == to_dim([1, 2], 2)
a = 50:75
@test reshape(a, 1, 1, 1, :) == to_dim(a, 4)
@test reshape([5], 1, 1, 1, 1) == to_dim(5, 4)


text:
using Mypackage
using Test

#3*x*y
@test my_f(3,3)==27

#import Mypackage
#Mypackage.my_f(2,1)

@testset "Mypackage.jl" begin
    # Write your tests here.
    @test my_f(3,3)==27
    @test my_f(2,3)==18
    @test my_f(1,3)==9
end


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 01-2018
#--------------------------------------------------------------------

function testcoordtrans(N::Int, M::Int)
    xp = chebgrid(N)
    xp2xg1 = Float64[coordtransL2G(M, 1, chebx(i, N)) for i in 1:N+1]
    xp2xg2 = Float64[coordtransL2G(M, 2, chebx(i, N)) for i in 1:N+1]
    xg2xp  = Float64[coordtransG2L(M, 2, xg) for xg in xp2xg2]
    (xp2xg1[end] == xp2xg2[1]) && (xg2xp ≈ xp)
end

function testshapeH2L2H(a4N::Array{Float64,4}, Nx::Int, Ny::Int)::Bool
    shapeL2H(shapeH2L(a4N), Nx-1, Ny-1) == a4N
end

function testshapeH2L2H(b2N::Array{Float64,2}, Nx::Int, Ny::Int)::Bool
    shapeL2H(shapeH2L(b2N), Nx-1, Ny-1) == b2N
end

function testchebint(Nx::Int, Ny::Int)::Float64
    X  = chebgrid(Nx)
    Y  = chebgrid(Ny)
    B  = Float64[x^2 - y^3 + x^4*y^2 for x in X, y in Y]
    wx = chebweights(Nx)
    wy = chebweights(Ny)
    numint = wx'*B*wy
    return abs.(numint - 8/5)
end

function testL1norm(Nx::Int, Ny::Int)::Float64
    X = chebgrid(Nx)
    Y = chebgrid(Ny)
    A = Float64[x^2 - y^3 for x in X, y in Y]
    B = Float64[x^2 - y^3 + x^4*y^2 for x in X, y in Y]
    errorL1 = L1norm(B, A, chebweights(Nx), chebweights(Ny))
    exactL1 = 4/15
    return abs(errorL1 - exactL1)
end

function testL2norm(Nx::Int, Ny::Int)::Float64
    X = chebgrid(Nx)
    Y = chebgrid(Ny)
    A = Float64[x^2 - y^3 for x in X, y in Y]
    B = Float64[x^2 - y^3 + x^4*y^2 for x in X, y in Y]
    errorL2 = L2norm(B, A, chebweights(Nx), chebweights(Ny))
    exactL2 = 2/sqrt(45)
    return abs(errorL2 - exactL2)
end

function testarray2dict2array(Nx, Ny, M)::Bool
    sPatch = rand((Nx+1)*M, (Ny+1)*M)
    return dict2array(array2dict(sPatch, Nx, Ny, M)) == sPatch
end

@test coordtransL2G(2, 1, -1.0) ≈ 0.0
@test coordtransL2G(2, 2,  1.0)  ≈ 0.0
@test testcoordtrans(8,4) == true
@test jacobian(12) == 1/12
@test testshapeH2L2H(randn(4,12,4,12), 4, 12) == true
@test testshapeH2L2H(randn(6,10), 6, 10) == true
@test testchebint(12,14)  < 1e-14
@test testL1norm(14, 11)  < 1e-14
@test testL2norm(11, 14)  < 1e-14
@test testarray2dict2array(4,5,3) == true


text:
module TestPatterns
using Test

using ReTest: and, or, not, interpolated, reachable, depth, pass, fail, iter
import ReTest

struct MockTestset
    id
    marks
    parent
    iter

    MockTestset() = new(rand(1:typemax(Int)), ReTest.Marks(), nothing, 1)
end

ReTest.tsdepth(::MockTestset) = 1

const basic_patterns = [and(), or(), not(0), interpolated, 0, r"", :label,
                        depth(2), pass, fail, iter(1)]
VERSION >= v"1.3" && push!(basic_patterns, reachable(1))

@testset "patterns: ==" begin
    for a = basic_patterns, b = basic_patterns
        if a === b
            @test a == b
            if !(a isa Regex || a isa Integer || a isa Symbol)
                @test -a == not(a)
            end
        else
            @test a != b
            for f in (and, or)
                @test f(a) == f(a)
                @test f(a) == f(deepcopy(a))
                @test f(a) != f(b)
                @test f(a, b) == f(a, b)
                @test f(a, b) == f(deepcopy(a), deepcopy(b))
                @test f(a, b) != f(b, a)
                @test f(and(a, b), or(b, a)) == f(and(a, b), or(b, a))
                @test f(and(a, b), or(a, b)) != f(and(a, b), or(b, a))
                @test f(not(a), and(a, b)) == f(not(a), and(a, b))
                @test f(and(a, b), not(a)) != f(not(a), and(a, b))
                @test f(not(a), r"a", and(a, b)) == f(not(a), r"a", and(a, b))
                @test f(not(a), r"a", and(a, b)) != f(not(a), r"b", and(a, b))
            end
            @test not(a) == not(a)
            @test not(a) == not(deepcopy(a))
            @test not(a) != not(b)
            @test not(not(a)) == not(deepcopy(not(a)))
            if VERSION >= v"1.3"
                @test reachable(a) == reachable(a)
                @test reachable(a) == reachable(deepcopy(a))
                @test reachable(a) != reachable(b)
                @test reachable(reachable(a)) == reachable(deepcopy(reachable(a)))
            end
        end
    end
end

@testset "patterns: interface" begin
    for a in basic_patterns
        @test ReTest.matches(a, "a", MockTestset()) isa Union{Missing, Bool}
        @test ReTest.matches(a, missing, MockTestset()) isa Union{Missing, Bool}
        @test ReTest.alwaysmatches(a, 1) isa Bool
        @test ReTest.has(a, Integer) isa Bool
        @test ReTest.has(a, ReTest.Iter) isa Bool
    end
end

@testset "patterns: not" begin
    pats = [or(1, 3), and(1, r"a"), not(1), interpolated, depth(3)]
    VERSION >= v"1.3" && push!(pats, reachable("c"))
    for p ∈ pats
        @test -p == not(p)
    end
end

end # TestPatterns

max_issues_repo_name: rfourquet/InlineTest

max_forks_repo_name: rfourquet/InlineTest


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: mediastore
using AWS.Compat
using AWS.UUIDs

"""
    create_container(container_name)
    create_container(container_name, params::Dict{String,<:Any})

Creates a storage container to hold objects. A container is similar to a bucket in the
Amazon S3 service.

# Arguments
- `container_name`: The name for the container. The name must be from 1 to 255 characters.
  Container names must be unique to your AWS account within a specific region. As an example,
  you could create a container named movies in every region, as long as you don’t have an
  existing container with that name.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Tags"`: An array of key:value pairs that you define. These values can be anything that
  you want. Typically, the tag key represents a category (such as \"environment\") and the
  tag value represents a specific value within that category (such as \"test,\"
  \"development,\" or \"production\"). You can add up to 50 tags to each container. For more
  information about tagging, including naming and usage conventions, see Tagging Resources in
  MediaStore.
"""
function create_container(ContainerName; aws_config::AbstractAWSConfig=global_aws_config())
    return mediastore(
        "CreateContainer",
        Dict{String,Any}("ContainerName" => ContainerName);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function create_container(
    ContainerName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "CreateContainer",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ContainerName" => ContainerName), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_container(container_name)
    delete_container(container_name, params::Dict{String,<:Any})

Deletes the specified container. Before you make a DeleteContainer request, delete any
objects in the container or in any folders in the container. You can delete only empty
containers.

# Arguments
- `container_name`: The name of the container to delete.

"""
function delete_container(ContainerName; aws_config::AbstractAWSConfig=global_aws_config())
    return mediastore(
        "DeleteContainer",
        Dict{String,Any}("ContainerName" => ContainerName);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_container(
    ContainerName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "DeleteContainer",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ContainerName" => ContainerName), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_container_policy(container_name)
    delete_container_policy(container_name, params::Dict{String,<:Any})

Deletes the access policy that is associated with the specified container.

# Arguments
- `container_name`: The name of the container that holds the policy.

"""
function delete_container_policy(
    ContainerName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mediastore(
        "DeleteContainerPolicy",
        Dict{String,Any}("ContainerName" => ContainerName);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_container_policy(
    ContainerName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "DeleteContainerPolicy",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ContainerName" => ContainerName), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_cors_policy(container_name)
    delete_cors_policy(container_name, params::Dict{String,<:Any})

Deletes the cross-origin resource sharing (CORS) configuration information that is set for
the container. To use this operation, you must have permission to perform the
MediaStore:DeleteCorsPolicy action. The container owner has this permission by default and
can grant this permission to others.

# Arguments
- `container_name`: The name of the container to remove the policy from.

"""
function delete_cors_policy(
    ContainerName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mediastore(
        "DeleteCorsPolicy",
        Dict{String,Any}("ContainerName" => ContainerName);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_cors_policy(
    ContainerName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "DeleteCorsPolicy",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ContainerName" => ContainerName), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_lifecycle_policy(container_name)
    delete_lifecycle_policy(container_name, params::Dict{String,<:Any})

Removes an object lifecycle policy from a container. It takes up to 20 minutes for the
change to take effect.

# Arguments
- `container_name`: The name of the container that holds the object lifecycle policy.

"""
function delete_lifecycle_policy(
    ContainerName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mediastore(
        "DeleteLifecyclePolicy",
        Dict{String,Any}("ContainerName" => ContainerName);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_lifecycle_policy(
    ContainerName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "DeleteLifecyclePolicy",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ContainerName" => ContainerName), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    delete_metric_policy(container_name)
    delete_metric_policy(container_name, params::Dict{String,<:Any})

Deletes the metric policy that is associated with the specified container. If there is no
metric policy associated with the container, MediaStore doesn't send metrics to CloudWatch.

# Arguments
- `container_name`: The name of the container that is associated with the metric policy
  that you want to delete.

"""
function delete_metric_policy(
    ContainerName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mediastore(
        "DeleteMetricPolicy",
        Dict{String,Any}("ContainerName" => ContainerName);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function delete_metric_policy(
    ContainerName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "DeleteMetricPolicy",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ContainerName" => ContainerName), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    describe_container()
    describe_container(params::Dict{String,<:Any})

Retrieves the properties of the requested container. This request is commonly used to
retrieve the endpoint of a container. An endpoint is a value assigned by the service when a
new container is created. A container's endpoint does not change after it has been
assigned. The DescribeContainer request returns a single Container object based on
ContainerName. To return all Container objects that are associated with a specified AWS
account, use ListContainers.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ContainerName"`: The name of the container to query.
"""
function describe_container(; aws_config::AbstractAWSConfig=global_aws_config())
    return mediastore(
        "DescribeContainer"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function describe_container(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mediastore(
        "DescribeContainer", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    get_container_policy(container_name)
    get_container_policy(container_name, params::Dict{String,<:Any})

Retrieves the access policy for the specified container. For information about the data
that is included in an access policy, see the AWS Identity and Access Management User Guide.

# Arguments
- `container_name`: The name of the container.

"""
function get_container_policy(
    ContainerName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mediastore(
        "GetContainerPolicy",
        Dict{String,Any}("ContainerName" => ContainerName);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_container_policy(
    ContainerName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "GetContainerPolicy",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ContainerName" => ContainerName), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_cors_policy(container_name)
    get_cors_policy(container_name, params::Dict{String,<:Any})

Returns the cross-origin resource sharing (CORS) configuration information that is set for
the container. To use this operation, you must have permission to perform the
MediaStore:GetCorsPolicy action. By default, the container owner has this permission and
can grant it to others.

# Arguments
- `container_name`: The name of the container that the policy is assigned to.

"""
function get_cors_policy(ContainerName; aws_config::AbstractAWSConfig=global_aws_config())
    return mediastore(
        "GetCorsPolicy",
        Dict{String,Any}("ContainerName" => ContainerName);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_cors_policy(
    ContainerName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "GetCorsPolicy",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ContainerName" => ContainerName), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_lifecycle_policy(container_name)
    get_lifecycle_policy(container_name, params::Dict{String,<:Any})

Retrieves the object lifecycle policy that is assigned to a container.

# Arguments
- `container_name`: The name of the container that the object lifecycle policy is assigned
  to.

"""
function get_lifecycle_policy(
    ContainerName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mediastore(
        "GetLifecyclePolicy",
        Dict{String,Any}("ContainerName" => ContainerName);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_lifecycle_policy(
    ContainerName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "GetLifecyclePolicy",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ContainerName" => ContainerName), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    get_metric_policy(container_name)
    get_metric_policy(container_name, params::Dict{String,<:Any})

Returns the metric policy for the specified container.

# Arguments
- `container_name`: The name of the container that is associated with the metric policy.

"""
function get_metric_policy(ContainerName; aws_config::AbstractAWSConfig=global_aws_config())
    return mediastore(
        "GetMetricPolicy",
        Dict{String,Any}("ContainerName" => ContainerName);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function get_metric_policy(
    ContainerName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "GetMetricPolicy",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ContainerName" => ContainerName), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    list_containers()
    list_containers(params::Dict{String,<:Any})

Lists the properties of all containers in AWS Elemental MediaStore.  You can query to
receive all the containers in one response. Or you can include the MaxResults parameter to
receive a limited number of containers in each response. In this case, the response
includes a token. To get the next set of containers, send the command again, this time with
the NextToken parameter (with the returned token as its value). The next set of responses
appears, with a token if there are still more containers to receive.  See also
DescribeContainer, which gets the properties of one container.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"MaxResults"`: Enter the maximum number of containers in the response. Use from 1 to 255
  characters.
- `"NextToken"`: Only if you used MaxResults in the first command, enter the token (which
  was included in the previous response) to obtain the next set of containers. This token is
  included in a response only if there actually are more containers to list.
"""
function list_containers(; aws_config::AbstractAWSConfig=global_aws_config())
    return mediastore(
        "ListContainers"; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end
function list_containers(
    params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mediastore(
        "ListContainers", params; aws_config=aws_config, feature_set=SERVICE_FEATURE_SET
    )
end

"""
    list_tags_for_resource(resource)
    list_tags_for_resource(resource, params::Dict{String,<:Any})

Returns a list of the tags assigned to the specified container.

# Arguments
- `resource`: The Amazon Resource Name (ARN) for the container.

"""
function list_tags_for_resource(Resource; aws_config::AbstractAWSConfig=global_aws_config())
    return mediastore(
        "ListTagsForResource",
        Dict{String,Any}("Resource" => Resource);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function list_tags_for_resource(
    Resource,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "ListTagsForResource",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("Resource" => Resource), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    put_container_policy(container_name, policy)
    put_container_policy(container_name, policy, params::Dict{String,<:Any})

Creates an access policy for the specified container to restrict the users and clients that
can access it. For information about the data that is included in an access policy, see the
AWS Identity and Access Management User Guide. For this release of the REST API, you can
create only one policy for a container. If you enter PutContainerPolicy twice, the second
command modifies the existing policy.

# Arguments
- `container_name`: The name of the container.
- `policy`: The contents of the policy, which includes the following:    One Version tag
  One Statement tag that contains the standard tags for the policy.

"""
function put_container_policy(
    ContainerName, Policy; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mediastore(
        "PutContainerPolicy",
        Dict{String,Any}("ContainerName" => ContainerName, "Policy" => Policy);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function put_container_policy(
    ContainerName,
    Policy,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "PutContainerPolicy",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("ContainerName" => ContainerName, "Policy" => Policy),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    put_cors_policy(container_name, cors_policy)
    put_cors_policy(container_name, cors_policy, params::Dict{String,<:Any})

Sets the cross-origin resource sharing (CORS) configuration on a container so that the
container can service cross-origin requests. For example, you might want to enable a
request whose origin is http://www.example.com to access your AWS Elemental MediaStore
container at my.example.container.com by using the browser's XMLHttpRequest capability. To
enable CORS on a container, you attach a CORS policy to the container. In the CORS policy,
you configure rules that identify origins and the HTTP methods that can be executed on your
container. The policy can contain up to 398,000 characters. You can add up to 100 rules to
a CORS policy. If more than one rule applies, the service uses the first applicable rule
listed. To learn more about CORS, see Cross-Origin Resource Sharing (CORS) in AWS Elemental
MediaStore.

# Arguments
- `container_name`: The name of the container that you want to assign the CORS policy to.
- `cors_policy`: The CORS policy to apply to the container.

"""
function put_cors_policy(
    ContainerName, CorsPolicy; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mediastore(
        "PutCorsPolicy",
        Dict{String,Any}("ContainerName" => ContainerName, "CorsPolicy" => CorsPolicy);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function put_cors_policy(
    ContainerName,
    CorsPolicy,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "PutCorsPolicy",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ContainerName" => ContainerName, "CorsPolicy" => CorsPolicy
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    put_lifecycle_policy(container_name, lifecycle_policy)
    put_lifecycle_policy(container_name, lifecycle_policy, params::Dict{String,<:Any})

Writes an object lifecycle policy to a container. If the container already has an object
lifecycle policy, the service replaces the existing policy with the new policy. It takes up
to 20 minutes for the change to take effect. For information about how to construct an
object lifecycle policy, see Components of an Object Lifecycle Policy.

# Arguments
- `container_name`: The name of the container that you want to assign the object lifecycle
  policy to.
- `lifecycle_policy`: The object lifecycle policy to apply to the container.

"""
function put_lifecycle_policy(
    ContainerName, LifecyclePolicy; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mediastore(
        "PutLifecyclePolicy",
        Dict{String,Any}(
            "ContainerName" => ContainerName, "LifecyclePolicy" => LifecyclePolicy
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function put_lifecycle_policy(
    ContainerName,
    LifecyclePolicy,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "PutLifecyclePolicy",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ContainerName" => ContainerName, "LifecyclePolicy" => LifecyclePolicy
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    put_metric_policy(container_name, metric_policy)
    put_metric_policy(container_name, metric_policy, params::Dict{String,<:Any})

The metric policy that you want to add to the container. A metric policy allows AWS
Elemental MediaStore to send metrics to Amazon CloudWatch. It takes up to 20 minutes for
the new policy to take effect.

# Arguments
- `container_name`: The name of the container that you want to add the metric policy to.
- `metric_policy`: The metric policy that you want to associate with the container. In the
  policy, you must indicate whether you want MediaStore to send container-level metrics. You
  can also include up to five rules to define groups of objects that you want MediaStore to
  send object-level metrics for. If you include rules in the policy, construct each rule with
  both of the following:   An object group that defines which objects to include in the
  group. The definition can be a path or a file name, but it can't have more than 900
  characters. Valid characters are: a-z, A-Z, 0-9, _ (underscore), = (equal), : (colon), .
  (period), - (hyphen), ~ (tilde), / (forward slash), and * (asterisk). Wildcards (*) are
  acceptable.   An object group name that allows you to refer to the object group. The name
  can't have more than 30 characters. Valid characters are: a-z, A-Z, 0-9, and _
  (underscore).

"""
function put_metric_policy(
    ContainerName, MetricPolicy; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mediastore(
        "PutMetricPolicy",
        Dict{String,Any}("ContainerName" => ContainerName, "MetricPolicy" => MetricPolicy);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function put_metric_policy(
    ContainerName,
    MetricPolicy,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "PutMetricPolicy",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}(
                    "ContainerName" => ContainerName, "MetricPolicy" => MetricPolicy
                ),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    start_access_logging(container_name)
    start_access_logging(container_name, params::Dict{String,<:Any})

Starts access logging on the specified container. When you enable access logging on a
container, MediaStore delivers access logs for objects stored in that container to Amazon
CloudWatch Logs.

# Arguments
- `container_name`: The name of the container that you want to start access logging on.

"""
function start_access_logging(
    ContainerName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mediastore(
        "StartAccessLogging",
        Dict{String,Any}("ContainerName" => ContainerName);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function start_access_logging(
    ContainerName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "StartAccessLogging",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ContainerName" => ContainerName), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    stop_access_logging(container_name)
    stop_access_logging(container_name, params::Dict{String,<:Any})

Stops access logging on the specified container. When you stop access logging on a
container, MediaStore stops sending access logs to Amazon CloudWatch Logs. These access
logs are not saved and are not retrievable.

# Arguments
- `container_name`: The name of the container that you want to stop access logging on.

"""
function stop_access_logging(
    ContainerName; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mediastore(
        "StopAccessLogging",
        Dict{String,Any}("ContainerName" => ContainerName);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function stop_access_logging(
    ContainerName,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "StopAccessLogging",
        Dict{String,Any}(
            mergewith(_merge, Dict{String,Any}("ContainerName" => ContainerName), params)
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    tag_resource(resource, tags)
    tag_resource(resource, tags, params::Dict{String,<:Any})

Adds tags to the specified AWS Elemental MediaStore container. Tags are key:value pairs
that you can associate with AWS resources. For example, the tag key might be \"customer\"
and the tag value might be \"companyA.\" You can specify one or more tags to add to each
container. You can add up to 50 tags to each container. For more information about tagging,
including naming and usage conventions, see Tagging Resources in MediaStore.

# Arguments
- `resource`: The Amazon Resource Name (ARN) for the container.
- `tags`: An array of key:value pairs that you want to add to the container. You need to
  specify only the tags that you want to add or update. For example, suppose a container
  already has two tags (customer:CompanyA and priority:High). You want to change the priority
  tag and also add a third tag (type:Contract). For TagResource, you specify the following
  tags: priority:Medium, type:Contract. The result is that your container has three tags:
  customer:CompanyA, priority:Medium, and type:Contract.

"""
function tag_resource(Resource, Tags; aws_config::AbstractAWSConfig=global_aws_config())
    return mediastore(
        "TagResource",
        Dict{String,Any}("Resource" => Resource, "Tags" => Tags);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function tag_resource(
    Resource,
    Tags,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "TagResource",
        Dict{String,Any}(
            mergewith(
                _merge, Dict{String,Any}("Resource" => Resource, "Tags" => Tags), params
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end

"""
    untag_resource(resource, tag_keys)
    untag_resource(resource, tag_keys, params::Dict{String,<:Any})

Removes tags from the specified container. You can specify one or more tags to remove.

# Arguments
- `resource`: The Amazon Resource Name (ARN) for the container.
- `tag_keys`: A comma-separated list of keys for tags that you want to remove from the
  container. For example, if your container has two tags (customer:CompanyA and
  priority:High) and you want to remove one of the tags (priority:High), you specify the key
  for the tag that you want to remove (priority).

"""
function untag_resource(
    Resource, TagKeys; aws_config::AbstractAWSConfig=global_aws_config()
)
    return mediastore(
        "UntagResource",
        Dict{String,Any}("Resource" => Resource, "TagKeys" => TagKeys);
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end
function untag_resource(
    Resource,
    TagKeys,
    params::AbstractDict{String};
    aws_config::AbstractAWSConfig=global_aws_config(),
)
    return mediastore(
        "UntagResource",
        Dict{String,Any}(
            mergewith(
                _merge,
                Dict{String,Any}("Resource" => Resource, "TagKeys" => TagKeys),
                params,
            ),
        );
        aws_config=aws_config,
        feature_set=SERVICE_FEATURE_SET,
    )
end


text:
### Useful utilities

# @dbg: Uncomment this for debugging
# DBG=false; dbg()=DBG; dbg(b::Bool)=(global DBG=b); macro dbg(x) :(DBG && $(esc(x))) end
# @dbg: Uncomment this for production
macro dbg(x) nothing end

# gpusync: Uncomment this for GPU profiling
#gpusync()=device_synchronize()
# gpusync: Uncomment this for production
gpusync()=nothing

# @date: Print date, expression; run and print elapsed time after execution
macro date(_x) :(println("$(now()) "*$(string(_x)));flush(STDOUT);@time $(esc(_x))) end


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

export InnerProductLayer

function FullyConnectedEnsemble(name::Symbol, net::Net, num_inputs::Int, num_outputs::Int; weight_init=xavier, bias_init=0)
    # Create a 1-D array of `num_outputs` WeightedNeurons
    neurons = Array(WeightedNeuron, num_outputs)

    weights = weight_init(Float32, num_inputs, num_outputs)
    ∇weights = zeros(Float32, num_inputs, num_outputs)

    bias = Array(Float32, 1, num_outputs)
    fill!(bias, bias_init)
    ∇bias = zeros(Float32, 1, num_outputs)

    # Instantiate each neuron
    for i in 1:num_outputs
        # One weight for each input neuron
        neurons[i] = WeightedNeuron(view(weights, :, i), view(∇weights, :, i),
                                    view(bias, :, i), view(∇bias, :, i))
    end
    # Construct the ensemble
    Ensemble(net, name, neurons, [Param(name,:weights, 1.0f0, 1.0f0),
                                  Param(name,:bias, 2.0f0, 0.0f0)])
end

function FullyConnectedEnsemble(net::Net, num_inputs::Int, num_outputs::Int; weight_init=xavier, bias_init=0)
    FullyConnectedEnsemble(gensym("ensemble"), net, num_inputs, num_outputs; weight_init=weight_init, bias_init=bias_init)
end

function InnerProductLayer(name::Symbol, net::Net,
                           input_ensemble::AbstractEnsemble,
                           num_outputs::Int; weight_init=xavier, bias_init::Real=0)
    ip = FullyConnectedEnsemble(name, net, length(input_ensemble),
                                num_outputs; weight_init=weight_init,
                                bias_init=bias_init)

    # Connect each neuron in input_ensemble to each neuron in ip
    add_connections(net, input_ensemble, ip,
                    (i) -> (tuple([Colon() for d in size(input_ensemble)]... )))
    ip
end

function InnerProductLayer(net::Net,
                           input_ensemble::AbstractEnsemble,
                           num_outputs::Int; weight_init=xavier, bias_init::Real=0)
    InnerProductLayer(gensym("ensemble"), net, input_ensemble, num_outputs; weight_init=weight_init, bias_init=bias_init)
end


text:
struct CVCapacitanceReport{U1,U2,U3,U4,U5,U6,U7,U8,U9,U10}
    Δt::U1
    ΔV::U2
    area::U3
    scan_rate::U4
    quadrant::Int
    C::U5
    C_specific::U6
    E::U7
    E_specific::U8
    P::U9
    P_specific::U10
end

function capacitance(df, scan_rate, quadrant, fixed_ΔV)
    Δt, ΔV, ∫IdV = discharge_area(df, quadrant, fixed_ΔV)

    C = ∫IdV / (ΔV * scan_rate)
    return Δt, ΔV, ∫IdV, C
end

function CVCapacitanceReport(datafile, df, quadrant, folder, setup)
    ext = ".dat"    # using processed data
    scan_rate = parse(Float64, filevalue(datafile, ext)) * uparse(datafile.legend_units)
    porosity = parse(Int, foldervalue(datafile))
    @unpack a, A, fixed_ΔV = setup

    Δt, ΔV, area, C = capacitance(df, scan_rate, quadrant, fixed_ΔV)
    C_specific = specific_capacitance(C, porosity, folder, a, A)
    E = energy(C, ΔV)
    E_specific = energy(C_specific, ΔV)
    P = power(E, Δt)
    P_specific = power(E_specific, Δt)

    CVCapacitanceReport(Δt, ΔV, area, scan_rate, quadrant, C, C_specific, E, E_specific, P, P_specific)
end

function add_aux!(df, cr::CVCapacitanceReport)
    push!(df[!, :Area], ustrip(u"V*A", cr.area))
    push!(df[!, :quadrant], cr.quadrant)
    push!(df[!, :Δt], ustrip(u"s", cr.Δt))
    push!(df[!, :ΔV], ustrip(u"V", cr.ΔV))
    push!(df[!, :scan_rate], ustrip(u"mV/s", cr.scan_rate))

    return nothing
end

function cv_report_units()
    join(["V A", "", "s", "V", "mV/s", replace_unicode("μF"), "F/g", "W*h", "Wh/kg", "W", "W/kg"], ',')
end

function cv_result()
    DataFrame(
        :Area=>Float64[],
        :quadrant=>Int[],
        :Δt=>Float64[],
        :ΔV=>Float64[],
        :scan_rate=>Float64[],
        :C=>Float64[],
        :C_specific=>Float64[],
        :E=>Float64[],
        :E_specific=>Float64[],
        :P=>Float64[],
        :P_specific=>Float64[])
end


text:
include("euler/euler.jl")
using .Calculus: hexagonal_index, get_hexagonal

function compute(n::Int)::Int
    i = hexagonal_index(n) + 1
    while ((√(48i ^ 2 - 24i + 1) + 1) / 6) % 1 ≠ 0 i += 1 end
    return get_hexagonal(i)
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler