text:
using PowerModelsACDC
using PowerModels
using Ipopt
#using CPLEX
using SCS
using Mosek
using MosekTools
using JuMP
files =
[
"./test/data/case5_2grids.m";
"./test/data/case3120sp_acdc.m";
# "./test/data/case5_lcc.m";
"./test/data/case5_b2bdc.m";
# "./test/data/case5_dc.m"; #Don't use normally
"./test/data/case5_dcgrid.m";
"./test/data/case5_dcgrid_b0.m";
"./test/data/case24_3zones_acdc.m";
"./test/data/case39_acdc.m";
"./test/data/case3120sp_acdc.m";
"./test/data/pglib_opf_case588_sdet_acdc.m"
]

ipopt = JuMP.with_optimizer(Ipopt.Optimizer, tol=1e-6, print_level=0)
mosek = JuMP.with_optimizer(Mosek.Optimizer)
scs = JuMP.with_optimizer(SCS.Optimizer)
#mosek = MosekSolver()
s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true)

objective = Dict{String, Any}()

function exctract_info(dict)
    return Dict("obj" => dict["objective"], "solve_time" => dict["solve_time"], "result" => dict, "u_error" => 0, "pf_error" => 0, "pfdc_error" => 0, "u_error_rel" => 0, "pf_error_rel" => 0, "pfdc_error_rel" => 0  )
end

function calc_gap(dict)
    for (filename, formulations) in dict
        for (formulation, results) in formulations
            print(formulation)
            ac = formulations["AC NLP"]["obj"]
            results["gap"] = (ac - results["obj"])/ac
        end
    end
end

function calc_errors(dict)
    for (filename, formulations) in dict
        for (formulation, results) in formulations
            print(formulation)
            for (bus_id, bus) in formulations["AC NLP"]["result"]["solution"]["bus"]
                results["u_error"] = results["u_error"] + (bus["vm"] - results["result"]["solution"]["bus"][bus_id]["vm"])^2
            end
            results["u_error"] = results["u_error"] / length(formulations["AC NLP"]["result"]["solution"]["bus"])
            for (branch_id, branch) in formulations["AC NLP"]["result"]["solution"]["branch"]
                results["pf_error"] = results["pf_error"] + (branch["pf"] - results["result"]["solution"]["branch"][branch_id]["pf"])^2
            end
            results["pf_error"] = results["pf_error"] / length(formulations["AC NLP"]["result"]["solution"]["branch"])
            for (branch_id, branch) in formulations["AC NLP"]["result"]["solution"]["branchdc"]
                results["pfdc_error"] = results["pfdc_error"] + (branch["pf"] - results["result"]["solution"]["branchdc"][branch_id]["pf"])^2
            end
            results["pfdc_error"] = results["pfdc_error"] / length(formulations["AC NLP"]["result"]["solution"]["branchdc"])
        end
    end
end

function calc_errors_rel(dict)
    for (filename, formulations) in dict
        for (formulation, results) in formulations
            print(formulation)
            for (bus_id, bus) in formulations["AC NLP"]["result"]["solution"]["bus"]
                results["u_error_rel"] = results["u_error_rel"] + abs(bus["vm"] - results["result"]["solution"]["bus"][bus_id]["vm"])/abs(bus["vm"]+1e-5)
            end
            results["u_error_rel"] = results["u_error_rel"] / length(formulations["AC NLP"]["result"]["solution"]["bus"])
            for (branch_id, branch) in formulations["AC NLP"]["result"]["solution"]["branch"]
                results["pf_error_rel"] = results["pf_error_rel"] + abs(branch["pf"] - results["result"]["solution"]["branch"][branch_id]["pf"])/abs(branch["pf"]+1e-5)
            end
            results["pf_error_rel"] = results["pf_error_rel"] / length(formulations["AC NLP"]["result"]["solution"]["branch"])
            for (branch_id, branch) in formulations["AC NLP"]["result"]["solution"]["branchdc"]
                results["pfdc_error_rel"] = results["pfdc_error_rel"] + abs(branch["pf"] - results["result"]["solution"]["branchdc"][branch_id]["pf"])/abs(branch["pf"]+1e-5)
            end
            results["pfdc_error_rel"] = results["pfdc_error_rel"] / length(formulations["AC NLP"]["result"]["solution"]["branchdc"])
        end
    end
end

#
function fix_things!(data)
    # fix number of poles
    # data["dcpol"] = 1
    #tap = 1
    for (i,conv) in data["convdc"]
        # remove transformer and phase reactor from relaxation
        # conv["reactor"] = 0
        #conv["transformer"] = 0
        #conv["filter"] = 0
        #conv["tm"] = 1
        #conv["LossA"] = 0
        #conv["LossB"] = 0
        #conv["LossC"] = 0
    end
end

for file in files
    data = _PM.parse_file(file)
    PowerModelsACDC.process_additional_data!(data)
    #fix_things!(data)

    case = Dict{String, Any}()
    er_dict = Dict{String, Any}()

    resultAC = run_acdcopf(data, ACPPowerModel, ipopt; setting = s)
    case["AC NLP"] = exctract_info(resultAC)

    resultQC = run_acdcopf(data, QCWRPowerModel, ipopt; setting = s)
    case["QC SOC"] = exctract_info(resultQC)
    # #
    # resultQCTri = run_acdcopf(data, QCWRTriPowerModel, ipopt; setting = s)
    # case["QCTri SOC"] = exctract_info(resultQCTri)
    # #
    resultSOCBIM = run_acdcopf(data, SOCWRPowerModel, ipopt; setting = s)
    case["BIM SOC"] = exctract_info(resultSOCBIM)
    # #
    resultSOCBFM = run_acdcopf_bf(data, SOCBFPowerModel, ipopt; setting = s)
    case["BFM SOC"] = exctract_info(resultSOCBFM)
    # #
    case["BIM SDP"] = exctract_info(resultAC)
    # resultSDP = run_acdcopf(data, SDPWRMPowerModel, mosek; setting = s)
    # case["BIM SDP"] = exctract_info(resultSDP)
    # # #
    resultDC = run_acdcopf(data, DCPPowerModel, ipopt; setting = s)
    case["DC LP"] = exctract_info(resultDC)
    #
    resultLPAC = run_acdcopf(data, LPACCPowerModel, ipopt; setting = s)
    case["LPAC"] = exctract_info(resultLPAC)

    objective[file] = case
end

calc_gap(objective)
calc_errors(objective)
calc_errors_rel(objective)

function print_table_opt_gap(dict)
    s = ""
    c = " & "
    s = s*"case"*c*"AC NLP"*c* "BIM SDP" *c*c*  "QC SOC" *c*c* "BIM SOC" *c*c* "BFM SOC" *c*c*"DC LP" *c*c*"LPAC"*c*raw"\ "[1]*raw"\ "* " \n"
    for (filename, ff) in dict
        s = s*filename
        l = 2
        l2 = 8
        p = 100

        s = s *c*string(round(ff["AC NLP"]["obj"], digits = l))
        s = s *c*string(round(ff["BIM SDP"]["obj"], digits = l))
        s = s *c*string(round(p*ff["BIM SDP"]["gap"], digits = l))
        s = s *c*string(round(ff["QC SOC"]["obj"], digits = l))
        s = s *c*string(round(p*ff["QC SOC"]["gap"], digits = l))
        # s = s *c*string(roundff["QCTri SOC"]["obj"])[1:l2]
        # s = s *c*string(roundp*ff["QCTri SOC"]["gap"])[1:l]
        s = s *c*string(round(ff["BIM SOC"]["obj"], digits = l))
        s = s *c*string(round(p*ff["BIM SOC"]["gap"], digits = l))
        s = s *c*string(round(ff["BFM SOC"]["obj"], digits = l))
        s = s *c*string(round(p*ff["BFM SOC"]["gap"], digits = l))
        s = s *c*string(round(ff["DC LP"]["obj"], digits = l))
        s = s *c*string(round(p*ff["DC LP"]["gap"], digits = l))
        s = s *c*string(round(ff["LPAC"]["obj"], digits = l))
        s = s *c*string(round(p*ff["LPAC"]["gap"], digits = l))
        s = s*raw"\ "[1]*raw"\ "*" \n"
    end
    return s
end
clearconsole()
stt = print_table_opt_gap(objective)
print(stt)
function print_table_errors(dict)
    s = ""
    c = " & "
    s = s*"case"*c*"BIM SDP"*c*c*c* "QC SOC"*c*c*c* "BIM SOC"*c*c*c* "BFM SOC"*c*c*c*raw"\ "[1]*raw"\ "* " \n"
    for (filename, ff) in dict
        s = s*filename
        l = 6
        l2 = 8
        p = 100
        s = s *c*string(round(ff["BIM SDP"]["u_error"], digits = l))
        s = s *c*string(round(p*ff["BIM SDP"]["pf_error"], digits = l))
        s = s *c*string(round(p*ff["BIM SDP"]["pfdc_error"], digits = l))
        s = s *c*string(round(ff["QC SOC"]["u_error"], digits = l))
        s = s *c*string(round(p*ff["QC SOC"]["pf_error"], digits = l))
        s = s *c*string(round(p*ff["QC SOC"]["pfdc_error"], digits = l))
        s = s *c*string(round(ff["BIM SOC"]["u_error"], digits = l))
        s = s *c*string(round(p*ff["BIM SOC"]["pf_error"], digits = l))
        s = s *c*string(round(p*ff["BIM SOC"]["pfdc_error"], digits = l))
        s = s *c*string(round(ff["BFM SOC"]["u_error"], digits = l))
        s = s *c*string(round(p*ff["BFM SOC"]["pf_error"], digits = l))
        s = s *c*string(round(p*ff["BFM SOC"]["pfdc_error"], digits = l))

        s = s*raw"\ "[1]*raw"\ "*" \n"
    end
    return s
end
stt = print_table_errors(objective)
print(stt)

function print_table_errors_rel(dict)
    s = ""
    c = " & "
    s = s*"case"*c*"BIM SDP"*c*c*c* "QC SOC"*c*c*c* "BIM SOC"*c*c*c* "BFM SOC"*c*c*c*raw"\ "[1]*raw"\ "* " \n"
    for (filename, ff) in dict
        s = s*filename
        l = 6
        l2 = 8
        p = 100
        s = s *c*string(round(p*ff["BIM SDP"]["u_error_rel"], digits = l))
        s = s *c*string(round(p*ff["BIM SDP"]["pf_error_rel"], digits = l))
        s = s *c*string(round(p*ff["BIM SDP"]["pfdc_error_rel"], digits = l))
        s = s *c*string(round(p*ff["QC SOC"]["u_error_rel"], digits = l))
        s = s *c*string(round(p*ff["QC SOC"]["pf_error_rel"], digits = l))
        s = s *c*string(round(p*ff["QC SOC"]["pfdc_error_rel"], digits = l))
        s = s *c*string(round(p*ff["BIM SOC"]["u_error_rel"], digits = l))
        s = s *c*string(round(p*ff["BIM SOC"]["pf_error_rel"], digits = l))
        s = s *c*string(round(p*ff["BIM SOC"]["pfdc_error_rel"], digits = l))
        s = s *c*string(round(p*ff["BFM SOC"]["u_error_rel"], digits = l))
        s = s *c*string(round(p*ff["BFM SOC"]["pf_error_rel"], digits = l))
        s = s *c*string(round(p*ff["BFM SOC"]["pfdc_error_rel"], digits = l))
        # s = s *c*string(round(p*ff["DC LP"]["u_error_rel"], digits = l))
        # s = s *c*string(round(p*ff["DC LP"]["pf_error_rel"], digits = l))
        # s = s *c*string(round(p*ff["DC LP"]["pfdc_error_rel"], digits = l))

        s = s*raw"\ "[1]*raw"\ "*" \n"
    end
    return s
end
stt = print_table_errors_rel(objective)
print(stt)
function print_table_comp_time(dict)
    s = ""
    c = " & "
    s = s*"case"*c*"BIM NLP"*c*"BIM SDP"*c* "QC SOC"*c*c*c* "BIM SOC"*c*c*c* "BFM SOC"*c*c*c*raw"\ "[1]*raw"\ "* " \n"
    for (filename, ff) in dict
        s = s*filename
        l = 6
        l2 = 8
        p = 100
        s = s *c*string(round(ff["AC NLP"]["solve_time"], digits = l))
        s = s *c*string(round(ff["BIM SDP"]["solve_time"], digits = l))
        s = s *c*string(round(ff["QC SOC"]["solve_time"], digits = l))
        s = s *c*string(round(ff["BIM SOC"]["solve_time"], digits = l))
        s = s *c*string(round(ff["BFM SOC"]["solve_time"], digits = l))
        s = s *c*string(round(ff["DC LP"]["solve_time"], digits = l))

        s = s*raw"\ "[1]*raw"\ "*" \n"
    end
    return s
