text:
"""
Reads in all the data stored in csv files
The general format for data is
    folder:
        gen.csv
        branch.csv
        bus.csv
        ..
        load.csv
Args:
    Path to folder with all the System data CSV files

Returns:
    Nested Data dictionary with key values as folder/file names and dataframes
    as values

"""
function read_csv_data(file_path::String, baseMVA::Float64)
    files = readdir(file_path)
    REGEX_DEVICE_TYPE = r"(.*?)\.csv"
    REGEX_IS_FOLDER = r"^[A-Za-z]+$"
    data =Dict{String,Any}()

    if length(files) == 0
        error("No files in the folder")
    else
        data["baseMVA"] = baseMVA
    end

    encountered_files = 0
    for d_file in files
        try
            if match(REGEX_IS_FOLDER, d_file) != nothing
                @info "Parsing csv files in $d_file ..."
                d_file_data = Dict{String,Any}()
                for file in readdir(joinpath(file_path,d_file))
                    if match(REGEX_DEVICE_TYPE, file) != nothing
                        @info "Parsing csv data in $file ..."
                        encountered_files += 1
                        fpath = joinpath(file_path,d_file,file)
                        raw_data = CSV.File(fpath) |> DataFrames.DataFrame
                        d_file_data[split(file,r"[.]")[1]] = raw_data
                    end
                end

                if length(d_file_data) > 0
                    data[d_file] = d_file_data
                    @info "Successfully parsed $d_file"
                end

            elseif match(REGEX_DEVICE_TYPE, d_file) != nothing
                @info "Parsing csv data in $d_file ..."
                encountered_files += 1
                fpath = joinpath(file_path,d_file)
                raw_data = CSV.File(fpath)|> DataFrames.DataFrame
                data[split(d_file,r"[.]")[1]] = raw_data
                @info "Successfully parsed $d_file"
            end
        catch ex
            @error "Error occurred while parsing $d_file" exception=ex
        end
    end
    if encountered_files == 0
        error("No csv files or folders in $file_path")
    end

    if "timeseries_pointers" in keys(data)
        @info "parsing timeseries data"
        tsp_raw = data["timeseries_pointers"]
        data["timeseries_data"] = Dict()
        if :Simulation in names(tsp_raw)
            for sim in unique(tsp_raw[:Simulation])
                data["timeseries_data"][String(sim)] = Dict()
            end
        end

        for r in eachrow(tsp_raw)
            fpath = joinpath(file_path,r[Symbol("Data File")])
            if isfile(fpath)
                # read data and insert into dict
                @info "parsing timeseries data in $fpath for $(r.Object)"
                raw_data = CSV.File(fpath) |> DataFrames.DataFrame |> read_datetime

                if length([c for c in names(raw_data) if String(c) == String(r.Object)]) == 1
                    raw_data = TimeSeries.TimeArray(raw_data[:DateTime],raw_data[Symbol(r.Object)])
                end

                if :Simulation in names(tsp_raw)
                    data["timeseries_data"][String(r.Simulation)][String(r.Object)] = raw_data
                else
                    data["timeseries_data"][String(r.Object)] = raw_data
                end
            else
                @warn "File referenced in timeseries_pointers.csv doesn't exist : $fpath"
            end
        end
    end

    return data
end

"""
Args:
    Dict with all the System data from CSV files ... see `read_csv_data()`'
Returns:
    A Power Systems Nested dictionary with keys as devices and values as data
    dictionary necessary to construct the device structs
    PS dictionary:
        "Bus" => Dict(bus_no => Dict("name" =>
                                     "number" => ... ) )
        "Generator" => Dict( "Thermal" => Dict( "name" =>
                                                "tech" => ...)
                             "Hydro"  => ..
                             "Renewable" => .. )
        "Branch" => ...
        "Load" => ...
        "LoadZones" => ...
        "BaseKV" => ..
        ...
"""
function csv2ps_dict(data::Dict{String,Any})
    ps_dict =Dict{String,Any}()

    if haskey(data,"baseMVA")
        ps_dict["baseMVA"] = data["baseMVA"]
    else
        @warn "Key error : key 'baseMVA' not found in PowerSystems dictionary, this will result in a ps_dict['baseMVA'] = 100.0"
        ps_dict["baseMVA"] = 100.0
    end

    if haskey(data,"bus")
        ps_dict["bus"] =  PowerSystems.bus_csv_parser(data["bus"])
        if :Area in names(data["bus"])
            ps_dict["loadzone"] =  PowerSystems.loadzone_csv_parser(data["bus"], ps_dict["bus"])
            ps_dict["load"] =  PowerSystems.load_csv_parser(data["bus"], ps_dict["bus"], ps_dict["loadzone"], ps_dict["baseMVA"], haskey(data,"load") ? data["load"] : nothing)
        else
            @warn "Missing Data : no 'Area' information for buses, cannot create loads based on areas"
            ps_dict["load"] =  PowerSystems.load_csv_parser(data["bus"], ps_dict["bus"], ps_dict["baseMVA"], haskey(data,"load") ? data["load"] : nothing)
        end
    else
        error("Key error : key 'bus' not found in PowerSystems dictionary, cannot construct any System Struct")
    end
    if haskey(data,"gen")
        ps_dict["gen"] =  PowerSystems.gen_csv_parser(data["gen"], ps_dict["bus"], ps_dict["baseMVA"])
    else
        @warn "Key error : key 'gen' not found in PowerSystems dictionary, this will result in an ps_dict['gen'] = nothing"
         ps_dict["gen"] = nothing
    end
    if haskey(data,"branch")
        ps_dict["branch"] =  PowerSystems.branch_csv_parser(data["branch"], ps_dict["bus"], ps_dict["baseMVA"])
    else
        @warn "Key error : key 'branch' not found in PowerSystems dictionary,
          \n This will result in an ps_dict['branch'] = nothing"
         ps_dict["branch"] = nothing
    end
    if haskey(data,"dc_branch")
        ps_dict["dcline"] =  PowerSystems.dc_branch_csv_parser(data["dc_branch"], ps_dict["bus"], ps_dict["baseMVA"])
    else
        @warn "Key error : key 'dc_branch' not found in PowerSystems dictionary,
          \n This will result in an ps_dict['dcline'] = nothing"
         ps_dict["dcline"] = nothing
    end
    if haskey(data,"reserves")
        ps_dict["services"] = PowerSystems.services_csv_parser(data["reserves"],data["gen"],data["bus"])
    else
        @warn "Key error : key 'reserves' not found in PowerSystems dictionary, this will result in a ps_dict['services'] =  nothing"
        ps_dict["services"] = nothing
    end

    if haskey(data,"timeseries_data")
        gen_map = _retrieve(ps_dict["gen"],"name",Dict())
        if haskey(data["timeseries_data"],"DAY_AHEAD")
            @info "adding DAY-AHEAD generator forcasats"
            _add_nested_dict!(ps_dict,["forecast","DA","gen"],_format_fcdict(data["timeseries_data"]["DAY_AHEAD"],gen_map))
            @info "adding DAY-AHEAD load forcasats"
            ps_dict["forecast"]["DA"]["load"] = data["timeseries_data"]["DAY_AHEAD"]["Load"]
        end
        if haskey(data["timeseries_data"],"REAL_TIME")
            @info "adding REAL-TIME generator forcasats"
            _add_nested_dict!(ps_dict,["forecast","RT","gen"],_format_fcdict(data["timeseries_data"]["REAL_TIME"],gen_map))
            @info "adding REAL-TIME load forcasats"
            ps_dict["forecast"]["RT"]["load"] = data["timeseries_data"]["REAL_TIME"]["Load"]
        end
    end

    return ps_dict
end

function _add_nested_dict!(d,keylist,value = Dict())
    if !haskey(d,keylist[1])
        d[keylist[1]] = length(keylist)==1 ? value : Dict{String,Any}()
    end
    if length(keylist) > 1
        _add_nested_dict!(d[keylist[1]],keylist[2:end],value)
    end
end

function _format_fcdict(fc,obj_map)
    paths = Dict()
    for (k,d) in fc
        if haskey(obj_map, k)
            _add_nested_dict!(paths,obj_map[k],d)
        end
    end
    return paths
end

"""
Args:
    Path to folder with all the System data CSV's files
Returns:
    A Power Systems Nested dictionary with keys as devices and values as data
    dictionary necessary to construct the device structs
    PS dictionary:
        "Bus" => Dict(bus_no => Dict("name" =>
                                     "number" => ... ) )
        "Generator" => Dict( "Thermal" => Dict( "name" =>
                                                "tech" => ...)
                             "Hydro"  => ..
                             "Renewable" => .. )
        "Branch" => ...
        "Load" => ...
        "LoadZones" => ...
        "BaseKV" => ..
        ...
"""
function csv2ps_dict(file_path::String, baseMVA::Float64)
    data =  read_csv_data(file_path, baseMVA)
    ps_dict = csv2ps_dict(data)
    return ps_dict
end

###########
#Bus data parser
###########

"""
Args:
    A DataFrame with the same column names as in RTS_GMLC bus.csv file

    "Bus ID" "Bus Name" "BaseKV" "Bus Type" "MW Load" "MVAR Load" "V Mag" "V
    Angle" "MW Shunt G" "MVAR Shunt B" "Area"

Returns:
    A Nested Dictionary with keys as Bus number and values as bus data
    dictionary with same keys as the device struct
"""
function bus_csv_parser(bus_raw::DataFrames.DataFrame,colnames = nothing)

    if colnames isa Nothing
        need_cols = ["Bus ID", "Bus Name", "BaseKV", "Bus Type", "V Mag", "V Angle"]
        tbl_cols = string.(names(bus_raw))
        colnames = Dict(zip(need_cols,[findall(tbl_cols.==c)[1] for c in need_cols]))
    end

    Buses_dict = Dict{Int64,Any}()
    for i in 1:DataFrames.nrow(bus_raw)
        Buses_dict[bus_raw[i,1]] = Dict{String,Any}("number" =>bus_raw[i,colnames["Bus ID"]] ,
                                                "name" => bus_raw[i,colnames["Bus Name"]],
                                                "bustype" => bus_raw[i,colnames["Bus Type"]],
                                                "angle" => bus_raw[i,colnames["V Angle"]],
                                                "voltage" => bus_raw[i,colnames["V Mag"]],
                                                "voltagelimits" => (min=0.95,max=1.05),
                                                "basevoltage" => bus_raw[i,colnames["BaseKV"]]
                                                )
    end
    return Buses_dict
end

function _get_value_or_nothing(value::Union{Real, String})::Union{Float64, Nothing}
    if value == "NA"
        return nothing
    end

    return Float64(value)
end

###########
#Generator data parser
###########