end
stt = print_table_comp_time(objective)
print(stt)


text:
using TerminalUserInterfaces
const TUI = TerminalUserInterfaces
using REPL
using InteractiveUtils

function main(max_count)
    TUI.initialize()

    y, x = 1, 1

    counter = 1
    t = TUI.Terminal()

    ADDITIVE = 1

    for _ in 1:max_count

        w, _ = TUI.terminal_size()

        r = TUI.Rect(x, y + 1, w, 2)

        b = TUI.Block(border = TUI.BorderNone, border_type = TUI.BorderTypeHeavy)
        pg = TUI.ProgressBar(b, counter / max_count)

        TUI.draw(t, pg, r)

        TUI.flush(t)

        counter += 1

    end

    println()

    TUI.cleanup()

end

main(1000)

ext: jl


text:
using Revise
using GeodesicTracer
using CarterBoyerLindquist
using StaticArrays
using Plots
using Random
using AccretionFormulae
using ComputedGeodesicEquations
using Dates
gr()

# setup
init_radii = 5:5
min_α = 3.0
max_α = 10.0
stp_α = 1
min_β = 0.0
max_β = 0.0
stp_β = 0.5

# values of α and β to search
α_vals = min_α:stp_α:max_α
β_vals = min_β:stp_β:max_β

total = length(α_vals)*length(β_vals)*length(init_radii) # number of combinations to test
found = 0   # number of solutions found so far
i = 1       # number of combinations tested
start = now() # timing start
for init_radius in init_radii
    for α in α_vals
        for β in β_vals
            println("Testing: $i / $total")

            # path details for the given parameters
            m = BoyerLindquist(M=1.0, a=1.0) #chosen metric
            R_isco = AccretionFormulae.r_isco(m.a, m.M)
            u = @SVector [0.0, init_radius, deg2rad(90.0), 0.0] #initial position (time, radius, theta, phi)
            v = map_impact_parameters(m, u, α, β) #impact parameters α, β to velocity vector at some position u in metric m
            sols = tracegeodesics(m, u, v, (0.0, 10000.0), μ=1.0;
                                    abstol=1e-14,
                                    reltol=1e-14,
                                    dtmax=1.0
                                    )

            # selecting only radius values, and only looking at the end
            new_sols = copy(selectdim(sols, 1, 6))
            new_sols = new_sols[Int32(round(length(sols)*0.75)):end]
            # finding the difference between its extremities
            min, max =  extrema(new_sols)
            diff = max - min

            # if the difference is small enough, and it doesn't crash into the ISCO
            # it is a valid solution, and is displayed
            if (diff <= 5) & (min > R_isco*1.1)
                found += 1
                plt = plot(sols,
                            vars=(4+4,2+4),
                            proj=:polar,
                            range=(0.0, 30.0),
                            legend=false,
                            titlefontsize=8
                            ) #plotting 4 (phi) and 2 (r), in specified radius range
                title!("Initial Radius=$init_radius, α=$α, β=$β; min=$min, max=$max")
                display(plt)
            end
            i+=1
        end
    end
end
duration = now() - start # timing finish
println("Done")
println("Found $found solutions in $duration.")
println("Success rate of $(found/total) .")


text:
# Depricated, tilsetMeta is the only one we use now

ATLAS_PATHS = Dict{String, String}(
    "tiles" => "Atlases/Gameplay/tilesets",
    "decals" => "Atlases/Gameplay/decals",
    "tilesetMeta" => "Content/Graphics"
)


text:
# Outer (tensor) product
function outer{N,a}(A::Array{N,a}, B::Vector{N})
    reshape(N[x*y for x in A, y in B], size(A)..., size(B)...)::Array{N,a+1}
end

# Covariance matrix with adjustable bias (Bessel's correction)
function _cov{T<:Real}(data::Matrix{T}; bias::Int=0, dense::Bool=true)
    num_samples, num_signals = size(data)
    @inbounds cntr = data .- mean(data, 1)

    # Dense matrix: all values are calculated, including duplicates
    if dense
        tensor = zeros(num_signals, num_signals)
        @inbounds for i = 1:num_samples
            c = cntr[i,:]
            tensor += c' * c
        end

    # Lower triangular matrix (duplicate values not calculated)
    # Convert to dense matrix using:
    #   tensor + tensor' - diagm(diag(tensor)))
    else
        tensor = zeros(num_signals, num_signals)
        @simd for i = 1:num_signals
            @simd for j = 1:i
                @inbounds tensor[i,j] = sum(cntr[:,i].*cntr[:,j])
            end
        end
    end
    tensor / (num_samples - bias)
end