"""
Args:
    A DataFrame with the same column names as in RTS_GMLC gen.csv file
    Parsed Bus PowerSystems dictionary
Returns:
    A Nested Dictionary with keys as generator types/names and values as
    generator data dictionary with same keys as the device struct
"""
function gen_csv_parser(gen_raw::DataFrames.DataFrame, Buses::Dict{Int64,Any}, baseMVA::Float64, colnames = nothing)
    Generators_dict = Dict{String,Any}()
    Generators_dict["Thermal"] = Dict{String,Any}()
    Generators_dict["Hydro"] = Dict{String,Any}()
    Generators_dict["Renewable"] = Dict{String,Any}()
    Generators_dict["Renewable"]["PV"]= Dict{String,Any}()
    Generators_dict["Renewable"]["RTPV"]= Dict{String,Any}()
    Generators_dict["Renewable"]["WIND"]= Dict{String,Any}()
    Generators_dict["Storage"] = Dict{String,Any}()

    if colnames isa Nothing
        need_cols = ["GEN UID", "Bus ID", "Fuel", "Unit Type",
                    "MW Inj", "MVAR Inj", "PMax MW", "PMin MW", "QMax MVAR", "QMin MVAR",
                    "Min Down Time Hr", "Min Up Time Hr", "Ramp Rate MW/Min",
                    "Start Heat Cold MBTU", "Non Fuel Start Cost \$", "Non Fuel Shutdown Cost \$", "Fuel Price \$/MMBTU",
                    "Output_pct_0", "Output_pct_1", "Output_pct_2", "Output_pct_3", "Output_pct_4", "HR_avg_0", "HR_incr_1", "HR_incr_2", "HR_incr_3", "HR_incr_4",
                    "Base MVA"]
        tbl_cols = string.(names(gen_raw))
        colnames = Dict(zip(need_cols,[findall(tbl_cols.==c)[1] for c in need_cols]))
    end

    cost_colnames = []
    for i in 0:length([n for n in names(gen_raw) if occursin("Output_pct_",String(n))])-1
        hr = [n for n in names(gen_raw) if !isa(match(Regex("HR_.*_$i"),String(n)),Nothing)][1]
        mw = Symbol("Output_pct_$i")
        push!(cost_colnames, (hr,mw))
    end

    pu_cols = ["PMin MW", "PMax MW", "MVAR Inj", "MW Inj", "QMin MVAR", "QMax MVAR", "Ramp Rate MW/Min"]
    [gen_raw[colnames[c]] = gen_raw[colnames[c]]./baseMVA for c in pu_cols] # P.U. conversion

    for gen in 1:DataFrames.nrow(gen_raw)
        pmax = _get_value_or_nothing(gen_raw[gen,colnames["PMax MW"]])

        if gen_raw[gen,colnames["Fuel"]] in ["Oil","Coal","NG","Nuclear"]

            fuel_cost = gen_raw[gen,colnames["Fuel Price \$/MMBTU"]]./1000

            var_cost = [(_get_value_or_nothing(gen_raw[gen,cn[1]]), _get_value_or_nothing(gen_raw[gen,cn[2]])) for cn in cost_colnames]
            var_cost = [(c[1]*c[2]*fuel_cost*baseMVA, c[2]).*pmax for c in var_cost if !in(nothing,c)]
            var_cost[2:end] = [(var_cost[i-1][1]+var_cost[i][1], var_cost[i][2]) for i in 2:length(var_cost)]

            bus_id =[Buses[i] for i in keys(Buses) if Buses[i]["number"] == gen_raw[gen,colnames["Bus ID"]]]

            Generators_dict["Thermal"][gen_raw[gen,colnames["GEN UID"]]] = Dict{String,Any}("name" => gen_raw[gen,colnames["GEN UID"]],
                                            "available" => true,
                                            "bus" => make_bus(bus_id[1]),
                                            "tech" => Dict{String,Any}("activepower" => gen_raw[gen,colnames["MW Inj"]],
                                                                        "activepowerlimits" => (min=_get_value_or_nothing(gen_raw[gen, colnames["PMin MW"]]), max=pmax),
                                                                        "reactivepower" => _get_value_or_nothing(gen_raw[gen, colnames["MVAR Inj"]]),
                                                                        "reactivepowerlimits" => (min=_get_value_or_nothing(gen_raw[gen, colnames["QMin MVAR"]]), max=_get_value_or_nothing(gen_raw[gen, colnames["QMax MVAR"]])),
                                                                        "ramplimits" => (up=gen_raw[gen,colnames["Ramp Rate MW/Min"]], down=gen_raw[gen, colnames["Ramp Rate MW/Min"]]),
                                                                        "timelimits" => (up=gen_raw[gen,colnames["Min Up Time Hr"]], down=gen_raw[gen, colnames["Min Down Time Hr"]])),
                                            "econ" => Dict{String,Any}("capacity" => pmax,
                                                                        "variablecost" => var_cost,
                                                                        "fixedcost" => 0.0,
                                                                        "startupcost" => gen_raw[gen,colnames["Start Heat Cold MBTU"]]*fuel_cost*1000,
                                                                        "shutdncost" => 0.0,
                                                                        "annualcapacityfactor" => nothing)
                                            )

        elseif gen_raw[gen,colnames["Fuel"]] in ["Hydro"]
            bus_id =[Buses[i] for i in keys(Buses) if Buses[i]["number"] == gen_raw[gen,colnames["Bus ID"]]]
            Generators_dict["Hydro"][gen_raw[gen,colnames["GEN UID"]]] = Dict{String,Any}("name" => gen_raw[gen,colnames["GEN UID"]],
                                            "available" => true, # change from staus to available
                                            "bus" => make_bus(bus_id[1]),
                                            "tech" => Dict{String,Any}( "installedcapacity" => pmax,
                                                                        "activepower" => gen_raw[gen, colnames["MW Inj"]],
                                                                        "activepowerlimits" => (min=_get_value_or_nothing(gen_raw[gen, colnames["PMin MW"]]), max=pmax),
                                                                        "reactivepower" => gen_raw[gen, colnames["MVAR Inj"]],
                                                                        "reactivepowerlimits" => (min=_get_value_or_nothing(gen_raw[gen, colnames["QMin MVAR"]]), max=_get_value_or_nothing(gen_raw[gen, colnames["QMax MVAR"]])),
                                                                        "ramplimits" => (up=gen_raw[gen, colnames["Ramp Rate MW/Min"]], down=gen_raw[gen, colnames["Ramp Rate MW/Min"]]),
                                                                        "timelimits" => (up=gen_raw[gen, colnames["Min Down Time Hr"]], down=gen_raw[gen, colnames["Min Down Time Hr"]])),
                                            "econ" => Dict{String,Any}("curtailcost" => 0.0,
                                                                        "interruptioncost" => nothing),
                                            "scalingfactor" => TimeSeries.TimeArray(collect(Dates.DateTime(Dates.today()):Dates.Hour(1):Dates.DateTime(Dates.today()+Dates.Day(1))), ones(25)) # TODO: connect scaling factors
                                            )

        elseif gen_raw[gen,colnames["Fuel"]] in ["Solar","Wind"]
            bus_id =[Buses[i] for i in keys(Buses) if Buses[i]["number"] == gen_raw[gen,colnames["Bus ID"]]]
            if gen_raw[gen,colnames["Unit Type"]] == "PV"
                Generators_dict["Renewable"]["PV"][gen_raw[gen,colnames["GEN UID"]]] = Dict{String, Any}("name" => gen_raw[gen,colnames["GEN UID"]],
                                                "available" => true, # change from staus to available
                                                "bus" => make_bus(bus_id[1]),
                                                "tech" => Dict{String, Any}("installedcapacity" => pmax,
                                                                            "reactivepowerlimits" => (min=_get_value_or_nothing(gen_raw[gen, colnames["QMin MVAR"]]), max=_get_value_or_nothing(gen_raw[gen, colnames["QMax MVAR"]])),
                                                                            "powerfactor" => 1),
                                                "econ" => Dict{String, Any}("curtailcost" => 0.0,
                                                                            "interruptioncost" => nothing),
                                                "scalingfactor" => TimeSeries.TimeArray(collect(Dates.DateTime(Dates.today()):Dates.Hour(1):Dates.DateTime(Dates.today()+Dates.Day(1))), ones(25)) # TODO: connect scaling factors
                                                )
            elseif gen_raw[gen,colnames["Unit Type"]] == "RTPV"
                Generators_dict["Renewable"]["RTPV"][gen_raw[gen,colnames["GEN UID"]]] = Dict{String,Any}("name" => gen_raw[gen,colnames["GEN UID"]],
                                                "available" => true, # change from staus to available
                                                "bus" => make_bus(bus_id[1]),
                                                "tech" => Dict{String, Any}("installedcapacity" => pmax,
                                                                            "reactivepowerlimits" => (min=_get_value_or_nothing(gen_raw[gen, colnames["QMin MVAR"]]), max=_get_value_or_nothing(gen_raw[gen, colnames["QMax MVAR"]])),
                                                                            "powerfactor" => 1),
                                                "econ" => Dict{String, Any}("curtailcost" => 0.0,
                                                                            "interruptioncost" => nothing),
                                                "scalingfactor" => TimeSeries.TimeArray(collect(Dates.DateTime(Dates.today()):Dates.Hour(1):Dates.DateTime(Dates.today()+Dates.Day(1))), ones(25)) # TODO: Connect scaling factors
                                                )
            elseif gen_raw[gen,colnames["Unit Type"]] == "WIND"
                Generators_dict["Renewable"]["WIND"][gen_raw[gen, colnames["GEN UID"]]] = Dict{String, Any}("name" => gen_raw[gen, colnames["GEN UID"]],
                                                "available" => true, # change from staus to available
                                                "bus" => make_bus(bus_id[1]),
                                                "tech" => Dict{String, Any}("installedcapacity" => pmax,
                                                                            "reactivepowerlimits" => (min=_get_value_or_nothing(gen_raw[gen, colnames["QMin MVAR"]]), max=_get_value_or_nothing(gen_raw[gen, colnames["QMax MVAR"]])),
                                                                            "powerfactor" => 1),
                                                "econ" => Dict{String, Any}("curtailcost" => 0.0,
                                                                            "interruptioncost" => nothing),
                                                "scalingfactor" => TimeSeries.TimeArray(collect(Dates.DateTime(Dates.today()):Dates.Hour(1):Dates.DateTime(Dates.today()+Dates.Day(1))), ones(25)) # TODO: connect scaling factors
                                                )
            end
        elseif gen_raw[gen,colnames["Fuel"]] in ["Storage"]
            bus_id =[Buses[i] for i in keys(Buses) if Buses[i]["number"] == gen_raw[gen,colnames["Bus ID"]]]
            Generators_dict["Storage"][gen_raw[gen,colnames["GEN UID"]]] = Dict{String,Any}("name" => gen_raw[gen,colnames["GEN UID"]],
                                            "available" => true, # change from staus to available
                                            "bus" => make_bus(bus_id[1]),
                                            "energy" => 0.0,
                                            "capacity" => (min=_get_value_or_nothing(gen_raw[gen,colnames["PMin MW"]]),max=pmax),
                                            "activepower" => gen_raw[gen,colnames["MW Inj"]],
                                            "inputactivepowerlimits" => (min=0.0,max=pmax),
                                            "outputactivepowerlimits" => (min=0.0,max=pmax),
                                            "efficiency" => (in= 0.0, out = 0.0),
                                            "reactivepower" => gen_raw[gen,colnames["MVAR Inj"]],
                                            "reactivepowerlimits" => (min = 0.0, max = 0.0),
                                            )
        end
    end
    return Generators_dict
end

###########
#Branch data parser
###########

"""
Args:
    A DataFrame with the same column names as in RTS_GMLC branch.csv file
    Parsed Bus PowerSystems dictionary
Returns:
    A Nested Dictionary with keys as branch types/names and values as
    line/transformer data dictionary with same keys as the device struct
"""
function branch_csv_parser(branch_raw::DataFrames.DataFrame, Buses::Dict, baseMVA::Float64, colnames=nothing)

    if colnames isa Nothing
        need_cols = ["From Bus", "To Bus", "Tr Ratio", "Cont Rating","UID","R","X","B"]
        tbl_cols = string.(names(branch_raw))
        colnames = Dict(zip(need_cols,[findall(tbl_cols.==c)[1] for c in need_cols]))
    end

    pu_cols = ["Cont Rating"]
    [branch_raw[colnames[c]] = branch_raw[colnames[c]]./baseMVA for c in pu_cols] # P.U. conversion

    Branches_dict = Dict{String,Any}()
    Branches_dict["Transformers"] = Dict{String,Any}()
    Branches_dict["Lines"] = Dict{String,Any}()
    for i in 1:DataFrames.nrow(branch_raw)
        bus_f = [Buses[f] for f in keys(Buses) if Buses[f]["number"] == branch_raw[i,colnames["From Bus"]]]
        bus_t = [Buses[t] for t in keys(Buses) if Buses[t]["number"] == branch_raw[i,colnames["To Bus"]]]
        if branch_raw[i,Symbol("Tr Ratio")] > 0.0
            Branches_dict["Transformers"][branch_raw[i,colnames["UID"]]] = Dict{String,Any}("name" => branch_raw[i,colnames["UID"]],
                                                        "available" => true,
                                                        "connectionpoints" => (from=make_bus(bus_f[1]),to=make_bus(bus_t[1])),
                                                        "r" => branch_raw[i,colnames["R"]],
                                                        "x" => branch_raw[i,colnames["X"]],
                                                        "primaryshunt" => branch_raw[i,colnames["B"]] ,  #TODO: add field in CSV
                                                        "alpha" => (branch_raw[i,colnames["B"]]/2) - (branch_raw[i,colnames["B"]]/2), #TODO: Phase-Shifting Transformer angle
                                                        "tap" => branch_raw[i,colnames["Tr Ratio"]],
                                                        "rate" => branch_raw[i,colnames["Cont Rating"]],
                                                        )
        else
            Branches_dict["Lines"][branch_raw[i,:UID]] = Dict{String,Any}("name" => branch_raw[i,:UID],
                                                        "available" => true,
                                                        "connectionpoints" => (from=make_bus(bus_f[1]),to=make_bus(bus_t[1])),
                                                        "r" => branch_raw[i,colnames["R"]],
                                                        "x" => branch_raw[i,colnames["X"]],
                                                        "b" => (from=(branch_raw[i,colnames["B"]]/2),to=(branch_raw[i,colnames["B"]]/2)),
                                                        "rate" =>  branch_raw[i,colnames["Cont Rating"]],
                                                        "anglelimits" => (min=-60.0,max =60.0) #TODO: add field in CSV
                                                        )

        end
    end
    return Branches_dict
end

###########
#DC Branch data parser
###########

"""
Args:
    A DataFrame with the same column names as in RTS_GMLC dc_branch.csv file
    Parsed Bus PowerSystems dictionary
Returns:
    A Nested Dictionary with keys as dc_branch types/names and values as
    dc_branch data dictionary with same keys as the device struct
"""
function dc_branch_csv_parser(dc_branch_raw::DataFrames.DataFrame, Buses::Dict, baseMVA::Float64, colnames=nothing)

    if colnames isa Nothing
        need_cols = ["UID","From Bus", "To Bus","From X Commutating", "From Tap Min", "From Tap Max","From Min Firing Angle","From Max Firing Angle", "To X Commutating", "To Tap Min", "To Tap Max","To Min Firing Angle","To Max Firing Angle", "MW Load"]
        tbl_cols = string.(names(dc_branch_raw))
        colnames = Dict(zip(need_cols,[findall(tbl_cols.==c)[1] for c in need_cols if c in tbl_cols]))
    end

    pu_cols = ["MW Load"]
    [dc_branch_raw[colnames[c]] = dc_branch_raw[colnames[c]]./baseMVA for c in pu_cols] # P.U. conversion

    DCBranches_dict = Dict{String,Any}()
    DCBranches_dict["HVDCLine"] = Dict{String,Any}()
    DCBranches_dict["VSCDCLine"] = Dict{String,Any}()

    for i in 1:DataFrames.nrow(dc_branch_raw)
        bus_f = [Buses[f] for f in keys(Buses) if Buses[f]["number"] == dc_branch_raw[i,colnames["From Bus"]]]
        bus_t = [Buses[t] for t in keys(Buses) if Buses[t]["number"] == dc_branch_raw[i,colnames["To Bus"]]]
        if (dc_branch_raw[i,colnames["From Max Firing Angle"]] + dc_branch_raw[i,colnames["To Max Firing Angle"]])/2 != 0.0 #TODO: Replace this with the correct conditional to create VSCDC or HVDC lines
            DCBranches_dict["VSCDCLine"][dc_branch_raw[i,colnames["UID"]]] = Dict{String,Any}("name" => dc_branch_raw[i,colnames["UID"]],
                                                        "available" => true,
                                                        "connectionpoints" => (from=make_bus(bus_f[1]),to=make_bus(bus_t[1])),
                                                        "rectifier_taplimits" => (min=dc_branch_raw[i,colnames["From Tap Min"]],max=dc_branch_raw[i,colnames["From Tap Max"]]),
                                                        "rectifier_xrc" => dc_branch_raw[i,colnames["From X Commutating"]], #TODO: What is this?
                                                        "rectifier_firingangle" => (min=dc_branch_raw[i,colnames["From Min Firing Angle"]],max=dc_branch_raw[i,colnames["From Max Firing Angle"]]),
                                                        "inverter_taplimits" => (min=dc_branch_raw[i,colnames["To Tap Min"]],max=dc_branch_raw[i,colnames["To Tap Max"]]),
                                                        "inverter_xrc" => dc_branch_raw[i,colnames["To X Commutating"]],  #TODO: What is this?
                                                        "inverter_firingangle" => (min=dc_branch_raw[i,colnames["To Min Firing Angle"]],max=dc_branch_raw[i,colnames["To Max Firing Angle"]])
                                                        )
        else
            DCBranches_dict["HVDCLine"][dc_branch_raw[i,colnames["UID"]]] = Dict{String,Any}("name" => dc_branch_raw[i,colnames["UID"]],
                                                        "available" => true,
                                                        "connectionpoints" => (from=make_bus(bus_f[1]),to=make_bus(bus_t[1])),
                                                        "activepowerlimits_from" => (min=-1*dc_branch_raw[i,colnames["MW Load"]], max=dc_branch_raw[i,colnames["MW Load"]]), #TODO: is there a better way to calculate this?
                                                        "activepowerlimits_to" => (min=-1*dc_branch_raw[i,colnames["Rating"]], max=dc_branch_raw[i,colnames["MW Load"]]), #TODO: is there a better way to calculate this?
                                                        "reactivepowerlimits_from" => (min=-1*dc_branch_raw[i,colnames["MW Load"]], max=dc_branch_raw[i,colnames["MW Load"]]), #TODO: is there a better way to calculate this?
                                                        "reactivepowerlimits_to" => (min=-1*dc_branch_raw[i,colnames["MW Load"]], max=dc_branch_raw[i,colnames["MW Load"]]), #TODO: is there a better way to calculate this?
                                                        "loss" => 0.0, #TODO: Can we infer this from the other data?
                                                        )

        end
    end
    return DCBranches_dict
end

###########
#Load data parser
###########

"""
Args:
    A DataFrame with the same column names as in RTS_GMLC bus.csv file
    Parsed Bus entry of PowerSystems dictionary
    Parsed LoadZone entry of PowerSystems dictionary
Optional Args:
    DataFrame of LoadZone timeseries data
    Dict of bus column names
    Dict of load LoadZone timeseries column names
Returns:
    A Nested Dictionary with keys as load names and values as load data
    dictionary with same keys as the device struct
"""
function load_csv_parser(bus_raw::DataFrames.DataFrame, Buses::Dict, LoadZone::Dict, baseMVA::Float64, load_raw=nothing,bus_colnames=nothing,load_colnames=nothing)
    Loads_dict = Dict{String,Any}()
    load_zone = nothing

    if bus_colnames isa Nothing
        need_cols = ["MW Load", "MVAR Load", "Bus ID"]
        tbl_cols = string.(names(bus_raw))
        bus_colnames = Dict(zip(need_cols,[findall(tbl_cols.==c)[1] for c in need_cols]))
    end

    if !isa(load_raw,Nothing)
        load_raw = read_datetime(load_raw)
        if load_colnames isa Nothing
            need_cols = ["DateTime",]
            tbl_cols = string.(names(load_raw))
            load_colnames = Dict(zip(need_cols,[findall(tbl_cols.==c)[1] for c in need_cols]))
        end
    end

    pu_cols = ["MW Load", "MVAR Load"]
    [bus_raw[bus_colnames[c]] = bus_raw[bus_colnames[c]]./baseMVA for c in pu_cols] # P.U. conversion

    for (k_b,b) in Buses
        for (k_l,l)  in LoadZone
            bus_numbers = [b.number for b in l["buses"] ]
            if b["number"] in bus_numbers
                load_zone = k_l
            end
        end
        p = [bus_raw[n,bus_colnames["MW Load"]] for n in 1:DataFrames.nrow(bus_raw) if bus_raw[n,bus_colnames["Bus ID"]] == b["number"]]
        q = [bus_raw[m,bus_colnames["MVAR Load"]] for m in 1:DataFrames.nrow(bus_raw) if bus_raw[m,bus_colnames["Bus ID"]] == b["number"]]
        ts = isa(load_raw,Nothing) ? nothing : TimeSeries.TimeArray(load_raw[load_colnames["DateTime"]],load_raw[load_zone]*(p[1]/LoadZone[load_zone]["maxactivepower"]))
        Loads_dict[b["name"]] = Dict{String,Any}("name" => b["name"],
                                            "available" => true,
                                            "bus" => make_bus(b),
                                            "maxactivepower" => p[1],
                                            "maxreactivepower" => q[1],
                                            "scalingfactor" => ts #TODO remove TS
                                            )
    end
    return Loads_dict