# Coskewness tensor (third-order)
function coskew{T<:Real}(data::Matrix{T};
                         standardize::Bool=false,
                         flatten::Bool=false,
                         bias::Int=0,
                         dense::Bool=true)

    # Center/whiten data
    cntr, num_samples, num_signals = center(data;
                                            standardize=standardize,
                                            bias=bias)

    # Flattened representation (i.e., unfolded to a matrix)
    if flatten

        # Dense: all values are calculated, including duplicates
        if dense
            tensor = zeros(num_signals, num_signals^2)
            @inbounds for i = 1:num_samples
                c = cntr[i,:]
                tensor += kron(c'*c, c)
            end

        # Triangular: duplicate values ignored
        else
            tensor = spzeros(num_signals, num_signals^2)
            @simd for i = 1:num_signals
                @simd for j = 1:i
                    @simd for k = 1:j
                        @inbounds tensor[j,(i-1)*num_signals+k] = sum(cntr[:,i].*cntr[:,j].*cntr[:,k])
                    end
                end
            end
        end

    # Cube: third-order tensor representation
    else
        tensor = zeros(num_signals, num_signals, num_signals)

        # Dense: all values are calculated, including duplicates
        if dense
            @inbounds for i = 1:num_samples
                c = cntr[i,:][:]
                tensor += outer(c*c', c)
            end

        # Triangular: duplicate values ignored
        else
            @simd for i = 1:num_signals
                @simd for j = 1:i
                    @simd for k = 1:j
                        @inbounds tensor[i,j,k] = sum(cntr[:,i].*cntr[:,j].*cntr[:,k])
                    end
                end
            end
        end
    end
    tensor / (num_samples - bias)
end

# Cokurtosis tensor (fourth-order)
function cokurt{T<:Real}(data::Matrix{T};
                         standardize::Bool=false,
                         flatten::Bool=false,
                         bias::Int=0,
                         dense::Bool=true)

    # Center/whiten data
    cntr, num_samples, num_signals = center(data;
                                            standardize=standardize,
                                            bias=bias)

    # Flattened representation (i.e., unfolded into a matrix)
    if flatten

        # Dense: all values are calculated, including duplicates
        if dense
            tensor = zeros(num_signals, num_signals^3)
            @inbounds for i = 1:num_samples
                c = cntr[i,:]
                tensor += kron(c'*c, c, c)
            end

        # Triangular: duplicate values ignored
        else
            tensor = spzeros(num_signals, num_signals^3)
            @simd for i = 1:num_signals
                @simd for j = 1:i
                    @simd for k = 1:j
                        @simd for l = 1:k
                            @inbounds tensor[k,(i-1)*num_signals^2 + (j-1)*num_signals + l] = sum(cntr[:,i].*cntr[:,j].*cntr[:,k].*cntr[:,l])
                        end
                    end
                end
            end
        end

    # Tesseract: fourth-order tensor representation
    else
        tensor = zeros(num_signals, num_signals, num_signals, num_signals)

        # Dense: all values are calculated, including duplicates
        if dense
            @inbounds for i = 1:num_samples
                c = cntr[i,:][:]
                tensor += outer(outer(c*c', c), c)
            end

        # Triangular: duplicate values ignored
        else
            @simd for i = 1:num_signals
                @simd for j = 1:i
                    @simd for k = 1:j
                        @simd for l = 1:k
                            @inbounds tensor[i,j,k,l] = sum(cntr[:,i].*cntr[:,j].*cntr[:,k].*cntr[:,l])
                        end
                    end
                end
            end
        end
    end
    tensor / (num_samples - bias)
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# 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 the copyright holder 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 HOLDER 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.
# *********************************************************************************
"""
    add_chp_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")

Adds the `CHP` results to the dictionary passed back from `run_reopt` using the solved model `m` and the `REoptInputs` for node `_n`.
Note: the node number is an empty string if evaluating a single `Site`.

CHP results:
- `size_kw` Power capacity size of the CHP system [kW]
- `size_supplemental_firing_kw` Power capacity of CHP supplementary firing system [kW]
- `year_one_fuel_used_mmbtu` Fuel consumed in year one [MMBtu]
- `year_one_electric_energy_produced_kwh` Electric energy produced in year one [kWh]
- `year_one_thermal_energy_produced_mmbtu` Thermal energy produced in year one [MMBtu]
- `year_one_electric_production_series_kw` Electric power production time-series array [kW]
- `year_one_to_grid_series_kw` Electric power exported time-series array [kW]
- `year_one_to_battery_series_kw` Electric power to charge the battery storage time-series array [kW]
- `year_one_to_load_series_kw` Electric power to serve the electric load time-series array [kW]
- `year_one_thermal_to_waste_series_mmbtu_per_hour` Thermal power wasted/unused/vented time-series array [MMBtu/hr]
- `year_one_thermal_to_load_series_mmbtu_per_hour` Thermal power to serve the heating load time-series array [MMBtu/hr]
- `year_one_chp_fuel_cost` Fuel cost from fuel consumed by the CHP system [\$]
- `lifecycle_chp_fuel_cost` Fuel cost from fuel consumed by the CHP system [\$]
"""
function add_chp_results(m::JuMP.AbstractModel, p::REoptInputs, d::Dict; _n="")
    r = Dict{String, Any}()
	r["size_kw"] = value(sum(m[Symbol("dvSize"*_n)][t] for t in p.techs.chp))
    r["size_supplemental_firing_kw"] = value(sum(m[Symbol("dvSupplementaryFiringSize"*_n)][t] for t in p.techs.chp))
	@expression(m, CHPFuelUsedKWH, sum(m[Symbol("dvFuelUsage"*_n)][t, ts] for t in p.techs.chp, ts in p.time_steps))
	r["year_one_fuel_used_mmbtu"] = round(value(CHPFuelUsedKWH) / MMBTU_TO_KWH, digits=3)
	@expression(m, Year1CHPElecProd,
		p.hours_per_timestep * sum(m[Symbol("dvRatedProduction"*_n)][t,ts] * p.production_factor[t, ts]
			for t in p.techs.chp, ts in p.time_steps))
	r["year_one_electric_energy_produced_kwh"] = round(value(Year1CHPElecProd), digits=3)
	@expression(m, Year1CHPThermalProdKWH,
		p.hours_per_timestep * sum(m[Symbol("dvThermalProduction"*_n)][t,ts] +
        m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] -
        m[Symbol("dvProductionToWaste"*_n)][t,ts]
            for t in p.techs.chp, ts in p.time_steps))
	r["year_one_thermal_energy_produced_mmbtu"] = round(value(Year1CHPThermalProdKWH) / MMBTU_TO_KWH, digits=3)
	@expression(m, CHPElecProdTotal[ts in p.time_steps],
		sum(m[Symbol("dvRatedProduction"*_n)][t,ts] * p.production_factor[t, ts] for t in p.techs.chp))
	r["year_one_electric_production_series_kw"] = round.(value.(CHPElecProdTotal), digits=3)
	@expression(m, CHPtoGrid[ts in p.time_steps], sum(m[Symbol("dvProductionToGrid"*_n)][t,u,ts]
			for t in p.techs.chp, u in p.export_bins_by_tech[t]))
	r["year_one_to_grid_series_kw"] = round.(value.(CHPtoGrid), digits=3)
	@expression(m, CHPtoBatt[ts in p.time_steps],
		sum(m[Symbol("dvProductionToStorage"*_n)][:elec,t,ts] for t in p.techs.chp))
	r["year_one_to_battery_series_kw"] = round.(value.(CHPtoBatt), digits=3)
	@expression(m, CHPtoLoad[ts in p.time_steps],
		sum(m[Symbol("dvRatedProduction"*_n)][t, ts] * p.production_factor[t, ts] * p.levelization_factor[t]
			for t in p.techs.chp) - CHPtoBatt[ts] - CHPtoGrid[ts])
	r["year_one_to_load_series_kw"] = round.(value.(CHPtoLoad), digits=3)
	# @expression(m, CHPtoHotTES[ts in p.time_steps],
	# 	sum(m[Symbol("dvProductionToStorage"*_n)]["HotTES",t,ts] for t in p.techs.chp))
	# r["year_one_thermal_to_tes_series_mmbtu_per_hour"] = round.(value.(CHPtoHotTES), digits=5)
    # if !isempty(p.SteamTurbineTechs)
    #     @expression(m, CHPToSteamTurbine[ts in p.time_steps], sum(m[Symbol("dvThermalToSteamTurbine"*_n)][t,ts] for t in p.techs.chp))
    #     r["year_one_thermal_to_steamturbine_series_mmbtu_per_hour"] = round.(value.(CHPToSteamTurbine), digits=3)
    # else
    #     CHPToSteamTurbine = zeros(p.TimeStepCount)
    #     r["year_one_thermal_to_steamturbine_series_mmbtu_per_hour"] = round.(CHPToSteamTurbine, digits=3)
    # end
	@expression(m, CHPThermalToWasteKWH[ts in p.time_steps],
		sum(m[Symbol("dvProductionToWaste"*_n)][t,ts] for t in p.techs.chp))
	r["year_one_thermal_to_waste_series_mmbtu_per_hour"] = round.(value.(CHPThermalToWasteKWH) / MMBTU_TO_KWH, digits=5)
	# @expression(m, CHPThermalToLoad[ts in p.time_steps],
	# 	sum(m[Symbol("dvThermalProduction"*_n)][t,ts] + m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts]
	# 		for t in p.techs.chp) - CHPtoHotTES[ts] - CHPToSteamTurbine[ts] - CHPThermalToWaste[ts])
    @expression(m, CHPThermalToLoadKWH[ts in p.time_steps],
        sum(m[Symbol("dvThermalProduction"*_n)][t,ts] for t in p.techs.chp) -
        CHPThermalToWasteKWH[ts])
	r["year_one_thermal_to_load_series_mmbtu_per_hour"] = round.(value.(CHPThermalToLoadKWH) / MMBTU_TO_KWH, digits=5)
    r["year_one_chp_fuel_cost"] = round(value(m[:TotalCHPFuelCosts] / p.pwf_fuel["CHP"]), digits=3)
	r["lifecycle_chp_fuel_cost"] = round(value(m[:TotalCHPFuelCosts]) * p.s.financial.offtaker_tax_pct, digits=3)
	#Standby charges and hourly O&M
	r["lifecycle_chp_standby_cost_us_dollars"] = round(value(m[Symbol("TotalCHPStandbyCharges")]) * p.s.financial.offtaker_tax_pct, digits=0)
	r["year_one_chp_standby_cost_us_dollars"] = round(value(m[Symbol("TotalCHPStandbyCharges")]) / p.pwf_e, digits=0)

    d["CHP"] = r
    nothing
end


text:
"""
Computation unit that uses shaders as part of a graphics or compute pipeline.

It exposes a program interface through its shader interfaces and its shader resources.
"""
@auto_hash_equals struct Program
    shaders::Dictionary{Vk.ShaderStageFlag, Shader}
end

function Program(cache::ShaderCache, shaders::ShaderSpecification...)
    shaders = map(shaders) do shader_spec
        shader_spec.stage => find_shader!(cache, shader_spec)
    end
    Program(dictionary(shaders))
end

function Program(device, shaders::ShaderSpecification...)
    Program(device.shader_cache, shaders...)
end

"""
Set of data required to call a program.
"""
struct ProgramInterface
    vbuffer::Buffer
    ibuffer::Optional{Buffer}
    descriptors::Vector{DescriptorSet}
    push_constants::Dictionary{Vk.PushConstantRange,Any}
end

"""
Program to be compiled into a pipeline with a specific state.
"""
@auto_hash_equals struct ProgramInstance
    program::Program
    state::DrawState
    targets::TargetAttachments
end


text:
# minimal working example for fitting a sem in parallel

using LinearAlgebra, Optim, BenchmarkTools

## we need to minimize a loss function that is quite similiar to the one below
# the loss function will hence be called a few hundreds to a few thousend times by an optimizer
# we will often find us in situation where we have several thousend (even millions) of these optimization problems
# but our focus for now is on parallelization of the a single loss function (which is a sum of independend matrix operations)

# number of parallelizable units
# in real life this value will usually be around 1-2000 and almost never > 10000
n = 100

parameters = [100.0, 300] # we get those from the optimizer;
                        # it decides which parameters are tested
specific_data = rand(n) # something that is unique between units

# pre-allocating some memory
# in real life those matrices will usually have dimensions around 10*10 to 50*50 and almost never > 200*200
pre = [rand(2,2) for i in 1:n]
for i in 1:n pre[i] = pre[i]*pre[i]' end
prod = copy(pre)
a = cholesky.(Hermitian.(pre))

F = zeros(n) # results of the loss function

function loss_perunit(pars, pre, a, specific_data)
    pre[1] = pars[2]*specific_data
    pre[2] = 2*pars[1]*specific_data
    pre[3] = 2*pars[1]*specific_data
    pre[4] = pars[2]*specific_data
    a = cholesky!(Hermitian(pre); check = false)
    if !isposdef(a) return Inf end
    ld = logdet(a)
    pre = LinearAlgebra.inv!(a)
    return ld + tr(pre)
end

### not parallelized
function loss(pars, pre, a, specific_data, n, F)
    for i in 1:n
        F[i] = loss_perunit(pars, pre[i], a[i], specific_data[i])
    end
    return sum(F)
end

@benchmark loss(parameters, pre, a, specific_data, n, F)

### parallelized
function loss_parallel(pars, pre, a, specific_data, n, F)
    Threads.@threads for i in 1:n
        F[i] = loss_perunit(pars, pre[i], a[i], specific_data[i])
    end
    return sum(F)
end

@benchmark loss_parallel(parameters, pre, a, specific_data, n, F)

# We have two questions:
# 1. The preallocated memory can be very big (e.g. 8 bytes * 50*50* 1000 * 2 = 40GB for example).
# However, in theory we believe only twice as many matrices have to be preallocated as there are threads,
# but we are not shure how to implement that.
# 2. As soon as one computation yields the value Inf; all other computations can be stopped; because the sum will be Inf in any case.


text:
# Warp Shuffle (B.14)

# TODO: does not work on sub-word (ie. Int16) or non-word divisible sized types

# TODO: should shfl_idx conform to 1-based indexing?

# TODO: these functions should dispatch based on the actual warp size
const ws = Int32(32)

# TODO: this functionality should throw <sm_30

# primitive intrinsics

# "two packed values specifying a mask for logically splitting warps into sub-segments
# and an upper bound for clamping the source lane index"
@inline pack(width::UInt32, mask::UInt32)::UInt32 = (convert(UInt32, ws - width) << 8) | mask

# NOTE: CUDA C disagrees with PTX on how shuffles are called
for (name, mode, mask) in (("_up",   :up,   UInt32(0x00)),
                           ("_down", :down, UInt32(0x1f)),
                           ("_xor",  :bfly, UInt32(0x1f)),
                           ("",      :idx,  UInt32(0x1f)))
    fname = Symbol("shfl$name")

    if cuda_driver_version >= v"9.0" && v"6.0" in ptx_support
        instruction = Symbol("shfl.sync.$mode.b32")
        fname_sync = Symbol("$(fname)_sync")

        # TODO: implement using LLVM intrinsics when we have D38090

        @eval begin
            export $fname_sync, $fname

            @inline $fname_sync(val::UInt32, src::UInt32, width::UInt32=$ws,
                                threadmask::UInt32=0xffffffff) =
                @asmcall($"$instruction \$0, \$1, \$2, \$3, \$4;", "=r,r,r,r,r", true,
                         UInt32, NTuple{4,UInt32},
                         val, src, pack(width, $mask), threadmask)

            # FIXME: replace this with a checked conversion once we have exceptions
            @inline $fname_sync(val::UInt32, src::Integer, width::Integer=$ws,
                                threadmask::UInt32=0xffffffff) =
                $fname_sync(val, unsafe_trunc(UInt32, src), unsafe_trunc(UInt32, width),
                            threadmask)

            @inline $fname(val::UInt32, src::Integer, width::Integer=$ws) =
                $fname_sync(val, src, width)
        end
    else
        intrinsic = Symbol("llvm.nvvm.shfl.$mode.i32")

        @eval begin
            export $fname
            @inline $fname(val::UInt32, src::UInt32, width::UInt32=$ws) =
                ccall($"$intrinsic", llvmcall, UInt32,
                      (UInt32, UInt32, UInt32),
                      val, src, pack(width, $mask))

            # FIXME: replace this with a checked conversion once we have exceptions
            @inline $fname(val::UInt32, src::Integer, width::Integer=$ws) =
                $fname(val, unsafe_trunc(UInt32, src), unsafe_trunc(UInt32, width))
        end
    end
end

# wide and aggregate intrinsics

for name in ["_up", "_down", "_xor", ""]
    fname = Symbol("shfl$name")
    @eval @inline $fname(src, args...) = recurse_value_invocation($fname, src, args...)

    fname_sync = Symbol("$(fname)_sync")
    @eval @inline $fname_sync(src, args...) = recurse_value_invocation($fname, src, args...)
end

# documentation

@doc """
    shfl(val, lane::Integer, width::Integer=32)

Shuffle a value from a directly indexed lane `lane`.
""" shfl

@doc """
    shfl_up(val, delta::Integer, width::Integer=32)

Shuffle a value from a lane with lower ID relative to caller.
""" shfl_up

@doc """
    shfl_down(val, delta::Integer, width::Integer=32)

Shuffle a value from a lane with higher ID relative to caller.
""" shfl_down

@doc """
    shfl_xor(val, mask::Integer, width::Integer=32)

Shuffle a value from a lane based on bitwise XOR of own lane ID with `mask`.
""" shfl_xor

@doc """
    shfl_sync(val, lane::Integer, width::Integer=32, threadmask::UInt32=0xffffffff)

Shuffle a value from a directly indexed lane `lane`. The default value for `threadmask`
performs the shuffle on all threads in the warp.
""" shfl_sync

@doc """
    shfl_up_sync(val, delta::Integer, width::Integer=32, threadmask::UInt32=0xffffffff)

Shuffle a value from a lane with lower ID relative to caller. The default value for
`threadmask` performs the shuffle on all threads in the warp.
""" shfl_up_sync

@doc """
    shfl_down_sync(val, delta::Integer, width::Integer=32, threadmask::UInt32=0xffffffff)

Shuffle a value from a lane with higher ID relative to caller. The default value for
`threadmask` performs the shuffle on all threads in the warp.
""" shfl_down_sync

@doc """
    shfl_xor_sync(val, mask::Integer, width::Integer=32, threadmask::UInt32=0xffffffff)

Shuffle a value from a lane based on bitwise XOR of own lane ID with `mask`. The default
value for `threadmask` performs the shuffle on all threads in the warp.
""" shfl_xor_sync


text:
include("../examples/baker.jl")
using LinearAlgebra
using JLD
using SharedArrays
using Distributed
function stable_sens(s,nSteps)
	vs = zeros(2)
	q = rand(2)
	q /= norm(q)
	le = 0.
	u = 2*pi*rand(2)
	u_trj = step(u, s, nSteps-1)
	du_trj = dstep(u_trj, s)
	x = pert(u_trj, 4)
	dJds = 0.
	for i = 1:nSteps
		vs .= du_trj[:,:,i]*vs + x[:,i]
		q .= du_trj[:,:,i]*q
		nm_q = norm(q)
		le += log(nm_q)/nSteps
		q ./= nm_q
		vs .-= dot(vs,q)*q
		x2 = (i==nSteps) ? step(u_trj[:,i],
			s, 1)[2,end] : u_trj[2,i+1]
		dJdu = [0., -4*sin(4*x2)]
		dJds += dot(dJdu, vs)/nSteps
	end
	#@show le, dJds
	return dJds
end
function get_stable_sens(s4)
	p = 4
	nSteps = 100
	s = zeros(p)
	# J = cos(4y)
	n_exps = size(s4)[1]
	dJds = zeros(n_exps)
	n_rep = 160000
	dJds_proc = SharedArray{Float64}(n_rep)
	for k=1:n_exps
		s[4] = s4[k]
		dJds_proc .= 0.
		t = @distributed for j=1:n_rep
			dJds_proc[j] = stable_sens(s,nSteps)/n_rep
		end
		wait(t)
		dJds[k] = sum(dJds_proc)
		@show dJds[k]
	end
	save("../data/stable_sens/dJds4.jld", "s4",
	     s4, "dJds", dJds)
end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
"""
    function mutilated()

A mutilated cow, from aspolito@csua.berkeley.edu

# Example
```jldoctest
julia> cowsay("I wouldn't have been so eager to make head cheese...", cow=Cowsay.mutilated)
 _______________________________________
/ I wouldn't have been so eager to make \\
\\ head cheese...                        /
 ---------------------------------------
       \\   \\_______
 v__v   \\  \\   O   )
 (oo)      ||----w |
 (__)      ||     ||  \\/\\
```
"""
function mutilated(;eyes="oo", tongue="  ", thoughts="\\")
the_cow =
"""
       $thoughts   \\_______
 v__v   $thoughts  \\   O   )
 ($eyes)      ||----w |
 (__)      ||     ||  \\/\\
  $tongue
"""
return the_cow
end


text:
using StaticArrays
using GeometryTypes

struct World{F<:Number,D<:Dogs,S<:Sheeps}
    boxsize::Base.RefValue{F}
    business_rate::F
    freedom_rate::F
    meanSleepyness::F
    meanSheepDiffusionTime::F
    dogs::D
    sheeps::S
end

function World(;
    v0 = 0,
    n_dogs = 0,
    boxsize::F = 0.0,
    business_rate::F = 20.0,
    freedom_rate::F = 20.0,
    meanSleepyness::F = 0.1,
    meanSheepDiffusionTime::F = 1000.0,
    motion::M = BrownianMotion(),
    sheeps = DenseSheeps{convert(Int,boxsize)}(10),
    ) where {F<:Number,M<:Motion}
    dogvec = SizedArray{Tuple{n_dogs},Dog{F}}(undef)
    @inbounds for i in 1:n_dogs
        ϕ = 2π * rand()
        dogvec[i] = Dog(
                        Point2{F}( boxsize .* rand(2) ),
                        Point2{F}( v0 .* [cos(ϕ),sin(ϕ)] ),
                     )
    end # for
    dogs = Dogs{M,n_dogs,Dog{F}}(dogvec, motion)
    return World{F,typeof(dogs),typeof(sheeps)}(Ref(boxsize),business_rate, freedom_rate, meanSleepyness, meanSheepDiffusionTime, dogs, sheeps)
end # function

function pbc( position::P, world::World ) where P
    p = position
    boxsize = world.boxsize[]
    tx = mod( p[1], boxsize )
    ty = mod( p[2], boxsize )
    while tx < 0; tx += boxsize; end
    while ty < 0; ty += boxsize; end

    return P( tx, ty )
end # function

function getSheepCoords( world::World{<:Number,<:Dogs,<:Sheeps{L}}, position ) where L
    x, y = position
    i = convert(typeof(L), (x * L) ÷ world.boxsize[]) + 1
    j = convert(typeof(L), (y * L) ÷ world.boxsize[]) + 1
    return i,j
end # function

function work!( world::World, dt )
    world.dogs.member .= work!.( Ref(world), world.dogs.member, dt )
    return nothing
end # function


text:
__precompile__()

module ArtBooster

using GR
using XGBoost
using Images

include("dataformats.jl")
include("boosting.jl")
include("rendering.jl")

export boostimage

end # module ArtBooster


text:
#   Unit "gcca.jl" of the Diagonalization.jl Package for Julia language
#
#   MIT License
#   Copyright (c) 2019-2021,
#   Marco Congedo, CNRS, Grenoble, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements the generalized maximum covariance analysis
#   and the generalized canonical correlation analysis.

"""
```julia
function gmca(𝐗::VecMat;
              covEst     :: StatsBase.CovarianceEstimator = SCM,
              dims       :: Into    = ○,
              meanX      :: Into    = 0,
          algorithm :: Symbol    = :OJoB,
          fullModel :: Bool      = false,
          sort      :: Bool      = true,
          init      :: VecMato   = ○,
          tol       :: Real      = 0.,
          maxiter   :: Int       = _maxiter(algorithm, eltype(𝐗[1])),
          verbose   :: Bool      = false,
          threaded  :: Bool      = true,
        eVar     :: TeVaro   = _minDim(𝐗),
        eVarMeth :: Function = searchsortedfirst,
        simple   :: Bool     = false)

```

Return a [LinearFilter](@ref) object.

**Generalized Maximum Covariance Analysis** of the set of ``m``
data matrices `𝐗`
using the given solving `algorithm` (*OJoB* by default).

If `fullModel` is true, the [gmca.3] problem here above is solved,
otherwise (default), the [gmca.2] problem here above is solved.

If `sort` is true (default), the column vectors of the matrices ``F_1,...,F_m``
are signed and permuted
as explained here above in [permutation for gMCA](@ref),
otherwise they will have arbitrary sign and will be in arbitrary order.

Regarding arguments `init`, `tol` and `maxiter`, see [Algorithms](@ref).

If `verbose` is true (false by default), the convergence attained
at each iteration will be printed in the REPL.

`eVar` and `eVarMeth` are used to define a
[subspace dimension](@ref) ``p`` using the accumulated regularized
eigenvalues of Eq. [gmca.7].

The default values are:
- `eVar` is set to the minimum dimension of the matrices in `𝐗`
- `eVarMeth=searchsortedfirst`

If `simple` is set to `true`, ``p`` is set equal to the dimension
of the covariance matrices that are computed on the matrices in `𝐗`,
which depends on the choice of `dims`,
and only the fields `.F` and `.iF`
are written in the constructed object.
This corresponds to the typical output of approximate diagonalization
algorithms.

if `threaded`=true (default) and the number of threads Julia is instructed
to use (the output of Threads.nthreads()), is higher than 1,
solving algorithms supporting multi-threading run in multi-threaded mode.
See [Algorithms](@ref) and
[these notes](https://marco-congedo.github.io/PosDefManifold.jl/dev/MainModule/#Threads-1)
on multi-threading.

**See also:** [MCA](@ref), [gCCA](@ref), [mAJD](@ref).

**Examples:**

```julia
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

####  Create data for testing the case k=1, m>1
# `t` is the number of samples,
# `m` is the number of datasets,
# `n` is the number of variables,
# `noise` must be smaller than 1.0. The smaller the noise,
#  the more data are correlated.
function getData(t, m, n, noise)
    # create m identical data matrices and rotate them by different
    # random orthogonal matrices V_1,...,V_m
    𝐕=[randU(n) for i=1:m] # random orthogonal matrices
    X=randn(n, t)  # data common to all subjects
    # each subject has this common part plus a random part
    𝐗=[𝐕[i]'*((1-noise)*X + noise*randn(n, t)) for i=1:m]
    return 𝐗
end

function getData(::Type{Complex{T}}, t, m, n, noise) where {T<:AbstractFloat}
    # create m identical data matrices and rotate them by different
    # random orthogonal matrices V_1,...,V_m
    𝐕=[randU(ComplexF64, n) for i=1:m] # random orthogonal matrices
    X=randn(ComplexF64, n, t)  # data common to all subjects
    # each subject has this common part plus a random part
    𝐗=[𝐕[i]'*((1-noise)*X + noise*randn(ComplexF64, n, t)) for i=1:m]
    return 𝐗
end

# REAL data: check that for the case m=2 gMCA gives the same result as MCA
t, m, n, noise = 20, 2, 6, 0.1
Xset=getData(t, m, n, noise)
Cx=(Xset[1]*Xset[1]')/t
Cy=(Xset[2]*Xset[2]')/t
Cxy=(Xset[1]*Xset[2]')/t

gm=gmca(Xset; simple=true)
m=mca(Cxy; simple=true)

@test (m.F[1]'*Cxy*m.F[2]) ≈ (gm.F[1]'*Cxy*gm.F[2])
# the following must be the identity matrix out of a possible sign ambiguity
@test abs.(m.F[1]'*gm.F[1]) ≈ I
@test abs.(m.F[2]'*gm.F[2]) ≈ I

# COMPLEX data: check that for the case m=2 gMCA gives the same result as MCA
t, m, n, noise = 20, 2, 6, 0.1
Xcset=getData(ComplexF64, t, m, n, noise)
Ccx=(Xcset[1]*Xcset[1]')/t
Ccy=(Xcset[2]*Xcset[2]')/t
Ccxy=(Xcset[1]*Xcset[2]')/t

gmc=gmca(Xcset; simple=true)
mc=mca(Ccxy; simple=true)

# for complex data just do a sanity check as the order of vectors
# is arbitrary
@test spForm(mc.F[1]'gmc.F[1])<0.01
@test spForm(mc.F[2]'gmc.F[2])<0.01

# REAL data: m>2 case
t, m, n, noise = 20, 4, 6, 0.1
Xset=getData(t, m, n, noise)

# ... selecting subspace dimension allowing an explained variance = 0.9
gm=gmca(Xset, eVar=0.9)

# name of the filter
gm.name

𝒞=Array{Matrix}(undef, 1, m, m)
for i=1:m, j=1:m 𝒞[1, i, j]=(Xset[i]*Xset[j]')/t end

using Plots
# plot regularized accumulated eigenvalues
plot(gm.arev)

# plot the original cross-covariance matrices and the rotated
# cross-covariance matrices

# Get all products 𝐔[i]' * 𝒞[l, i, j] * 𝐔[j]
function _rotate_crossCov(𝐔, 𝒞, m, k)
    𝒮=Array{Matrix}(undef, k, m, m)
    @inbounds for l=1:k, i=1:m, j=1:m 𝒮[l, i, j]=𝐔[i]'*𝒞[l, i, j]*𝐔[j] end
    return 𝒮
end

# Put all cross-covariances in a single matrix of dimension m*n x m*n for visualization
function 𝒞2Mat(𝒞::AbstractArray, m, k)
    n=size(𝒞[1, 1, 1], 1)
    C=Matrix{Float64}(undef, m*n, m*n)
    for i=1:m, j=1:m, x=1:n, y=1:n C[i*n-n+x, j*n-n+y]=𝒞[k, i, j][x, y] end
    return C
end

 C=𝒞2Mat(𝒞, m, 1)
 Cmax=maximum(abs.(C));
 h1 = heatmap(C, clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-covariances")
 𝒮=_rotate_crossCov(gm.F, 𝒞, m, 1)
 S=𝒞2Mat(𝒮, m, 1)
 Smax=maximum(abs.(S));
 h2 = heatmap(S, clim=(0, Smax), yflip=true, c=:amp, title="all rotated cross-covariances")
 📈=plot(h1, h2, size=(700,300))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FiggMCA.png")

```

 ![Figure gMCA](assets/FiggMCA.png)

 In the figure here above, the rotated cross-covariance matrices have the expected
 *strip-diagonal* form, that is, each block ``F_i^T\\frac{1}{t}(X_iX_j^T)F_j``,
 for ``i,j∈[1,...,m]``, is approximately diagonal. Each block is ``5⋅5`` because
 setting `eVar=0.9` the subspace dimension has been set to 5.

```julia
# COMPLEX data: m>2 case
t, m, n, noise = 20, 4, 6, 0.1
Xcset=getData(ComplexF64, t, m, n, noise)

# ... selecting subspace dimension allowing an explained variance = 0.9
gmc=gmca(Xcset, eVar=0.9)
```

"""
function gmca(𝐗::VecMat;
              covEst     :: StatsBase.CovarianceEstimator = SCM,
              dims       :: Into    = ○,
              meanX      :: Into    = 0,
          algorithm :: Symbol    = :OJoB,
          fullModel :: Bool      = false,
          sort      :: Bool      = true,
          init      :: VecMato   = ○,
          tol       :: Real      = 0.,
          maxiter   :: Int       = _maxiter(algorithm, eltype(𝐗[1])),
          verbose   :: Bool      = false,
          threaded  :: Bool      = true,
        eVar     :: TeVaro   = _minDim(𝐗),
        eVarMeth :: Function = searchsortedfirst,
        simple   :: Bool     = false)

   dims===○ && (dims=_set_dims(𝐗))
   _check_data(𝐗, 1, covEst, meanX, ○)===○ && return
   _check_data(𝐗, 2, covEst, meanX, ○)===○ && return
   (n, t)=dims==1 ? reverse(size(𝐗[1])) : size(𝐗[1])
   m=length(𝐗)
   args=("generalized Maximum Covariance Analysis", false)

   if algorithm ∈(:OJoB, :NoJoB)
      𝐔, 𝐕, λ, iter, conv=JoB(𝐗, m, 1, :d, algorithm, eltype(𝐗[1]);
               covEst=covEst, dims=dims, meanX=meanX,
               fullModel=fullModel, sort=sort, init=init,
                  tol=tol, maxiter=maxiter, verbose=verbose, threaded=threaded,
               eVar=eVar, eVarMeth=eVarMeth)
   # elseif...
   else
      throw(ArgumentError(📌*", gmca constructor: invalid `algorithm` argument"))
   end

   λ = _checkλ(λ) # make sure no imaginary noise is present (for complex data)

   simple ? LF(𝐔, 𝐕, Diagonal(λ), ○, ○, ○, args...) :
   begin
      p, arev = _getssd!(eVar, λ, n, eVarMeth) # find subspace
      LF([𝐔[i][:, 1:p] for i=1:m], [𝐕[i][1:p, :] for i=1:m], Diagonal(λ[1:p]), arev[p], λ, arev, args...)
   end
end

"""
```julia
function gcca(𝐗::VecMat;
              covEst     :: StatsBase.CovarianceEstimator = SCM,
              dims       :: Into    = ○,
              meanX      :: Into    = 0,
          algorithm :: Symbol    = :OJoB,
          sort      :: Bool      = true,
          init      :: VecMato   = ○,
          tol       :: Real      = 0.,
          maxiter   :: Int       = _maxiter(algorithm, eltype(𝐗[1])),
          verbose   :: Bool      = false,
          threaded  :: Bool      = true,
        eVar     :: TeVaro   = _minDim(𝐗),
        eVarMeth :: Function = searchsortedfirst,
        simple   :: Bool     = false)

```

Return a [LinearFilter](@ref) object.

**Generalized Canonical Correlation Analysis** of the set of ``m``
data matrices `𝐗`
using the given solving `algorithm` (*OJoB* by default).

If `sort` is true (default), the column vectors of the matrices ``F_1,...,F_m``
are signed and permuted
as explained here above in [permutation for gCCA](@ref),
otherwise they will have arbitrary sign and will be in arbitrary order.

Regarding arguments `init`, `tol` and `maxiter`, see [Algorithms](@ref).

If `verbose` is true (false by default), the convergence attained
at each iteration will be printed in the REPL.

`eVar` and `eVarMeth` are used to define a
[subspace dimension](@ref) ``p`` using the accumulated regularized
eigenvalues in Eq. [gcca.7]

The default values are:
- `eVar` is set to the minimum dimension of the matrices in `𝐗`
- `eVarMeth=searchsortedfirst`

If `simple` is set to `true`, ``p`` is set equal to the dimension
of the covariance matrices that are computed on the matrices in `𝐗`,
which depends on the choice of `dims`,
and only the fields `.F` and `.iF`
are written in the constructed object.
This corresponds to the typical output of approximate diagonalization
algorithms.

if `threaded`=true (default) and the number of threads Julia is instructed
to use (the output of Threads.nthreads()), is higher than 1,
solving algorithms supporting multi-threading run in multi-threaded mode.
See [Algorithms](@ref) and
[these notes](https://marco-congedo.github.io/PosDefManifold.jl/dev/MainModule/#Threads-1)
on multi-threading.

**See also:** [MCA](@ref), [gCCA](@ref), [mAJD](@ref).

**Examples:**

```julia
using Diagonalizations, LinearAlgebra, PosDefManifold, Test

####  Create data for testing the case k=1, m>1
# `t` is the number of samples,
# `m` is the number of datasets,
# `n` is the number of variables,
# `noise` must be smaller than 1.0. The smaller the noise,
#  the more data are correlated.
function getData(t, m, n, noise)
    # create m identical data matrices and rotate them by different
    # random orthogonal matrices V_1,...,V_m
    𝐕=[randU(n) for i=1:m] # random orthogonal matrices
    X=randn(n, t)  # data common to all subjects
    # each subject has this common part plus a random part
    𝐗=[𝐕[i]'*((1-noise)*X + noise*randn(n, t)) for i=1:m]
    return 𝐗
end

function getData(::Type{Complex{T}}, t, m, n, noise) where {T<:AbstractFloat}
    # create m identical data matrices and rotate them by different
    # random orthogonal matrices V_1,...,V_m
    𝐕=[randU(ComplexF64, n) for i=1:m] # random orthogonal matrices
    X=randn(ComplexF64, n, t)  # data common to all subjects
    # each subject has this common part plus a random part
    𝐗=[𝐕[i]'*((1-noise)*X + noise*randn(ComplexF64, n, t)) for i=1:m]
    return 𝐗
end

# REAL data: check that for the case m=2 gCCA gives the same result as CCA
t, m, n, noise = 20, 2, 6, 0.1
Xset=getData(t, m, n, noise)
Cx=(Xset[1]*Xset[1]')/t
Cy=(Xset[2]*Xset[2]')/t
Cxy=(Xset[1]*Xset[2]')/t

gc=gcca(Xset; simple=true)

c=cca(Hermitian(Cx), Hermitian(Cy), Cxy; simple=true)

@test (c.F[1]'*Cxy*c.F[2]) ≈ (gc.F[1]'*Cxy*gc.F[2])
@test gc.F[1]'*Cx*gc.F[1]≈I
@test gc.F[2]'*Cy*gc.F[2]≈I
D=gc.F[1]'*Cxy*gc.F[2]
@test norm(D-Diagonal(D))+1≈1.

# COMPLEX data: check that for the case m=2 gCCA gives the same result as CCA
t, m, n, noise = 20, 2, 6, 0.1
Xcset=getData(ComplexF64, t, m, n, noise)
Ccx=(Xcset[1]*Xcset[1]')/t
Ccy=(Xcset[2]*Xcset[2]')/t
Ccxy=(Xcset[1]*Xcset[2]')/t

gcc=gcca(Xcset; simple=true)
cc=cca(Hermitian(Ccx), Hermitian(Ccy), Ccxy; simple=true)

# for complex data just do a sanity check as the order of vectors
# is arbitrary. The following two tests currently fail
# @test spForm(cc.F[1]'gcc.F[1])<0.001
# @test spForm(cc.F[2]'gcc.F[2])<0.001

@test gcc.F[1]'*Ccx*gcc.F[1]≈I
@test gcc.F[2]'*Ccy*gcc.F[2]≈I
# sanity check only as there is noise in the complex case
D=gcc.F[1]'*Ccxy*gcc.F[2]
@test norm(D-Diagonal(D))/(n^2-n)<0.001

# REAL data: m>2 case
t, m, n, noise = 20, 4, 6, 0.1
Xset=getData(t, m, n, noise)

# ... selecting subspace dimension allowing an explained variance = 0.9
gc=gcca(Xset, eVar=0.9)

# name of the filter
gc.name

𝒞=Array{Matrix}(undef, 1, m, m)
for i=1:m, j=1:m 𝒞[1, i, j]=(Xset[i]*Xset[j]')/t end

using Plots
# plot regularized accumulated eigenvalues
plot(gc.arev)

# plot the original cross-covariance matrices and the rotated
# cross-covariance matrices

# Get all products 𝐔[i]' * 𝒞[l, i, j] * 𝐔[j]
function _rotate_crossCov(𝐔, 𝒞, m, k)
    𝒮=Array{Matrix}(undef, k, m, m)
    @inbounds for l=1:k, i=1:m, j=1:m 𝒮[l, i, j]=𝐔[i]'*𝒞[l, i, j]*𝐔[j] end
    return 𝒮
end

# Put all cross-covariances in a single matrix of dimension m*n x m*n for visualization
function 𝒞2Mat(𝒞::AbstractArray, m, k)
    n=size(𝒞[1, 1, 1], 1)
    C=Matrix{Float64}(undef, m*n, m*n)
    for i=1:m, j=1:m, x=1:n, y=1:n C[i*n-n+x, j*n-n+y]=𝒞[k, i, j][x, y] end
    return C
end

 C=𝒞2Mat(𝒞, m, 1)
 Cmax=maximum(abs.(C));
 h1 = heatmap(C, clim=(-Cmax, Cmax), yflip=true, c=:bluesreds, title="all cross-covariances")
 𝒮=_rotate_crossCov(gc.F, 𝒞, m, 1)
 S=𝒞2Mat(𝒮, m, 1)
 h2 = heatmap(S, clim=(0, 1), yflip=true, c=:amp, title="all rotated cross-covariances")
 📈=plot(h1, h2, size=(700,300))
# savefig(📈, homedir()*"\\Documents\\Code\\julia\\Diagonalizations\\docs\\src\\assets\\FiggCCA.png")

```

  ![Figure gCCA](assets/FiggCCA.png)

  In the figure here above, the rotated cross-covariance matrices have the expected
  *strip-diagonal* form, that is, each block ``F_i^T\\frac{1}{t}(X_iX_j^T)F_j``,
  for ``i,j∈[1,...,m]``, is approximately diagonal. Each block is ``6⋅6`` because
  setting `eVar=0.9` has not reduced the original dimension.
  The solution is similar to the [gMCA](@ref), but here the diagonal
  of the rotated block matrix is the identity.

```julia
# COMPLEX data: m>2 case
t, m, n, noise = 20, 4, 6, 0.1
Xcset=getData(ComplexF64, t, m, n, noise)

# ... selecting subspace dimension allowing an explained variance = 0.9
gcc=gcca(Xcset, eVar=0.9)
```

"""
function gcca(𝐗::VecMat;
              covEst     :: StatsBase.CovarianceEstimator = SCM,
              dims       :: Into    = ○,
              meanX      :: Into    = 0,
          algorithm :: Symbol    = :OJoB,
          sort      :: Bool      = true,
          init      :: VecMato   = ○,
          tol       :: Real      = 0.,
          maxiter   :: Int       = _maxiter(algorithm, eltype(𝐗[1])),
          verbose   :: Bool      = false,
          threaded  :: Bool      = true,
        eVar     :: TeVaro   = _minDim(𝐗),
        eVarMeth :: Function = searchsortedfirst,
        simple   :: Bool     = false)

   dims===○ && (dims=_set_dims(𝐗))
   _check_data(𝐗, 1, covEst, meanX, ○)===○ && return
   _check_data(𝐗, 2, covEst, meanX, ○)===○ && return
   (n, t)=dims==1 ? reverse(size(𝐗[1])) : size(𝐗[1])
   m=length(𝐗)
   args=("generalized Canonical Correlation Analysis", false)

   if algorithm ==:OJoB #∈(:OJoB)
      𝐔, 𝐕, λ, iter, conv=JoB(𝐗, m, 1, :d, algorithm, eltype(𝐗[1]);
               covEst=covEst, dims=dims, meanX=meanX,
               fullModel=false, preWhite=true, sort=sort, init=init,
                  tol=tol, maxiter=maxiter, verbose=verbose, threaded=threaded,
               eVar=eVar, eVarMeth=eVarMeth)
   # elseif...
   else
      algorithm == :NoJoB && @warn "The NoJoB algorithm does not suit gCCA."
      throw(ArgumentError(📌*", gcca constructor: invalid `algorithm` argument"))
   end

   λ = _checkλ(λ) # make sure no imaginary noise is present (for complex data)

   simple ? LF(𝐔, 𝐕, Diagonal(λ), ○, ○, ○, args...) :
   begin
      p, arev = _getssd!(eVar, λ, n, eVarMeth) # find subspace
      LF([𝐔[i][:, 1:p] for i=1:m], [𝐕[i][1:p, :] for i=1:m], Diagonal(λ[1:p]), arev[p], λ, arev, args...)
   end
end


text:
# Objects
str2obj = Dict("all" => 0,
               "bus" => 1,
               "generator" => 2,
               "branch" => 3,
               "shunt" => 4,
               "load" => 5,
               "variable generator" => 6,
               "battery" => 7,
               "unknown" => 8)

# Flags
str2flag = Dict("variable" => 0x01,
                "fixed" => 0x02,
                "bounded" => 0x04,
                "sparse" => 0x08)

# Quantities
str2q_bus = Dict("all" => 0xFF,
                 "voltage magnitude" => 0x01,
                 "voltage angle" => 0x02)

str2q_branch = Dict("all" => 0xFF,
                    "tap ratio" => 0x01,
                    "phase shift" => 0x02)

str2q_gen = Dict("all" => 0xFF,
                 "active power" => 0x01,
                 "reactive power" => 0x02)

str2q_shunt = Dict("all" => 0xFF,
                   "susceptance" => 0x01)

str2q_load = Dict("all" => 0xFF,
                  "active power" => 0x01,
                  "reactive power" => 0x02)

str2q_vargen = Dict("all" => 0xFF,
                    "active power" => 0x01,
                    "reactive power" => 0x02)

str2q_bat = Dict("all" => 0xFF,
                 "charging power" => 0x01,
                 "energy level" => 0x02)

str2q = Dict("all" => Dict("all" => 0),
             "bus" => str2q_bus,
             "branch" => str2q_branch,
             "generator" => str2q_gen,
             "shunt" => str2q_shunt,
             "load" => str2q_load,
             "variable generator" => str2q_vargen,
             "battery" => str2q_bat)

# Properties
str2prop_bus = Dict("any" => 0x00,
                    "slack" => 0x01,
                    "regulated by generator" => 0x02,
                    "regulated by transformer" => 0x04,
                    "regulated by shunt" => 0x08,
                    "not regulated by generator" => 0x10,
                    "not slack" => 0x20)

str2prop_branch = Dict("any" => 0x00,
                       "tap changer" => 0x01,
                       "tap changer - v" => 0x02,
                       "tap changer - Q" => 0x04,
                       "phase shifter" => 0x08,
                       "not on outage" => 0x10)

str2prop_gen = Dict("any" => 0x00,
                    "slack" => 0x01,
                    "regulator" => 0x02,
                    "not regulator" => 0x04,
                    "not slack" => 0x08,
                    "not on outage" => 0x10,
                    "adjustable active power" => 0x20)

str2prop_shunt =  Dict("any" => 0x00,
                       "switching - v" => 0x01)

str2prop_load = Dict("any" => 0x00,
                     "adjustable active power" => 0x01)

str2prop_vargen = Dict("any" => 0x00)

str2prop_bat = Dict("any" => 0x00)

str2prop = Dict("all" => Dict("all" => 0),
                "bus" => str2prop_bus,
                "branch" => str2prop_branch,
                "generator" => str2prop_gen,
                "shunt" => str2prop_shunt,
                "load" => str2prop_load,
                "variable generator" => str2prop_vargen,
                "battery" => str2prop_bat)


text:
#_______________________________________________________________________________
# Input variables

MaxCareerPathLength = 40
# GuyCareerPaths , InitManpower , MPObjectives

#_______________________________________________________________________________
# Variables number

    # Compute which careers could each Initial manpower subset go with

InitMPPartCP = Array{Vector{Int}}(length(InitManpower))
InitMPDivisionNb = 0
for i in 1:length(InitManpower)
    InitIpartCP = Vector{Int}()
    InitManpower[i].Seniority += 1
    while length(InitIpartCP) == 0 && InitManpower[i].Seniority > 0
        InitManpower[i].Seniority -= 1
        for j in 1:length(GuyCareerPaths)
            index = get_status_index_after_duration(GuyCareerPaths[j], InitManpower[i].Seniority + 1)
            if index == -1
                continue
            end
            if contains(get_Name_Level(GuyCareerPaths[j].Path[index], AcademicLevel), InitManpower[i].Academiclvl) && contains(get_Name_Level(GuyCareerPaths[j].Path[index], Affiliation), InitManpower[i].Affiliation)
                push!(InitIpartCP, j)
            end
        end
    end
    if length(InitIpartCP) == 0
        InitManpower[i].Nb = 0
        warn("Initial subpopulation \"Academ:$(InitManpower[i].Academiclvl),Seniority:$(InitManpower[i].ActualSeniority)\" is not considered.")
    end
    InitMPPartCP[i] = InitIpartCP
    InitMPDivisionNb += length(InitIpartCP)
end

AnnualRecDivNb = length(GuyCareerPaths) * NBYears

DeviationVarNb = length(MPObjectives) * NBYears * 2

DependDeviationVar = NBConstDepend * NBYears * 2

VarNb = InitMPDivisionNb + AnnualRecDivNb + DeviationVarNb + DependDeviationVar

#_______________________________________________________________________________
# Constraints number

InitConst = length(InitManpower)

RecruitmentConst = NBYears * 2

GoalsConst = (NBYears * length(MPObjectives))*2

if AllowableDeviation < 0
    RecruitmentDeviation = 0
else
    RecruitmentDeviation = (NBYears - 1) * 2
end

DependenciesConst = NBYears * NBConstDepend

ConstNb = InitConst + RecruitmentConst + GoalsConst + RecruitmentDeviation + DependenciesConst

#_______________________________________________________________________________
# MIP construction
    # Definition of Arrays

Cost = Array{Float64}(VarNb) # Done

A = zeros(Float64, ConstNb, VarNb) # Init Done, Rec Done, Goal Done

sense = Array{Char}(ConstNb) # Done

b = zeros(Float64, ConstNb)#Array{Float64}(ConstNb) # Done

vartypes = Array{Symbol}(VarNb) # Done

    # Fill Cost
Cost[1:(InitMPDivisionNb)] = 0
Cost[1 + InitMPDivisionNb:(InitMPDivisionNb + AnnualRecDivNb)] = 0
for i in 1:length(MPObjectives)
    for j in 1:NBYears
        Cost[(InitMPDivisionNb + AnnualRecDivNb) +                                  (j - 1)*length(MPObjectives) + i ] = 10^(MPObjectives[i].priority)
        Cost[(InitMPDivisionNb + AnnualRecDivNb) + (length(MPObjectives)*NBYears) + (j - 1)*length(MPObjectives) + i ] = 10^(MPObjectives[i].priority)
    end
end
tempoDepLine = 0
for i in 1:length(Dependencies)
    for j in 1:NBYears
        for k in 1:length(Dependencies[i].Sets)
            Cost[(InitMPDivisionNb + AnnualRecDivNb + DeviationVarNb) +                             (j - 1)*NBConstDepend + tempoDepLine + k ] = 10^(Dependencies[i].Priority)
            Cost[(InitMPDivisionNb + AnnualRecDivNb + DeviationVarNb) + (NBConstDepend * NBYears) + (j - 1)*NBConstDepend + tempoDepLine + k ] = 10^(Dependencies[i].Priority)
        end
    end
    tempoDepLine += length(Dependencies[i].Sets)
end
    # Fill A
        # Initial Manpower Constraints
AdvacmentY = 0
for i in 1:InitConst
    for j in 1:length(InitMPPartCP[i])
        A[i, AdvacmentY + j] = 1
    end
    AdvacmentY += length(InitMPPartCP[i])
end

        # Annual Recruitment
            # Max recruitment
for i in 1:Int(RecruitmentConst/2)
    for j in 1:length(GuyCareerPaths)
        A[InitConst + i, InitMPDivisionNb + (i-1)*length(GuyCareerPaths) + j] = 1
    end
end
            # Min recruitment   ________________________________________________ To verify
for i in 1:Int(RecruitmentConst/2)
    for j in 1:length(GuyCareerPaths)
        A[InitConst + NBYears + i, InitMPDivisionNb + (i-1)*length(GuyCareerPaths) + j] = 1
    end
end

            # compute which careers go with each goal

SetsPartCPinGoals = Vector{Vector{Vector{Int}}}()

for i in 1:length(GuyCareerPaths)
    AllParticipations = Vector{Vector{Int}}()
    for j in 1:length(GuyCareerPaths[i].Path)
        Participation = Vector{Int}()

        for k in 1:length(MPObjectives)
            isparticipating = true
            for subObj in 1:length(MPObjectives[k].targetLevel)
                if !contains(get_Name_Level(GuyCareerPaths[i].Path[j], MPObjectives[k].targetLevel[subObj]), MPObjectives[k].Objectives[subObj])#_______________________________________________________________________________
                    isparticipating = false#push!(Participation, k)
                    break
                end
            end

            if isparticipating
                push!(Participation, k)
            end
        end
        push!(AllParticipations, Participation)
    end
    push!(SetsPartCPinGoals, AllParticipations)
end

            # Compute Career Paths Length

for i in 1:length(GuyCareerPaths)
    compute_Length(GuyCareerPaths[i])
end

        # Goals___________________________________________________ (upper bound)

            # Initial Manpower participation

TempActualVar = 1
for imp in 1:length(InitManpower)
    for i in 1:length(InitMPPartCP[imp])
        for j in (1 + InitManpower[imp].Seniority):GuyCareerPaths[InitMPPartCP[imp][i]].Length
            indx = get_status_index_after_duration(GuyCareerPaths[InitMPPartCP[imp][i]], j)
            if indx == -1
                break
            end
            for k in 1:length(SetsPartCPinGoals[InitMPPartCP[imp][i]][indx])
                if (j - InitManpower[imp].Seniority - 1) <= NBYears
                    A[InitConst + RecruitmentConst + (j - InitManpower[imp].Seniority - 1)*length(MPObjectives) + SetsPartCPinGoals[InitMPPartCP[imp][i]][indx][k], TempActualVar] = 1
                end
            end
        end
        TempActualVar += 1
    end
end

            # Annual Recruitment participation

for y in 1:NBYears
    for i in 1:length(GuyCareerPaths)
        for j in 1:GuyCareerPaths[i].Length # j => year for const
            indx = get_status_index_after_duration(GuyCareerPaths[i], j)
            if indx == -1
                break
            end
            for k in 1:length(SetsPartCPinGoals[i][indx])
                if (y + j - 1) <= NBYears
                    A[InitConst + RecruitmentConst + (y - 1)*length(MPObjectives) + (j - 1)*length(MPObjectives) + SetsPartCPinGoals[i][indx][k], InitMPDivisionNb + (y - 1)*length(GuyCareerPaths) + i] = 1 - AttritionAtYear(GuyCareerPaths[i], j)
                end
            end
        end
    end
end

            # deviation variables

for y in 1:NBYears
    for eq in 1:length(MPObjectives)

        A[InitConst + RecruitmentConst + (y - 1)*length(MPObjectives) + eq, InitMPDivisionNb + AnnualRecDivNb +                                  (y - 1)*length(MPObjectives) + eq] =  1
        A[InitConst + RecruitmentConst + (y - 1)*length(MPObjectives) + eq, InitMPDivisionNb + AnnualRecDivNb + length(MPObjectives) * NBYears + (y - 1)*length(MPObjectives) + eq] = -1

    end
end

    # Goals___________________________________________________ (lower bound)

    # Initial Manpower participation

TempActualVar = 1
for imp in 1:length(InitManpower)
    for i in 1:length(InitMPPartCP[imp])
        for j in (1 + InitManpower[imp].Seniority):GuyCareerPaths[InitMPPartCP[imp][i]].Length
            indx = get_status_index_after_duration(GuyCareerPaths[InitMPPartCP[imp][i]], j)
            if indx == -1
                break
            end
            for k in 1:length(SetsPartCPinGoals[InitMPPartCP[imp][i]][indx])
                if (j - InitManpower[imp].Seniority - 1) <= NBYears
                    A[InitConst + RecruitmentConst + Int(GoalsConst/2) + (j - InitManpower[imp].Seniority - 1)*length(MPObjectives) + SetsPartCPinGoals[InitMPPartCP[imp][i]][indx][k], TempActualVar] = 1
                end
            end
        end
        TempActualVar += 1
    end
end

    # Annual Recruitment participation

for y in 1:NBYears
    for i in 1:length(GuyCareerPaths)
        for j in 1:GuyCareerPaths[i].Length # j => year for const
            indx = get_status_index_after_duration(GuyCareerPaths[i], j)
            if indx == -1
                break
            end
            for k in 1:length(SetsPartCPinGoals[i][indx])
                if (y + j - 1) <= NBYears
                    A[InitConst + RecruitmentConst + Int(GoalsConst/2) + (y - 1)*length(MPObjectives) + (j - 1)*length(MPObjectives) + SetsPartCPinGoals[i][indx][k], InitMPDivisionNb + (y - 1)*length(GuyCareerPaths) + i] = 1 - AttritionAtYear(GuyCareerPaths[i], j)
                end
            end
        end
    end
end

    # deviation variables
for y in 1:NBYears
    for eq in 1:length(MPObjectives)

        A[InitConst + RecruitmentConst + Int(GoalsConst/2) + (y - 1)*length(MPObjectives) + eq, InitMPDivisionNb + AnnualRecDivNb +                                  (y - 1)*length(MPObjectives) + eq] =  1
        A[InitConst + RecruitmentConst + Int(GoalsConst/2) + (y - 1)*length(MPObjectives) + eq, InitMPDivisionNb + AnnualRecDivNb + length(MPObjectives) * NBYears + (y - 1)*length(MPObjectives) + eq] = -1

    end
end

            # Recruitment deviation *******************************************
if AllowableDeviation >= 0
    for i in 1:NBYears-1
        for j in 1:length(GuyCareerPaths)
            A[InitConst + RecruitmentConst + GoalsConst + i, InitMPDivisionNb + (i)*length(GuyCareerPaths) + j] = 1
            A[InitConst + RecruitmentConst + GoalsConst + i, InitMPDivisionNb + (i-1)*length(GuyCareerPaths) + j] = -( 1 +  AllowableDeviation)
        end
    end

    for i in 1:NBYears-1
        for j in 1:length(GuyCareerPaths)
            A[InitConst + RecruitmentConst + GoalsConst + (NBYears - 1) + i, InitMPDivisionNb + (i)*length(GuyCareerPaths) + j] = 1
            A[InitConst + RecruitmentConst + GoalsConst + (NBYears - 1) + i, InitMPDivisionNb + (i-1)*length(GuyCareerPaths) + j] = -(1 - AllowableDeviation)
        end
    end
end

            # Dependencies Equations *******************************************
tempoDepLine = 0
for i in 1:length(Dependencies)
    for j in 1:NBYears
        for k in 1:length(Dependencies[i].Sets)
            # Career Paths Variables
            for l in 1:length(Dependencies[i].Sets)
                if l == k
                    for m in 1:length(Dependencies[i].Sets[k])
                        A[InitConst + RecruitmentConst + GoalsConst + RecruitmentDeviation + (j-1)*NBConstDepend + tempoDepLine + k, InitMPDivisionNb + length(GuyCareerPaths)*(j - 1) + Dependencies[i].Sets[l][m]] = 1 - Dependencies[i].Percentages[k]
                    end
                else
                    for m in 1:length(Dependencies[i].Sets[k])
                        A[InitConst + RecruitmentConst + GoalsConst + RecruitmentDeviation + (j-1)*NBConstDepend + tempoDepLine + k, InitMPDivisionNb + length(GuyCareerPaths)*(j - 1) + Dependencies[i].Sets[l][m]] = Dependencies[i].Percentages[k]
                    end
                end
            end

            # Deviation Variables
            A[(InitMPDivisionNb + AnnualRecDivNb + DeviationVarNb) +                             (j - 1)*NBConstDepend + tempoDepLine + k ] = 1
            A[(InitMPDivisionNb + AnnualRecDivNb + DeviationVarNb) + (NBConstDepend * NBYears) + (j - 1)*NBConstDepend + tempoDepLine + k ] = -1
        end
    end
    tempoDepLine += length(Dependencies[i].Sets)
end

    # Fill sense

sense[1:InitConst] = '='
sense[(InitConst + 1):(InitConst + Int(RecruitmentConst/2))] = '<'
sense[(InitConst + Int(RecruitmentConst/2) + 1):(InitConst + RecruitmentConst)] = '>'
sense[(InitConst + RecruitmentConst + 1):(InitConst + RecruitmentConst + Int(GoalsConst/2))] = '<'
sense[(InitConst + RecruitmentConst + Int(GoalsConst/2) + 1):(InitConst + RecruitmentConst + GoalsConst)] = '>'
if AllowableDeviation >= 0
    sense[(InitConst + RecruitmentConst + GoalsConst + 1):(InitConst + RecruitmentConst + GoalsConst + (NBYears - 1))] = '<'
    sense[(InitConst + RecruitmentConst + GoalsConst + (NBYears - 1) + 1):(InitConst + RecruitmentConst + GoalsConst + (NBYears - 1) + (NBYears - 1))] = '>'
end
sense[(InitConst + RecruitmentConst + GoalsConst + RecruitmentDeviation + 1):(InitConst + RecruitmentConst + GoalsConst + RecruitmentDeviation + DependenciesConst)] = '='

    # Fill b
for i in 1:InitConst
    b[i] = InitManpower[i].Nb
end
b[(InitConst +                           1):(InitConst + Int(RecruitmentConst/2))] = MaxRecruitment
b[(InitConst + Int(RecruitmentConst/2) + 1):(InitConst +  RecruitmentConst   )] = MinRecruitment
# Goals second part
index = InitConst + RecruitmentConst + 1
for y in 1:NBYears
    for i in 1:length(MPObjectives)
        b[index] = MPObjectives[i].Number*(1 + (MPObjectives[i].InitTolerance - (MPObjectives[i].Alfa)*(y-1)) )
        index += 1
    end
end
for y in 1:NBYears
    for i in 1:length(MPObjectives)
        b[index] = MPObjectives[i].Number*(1 - (MPObjectives[i].InitTolerance - (MPObjectives[i].Alfa)*(y-1)))
        index += 1
    end
end

    # Fill vartypes

if IntegerSolution
    vartypes[1:(InitMPDivisionNb + AnnualRecDivNb)] = :Int
    vartypes[(InitMPDivisionNb + AnnualRecDivNb + 1):end] = :Cont
else
    vartypes[1:end] = :Cont
end

writedlm("matrix.txt", A)

#_______________________________________________________________________________
# MIP resolution
println("Solving LP...")
using MathProgBase, CPLEX#, Gurobi
stattoendStart = now()
if IntegerSolution
    sol = mixintprog(Cost, A, sense, b, vartypes, 0, Inf, CplexSolver(CPXPARAM_MIP_Tolerances_MIPGap=Tolerances_MIPGap))#CbcSolver(allowableGap=0.8)) #GurobiSolver(Presolve=0)
else
    sol = linprog(Cost, A, sense, b, 0, Inf, CplexSolver(CPX_PARAM_EPINT = 0.5))
end
stattoendEnd = now()
println( "ended with: $(sol.status). Elapsed time: $(Dates.canonicalize(Dates.CompoundPeriod(Dates.Millisecond(stattoendEnd - stattoendStart))))." )
writedlm("solution.txt", sol.sol)

#_______________________________________________________________________________
# Computing results

RequirementFulfillment = zeros(Float64, NBYears, length(MPObjectives))
sumi = 0
for i in 1:NBYears
    for j in 1:length(MPObjectives)
        for k in 1:(InitMPDivisionNb + AnnualRecDivNb)
            if A[InitConst + RecruitmentConst + (i - 1)*length(MPObjectives) + j, k] != 0
                RequirementFulfillment[i, j] +=  sol.sol[k]*A[InitConst + RecruitmentConst + (i - 1)*length(MPObjectives) + j, k]#round(Int,)
            end
        end
    end
end

YearlyRecruitment = zeros(Float64, NBYears)

for i in 1:NBYears
    for j in 1:length(GuyCareerPaths)
        YearlyRecruitment[i] += sol.sol[InitMPDivisionNb + (i - 1)*length(GuyCareerPaths) + j]
    end
end

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
"Market state in goldencross strategy."
function goldencrossmktstate(mafast::Float64, maslow::Float64)
  if isnan(mafast) || isnan(maslow)
    return :undefined
  end
  if mafast > maslow
    return :trendup
  end
  if mafast < maslow
    return :trenddown
  end
  # this would be rare
  return :undefined
end

"""
Target position for goldencross strategy.
This simplest form involves only market orders, long-side enter.
...
Returns `(poschg::Int64, Vector[limitprice, stopprice]`.
"""
function goldencrossposlogic(mktstate::Symbol,
                             targetqty::Int64,
                             position_actual_mut::Vector{Int64})
  posact = position_actual_mut[1]
  if posact < 0
    # should not happen: long-only strategy
    return (-posact, Array(Float64, 0))
  end
  if posact > targetqty && mktstate == :trendup
    # should not happen if only this strategy is active
    # sell excess (TODO do we need this? keeping for completeness only)
    return (-(posact-targetqty), Array(Float64, 0))
  end

  @match (posact < targetqty, mktstate) begin
    # enter long position: market buy
    (true, :trendup) => (targetqty - posact, Array(Float64, 0))
    # exit from long position: market sell
    (_, :trenddown), if posact > 0 end => (-posact, Array(Float64, 0))
    # position larger than target (should not normally happen)
    (false, :trenddown), if posact > 0 end => (-posact, Array(Float64, 0))
    # no match among action-triggering contidions above
    othercase => (0, Array(Float64, 0))
  end
end

"Target signal for goldencross strategy."
function goldencrosstarget(s_ohlc::Signal{OHLC},
                           ohlc_inds::Dict{Symbol,Int64},
                           position_actual_mut::Vector{Int64},
                           targetqty::Int64,
                           nsma_fast::Int64 = 50,
                           nsma_slow::Int64 = 200)
  posact = position_actual_mut[1]
  targetqty > 0 || error("target quantity must be positive")

  # signals updating at each timestep (aka OHLC bar)
  buffpclose!(buff::Vector{Float64}, tohlc::OHLC) =
    sighistbuffer!(buff, tohlc[2][ohlc_inds[:close]])
  pcloseinit = s_ohlc.value[2][ohlc_inds[:close]]
  s_sma_fast = map(mean,
                   foldp(buffpclose!, initbuff(nsma_fast, pcloseinit), s_ohlc),
                   typ=Float64)
  s_sma_slow = map(mean,
                   foldp(buffpclose!, initbuff(nsma_slow, pcloseinit), s_ohlc),
                   typ=Float64)

  # market state signal
  s_mktstate = map(goldencrossmktstate, s_sma_fast, s_sma_slow, typ=Symbol)

  # target position updates only when market state input changes
  s_target = map(mks -> goldencrossposlogic(mks,
                                            targetqty,
                                            position_actual_mut),
                 s_mktstate)
  # targeting tuple to pass to the order processing function
  return s_target
end


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta
import ProtoBuf.google.protobuf

mutable struct SignCertificateRequest <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SignCertificateRequest(; kwargs...)
        obj = new(meta(SignCertificateRequest), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct SignCertificateRequest
const __meta_SignCertificateRequest = Ref{ProtoMeta}()
function meta(::Type{SignCertificateRequest})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SignCertificateRequest)
            __meta_SignCertificateRequest[] = target = ProtoMeta(SignCertificateRequest)
            allflds = Pair{Symbol,Union{Type,String}}[:id => AbstractString, :token => AbstractString, :trust_domain => AbstractString, :namespace => AbstractString, :certificate_signing_request => Vector{UInt8}]
            meta(target, SignCertificateRequest, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SignCertificateRequest[]
    end
end
function Base.getproperty(obj::SignCertificateRequest, name::Symbol)
    if name === :id
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :token
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :trust_domain
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :namespace
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :certificate_signing_request
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    else
        getfield(obj, name)
    end
end

mutable struct SignCertificateResponse <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function SignCertificateResponse(; kwargs...)
        obj = new(meta(SignCertificateResponse), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct SignCertificateResponse
const __meta_SignCertificateResponse = Ref{ProtoMeta}()
function meta(::Type{SignCertificateResponse})
    ProtoBuf.metalock() do
        if !isassigned(__meta_SignCertificateResponse)
            __meta_SignCertificateResponse[] = target = ProtoMeta(SignCertificateResponse)
            allflds = Pair{Symbol,Union{Type,String}}[:workload_certificate => Vector{UInt8}, :trust_chain_certificates => Base.Vector{Vector{UInt8}}, :valid_until => ProtoBuf.google.protobuf.Timestamp]
            meta(target, SignCertificateResponse, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_SignCertificateResponse[]
    end
end
function Base.getproperty(obj::SignCertificateResponse, name::Symbol)
    if name === :workload_certificate
        return (obj.__protobuf_jl_internal_values[name])::Vector{UInt8}
    elseif name === :trust_chain_certificates
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{Vector{UInt8}}
    elseif name === :valid_until
        return (obj.__protobuf_jl_internal_values[name])::ProtoBuf.google.protobuf.Timestamp
    else
        getfield(obj, name)
    end
end

# service methods for CA
const _CA_methods = MethodDescriptor[
        MethodDescriptor("SignCertificate", 1, SignCertificateRequest, SignCertificateResponse)
    ] # const _CA_methods
const _CA_desc = ServiceDescriptor("dapr.proto.sentry.v1.CA", 1, _CA_methods)

CA(impl::Module) = ProtoService(_CA_desc, impl)

mutable struct CAStub <: AbstractProtoServiceStub{false}
    impl::ProtoServiceStub
    CAStub(channel::ProtoRpcChannel) = new(ProtoServiceStub(_CA_desc, channel))
end # mutable struct CAStub

mutable struct CABlockingStub <: AbstractProtoServiceStub{true}
    impl::ProtoServiceBlockingStub
    CABlockingStub(channel::ProtoRpcChannel) = new(ProtoServiceBlockingStub(_CA_desc, channel))
end # mutable struct CABlockingStub

SignCertificate(stub::CAStub, controller::ProtoRpcController, inp::SignCertificateRequest, done::Function) = call_method(stub.impl, _CA_methods[1], controller, inp, done)
SignCertificate(stub::CABlockingStub, controller::ProtoRpcController, inp::SignCertificateRequest) = call_method(stub.impl, _CA_methods[1], controller, inp)

export SignCertificateRequest, SignCertificateResponse, CA, CAStub, CABlockingStub, SignCertificate


text:
__precompile__()

module HermiteInterpolation

using   #Reduce,
        SpecialFunctions,
        Base.Cartesian,
   #     Sobol,
   #     GPUArrays,
        Compat,
        HermiteQuadratureRules,
        StaticArrays

# import HermiteQuadratureRules:
#         hermite_nodes,
#         hermite_nodes_diffs,
#         hermite_nodes_posi_diffs,
#         hermite_weights

export Marginal,
        Marginalize,
        _marginalize,
        SparseGrid

# include("sobol.jl")
# include("dfttools.jl")
include("sparse_grids/sparse_grids.jl")
include("polynomial/polynomials.jl")
include("polynomial/solution_matrix.jl")
# include("high_degree_polynomial/polynomials.jl")

# export  sobol_vec,
#         max_degree_poly,
#         lnorm_poly,
#         create_linear_map,
#         linear_map_stable

#function __init__()
#    Reduce.Rational(false)
#end

end # module


text:
############################################################
# judiExtendedSource ## ####################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017

export judiExtendedSource, judiExtendedSourceException

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

struct judiExtendedSource{vDT<:Number} <: joAbstractLinearOperator{vDT,vDT}
    name::String
    m::Integer
    n::Integer
    info::Info
    wavelet
    weights
end

mutable struct judiExtendedSourceException <: Exception
    msg :: String
end

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

## outer constructors

"""
    judiExtendedSource
        name::String
        m::Integer
        n::Integer
        info::Info
        wavelet
        weights
Abstract sparse vector for an extended source. The `judiRHS` vector has the\\
dimensions of the full time history of the wavefields, but contains only the wavelet, \\
as well as the weights of the extended source.
Constructor
==========
    judiExtendedSource(info, wavelet, weights)
Examples
========
Assuming `Pw` ia a projection operators of type `judiLRWF` and `w` is a seismic weight \\
vector of type `judiWeights`, then a `judiExtendedSource` vector can be created as follows:
    q_ext = adjoint(Pw)*w    # abstract extended source
"""
function judiExtendedSource(info,wavelet,weights;vDT::DataType=Float32)
    vDT == Float32 || throw(judiExtendedSourceException("Domain type not supported"))
    # length of vector
    m = info.n * sum(info.nt)
    n = 1
    return judiExtendedSource{Float32}("judiExtendedSource",m,n,info,wavelet,weights)
end

####################################################################
## overloaded Base functions

# conj(jo)
conj(A::judiExtendedSource{vDT}) where vDT =
    judiExtendedSource{vDT}("conj("*A.name*")",A.m,A.n,A.info,A.wavelet,A.weights)

# transpose(jo)
transpose(A::judiExtendedSource{vDT}) where vDT =
    judiExtendedSource{vDT}(""*A.name*".'",A.n,A.m,A.info,A.wavelet,A.weights)

# adjoint(jo)
adjoint(A::judiExtendedSource{vDT}) where vDT =
    judiExtendedSource{vDT}(""*A.name*"'",A.n,A.m,A.info,A.wavelet,A.weights)

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

# +(judiExtendedSource,judiExtendedSource)
function +(A::judiExtendedSource{avDT}, B::judiExtendedSource{bvDT}) where {avDT, bvDT}
    # Error checking
    size(A) == size(B) || throw(judiExtendedSourceException("Shape mismatch: A:$(size(A)), v: $(size(B))"))
    compareInfo(A.info, B.info) == true || throw(judiExtendedSourceException("info mismatch"))

    # Size
    m = A.info.n * sum(A.info.nt)
    n = 1

    # wavelet and weights
    A.wavelet == B.wavelet ? nothing : throw(judiExtendedSourceException("Can only add two extended
                                                                          sources with same wavelet"))
    wavelet = A.wavelet
    weights = A.weights .+ B.weights

    nvDT = promote_type(avDT,bvDT)

    return judiExtendedSource{nvDT}("judiRHS",m,n,A.info, wavelet, weights)
end

# -(judiExtendedSource,judiExtendedSource)
function -(A::judiExtendedSource{avDT}, B::judiExtendedSource{bvDT}) where {avDT, bvDT}

        # Error checking
        size(A) == size(B) || throw(judiExtendedSourceException("Shape mismatch: A:$(size(A)), v: $(size(B))"))
        compareInfo(A.info, B.info) == true || throw(judiExtendedSourceException("info mismatch"))

        # Size
        m = A.info.n * sum(A.info.nt)
        n = 1

        # wavelet and weights
        A.wavelet == B.wavelet ? nothing : throw(judiExtendedSourceException("Can only add two extended
                                                                            sources with same wavelet"))
        wavelet = A.wavelet
        weights = A.weights .- B.weights

        nvDT = promote_type(avDT,bvDT)

        return judiExtendedSource{nvDT}("judiRHS",m,n,A.info, wavelet, weights)
    end

function subsample(a::judiExtendedSource{avDT},srcnum) where avDT
    info = Info(a.info.n, length(srcnum), a.info.nt[srcnum])
    return judiExtendedSource(info,a.wavelet[srcnum],a.weights[srcnum];vDT=avDT)
end

getindex(x::judiExtendedSource,a) = subsample(x,a)

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
using CSV
using Statistics
using DataFrames
using Printf

function process(input)
    try
        df = CSV.read(input, DataFrame)
        p = df."Power(mAh)"
        println("power: $p")
        m = Statistics.mean(p)
        d = Statistics.stdm(p, m)
        @printf("power: %.2f ± %.2f\n", m, d)
    catch e
    end
end

@time process("all.csv")

max_stars_repo_name: alsam/scriptics-samples

max_issues_repo_name: alsam/scriptics-samples

max_forks_repo_name: alsam/scriptics-samples


text:
using BLP
using Test
import Random

include("iteration.jl")
include("inversion.jl")


text:
# Binary heap (non-mutable)

#################################################
#
#   core implementation
#
#################################################

function _heap_bubble_up!{Comp,T}(comp::Comp, valtree::Array{T}, i::Int)
    i0::Int = i
    @inbounds v = valtree[i]

    while i > 1  # nd is not root
        p = i >> 1
        @inbounds vp = valtree[p]

        if compare(comp, v, vp)
            # move parent downward
            @inbounds valtree[i] = vp
            i = p
        else
            break
        end
    end

    if i != i0
        @inbounds valtree[i] = v
    end
end

function _heap_bubble_down!{Comp,T}(comp::Comp, valtree::Array{T}, i::Int)
    @inbounds v::T = valtree[i]
    swapped = true
    n = length(valtree)
    last_parent = n >> 1

    while swapped && i <= last_parent
        lc = i << 1
        if lc < n   # contains both left and right children
            rc = lc + 1
            @inbounds lv = valtree[lc]
            @inbounds rv = valtree[rc]
            if compare(comp, rv, lv)
                if compare(comp, rv, v)
                    @inbounds valtree[i] = rv
                    i = rc
                else
                    swapped = false
                end
            else
                if compare(comp, lv, v)
                    @inbounds valtree[i] = lv
                    i = lc
                else
                    swapped = false
                end
            end
        else        # contains only left child
            @inbounds lv = valtree[lc]
            if compare(comp, lv, v)
                @inbounds valtree[i] = lv
                i = lc
            else
                swapped = false
            end
        end
    end

    valtree[i] = v
end

function _binary_heap_pop!{Comp,T}(comp::Comp, valtree::Array{T})
    # extract root
    v = valtree[1]

    if length(valtree) == 1
        empty!(valtree)
    else
        valtree[1] = pop!(valtree)
        if length(valtree) > 1
            _heap_bubble_down!(comp, valtree, 1)
        end
    end
    v
end

function _make_binary_heap{Comp,T}(comp::Comp, ty::Type{T}, xs)
    n = length(xs)
    valtree = copy(xs)
    for i = 2 : n
        _heap_bubble_up!(comp, valtree, i)
    end
    valtree
end

#################################################
#
#   heap type and constructors
#
#################################################

type BinaryHeap{T,Comp} <: AbstractHeap{T}
    comparer::Comp
    valtree::Array{T}

    function BinaryHeap(comp::Comp)
        new(comp, Array(T,0))
    end

    function BinaryHeap(comp::Comp, xs)  # xs is an iterable collection of values
        valtree = _make_binary_heap(comp, T, xs)
        new(comp, valtree)
    end
end

function binary_minheap{T}(ty::Type{T})
    BinaryHeap{T,LessThan}(LessThan())
end

binary_maxheap{T}(ty::Type{T}) = BinaryHeap{T,GreaterThan}(GreaterThan())
binary_minheap{T}(xs::AbstractVector{T}) = BinaryHeap{T,LessThan}(LessThan(), xs)
binary_maxheap{T}(xs::AbstractVector{T}) = BinaryHeap{T,GreaterThan}(GreaterThan(), xs)

#################################################
#
#   interfaces
#
#################################################

length(h::BinaryHeap) = length(h.valtree)

isempty(h::BinaryHeap) = isempty(h.valtree)

function push!{T}(h::BinaryHeap{T}, v::T)
    valtree = h.valtree
    push!(valtree, v)
    _heap_bubble_up!(h.comparer, valtree, length(valtree))
end

top(h::BinaryHeap) = h.valtree[1]

pop!{T}(h::BinaryHeap{T}) = _binary_heap_pop!(h.comparer, h.valtree)


text:
using HPAT

@acc hpat function cumsum_test(file_name)
    df = DataSource(DataTable{:id=Int64,:x=Float64,:y=Float64}, HDF5, file_name)
    df1 = df[:id<100]
    return sum(df1[:x])
end

cumsum_test(ENV["HOME"]*"/tmp/udf.hdf5")

max_stars_repo_name: ehsantn/data-analytics-benchmarks

max_issues_repo_name: ehsantn/data-analytics-benchmarks

max_forks_repo_name: ehsantn/data-analytics-benchmarks


text:
using Base.Test
using AbstractTrees

function runtest_leaves(f, masks = false)
    root = IOTreeRoot()
    t = GenerativeIOTree(root, AbstractTrees.firststate(Leaves(root)), 1)
    while true
        show(STDOUT, Tree(t.tree))
        next!(t, f(t)) || break
    end
    show(STDOUT, Tree(t.tree))
    !masks ? collect(map(x->x.value, Leaves(t.tree))) :
            collect(map(x->x.mask=>x.value, Leaves(t.tree)))
end

# Simple splitting on the third byte
function parse_simple(io)
    a = read(io, UInt8)
    b = read(io, UInt8)
    c = read(io, UInt8)
    return (c & (UInt8(1) << 4)) == 0
end
@test runtest_leaves(parse_simple, true) == [bitwise_masks[5]=>false, .~(bitwise_masks[5])=>true]

function parse_multiple(io)
    a = read(io, UInt8)
    b = read(io, UInt8)
    c = read(io, UInt8)
    return c == 0x42
end
z = fill!(BitVector(256), false)
z[0x42+1] = 1
@test runtest_leaves(parse_multiple, true) == [ChildMask(.~(z))=>false, ChildMask(z)=>true]

# Simple splitting, but on an earlier byte
function parse_early(io)
    a = read(io, UInt8)
    b = read(io, UInt8)
    c = read(io, UInt8)
    return (c & (UInt8(1) << 4)) == 0 && (a & (UInt8(1) << 2)) == 0
end
runtest_leaves(parse_early, true)

# Just reading the first byte
function parse_first(io)
    a = read(io, UInt8)
    return false
end
runtest_leaves(parse_first, true)

# Just reading the first byte, but actually using it in the decision
function parse_first_decision(io)
    a = read(io, UInt8)
    return (a & (UInt8(1) << 4)) == 0 && (a & (UInt8(1) << 2)) == 0
end
runtest_leaves(parse_first_decision, true)