end

"""
Args:
    A DataFrame with the same column names as in RTS_GMLC bus.csv file
    Parsed Bus entry of PowerSystems dictionary
Optional Args:
    DataFrame of LoadZone timeseries data
    Dict of bus column names
    Dict of load LoadZone timeseries column names
Returns:
    A Nested Dictionary with keys as load names and values as load data
    dictionary with same keys as the device struct
"""
function load_csv_parser(bus_raw::DataFrames.DataFrame, Buses::Dict, baseMVA::Float64, load_raw=nothing, bus_colnames=nothing, load_colnames=nothing)
    Loads_dict = Dict{String,Any}()
    load_zone = nothing

    if bus_colnames isa Nothing
        need_cols = ["MW Load", "MVAR Load", "Bus ID"]
        tbl_cols = string.(names(bus_raw))
        bus_colnames = Dict(zip(need_cols,[findall(tbl_cols.==c)[1] for c in need_cols]))
    end

    if !isa(load_raw,nothing)
        load_raw = read_datetime(load_raw)
        if load_colnames isa Nothing
            need_cols = ["DateTime"]
            tbl_cols = string.(names(load_raw))
            load_colnames = Dict(zip(need_cols,[findall(tbl_cols.==c)[1] for c in need_cols]))
        end
    end

    pu_cols = ["MW Load", "MVAR Load"]
    [bus_raw[colnames[c]] = bus_raw[colnames[c]]./baseMVA for c in pu_cols] # P.U. conversion

    for (k_b,b) in Buses
        p = [bus_raw[n,bus_colnames["MW Load"]] for n in 1:DataFrames.nrow(bus_raw) if bus_raw[n,bus_colnames["Bus ID"]] == b["number"]]
        q = [bus_raw[m,bus_colnames["MVAR Load"]] for m in 1:DataFrames.nrow(bus_raw) if bus_raw[m,bus_colnames["Bus ID"]] == b["number"]]
        ts = isa(load_raw,nothing) ? nothing : TimeSeries.TimeArray(load_raw[load_colnames["DateTime"]],load_raw[b["number"]]/(p[1]*baseMVA))
        Loads_dict[b["name"]] = Dict{String,Any}("name" => b["name"],
                                            "available" => true,
                                            "bus" => make_bus(b),
                                            "maxactivepower" => p[1],
                                            "maxreactivepower" => q[1],
                                            "scalingfactor" => ts #TODO remove TS
                                            )
    end
    return Loads_dict
end

###########
#LoadZone data parser
###########

"""
Args:
    A DataFrame with the same column names as in RTS_GMLC bus.csv file
    Parsed Bus PowerSystems dictionary
Returns:
    A Nested Dictionary with keys as loadzone names and values as loadzone data
    dictionary with same keys as the device struct
"""
function loadzone_csv_parser(bus_raw::DataFrames.DataFrame, Buses::Dict, colnames=nothing)
    if colnames isa Nothing
        need_cols = ["MW Load", "MVAR Load", "Bus ID", "Area"]
        tbl_cols = string.(names(bus_raw))
        colnames = Dict(zip(need_cols,[findall(tbl_cols.==c)[1] for c in need_cols]))
    end

    LoadZone_dict = Dict{Int64,Any}()
    lbs = zip(unique(bus_raw[colnames["Area"]]),[sum(bus_raw[colnames["Area"]].==a) for a in unique(bus_raw[colnames["Area"]])])
    for (zone,count) in lbs
        b_numbers = [bus_raw[b,colnames["Bus ID"]] for b in 1:DataFrames.nrow(bus_raw) if bus_raw[b,colnames["Area"]] == zone ]
        buses = [make_bus(Buses[i]) for i in keys(Buses) if Buses[i]["number"] in  b_numbers]
        activepower = [bus_raw[b,colnames["MW Load"]] for b in 1:DataFrames.nrow(bus_raw) if bus_raw[b,colnames["Area"]] == zone]
        reactivepower = [bus_raw[b,colnames["MVAR Load"]] for b in 1:DataFrames.nrow(bus_raw) if bus_raw[b,colnames["Area"]] == zone]
        LoadZone_dict[zone] = Dict{String,Any}("number" => zone,
                                                        "name" => zone ,
                                                        "buses" => buses,
                                                        "maxactivepower" => sum(activepower),
                                                        "maxreactivepower" => sum(reactivepower)
                                                        )
    end
    return LoadZone_dict
end

###########
#reserves data parser
###########

"""
Args:
    A DataFrame with the same column names as in RTS_GMLC reserves.csv file
    A DataFrame with the same column names as in RTS_GMLC gen.csv file
    A DataFrame with the same column names as in RTS_GMLC bus.csv file
Returns:
    A Nested Dictionary with keys as loadzone names and values as loadzone data
    dictionary with same keys as the device struct
"""
function services_csv_parser(rsv_raw::DataFrames.DataFrame,gen_raw::DataFrames.DataFrame,bus_raw::DataFrames.DataFrame,colnames=nothing, gen_colnames = nothing,bus_colnames = nothing)
    if colnames isa Nothing
        need_cols = ["Reserve Product", "Timeframe (sec)", "Eligible Regions", "Eligible Gen Categories"]
        tbl_cols = string.(names(rsv_raw))
        colnames = Dict(zip(need_cols,[findall(tbl_cols.==c)[1] for c in need_cols]))
    end

    if gen_colnames isa Nothing
        need_cols = ["GEN UID", "Bus ID", "Category"]
        tbl_cols = string.(names(gen_raw))
        gen_colnames = Dict(zip(need_cols,[findall(tbl_cols.==c)[1] for c in need_cols]))
    end

    if bus_colnames isa Nothing
        need_cols = ["Bus ID", "Area"]
        tbl_cols = string.(names(bus_raw))
        bus_colnames = Dict(zip(need_cols,[findall(tbl_cols.==c)[1] for c in need_cols]))
    end

    services_dict = Dict{String,Any}()
    for r in 1:DataFrames.nrow(rsv_raw)
        gen_cats = strip(rsv_raw[r,colnames["Eligible Gen Categories"]],[i for i in "()"]) |> (x->split(x, ","))
        regions = strip(rsv_raw[r,colnames["Eligible Regions"]],[i for i in "()"]) |> (x->split(x,","))
        contributingdevices = []
        [push!(contributingdevices,gen_raw[g,gen_colnames["GEN UID"]]) for g in 1:DataFrames.nrow(gen_raw) if
                    (gen_raw[g,gen_colnames["Category"]] in gen_cats) &
                    (string(bus_raw[bus_raw[bus_colnames["Bus ID"]] .== gen_raw[g,gen_colnames["Bus ID"]],bus_colnames["Area"]][1]) in regions)];

        services_dict[rsv_raw[r,colnames["Reserve Product"]]] = Dict{String,Any}("name" => rsv_raw[r,colnames["Reserve Product"]],
                                                                "contributingdevices" => contributingdevices,
                                                                "timeframe" => rsv_raw[r,colnames["Timeframe (sec)"]])
    end

    return services_dict
end

# Remove missing values form dataframes
#TODO : Remove "NA" Strings from the data created by CSV.read()
"""
Arg:
    Any DataFrame with Missing values / "NA" strings that are either created by
    readtable() or CSV.read()
Returns:
    DataFrame with missing values replaced by 0
"""
function remove_missing(df)
    for col in names(df)
        df[ismissing.(df[col]), col] = 0
    end
    return df
end


text:
using NeuralQuantum, Test
using NeuralQuantum: set_index!
num_types = [Float32, Float64]
atol_types  = [1e-5, 1e-8]

machines = Dict()

ma = (T, hi) -> RBMSplit(T, hi, 2)
machines["RBMSplit"] = ma

ma = (T, hi) -> RBM(T, hi, 2, NeuralQuantum.logℒ)
machines["RBM_softplus"] = ma

ma = (T, hi) -> RBM(T, hi, 2, NeuralQuantum.logℒ2)
machines["RBM_cosh"] = ma

ma = (T, hi) -> NDM(T, hi, 1, 2, NeuralQuantum.logℒ)
machines["NDM_softplus"] = ma

ma = (T, hi) -> NDM(T, hi, 1, 2, NeuralQuantum.logℒ2)
machines["NDM_cosh"] = ma

graph = HyperCube([N], periodic=true)
symm  = translational_symm_table(graph)
ma = (T, hi) -> NDMSymm(T, hi, 1, 2, symm, NeuralQuantum.logℒ2)
machines["NDMSymm_cosh"] = ma

ma = (T, hi) -> begin
    N = nsites(hi)
    ch = Chain(Dense(T, N, 2*N, af_softplus), sum_autobatch)
    return PureStateAnsatz(ch, hi)
end
machines["chain_pure_softplus"] = ma

ma = (T, hi) -> begin
    N = nsites(hi)
    ch = Chain(Dense(T, N, 2*N, af_logcosh), sum_autobatch)
    return PureStateAnsatz(ch, N)
end
machines["chain_pure_cosh"] = ma

ma = (T, hi) -> begin
    N = nsites(hi)
    ch = Chain(Dense(T, N, 2*N, af_logcosh), Dense(T, 2*N, N, af_logcosh), sum_autobatch)
    return PureStateAnsatz(ch, N)
end
machines["chain_pure_cosh_2"] = ma

ma = (T, hi) -> begin
    N = nsites(hi)
    ch = Chain(DenseSplit(T, N, 2*N, af_softplus), sum_autobatch)
    return MixedStateAnsatz(ch, N)
end
#machines["chain_mixed_softplus"] = ma

N = 4
T = Float32

@testset "test cached dispatch - values: $name" for name=keys(machines)
    hilb = HomogeneousFock(N, 2)

    net = machines[name](T,hilb)
    if net isa NeuralQuantum.MatrixNet
        hilb = SuperOpSpace(hilb)
    end
    cnet = cached(net)
    v = state(T, hilb, net)

    @test net(v) ≈ cnet(v)
    if !(v isa AbstractVector)
        @test net(v...) ≈ cnet(v...)
    end
end

@testset "test cached dispatch - allocating gradients: $name" for name=keys(machines)
    hilb = HomogeneousFock(N, 2)

    net = machines[name](T,hilb)

    if net isa NeuralQuantum.MatrixNet
        hilb = SuperOpSpace(hilb)
    end
    cnet = cached(net)
    v = state(T, hilb, net)

    @test ∇logψ(net, v) ≈ ∇logψ(cnet, v)
end

@testset "test cached dispatch - inplace gradients: $name" for name=keys(machines)
    hilb = HomogeneousFock(N, 2)

    net = machines[name](T,hilb)

    if net isa NeuralQuantum.MatrixNet
        hilb = SuperOpSpace(hilb)
    end
    cnet = cached(net)
    v = state(T, hilb, net)

    g1 = grad_cache(net)
    g2 = grad_cache(net)

    v1, gg1 = logψ_and_∇logψ!(g1, net, v)
    v2, gg2 = logψ_and_∇logψ!(g2, cnet, v)
    @test v1 ≈ v2
    @test gg1 ≈ gg2
    @test g1 === gg1 && g2 === gg2

    v1, gg1 = logψ_and_∇logψ!(g1, net, v)
    v2, gg2 = logψ_and_∇logψ!(g2, cnet, v)
    @test v1 ≈ v2
    @test gg1 ≈ gg2
    @test g1 === gg1 && g2 === gg2
end

ext: jl


text:
# ==================================================== #
# ===                                              === #
# === Methods to create, resample sequence events. === #
# ===                                              === #
# ==================================================== #

"""Constructs empty sequence event."""
function SeqEvent(
        num_sequence_types::Int64,
        num_warp_values::Int64
    )

    SeqEvent(
        0,
        zeros(num_sequence_types, num_warp_values),  # summed potentials
        zeros(num_sequence_types, num_warp_values),  # summed precisions
        zeros(num_sequence_types, num_warp_values),  # summed log Z
        zeros(num_sequence_types, num_warp_values),  # posterior on sequence type.
        -1,   # sampled_type is ignored until spike_count > 0
        -1,   # sampled_warp is ignored until spike_count > 0
        0.0,  # sampled_timestamp is ignored until spike_count > 0
        0.0   # sampled_amplitude is ignored until spike_count > 0
    )
end

"""Resets event to be empty."""
function reset!(e::SeqEvent)
    e.spike_count = 0
    fill!(e.summed_potentials, 0)
    fill!(e.summed_precisions, 0)
    fill!(e.summed_logZ, 0)
    fill!(e.seq_type_posterior, 0)
end

# =============================================================== #
# ===                                                         === #
# === Methods to maintain an array of sequence event structs. === #
# ===                                                         === #
# =============================================================== #

# Initialize SeqEventList
function SeqEventList(
        num_sequence_types::Int64,
        num_warp_values::Int64
    )
    events = [SeqEvent(num_sequence_types, num_warp_values)]  # one empty event
    indices = Int64[]
    SeqEventList(events, indices)
end

# Makes SeqEventList a read-only array.
getindex(ev::SeqEventList, i::Int64) = ev.events[i]

# Length of the event list is the number of occupied clusters.
length(ev::SeqEventList) = length(ev.indices)

# Define interface for iteration. Enables "for event in sequences"
# to iterate over the occupied sequence events.
iterate(ev::SeqEventList) = (
    isempty(ev.indices) ? nothing : (ev.events[ev.indices[1]], 2)
)

iterate(ev::SeqEventList, i::Int64) = (
    (i > length(ev)) ? nothing : (ev.events[ev.indices[i]], i + 1)
)

"""
Finds either an empty SeqEvent struct, or creates a new
one. Returns the index of the new cluster.
"""
function add_event!(ev::SeqEventList)

    # Check if any indices are skipped. If so, use the smallest skipped
    # integer as the index for the new event.
    i = 1
    for j in ev.indices

        # We have j == ev.indices[i].

        # If (indices[i] != i) then events[i] is empty.
        if i != j
            insert!(ev.indices, i, i)  # mark events[i] as occupied
            return i                   # i is the index for the new event.
        end

        # Increment to check if (i + 1) is empty.
        i += 1
    end

    # If we reached here without returning, then indices is a vector
    # [1, 2, ..., K] without any skipped integers. So we'll use K + 1
    # as the new integer index.
    push!(ev.indices, length(ev.indices) + 1)

    # Create a new SeqEvent object if necessary.
    if length(ev.events) < length(ev.indices)
        R, W = size(ev[1].summed_potentials)
        push!(ev.events, SeqEvent(R, W))
    end

    # Return index of the empty SeqEvent struct.
    return ev.indices[end]
end

"""
Marks a SeqEvent struct as empty and resets its sufficient
statistics. This does not delete the SeqEvent.
"""
function remove_event!(ev::SeqEventList, index::Int64)
    reset!(ev.events[index])
    deleteat!(ev.indices, searchsorted(ev.indices, index))
end

"""
Recompute sufficient statistics for all sequence events.
"""
function recompute!(
        model::SeqModel,
        spikes::Vector{Spike},
        assignments::AbstractVector{Int64}
    )

    # Grab sequence event list.
    ev = model.sequence_events

    # Reset all events to zero spikes.
    for k in ev.indices
        reset!(ev.events[k])
    end
    empty!(ev.indices)

    # Add spikes back to their previously assigned event.
    for (s, k) in zip(spikes, assignments)

        # Skip background spikes.
        (k < 0) && continue

        # Check that event k exists.
        while k > length(ev.events)
            push!(
                ev.events,
                SeqEvent(
                    num_sequence_types(model),
                    num_warp_values(model)
                )
            )
        end

        # Add datapoint to event k.
        add_datapoint!(model, s, k, recompute_posterior=false)

        # Make sure that event k is marked as non-empty.
        j = searchsortedfirst(ev.indices, k)
        if (j > length(ev.indices)) || (ev.indices[j] != k)
            insert!(ev.indices, j, k)
        end
    end

    # Set the posterior, since we didn't do so when adding datapoints
    for k in ev.indices
        set_posterior!(model, k)
    end
end

"""
Returns a vector of EventSummaryInfo structs
summarizing latent events and throwing away
sufficient statistics.
"""
function event_list_summary(model::SeqModel)
    infos = EventSummaryInfo[]
    ev = model.sequence_events
    warp_vals = model.priors.warp_values

    for ind in ev.indices
        push!(
            infos,
            EventSummaryInfo(
                ind,
                ev[ind].sampled_timestamp,
                ev[ind].sampled_type,
                warp_vals[ev[ind].sampled_warp],
                ev[ind].sampled_amplitude
            )
        )
    end

    return infos
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

## types ##

typealias Dims{N} NTuple{N,Int}
typealias DimsInteger{N} NTuple{N,Integer}
typealias Indices{N} NTuple{N,AbstractUnitRange}

"""
    <:(T1, T2)

Subtype operator, equivalent to `issubtype(T1,T2)`.

```jldoctest
julia> Float64 <: AbstractFloat
true

julia> Vector{Int} <: AbstractArray
true

julia> Matrix{Float64} <: Matrix{AbstractFloat}
false
```
"""
const (<:) = issubtype

"""
    supertype(T::DataType)

Return the supertype of DataType `T`.

```jldoctest
julia> supertype(Int32)
Signed
```
"""
function supertype(T::DataType)
    @_pure_meta
    T.super
end

function supertype(T::UnionAll)
    @_pure_meta
    UnionAll(T.var, supertype(T.body))
end

## generic comparison ##

==(x, y) = x === y

"""
    isequal(x, y)

Similar to `==`, except treats all floating-point `NaN` values as equal to each other, and
treats `-0.0` as unequal to `0.0`. The default implementation of `isequal` calls `==`, so if
you have a type that doesn't have these floating-point subtleties then you probably only
need to define `==`.

`isequal` is the comparison function used by hash tables (`Dict`). `isequal(x,y)` must imply
that `hash(x) == hash(y)`.

This typically means that if you define your own `==` function then you must define a
corresponding `hash` (and vice versa). Collections typically implement `isequal` by calling
`isequal` recursively on all contents.

Scalar types generally do not need to implement `isequal` separate from `==`, unless they
represent floating-point numbers amenable to a more efficient implementation than that
provided as a generic fallback (based on `isnan`, `signbit`, and `==`).

```jldoctest
julia> isequal([1., NaN], [1., NaN])
true

julia> [1., NaN] == [1., NaN]
false

julia> 0.0 == -0.0
true

julia> isequal(0.0, -0.0)
false
```
"""
isequal(x, y) = x == y

isequal(x::AbstractFloat, y::AbstractFloat) = (isnan(x) & isnan(y)) | (signbit(x) == signbit(y)) & (x == y)
isequal(x::Real,          y::AbstractFloat) = (isnan(x) & isnan(y)) | (signbit(x) == signbit(y)) & (x == y)
isequal(x::AbstractFloat, y::Real         ) = (isnan(x) & isnan(y)) | (signbit(x) == signbit(y)) & (x == y)

isless(x::AbstractFloat, y::AbstractFloat) = (!isnan(x) & isnan(y)) | (signbit(x) & !signbit(y)) | (x < y)
isless(x::Real,          y::AbstractFloat) = (!isnan(x) & isnan(y)) | (signbit(x) & !signbit(y)) | (x < y)
isless(x::AbstractFloat, y::Real         ) = (!isnan(x) & isnan(y)) | (signbit(x) & !signbit(y)) | (x < y)

function ==(T::Type, S::Type)
    @_pure_meta
    typeseq(T, S)
end
function !=(T::Type, S::Type)
    @_pure_meta
    !(T == S)
end
==(T::TypeVar, S::Type) = false
==(T::Type, S::TypeVar) = false

## comparison fallbacks ##

"""
    !=(x, y)
    ≠(x,y)

Not-equals comparison operator. Always gives the opposite answer as `==`. New types should
generally not implement this, and rely on the fallback definition `!=(x,y) = !(x==y)` instead.

```jldoctest
julia> 3 != 2
true

julia> "foo" ≠ "foo"
false
```
"""
!=(x, y) = !(x == y)
const ≠ = !=

"""
    ===(x,y) -> Bool
    ≡(x,y) -> Bool

Determine whether `x` and `y` are identical, in the sense that no program could distinguish
them. Compares mutable objects by address in memory, and compares immutable objects (such as
numbers) by contents at the bit level. This function is sometimes called `egal`.

```jldoctest
julia> a = [1, 2]; b = [1, 2];

julia> a == b
true

julia> a === b
false

julia> a === a
true
```
"""
===
const ≡ = ===

"""
    !==(x, y)
    ≢(x,y)

Equivalent to `!(x === y)`.

```jldoctest
julia> a = [1, 2]; b = [1, 2];

julia> a ≢ b
true

julia> a ≢ a
false
```
"""
!==(x, y) = !(x === y)
const ≢ = !==

"""
    <(x, y)

Less-than comparison operator. New numeric types should implement this function for two
arguments of the new type. Because of the behavior of floating-point NaN values, `<`
implements a partial order. Types with a canonical partial order should implement `<`, and
types with a canonical total order should implement `isless`.

```jldoctest
julia> 'a' < 'b'
true

julia> "abc" < "abd"
true

julia> 5 < 3
false
```
"""
<(x, y) = isless(x, y)

"""
    >(x, y)

Greater-than comparison operator. Generally, new types should implement `<` instead of this
function, and rely on the fallback definition `>(x, y) = y < x`.

```jldoctest
julia> 'a' > 'b'
false

julia> 7 > 3 > 1
true

julia> "abc" > "abd"
false

julia> 5 > 3
true
```
"""
>(x, y) = y < x

"""
    <=(x, y)
    ≤(x,y)

Less-than-or-equals comparison operator.

```jldoctest
julia> 'a' <= 'b'
true

julia> 7 ≤ 7 ≤ 9
true

julia> "abc" ≤ "abc"
true

julia> 5 <= 3
false
```
"""
<=(x, y) = !(y < x)
const ≤ = <=

"""
    >=(x, y)
    ≥(x,y)

Greater-than-or-equals comparison operator.

```jldoctest
julia> 'a' >= 'b'
false

julia> 7 ≥ 7 ≥ 3
true

julia> "abc" ≥ "abc"
true

julia> 5 >= 3
true
```
"""
>=(x, y) = (y <= x)
const ≥ = >=

# this definition allows Number types to implement < instead of isless,
# which is more idiomatic:
isless(x::Real, y::Real) = x<y
lexcmp(x::Real, y::Real) = isless(x,y) ? -1 : ifelse(isless(y,x), 1, 0)

"""
    ifelse(condition::Bool, x, y)

Return `x` if `condition` is `true`, otherwise return `y`. This differs from `?` or `if` in
that it is an ordinary function, so all the arguments are evaluated first. In some cases,
using `ifelse` instead of an `if` statement can eliminate the branch in generated code and
provide higher performance in tight loops.

```jldoctest
julia> ifelse(1 > 2, 1, 2)
2
```
"""
ifelse(c::Bool, x, y) = select_value(c, x, y)

"""
    cmp(x,y)

Return -1, 0, or 1 depending on whether `x` is less than, equal to, or greater than `y`,
respectively. Uses the total order implemented by `isless`. For floating-point numbers, uses `<`
but throws an error for unordered arguments.

```jldoctest
julia> cmp(1, 2)
-1

julia> cmp(2, 1)
1

julia> cmp(2+im, 3-im)
ERROR: MethodError: no method matching isless(::Complex{Int64}, ::Complex{Int64})
[...]
```
"""
cmp(x, y) = isless(x, y) ? -1 : ifelse(isless(y, x), 1, 0)

"""
    lexcmp(x, y)

Compare `x` and `y` lexicographically and return -1, 0, or 1 depending on whether `x` is
less than, equal to, or greater than `y`, respectively. This function should be defined for
lexicographically comparable types, and `lexless` will call `lexcmp` by default.

```jldoctest
julia> lexcmp("abc", "abd")
-1

julia> lexcmp("abc", "abc")
0
```
"""
lexcmp(x, y) = cmp(x, y)

"""
    lexless(x, y)

Determine whether `x` is lexicographically less than `y`.

```jldoctest
julia> lexless("abc", "abd")
true
```
"""
lexless(x, y) = lexcmp(x,y) < 0

# cmp returns -1, 0, +1 indicating ordering
cmp(x::Integer, y::Integer) = ifelse(isless(x, y), -1, ifelse(isless(y, x), 1, 0))

"""
    max(x, y, ...)

Return the maximum of the arguments. See also the [`maximum`](@ref) function
to take the maximum element from a collection.

```jldoctest
julia> max(2, 5, 1)
5
```
"""
max(x, y) = ifelse(y < x, x, y)

"""
    min(x, y, ...)

Return the minimum of the arguments. See also the [`minimum`](@ref) function
to take the minimum element from a collection.

```jldoctest
julia> min(2, 5, 1)
1
```
"""
min(x,y) = ifelse(y < x, y, x)

"""
    minmax(x, y)

Return `(min(x,y), max(x,y))`. See also: [`extrema`](@ref) that returns `(minimum(x), maximum(x))`.

```jldoctest
julia> minmax('c','b')
('b', 'c')
```
"""
minmax(x,y) = y < x ? (y, x) : (x, y)

scalarmax(x,y) = max(x,y)
scalarmax(x::AbstractArray, y::AbstractArray) = throw(ArgumentError("ordering is not well-defined for arrays"))
scalarmax(x               , y::AbstractArray) = throw(ArgumentError("ordering is not well-defined for arrays"))
scalarmax(x::AbstractArray, y               ) = throw(ArgumentError("ordering is not well-defined for arrays"))

scalarmin(x,y) = min(x,y)
scalarmin(x::AbstractArray, y::AbstractArray) = throw(ArgumentError("ordering is not well-defined for arrays"))
scalarmin(x               , y::AbstractArray) = throw(ArgumentError("ordering is not well-defined for arrays"))
scalarmin(x::AbstractArray, y               ) = throw(ArgumentError("ordering is not well-defined for arrays"))

## definitions providing basic traits of arithmetic operators ##

"""
    identity(x)

The identity function. Returns its argument.

```jldoctest
julia> identity("Well, what did you expect?")
"Well, what did you expect?"
```
"""
identity(x) = x

+(x::Number) = x
*(x::Number) = x
(&)(x::Integer) = x
(|)(x::Integer) = x
xor(x::Integer) = x

const ⊻ = xor

# foldl for argument lists. expand recursively up to a point, then
# switch to a loop. this allows small cases like `a+b+c+d` to be inlined
# efficiently, without a major slowdown for `+(x...)` when `x` is big.
afoldl(op,a) = a
afoldl(op,a,b) = op(a,b)
afoldl(op,a,b,c...) = afoldl(op, op(a,b), c...)
function afoldl(op,a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,qs...)
    y = op(op(op(op(op(op(op(op(op(op(op(op(op(op(op(a,b),c),d),e),f),g),h),i),j),k),l),m),n),o),p)
    for x in qs; y = op(y,x); end
    y
end

for op in (:+, :*, :&, :|, :xor, :min, :max, :kron)
    @eval begin
        # note: these definitions must not cause a dispatch loop when +(a,b) is
        # not defined, and must only try to call 2-argument definitions, so
        # that defining +(a,b) is sufficient for full functionality.
        ($op)(a, b, c, xs...) = afoldl($op, ($op)(($op)(a,b),c), xs...)
        # a further concern is that it's easy for a type like (Int,Int...)
        # to match many definitions, so we need to keep the number of
        # definitions down to avoid losing type information.
    end
end

"""
    \\(x, y)

Left division operator: multiplication of `y` by the inverse of `x` on the left. Gives
floating-point results for integer arguments.

```jldoctest
julia> 3 \\ 6
2.0

julia> inv(3) * 6
2.0

julia> A = [1 2; 3 4]; x = [5, 6];

julia> A \\ x
2-element Array{Float64,1}:
 -4.0
  4.5

julia> inv(A) * x
2-element Array{Float64,1}:
 -4.0
  4.5
```
"""
\(x,y) = (y'/x')'

# Core <<, >>, and >>> take either Int or UInt as second arg. Signed shift
# counts can shift in either direction, and are translated here to unsigned
# counts. Integer datatypes only need to implement the unsigned version.

"""
    <<(x, n)

Left bit shift operator, `x << n`. For `n >= 0`, the result is `x` shifted left
by `n` bits, filling with `0`s. This is equivalent to `x * 2^n`. For `n < 0`,
this is equivalent to `x >> -n`.

```jldoctest
julia> Int8(3) << 2
12

julia> bits(Int8(3))
"00000011"

julia> bits(Int8(12))
"00001100"
```
See also [`>>`](@ref), [`>>>`](@ref).
"""
function <<(x::Integer, c::Integer)
    typemin(Int) <= c <= typemax(Int) && return x << (c % Int)
    (x >= 0 || c >= 0) && return zero(x)
    oftype(x, -1)
end
<<(x::Integer, c::Unsigned) = c <= typemax(UInt) ? x << (c % UInt) : zero(x)
<<(x::Integer, c::Int) = c >= 0 ? x << unsigned(c) : x >> unsigned(-c)

"""
    >>(x, n)

Right bit shift operator, `x >> n`. For `n >= 0`, the result is `x` shifted
right by `n` bits, where `n >= 0`, filling with `0`s if `x >= 0`, `1`s if `x <
0`, preserving the sign of `x`. This is equivalent to `fld(x, 2^n)`. For `n <
0`, this is equivalent to `x << -n`.

```jldoctest
julia> Int8(13) >> 2
3

julia> bits(Int8(13))
"00001101"

julia> bits(Int8(3))
"00000011"

julia> Int8(-14) >> 2
-4

julia> bits(Int8(-14))
"11110010"

julia> bits(Int8(-4))
"11111100"
```
See also [`>>>`](@ref), [`<<`](@ref).
"""
function >>(x::Integer, c::Integer)
    typemin(Int) <= c <= typemax(Int) && return x >> (c % Int)
    (x >= 0 || c < 0) && return zero(x)
    oftype(x, -1)
end
>>(x::Integer, c::Unsigned) = c <= typemax(UInt) ? x >> (c % UInt) : zero(x)
>>(x::Integer, c::Int) = c >= 0 ? x >> unsigned(c) : x << unsigned(-c)

"""
    >>>(x, n)

Unsigned right bit shift operator, `x >>> n`. For `n >= 0`, the result is `x`
shifted right by `n` bits, where `n >= 0`, filling with `0`s. For `n < 0`, this
is equivalent to `x << -n`.

For `Unsigned` integer types, this is equivalent to [`>>`](@ref). For
`Signed` integer types, this is equivalent to `signed(unsigned(x) >> n)`.

```jldoctest
julia> Int8(-14) >>> 2
60

julia> bits(Int8(-14))
"11110010"

julia> bits(Int8(60))
"00111100"
```
`BigInt`s are treated as if having infinite size, so no filling is required and this
is equivalent to [`>>`](@ref).

See also [`>>`](@ref), [`<<`](@ref).
"""
>>>(x::Integer, c::Integer) =
    typemin(Int) <= c <= typemax(Int) ? x >>> (c % Int) : zero(x)
>>>(x::Integer, c::Unsigned) = c <= typemax(UInt) ? x >>> (c % UInt) : zero(x)
>>>(x::Integer, c::Int) = c >= 0 ? x >>> unsigned(c) : x << unsigned(-c)

# fallback div, fld, and cld implementations
# NOTE: C89 fmod() and x87 FPREM implicitly provide truncating float division,
# so it is used here as the basis of float div().
div{T<:Real}(x::T, y::T) = convert(T,round((x-rem(x,y))/y))

"""
    fld(x, y)

Largest integer less than or equal to `x/y`.

```jldoctest
julia> fld(7.3,5.5)
1.0
```
"""
fld{T<:Real}(x::T, y::T) = convert(T,round((x-mod(x,y))/y))

"""
    cld(x, y)

Smallest integer larger than or equal to `x/y`.
```jldoctest
julia> cld(5.5,2.2)
3.0
```
"""
cld{T<:Real}(x::T, y::T) = convert(T,round((x-modCeil(x,y))/y))
#rem{T<:Real}(x::T, y::T) = convert(T,x-y*trunc(x/y))
#mod{T<:Real}(x::T, y::T) = convert(T,x-y*floor(x/y))
modCeil{T<:Real}(x::T, y::T) = convert(T,x-y*ceil(x/y))

# operator alias

"""
    rem(x, y)
    %(x, y)

Remainder from Euclidean division, returning a value of the same sign as `x`, and smaller in
magnitude than `y`. This value is always exact.

```jldoctest
julia> x = 15; y = 4;

julia> x % y
3

julia> x == div(x, y) * y + rem(x, y)
true
```
"""
rem
const % = rem

"""
    div(x, y)
    ÷(x, y)

The quotient from Euclidean division. Computes `x/y`, truncated to an integer.

```jldoctest
julia> 9 ÷ 4
2

julia> -5 ÷ 3
-1
```
"""
div
const ÷ = div

"""
    mod1(x, y)

Modulus after flooring division, returning a value `r` such that `mod(r, y) == mod(x, y)`
in the range ``(0, y]`` for positive `y` and in the range ``[y,0)`` for negative `y`.

```jldoctest
julia> mod1(4, 2)
2

julia> mod1(4, 3)
1
```
"""
mod1{T<:Real}(x::T, y::T) = (m = mod(x, y); ifelse(m == 0, y, m))
# efficient version for integers
mod1{T<:Integer}(x::T, y::T) = mod(x + y - T(1), y) + T(1)

"""
    fld1(x, y)

Flooring division, returning a value consistent with `mod1(x,y)`

See also: [`mod1`](@ref).

```jldoctest
julia> x = 15; y = 4;

julia> fld1(x, y)
4

julia> x == fld(x, y) * y + mod(x, y)
true

julia> x == (fld1(x, y) - 1) * y + mod1(x, y)
true
```
"""
fld1{T<:Real}(x::T, y::T) = (m=mod(x,y); fld(x-m,y))
# efficient version for integers
fld1{T<:Integer}(x::T, y::T) = fld(x+y-T(1),y)

"""
    fldmod1(x, y)

Return `(fld1(x,y), mod1(x,y))`.

See also: [`fld1`](@ref), [`mod1`](@ref).
"""
fldmod1{T<:Real}(x::T, y::T) = (fld1(x,y), mod1(x,y))
# efficient version for integers
fldmod1{T<:Integer}(x::T, y::T) = (fld1(x,y), mod1(x,y))

# transpose

"""
    ctranspose(A)

The conjugate transposition operator (`'`).

# Example

```jldoctest
julia> A =  [3+2im 9+2im; 8+7im  4+6im]
2×2 Array{Complex{Int64},2}:
 3+2im  9+2im
 8+7im  4+6im

julia> ctranspose(A)
2×2 Array{Complex{Int64},2}:
 3-2im  8-7im
 9-2im  4-6im
```
"""
ctranspose(x) = conj(transpose(x))
conj(x) = x

# transposed multiply

"""
    Ac_mul_B(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᴴ⋅B``.
"""
Ac_mul_B(a,b)  = ctranspose(a)*b

"""
    A_mul_Bc(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``A⋅Bᴴ``.
"""
A_mul_Bc(a,b)  = a*ctranspose(b)

"""
    Ac_mul_Bc(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᴴ Bᴴ``.
"""
Ac_mul_Bc(a,b) = ctranspose(a)*ctranspose(b)

"""
    At_mul_B(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᵀ⋅B``.
"""
At_mul_B(a,b)  = transpose(a)*b

"""
    A_mul_Bt(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``A⋅Bᵀ``.
"""
A_mul_Bt(a,b)  = a*transpose(b)

"""
    At_mul_Bt(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᵀ⋅Bᵀ``.
"""
At_mul_Bt(a,b) = transpose(a)*transpose(b)

# transposed divide

"""
    Ac_rdiv_B(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᴴ / B``.
"""
Ac_rdiv_B(a,b)  = ctranspose(a)/b

"""
    A_rdiv_Bc(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``A / Bᴴ``.
"""
A_rdiv_Bc(a,b)  = a/ctranspose(b)

"""
    Ac_rdiv_Bc(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᴴ / Bᴴ``.
"""
Ac_rdiv_Bc(a,b) = ctranspose(a)/ctranspose(b)

"""
    At_rdiv_B(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᵀ / B``.
"""
At_rdiv_B(a,b)  = transpose(a)/b

"""
    A_rdiv_Bt(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``A / Bᵀ``.
"""
A_rdiv_Bt(a,b)  = a/transpose(b)

"""
    At_rdiv_Bt(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᵀ / Bᵀ``.
"""
At_rdiv_Bt(a,b) = transpose(a)/transpose(b)

"""
    Ac_ldiv_B(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᴴ`` \\ ``B``.
"""
Ac_ldiv_B(a,b)  = ctranspose(a)\b

"""
    A_ldiv_Bc(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``A`` \\ ``Bᴴ``.
"""
A_ldiv_Bc(a,b)  = a\ctranspose(b)

"""
    Ac_ldiv_Bc(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᴴ`` \\ ``Bᴴ``.
"""
Ac_ldiv_Bc(a,b) = ctranspose(a)\ctranspose(b)

"""
    At_ldiv_B(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᵀ`` \\ ``B``.
"""
At_ldiv_B(a,b)  = transpose(a)\b

"""
    A_ldiv_Bt(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``A`` \\ ``Bᵀ``.
"""
A_ldiv_Bt(a,b)  = a\transpose(b)

"""
    At_ldiv_Bt(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᵀ`` \\ ``Bᵀ``.
"""
At_ldiv_Bt(a,b) = At_ldiv_B(a,transpose(b))

"""
    Ac_ldiv_Bt(A, B)

For matrices or vectors ``A`` and ``B``, calculates ``Aᴴ`` \\ ``Bᵀ``.
"""
Ac_ldiv_Bt(a,b) = Ac_ldiv_B(a,transpose(b))

widen{T<:Number}(x::T) = convert(widen(T), x)

"""
    eltype(type)

Determine the type of the elements generated by iterating a collection of the given `type`.
For associative collection types, this will be a `Pair{KeyType,ValType}`. The definition
`eltype(x) = eltype(typeof(x))` is provided for convenience so that instances can be passed
instead of types. However the form that accepts a type argument should be defined for new
types.

```jldoctest
julia> eltype(ones(Float32,2,2))
Float32

julia> eltype(ones(Int8,2,2))
Int8
```
"""
eltype(::Type) = Any
eltype(::Type{Any}) = Any
eltype(t::DataType) = eltype(supertype(t))
eltype(x) = eltype(typeof(x))

# function pipelining

"""
    |>(x, f)

Applies a function to the preceding argument. This allows for easy function chaining.

```jldoctest
julia> [1:5;] |> x->x.^2 |> sum |> inv
0.01818181818181818
```
"""
|>(x, f) = f(x)

# function composition

"""
    f ∘ g

Compose functions: i.e. `(f ∘ g)(args...)` means `f(g(args...))`. The `∘` symbol can be
entered in the Julia REPL (and most editors, appropriately configured) by typing `\\circ<tab>`.
Example:

```jldoctest
julia> map(uppercase∘hex, 250:255)
6-element Array{String,1}:
 "FA"
 "FB"
 "FC"
 "FD"
 "FE"
 "FF"
```
"""
∘(f, g) = (x...)->f(g(x...))

"""
    !f::Function

Predicate function negation: when the argument of `!` is a function, it returns a
function which computes the boolean negation of `f`. Example:

```jldoctest
julia> str = "∀ ε > 0, ∃ δ > 0: |x-y| < δ ⇒ |f(x)-f(y)| < ε"
"∀ ε > 0, ∃ δ > 0: |x-y| < δ ⇒ |f(x)-f(y)| < ε"

julia> filter(isalpha, str)
"εδxyδfxfyε"

julia> filter(!isalpha, str)
"∀  > 0, ∃  > 0: |-| <  ⇒ |()-()| < "
```
"""
!(f::Function) = (x...)->!f(x...)

# array shape rules

promote_shape(::Tuple{}, ::Tuple{}) = ()

function promote_shape(a::Tuple{Int,}, b::Tuple{Int,})
    if a[1] != b[1]
        throw(DimensionMismatch("dimensions must match"))
    end
    return a
end

function promote_shape(a::Tuple{Int,Int}, b::Tuple{Int,})
    if a[1] != b[1] || a[2] != 1
        throw(DimensionMismatch("dimensions must match"))
    end
    return a
end

promote_shape(a::Tuple{Int,}, b::Tuple{Int,Int}) = promote_shape(b, a)

function promote_shape(a::Tuple{Int, Int}, b::Tuple{Int, Int})
    if a[1] != b[1] || a[2] != b[2]
        throw(DimensionMismatch("dimensions must match"))
    end
    return a
end

"""
    promote_shape(s1, s2)

Check two array shapes for compatibility, allowing trailing singleton dimensions, and return
whichever shape has more dimensions.

```jldoctest
julia> a = ones(3,4,1,1,1);

julia> b = ones(3,4);

julia> promote_shape(a,b)
(Base.OneTo(3), Base.OneTo(4), Base.OneTo(1), Base.OneTo(1), Base.OneTo(1))

julia> promote_shape((2,3,1,4), (2, 3, 1, 4, 1))
(2, 3, 1, 4, 1)
```
"""
function promote_shape(a::Dims, b::Dims)
    if length(a) < length(b)
        return promote_shape(b, a)
    end
    for i=1:length(b)
        if a[i] != b[i]
            throw(DimensionMismatch("dimensions must match"))
        end
    end
    for i=length(b)+1:length(a)
        if a[i] != 1
            throw(DimensionMismatch("dimensions must match"))
        end
    end
    return a
end

function promote_shape(a::AbstractArray, b::AbstractArray)
    promote_shape(indices(a), indices(b))
end

function promote_shape(a::Indices, b::Indices)
    if length(a) < length(b)
        return promote_shape(b, a)
    end
    for i=1:length(b)
        if a[i] != b[i]
            throw(DimensionMismatch("dimensions must match"))
        end
    end
    for i=length(b)+1:length(a)
        if a[i] != 1:1
            throw(DimensionMismatch("dimensions must match"))
        end
    end
    return a
end

function throw_setindex_mismatch(X, I)
    if length(I) == 1
        throw(DimensionMismatch("tried to assign $(length(X)) elements to $(I[1]) destinations"))
    else
        throw(DimensionMismatch("tried to assign $(dims2string(size(X))) array to $(dims2string(I)) destination"))
    end
end

# check for valid sizes in A[I...] = X where X <: AbstractArray
# we want to allow dimensions that are equal up to permutation, but only
# for permutations that leave array elements in the same linear order.
# those are the permutations that preserve the order of the non-singleton
# dimensions.
function setindex_shape_check(X::AbstractArray, I::Integer...)
    li = ndims(X)
    lj = length(I)
    i = j = 1
    while true
        ii = length(indices(X,i))
        jj = I[j]
        if i == li || j == lj
            while i < li
                i += 1
                ii *= length(indices(X,i))
            end
            while j < lj
                j += 1
                jj *= I[j]
            end
            if ii != jj
                throw_setindex_mismatch(X, I)
            end
            return
        end
        if ii == jj
            i += 1
            j += 1
        elseif ii == 1
            i += 1
        elseif jj == 1
            j += 1
        else
            throw_setindex_mismatch(X, I)
        end
    end
end

setindex_shape_check(X::AbstractArray) =
    (_length(X)==1 || throw_setindex_mismatch(X,()))

setindex_shape_check(X::AbstractArray, i::Integer) =
    (_length(X)==i || throw_setindex_mismatch(X, (i,)))

setindex_shape_check{T}(X::AbstractArray{T,1}, i::Integer) =
    (_length(X)==i || throw_setindex_mismatch(X, (i,)))

setindex_shape_check{T}(X::AbstractArray{T,1}, i::Integer, j::Integer) =
    (_length(X)==i*j || throw_setindex_mismatch(X, (i,j)))

function setindex_shape_check{T}(X::AbstractArray{T,2}, i::Integer, j::Integer)
    if length(X) != i*j
        throw_setindex_mismatch(X, (i,j))
    end
    sx1 = length(indices(X,1))
    if !(i == 1 || i == sx1 || sx1 == 1)
        throw_setindex_mismatch(X, (i,j))
    end
end
setindex_shape_check(X, I...) = nothing # Non-arrays broadcast to all idxs

# convert to a supported index type (array or Int)
"""
    to_index(A, i)

Convert index `i` to an `Int` or array of indices to be used as an index into array `A`.

Custom array types may specialize `to_index(::CustomArray, i)` to provide
special indexing behaviors. Note that some index types (like `Colon`) require
more context in order to transform them into an array of indices; those get
converted in the more complicated `to_indices` function. By default, this
simply calls the generic `to_index(i)`. This must return either an `Int` or an
`AbstractArray` of scalar indices that are supported by `A`.
"""
to_index(A, i) = to_index(i)

"""
    to_index(i)

Convert index `i` to an `Int` or array of `Int`s to be used as an index for all arrays.

Custom index types may specialize `to_index(::CustomIndex)` to provide special
indexing behaviors. This must return either an `Int` or an `AbstractArray` of
`Int`s.
"""
to_index(i::Integer) = convert(Int,i)::Int
to_index(I::AbstractArray{Bool}) = LogicalIndex(I)
to_index(I::AbstractArray) = I
to_index{T<:Union{AbstractArray, Colon}}(I::AbstractArray{T}) = throw(ArgumentError("invalid index: $I"))
to_index(::Colon) = throw(ArgumentError("colons must be converted by to_indices(...)"))
to_index(i) = throw(ArgumentError("invalid index: $i"))

# The general to_indices is mostly defined in multidimensional.jl, but this
# definition is required for bootstrap:
"""
    to_indices(A, I::Tuple)

Convert the tuple `I` to a tuple of indices for use in indexing into array `A`.

The returned tuple must only contain either `Int`s or `AbstractArray`s of
scalar indices that are supported by array `A`. It will error upon encountering
a novel index type that it does not know how to process.

For simple index types, it defers to the unexported `Base.to_index(A, i)` to
process each index `i`. While this internal function is not intended to be
called directly, `Base.to_index` may be extended by custom array or index types
to provide custom indexing behaviors.

More complicated index types may require more context about the dimension into
which they index. To support those cases, `to_indices(A, I)` calls
`to_indices(A, indices(A), I)`, which then recursively walks through both the
given tuple of indices and the dimensional indices of `A` in tandem. As such,
not all index types are guaranteed to propagate to `Base.to_index`.
"""
to_indices(A, I::Tuple) = (@_inline_meta; to_indices(A, indices(A), I))
to_indices(A, inds, ::Tuple{}) = ()
to_indices(A, inds, I::Tuple{Any, Vararg{Any}}) =
    (@_inline_meta; (to_index(A, I[1]), to_indices(A, _maybetail(inds), tail(I))...))

_maybetail(::Tuple{}) = ()
_maybetail(t::Tuple) = tail(t)

"""
   Slice(indices)

Represent an AbstractUnitRange of indices as a vector of the indices themselves.

Upon calling `to_indices()`, Colons are converted to Slice objects to represent
the indices over which the Colon spans. Slice objects are themselves unit
ranges with the same indices as those they wrap. This means that indexing into
Slice objects with an integer always returns that exact integer, and they
iterate over all the wrapped indices, even supporting offset indices.
"""
immutable Slice{T<:AbstractUnitRange} <: AbstractUnitRange{Int}
    indices::T
end
indices(S::Slice) = (S.indices,)
unsafe_indices(S::Slice) = (S.indices,)
indices1(S::Slice) = S.indices
first(S::Slice) = first(S.indices)
last(S::Slice) = last(S.indices)
errmsg(A) = error("size not supported for arrays with indices $(indices(A)); see http://docs.julialang.org/en/latest/devdocs/offset-arrays/")
size(S::Slice) = first(S.indices) == 1 ? (length(S.indices),) : errmsg(S)
length(S::Slice) = first(S.indices) == 1 ? length(S.indices) : errmsg(S)
unsafe_length(S::Slice) = first(S.indices) == 1 ? unsafe_length(S.indices) : errmsg(S)
getindex(S::Slice, i::Int) = (@_inline_meta; @boundscheck checkbounds(S, i); i)
show(io::IO, r::Slice) = print(io, "Base.Slice(", r.indices, ")")
start(S::Slice) = start(S.indices)
next(S::Slice, s) = next(S.indices, s)
done(S::Slice, s) = done(S.indices, s)

# Addition/subtraction of ranges
for f in (:+, :-)
    @eval begin
        function $f(r1::OrdinalRange, r2::OrdinalRange)
            r1l = length(r1)
            (r1l == length(r2) ||
             throw(DimensionMismatch("argument dimensions must match")))
            range($f(first(r1),first(r2)), $f(step(r1),step(r2)), r1l)
        end

        function $f{T}(r1::LinSpace{T}, r2::LinSpace{T})
            len = r1.len
            (len == r2.len ||
             throw(DimensionMismatch("argument dimensions must match")))
            linspace(convert(T, $f(first(r1), first(r2))),
                     convert(T, $f(last(r1), last(r2))), len)
        end

        $f(r1::Union{StepRangeLen, OrdinalRange, LinSpace},
           r2::Union{StepRangeLen, OrdinalRange, LinSpace}) =
               $f(promote_noncircular(r1, r2)...)
    end
end

function +{T,S}(r1::StepRangeLen{T,S}, r2::StepRangeLen{T,S})
    len = length(r1)
    (len == length(r2) ||
        throw(DimensionMismatch("argument dimensions must match")))
    StepRangeLen(first(r1)+first(r2), step(r1)+step(r2), len)
end

-(r1::StepRangeLen, r2::StepRangeLen) = +(r1, -r2)

# Pair

immutable Pair{A,B}
    first::A
    second::B
end

const => = Pair

start(p::Pair) = 1
done(p::Pair, i) = i>2
next(p::Pair, i) = (getfield(p,i), i+1)
eltype{A,B}(p::Pair{A,B}) = Union{A,B}

indexed_next(p::Pair, i::Int, state) = (getfield(p,i), i+1)

hash(p::Pair, h::UInt) = hash(p.second, hash(p.first, h))

==(p::Pair, q::Pair) = (p.first==q.first) & (p.second==q.second)
isequal(p::Pair, q::Pair) = isequal(p.first,q.first) & isequal(p.second,q.second)

isless(p::Pair, q::Pair) = ifelse(!isequal(p.first,q.first), isless(p.first,q.first),
                                                             isless(p.second,q.second))
getindex(p::Pair,i::Int) = getfield(p,i)
getindex(p::Pair,i::Real) = getfield(p, convert(Int, i))
reverse{A,B}(p::Pair{A,B}) = Pair{B,A}(p.second, p.first)

endof(p::Pair) = 2
length(p::Pair) = 2
first(p::Pair) = p.first
last(p::Pair) = p.second

convert{A,B}(::Type{Pair{A,B}}, x::Pair{A,B}) = x
function convert{A,B}(::Type{Pair{A,B}}, x::Pair)
    Pair{A, B}(convert(A, x[1]), convert(B, x[2]))
end

promote_rule{A1, B1, A2, B2}(::Type{Pair{A1, B1}}, ::Type{Pair{A2, B2}}) =
    Pair{promote_type(A1, A2), promote_type(B1, B2)}

# some operators not defined yet
global //, >:, <|, hcat, hvcat, ⋅, ×, ∈, ∉, ∋, ∌, ⊆, ⊈, ⊊, ∩, ∪, √, ∛

this_module = current_module()
baremodule Operators

export
    !,
    !=,
    !==,
    ===,
    xor,
    %,
    ÷,
    &,
    *,
    +,
    -,
    /,
    //,
    <,
    <:,
    >:,
    <<,
    <=,
    ==,
    >,
    >=,
    ≥,
    ≤,
    ≠,
    >>,
    >>>,
    \,
    ^,
    |,
    |>,
    <|,
    ~,
    ⋅,
    ×,
    ∈,
    ∉,
    ∋,
    ∌,
    ⊆,
    ⊈,
    ⊊,
    ∩,
    ∪,
    √,
    ∛,
    ⊻,
    ∘,
    colon,
    hcat,
    vcat,
    hvcat,
    getindex,
    setindex!,
    transpose,
    ctranspose

import ..this_module: !, !=, xor, %, ÷, &, *, +, -,
    /, //, <, <:, <<, <=, ==, >, >=, >>, >>>,
    <|, |>, \, ^, |, ~, !==, ===, >:, colon, hcat, vcat, hvcat, getindex, setindex!,
    transpose, ctranspose,
    ≥, ≤, ≠, ⋅, ×, ∈, ∉, ∋, ∌, ⊆, ⊈, ⊊, ∩, ∪, √, ∛, ⊻, ∘

end