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.
# *********************************************************************************
"""
    DomesticHotWaterLoad

There are many ways in which a DomesticHotWaterLoad can be defined:
1. When using either `doe_reference_name` or `blended_doe_reference_names` in an `ElectricLoad` one only needs to provide the input key "DomesticHotWaterLoad" in the `Scenario` (JSON or Dict). In this case the values from DoE reference names from the `ElectricLoad` will be used to define the `DomesticHotWaterLoad`.
2. One can provide the `doe_reference_name` or `blended_doe_reference_names` directly in the `DomesticHotWaterLoad` key within the `Scenario`. These values can be combined with the `annual_mmbtu` or `monthly_mmbtu` inputs to scale the DoE reference profile(s).
3. One can provide the `fuel_loads_mmbtu_per_hour` value in the `DomesticHotWaterLoad` key within the `Scenario`.

```julia
function DomesticHotWaterLoad(;
    doe_reference_name::String = "",
    city::String = "",
    blended_doe_reference_names::Array{String, 1} = String[],
    blended_doe_reference_percents::Array{<:Real,1} = Real[],
    annual_mmbtu::Union{Real, Nothing} = nothing,
    monthly_mmbtu::Array{<:Real,1} = Real[],
    fuel_loads_mmbtu_per_hour::Array{<:Real,1} = Real[],
    time_steps_per_hour::Int = 1,
    latitude::Real=0.0,
    longitude::Real=0.0
)
```
"""
struct DomesticHotWaterLoad
    loads_kw::Array{Real, 1}

    function DomesticHotWaterLoad(;
        doe_reference_name::String = "",
        city::String = "",
        blended_doe_reference_names::Array{String, 1} = String[],
        blended_doe_reference_percents::Array{<:Real,1} = Real[],
        annual_mmbtu::Union{Real, Nothing} = nothing,
        monthly_mmbtu::Array{<:Real,1} = Real[],
        # addressable_load_fraction,  # TODO
        fuel_loads_mmbtu_per_hour::Array{<:Real,1} = Real[],  # TODO shouldn't this be per_time_step?
        time_steps_per_hour::Int = 1,
        latitude::Real=0.0,
        longitude::Real=0.0
    )
        if length(fuel_loads_mmbtu_per_hour) > 0

            if !(length(fuel_loads_mmbtu_per_hour) / time_steps_per_hour ≈ 8760)
                throw(@error "Provided domestic hot water load does not match the time_steps_per_hour.")
            end

            loads_kw = fuel_loads_mmbtu_per_hour .* (KWH_PER_MMBTU * EXISTING_BOILER_EFFICIENCY)

        elseif !isempty(doe_reference_name)
            loads_kw = BuiltInDomesticHotWaterLoad(city, doe_reference_name, latitude, longitude, 2017, annual_mmbtu, monthly_mmbtu)

        elseif length(blended_doe_reference_names) > 1 &&
            length(blended_doe_reference_names) == length(blended_doe_reference_percents)
            loads_kw = blend_and_scale_doe_profiles(BuiltInDomesticHotWaterLoad, latitude, longitude, 2017,
                                                    blended_doe_reference_names, blended_doe_reference_percents, city,
                                                    annual_mmbtu, monthly_mmbtu)
        else
            error("Cannot construct DomesticHotWaterLoad. You must provide either [fuel_loads_mmbtu_per_hour],
                [doe_reference_name, city],
                or [blended_doe_reference_names, blended_doe_reference_percents, city].")
        end

        if length(loads_kw) < 8760*time_steps_per_hour
            loads_kw = repeat(loads_kw, inner=Int(time_steps_per_hour / (length(loads_kw)/8760)))
            @info "Repeating domestic hot water loads in each hour to match the time_steps_per_hour."
        end

        new(
            loads_kw
        )
    end
end

"""
    function SpaceHeatingLoad(;
        doe_reference_name::String = "",
        city::String = "",
        blended_doe_reference_names::Array{String, 1} = String[],
        blended_doe_reference_percents::Array{<:Real,1} = Real[],
        annual_mmbtu::Union{Real, Nothing} = nothing,
        monthly_mmbtu::Array{<:Real,1} = Real[],
        fuel_loads_mmbtu_per_hour::Array{<:Real,1} = Real[]
    )

There are many ways to define a `SpaceHeatingLoad`:
1. a time-series via the `fuel_loads_mmbtu_per_hour`,
2. scaling a DoE Commercial Reference Building (CRB) profile or a blend of CRB profiles to either the `annual_mmbtu` or `monthly_mmbtu` values;
3. or using the `doe_reference_name` or `blended_doe_reference_names` from the `ElectricLoad`.

When using an `ElectricLoad` defined from a `doe_reference_name` or `blended_doe_reference_names`
one only needs to provide an empty Dict in the scenario JSON to add a `SpaceHeatingLoad` to a
`Scenario`, i.e.:
```json
...
"ElectricLoad": {"doe_reference_name": "MidriseApartment"},
"SpaceHeatingLoad" : {},
...
```
In this case the values provided for `doe_reference_name`, or  `blended_doe_reference_names` and
`blended_doe_reference_percents` are copied from the `ElectricLoad` to the `SpaceHeatingLoad`.
"""
struct SpaceHeatingLoad
    loads_kw::Array{Real, 1}

    function SpaceHeatingLoad(;
        doe_reference_name::String = "",
        city::String = "",
        blended_doe_reference_names::Array{String, 1} = String[],
        blended_doe_reference_percents::Array{<:Real,1} = Real[],
        annual_mmbtu::Union{Real, Nothing} = nothing,
        monthly_mmbtu::Array{<:Real,1} = Real[],
        # addressable_load_fraction,  # TODO
        fuel_loads_mmbtu_per_hour::Array{<:Real,1} = Real[],
        time_steps_per_hour::Int = 1,
        latitude::Real=0.0,
        longitude::Real=0.0
    )
        if length(fuel_loads_mmbtu_per_hour) > 0

            if !(length(fuel_loads_mmbtu_per_hour) / time_steps_per_hour ≈ 8760)
                throw(@error "Provided space heating load does not match the time_steps_per_hour.")
            end

            loads_kw = fuel_loads_mmbtu_per_hour .* (KWH_PER_MMBTU * EXISTING_BOILER_EFFICIENCY)

        elseif !isempty(doe_reference_name)
            loads_kw = BuiltInSpaceHeatingLoad(city, doe_reference_name, latitude, longitude, 2017, annual_mmbtu, monthly_mmbtu)

        elseif length(blended_doe_reference_names) > 1 &&
            length(blended_doe_reference_names) == length(blended_doe_reference_percents)
            loads_kw = blend_and_scale_doe_profiles(BuiltInSpaceHeatingLoad, latitude, longitude, 2017,
                                                    blended_doe_reference_names, blended_doe_reference_percents, city,
                                                    annual_mmbtu, monthly_mmbtu)
        else
            error("Cannot construct BuiltInSpaceHeatingLoad. You must provide either [fuel_loads_mmbtu_per_hour],
                [doe_reference_name, city],
                or [blended_doe_reference_names, blended_doe_reference_percents, city].")
        end

        if length(loads_kw) < 8760*time_steps_per_hour
            loads_kw = repeat(loads_kw, inner=Int(time_steps_per_hour / (length(loads_kw)/8760)))
            @info "Repeating space heating loads in each hour to match the time_steps_per_hour."
        end

        new(
            loads_kw
        )
    end
end

"""
function get_existing_chiller_default_cop(; existing_chiller_max_thermal_factor_on_peak_load=nothing,
                                            loads_kw=nothing,
                                            loads_kw_thermal=nothing)
This function returns the default value for ExistingChiller.cop based on:
    1. No information about load, returns average of lower and higher cop default values (`cop_unknown_thermal`)
    2. If the cooling electric `loads_kw` is known, we first guess the thermal load profile using `cop_unknown_thermal`,
        and then we use the default logic to determine the `existing_chiller_cop` based on the peak thermal load with a thermal factor multiplier.
    3. If the cooling thermal `loads_kw_thermal` is known, same as 2. but we don't have to guess the cop to convert electric to thermal load first.
"""
function get_existing_chiller_default_cop(; existing_chiller_max_thermal_factor_on_peak_load=nothing, loads_kw=nothing, loads_kw_thermal=nothing)
    cop_less_than_100_ton = 4.40
    cop_more_than_100_ton = 4.69
    cop_unknown_thermal = (cop_less_than_100_ton + cop_more_than_100_ton) / 2.0
    max_cooling_load_ton = nothing
    if !isnothing(loads_kw_thermal)
        max_cooling_load_ton = maximum(loads_kw_thermal) / KWH_THERMAL_PER_TONHOUR
    elseif !isnothing(loads_kw)
        max_cooling_load_ton = maximum(loads_kw) / KWH_THERMAL_PER_TONHOUR * cop_unknown_thermal
    end
    if isnothing(max_cooling_load_ton) || isnothing(existing_chiller_max_thermal_factor_on_peak_load)
        return cop_unknown_thermal
    elseif max_cooling_load_ton * existing_chiller_max_thermal_factor_on_peak_load < 100.0
        return cop_less_than_100_ton
    else
        return cop_more_than_100_ton
    end
end

"""
    function CoolingLoad(;
        doe_reference_name::String = "",
        city::String = "",
        blended_doe_reference_names::Array{String, 1} = String[],
        blended_doe_reference_percents::Array{<:Real,1} = Real[],
        annual_tonhour::Union{Real, Nothing} = nothing,
        monthly_tonhour::Array{<:Real,1} = Real[],
        thermal_loads_ton::Array{<:Real,1} = Real[],
        annual_fraction_of_electric_load::Union{Real, Nothing} = nothing,
        monthly_fractions_of_electric_load::Array{<:Real,1} = Real[],
        per_time_step_fractions_of_electric_load::Array{<:Real,1} = Real[],
        existing_chiller_cop::Real = nothing,
        existing_chiller_max_thermal_factor_on_peak_load::Real= nothing
    )

There are many ways to define a `CoolingLoad`:
1. a time-series via the `thermal_loads_ton`,
2. DoE Commercial Reference Building (CRB) profile or a blend of CRB profiles which uses the buildings' fraction of total electric for cooling profile applied to the `ElectricLoad`
3. scaling a DoE Commercial Reference Building (CRB) profile or a blend of CRB profiles using `annual_tonhour` or `monthly_tonhour`
4. the `annual_fraction_of_electric_load` or `monthly_fractions_of_electric_load` values which gets applied to the `ElectricLoad` to determine the cooling electric load;
5. or using the `doe_reference_name` or `blended_doe_reference_names` from the `ElectricLoad`.

The electric-based `loads_kw` of the `CoolingLoad` is a _subset_ of the total electric load `ElectricLoad`, so `CoolingLoad.loads_kw` for the BAU/conventional electric consumption of the `existing_chiller` is subtracted from the `ElectricLoad` for the non-cooling electric load balance constraint in the model.

When using an `ElectricLoad` defined from a `doe_reference_name` or `blended_doe_reference_names`
one only needs to provide an empty Dict in the scenario JSON to add a `CoolingLoad` to a
`Scenario`, i.e.:
```json
...
"ElectricLoad": {"doe_reference_name": "MidriseApartment"},
"CoolingLoad" : {},
...
```
In this case the values provided for `doe_reference_name`, or  `blended_doe_reference_names` and
`blended_doe_reference_percents` are copied from the `ElectricLoad` to the `CoolingLoad`.
"""

struct CoolingLoad
    loads_kw_thermal::Array{Real, 1}
    existing_chiller_cop::Real

    function CoolingLoad(;
        doe_reference_name::String = "",
        city::String = "",
        blended_doe_reference_names::Array{String, 1} = String[],
        blended_doe_reference_percents::Array{<:Real,1} = Real[],
        annual_tonhour::Union{Real, Nothing} = nothing,
        monthly_tonhour::Array{<:Real,1} = Real[],
        thermal_loads_ton::Array{<:Real,1} = Real[],
        annual_fraction_of_electric_load::Union{Real, Nothing} = nothing,
        monthly_fractions_of_electric_load::Array{<:Real,1} = Real[],
        per_time_step_fractions_of_electric_load::Array{<:Real,1} = Real[],
        site_electric_load_profile = Real[],
        time_steps_per_hour::Int = 1,
        latitude::Float64=0.0,
        longitude::Float64=0.0,
        existing_chiller_cop::Union{Real, Nothing} = nothing,
        existing_chiller_max_thermal_factor_on_peak_load::Union{Real, Nothing}= nothing
    )
        # determine the timeseries of loads_kw_thermal
        loads_kw_thermal = nothing
        loads_kw = nothing
        if length(thermal_loads_ton) > 0
            if !(length(thermal_loads_ton) / time_steps_per_hour ≈ 8760)
                throw(@error "Provided cooling load does not match the time_steps_per_hour.")
            end
            loads_kw_thermal = thermal_loads_ton .* KWH_THERMAL_PER_TONHOUR

        elseif !isempty(per_time_step_fractions_of_electric_load) && (length(site_electric_load_profile) / time_steps_per_hour ≈ 8760)
            if !(length(per_time_step_fractions_of_electric_load) / time_steps_per_hour ≈ 8760)
                @error "Provided cooling per_time_step_fractions_of_electric_load array does not match the time_steps_per_hour."
            end
            loads_kw = per_time_step_fractions_of_electric_load .* site_electric_load_profile

        elseif !isempty(monthly_fractions_of_electric_load) && (length(site_electric_load_profile) / time_steps_per_hour ≈ 8760)
            if !(length(monthly_fractions_of_electric_load) ≈ 12)
                @error "Provided cooling monthly_fractions_of_electric_load array does not have 12 values."
            end
            timeseries = collect(DateTime(2017,1,1) : Minute(60/time_steps_per_hour) :
                                 DateTime(2017,1,1) + Minute(8760*60 - 60/time_steps_per_hour))
            loads_kw = [monthly_fractions_of_electric_load[month(dt)] * site_electric_load_profile[ts] for (ts, dt)
                        in enumerate(timeseries)]

        elseif !isnothing(annual_fraction_of_electric_load) && (length(site_electric_load_profile) / time_steps_per_hour ≈ 8760)
            loads_kw = annual_fraction_of_electric_load * site_electric_load_profile

        elseif !isempty(doe_reference_name)
            if isnothing(annual_tonhour) && isempty(monthly_tonhour)
                loads_kw = get_default_fraction_of_total_electric(city, doe_reference_name,
                                            latitude, longitude, 2017) .* site_electric_load_profile
            else
                loads_kw = BuiltInCoolingLoad(city, doe_reference_name, latitude, longitude, 2017,
                                          annual_tonhour, monthly_tonhour, existing_chiller_cop)
            end

        elseif length(blended_doe_reference_names) > 1 &&
            length(blended_doe_reference_names) == length(blended_doe_reference_percents)
            if isnothing(annual_tonhour) && isempty(monthly_tonhour)
                loads_kw = zeros(Int(8760/time_steps_per_hour))
                for (i, building) in enumerate(blended_doe_reference_names)
                    default_fraction = get_default_fraction_of_total_electric(city, building, latitude, longitude, 2017)
                    modified_fraction = default_fraction * blended_doe_reference_percents[i]
                    if length(site_electric_load_profile) > 8784
                        modified_fraction = repeat(modified_fraction, inner=time_steps_per_hour / (length(site_electric_load_profile)/8760))
                        @info "Repeating cooling electric load in each hour to match the time_steps_per_hour."
                    end
                    loads_kw += site_electric_load_profile .* modified_fraction
                end
            else
                loads_kw = blend_and_scale_doe_profiles(BuiltInCoolingLoad, latitude, longitude, 2017,
                                                        blended_doe_reference_names,
                                                        blended_doe_reference_percents, city,
                                                        annual_tonhour, monthly_tonhour)
            end

        else
            error("Cannot construct BuiltInCoolingLoad. You must provide either [thermal_loads_ton],
                [doe_reference_name, city], [blended_doe_reference_names, blended_doe_reference_percents, city],
                or the site_electric_load_profile along with one of [per_time_step_fractions_of_electric_load, monthly_fractions_of_electric_load, annual_fraction_of_electric_load].")
        end

        # Now that either the cooling electric loads_kw or the cooling thermal loads_kw_thermal is known, update existing_chiller_cop if it was not input
        if isnothing(existing_chiller_cop)
            existing_chiller_cop = get_existing_chiller_default_cop(;existing_chiller_max_thermal_factor_on_peak_load=existing_chiller_max_thermal_factor_on_peak_load,
                                loads_kw=loads_kw, loads_kw_thermal=loads_kw_thermal)
        end

        if isnothing(loads_kw_thermal)  # have to convert electric loads_kw to thermal load
            if !isnothing(annual_tonhour) || !isempty(monthly_tonhour)
                # cop_unknown_thermal (4.55) was used to convert thermal to electric in BuiltInCoolingLoad, so need to use the same here to convert back
                #  in order to preserve the input tonhour amounts - however, the updated/actual existing_chiller_cop is still assigned based on user input or actual max ton load conditional defaults
                chiller_cop = get_existing_chiller_default_cop()
            else
                chiller_cop = existing_chiller_cop
            end
            loads_kw_thermal = chiller_cop * loads_kw
        end

        if length(loads_kw_thermal) < 8760*time_steps_per_hour
            loads_kw_thermal = repeat(loads_kw_thermal, inner=Int(time_steps_per_hour /
                               (length(loads_kw_thermal)/8760)))
            @info "Repeating cooling loads in each hour to match the time_steps_per_hour."
        end

        new(
            loads_kw_thermal,
            existing_chiller_cop
        )
    end
end

function get_default_fraction_of_total_electric(city, doe_reference_name, latitude, longitude, year)
    crb_total_elec_loads_kw = BuiltInElectricLoad(city, doe_reference_name, latitude, longitude, year)

    crb_cooling_elec_loads_kw = BuiltInCoolingLoad(city, doe_reference_name, latitude, longitude, year, nothing, Real[], nothing)

    default_fraction_of_total_electric_profile = crb_cooling_elec_loads_kw ./
                                                    max.(crb_total_elec_loads_kw, repeat([1.0E-6], length(crb_total_elec_loads_kw)))
    return default_fraction_of_total_electric_profile
end

function BuiltInDomesticHotWaterLoad(
    city::String,
    buildingtype::String,
    latitude::Real,
    longitude::Real,
    year::Int,
    annual_mmbtu::Union{<:Real, Nothing}=nothing,
    monthly_mmbtu::Vector{<:Real}=Real[]
    )
    dhw_annual_mmbtu = Dict(
        "Miami" => Dict(
            "FastFoodRest" => 53.47209411,
            "FullServiceRest" => 158.0518043,
            "Hospital" => 442.7295435,
            "LargeHotel" => 3713.248373,
            "LargeOffice" => 127.9412792,
            "MediumOffice" => 22.09603477,
            "MidriseApartment" => 158.0580017,
            "Outpatient" => 27.60091429,
            "PrimarySchool" => 105.3179165,
            "RetailStore" => 0.0,
            "SecondarySchool" => 250.1299246,
            "SmallHotel" => 242.232695,
            "SmallOffice" => 9.891779415,
            "StripMall" => 0.0,
            "Supermarket" => 17.94985187,
            "Warehouse" => 0.0,
            "FlatLoad" => 333.0450133
        ),
        "Houston" => Dict(
            "FastFoodRest" => 62.56835989,
            "FullServiceRest" => 188.292814,
            "Hospital" => 530.6352726,
            "LargeHotel" => 4685.666667,
            "LargeOffice" => 160.8917808,
            "MediumOffice" => 25.9266894,
            "MidriseApartment" => 199.2544784,
            "Outpatient" => 32.87691943,
            "PrimarySchool" => 128.0705362,
            "RetailStore" => 0.0,
            "SecondarySchool" => 314.7654465,
            "SmallHotel" => 290.6293673,
            "SmallOffice" => 10.27839347,
            "StripMall" => 0.0,
            "Supermarket" => 19.86608717,
            "Warehouse" => 0.0,
            "FlatLoad" => 415.6076756
        ),
        "Phoenix" => Dict(
            "FastFoodRest" => 57.34025418,
            "FullServiceRest" => 170.9086319,
            "Hospital" => 480.098265,
            "LargeHotel" => 4127.191046,
            "LargeOffice" => 141.8507451,
            "MediumOffice" => 23.71397275,
            "MidriseApartment" => 175.5949563,
            "Outpatient" => 29.83372212,
            "PrimarySchool" => 116.811664,
            "RetailStore" => 0.0,
            "SecondarySchool" => 285.2339344,
            "SmallHotel" => 262.8487714,
            "SmallOffice" => 10.05557471,
            "StripMall" => 0.0,
            "Supermarket" => 18.7637822,
            "Warehouse" => 0.0,
            "FlatLoad" => 368.7653325
        ),
        "Atlanta" => Dict(
            "FastFoodRest" => 71.33170579,
            "FullServiceRest" => 217.4332205,
            "Hospital" => 615.3498557,
            "LargeHotel" => 5622.340656,
            "LargeOffice" => 192.7164525,
            "MediumOffice" => 29.62182675,
            "MidriseApartment" => 238.9315749,
            "Outpatient" => 37.9759973,
            "PrimarySchool" => 148.8362119,
            "RetailStore" => 0.0,
            "SecondarySchool" => 372.2083434,
            "SmallHotel" => 337.2688069,
            "SmallOffice" => 10.65138846,
            "StripMall" => 0.0,
            "Supermarket" => 21.71038069,
            "Warehouse" => 0.0,
            "FlatLoad" => 494.7735263
        ),
        "LasVegas" => Dict(
            "FastFoodRest" => 63.63848459,
            "FullServiceRest" => 191.8494897,
            "Hospital" => 540.9697668,
            "LargeHotel" => 4800.331564,
            "LargeOffice" => 164.7154124,
            "MediumOffice" => 26.36796732,
            "MidriseApartment" => 204.1120165,
            "Outpatient" => 33.48190098,
            "PrimarySchool" => 131.9651451,
            "RetailStore" => 0.0,
            "SecondarySchool" => 327.441087,
            "SmallHotel" => 296.3578765,
            "SmallOffice" => 10.32392915,
            "StripMall" => 0.0,
            "Supermarket" => 20.08676069,
            "Warehouse" => 0.0,
            "FlatLoad" => 425.7275876
        ),
        "LosAngeles" => Dict(
            "FastFoodRest" => 69.63212501,
            "FullServiceRest" => 211.7827529,
            "Hospital" => 598.9350422,
            "LargeHotel" => 5440.174033,
            "LargeOffice" => 186.6199083,
            "MediumOffice" => 28.91483286,
            "MidriseApartment" => 231.215325,
            "Outpatient" => 37.00823296,
            "PrimarySchool" => 142.8059487,
            "RetailStore" => 0.0,
            "SecondarySchool" => 352.7467563,
            "SmallHotel" => 328.1935523,
            "SmallOffice" => 10.58011717,
            "StripMall" => 0.0,
            "Supermarket" => 21.35337379,
            "Warehouse" => 0.0,
            "FlatLoad" => 478.7476251
        ),
        "SanFrancisco" => Dict(
            "FastFoodRest" => 77.13092952,
            "FullServiceRest" => 236.7180594,
            "Hospital" => 671.40531,
            "LargeHotel" => 6241.842643,
            "LargeOffice" => 213.8445094,
            "MediumOffice" => 32.07909301,
            "MidriseApartment" => 265.1697301,
            "Outpatient" => 41.35500136,
            "PrimarySchool" => 160.4507431,
            "RetailStore" => 0.0,
            "SecondarySchool" => 401.395655,
            "SmallHotel" => 368.0979112,
            "SmallOffice" => 10.90004379,
            "StripMall" => 0.0,
            "Supermarket" => 22.9292287,
            "Warehouse" => 0.0,
            "FlatLoad" => 546.4574286
        ),
        "Baltimore" => Dict(
            "FastFoodRest" => 78.2191761,
            "FullServiceRest" => 240.338156,
            "Hospital" => 681.9322322,
            "LargeHotel" => 6358.710286,
            "LargeOffice" => 217.7306132,
            "MediumOffice" => 32.52815422,
            "MidriseApartment" => 270.1195541,
            "Outpatient" => 41.96148216,
            "PrimarySchool" => 165.3116185,
            "RetailStore" => 0.0,
            "SecondarySchool" => 417.9512972,
            "SmallHotel" => 373.906416,
            "SmallOffice" => 10.94554028,
            "StripMall" => 0.0,
            "Supermarket" => 23.15795696,
            "Warehouse" => 0.0,
            "FlatLoad" => 557.0507802
        ),
        "Albuquerque" => Dict(
            "FastFoodRest" => 76.9149868,
            "FullServiceRest" => 235.9992545,
            "Hospital" => 669.3128607,
            "LargeHotel" => 6219.08303,
            "LargeOffice" => 212.9944774,
            "MediumOffice" => 31.97726287,
            "MidriseApartment" => 264.2063457,
            "Outpatient" => 41.20639013,
            "PrimarySchool" => 162.1556119,
            "RetailStore" => 0.0,
            "SecondarySchool" => 409.1649863,
            "SmallHotel" => 366.9712928,
            "SmallOffice" => 10.88949351,
            "StripMall" => 0.0,
            "Supermarket" => 22.88525618,
            "Warehouse" => 0.0,
            "FlatLoad" => 545.235078
        ),
        "Seattle" => Dict(
            "FastFoodRest" => 81.80231236,
            "FullServiceRest" => 252.2609525,
            "Hospital" => 716.6111323,
            "LargeHotel" => 6741.736717,
            "LargeOffice" => 230.8057849,
            "MediumOffice" => 34.04746055,
            "MidriseApartment" => 286.3412104,
            "Outpatient" => 44.07342164,
            "PrimarySchool" => 172.0233322,
            "RetailStore" => 0.0,
            "SecondarySchool" => 434.0806311,
            "SmallHotel" => 392.968915,
            "SmallOffice" => 11.09863592,
            "StripMall" => 0.0,
            "Supermarket" => 23.91178737,
            "Warehouse" => 0.0,
            "FlatLoad" => 588.8601433
        ),
        "Chicago" => Dict(
            "FastFoodRest" => 84.2645196,
            "FullServiceRest" => 260.4454844,
            "Hospital" => 740.4172516,
            "LargeHotel" => 7005.083356,
            "LargeOffice" => 239.7065959,
            "MediumOffice" => 35.08184587,
            "MidriseApartment" => 297.4938584,
            "Outpatient" => 45.49600079,
            "PrimarySchool" => 179.4639347,
            "RetailStore" => 0.0,
            "SecondarySchool" => 456.8817409,
            "SmallHotel" => 406.0751832,
            "SmallOffice" => 11.2033023,
            "StripMall" => 0.0,
            "Supermarket" => 24.4292392,
            "Warehouse" => 0.0,
            "FlatLoad" => 611.6276445
        ),
        "Boulder" => Dict(
            "FastFoodRest" => 83.95201542,
            "FullServiceRest" => 259.3997752,
            "Hospital" => 737.372005,
            "LargeHotel" => 6971.32924,
            "LargeOffice" => 238.572519,
            "MediumOffice" => 34.9486709,
            "MidriseApartment" => 296.06471,
            "Outpatient" => 45.31437164,
            "PrimarySchool" => 178.3378526,
            "RetailStore" => 0.0,
            "SecondarySchool" => 453.228537,
            "SmallHotel" => 404.4154946,
            "SmallOffice" => 11.18970855,
            "StripMall" => 0.0,
            "Supermarket" => 24.36505320,
            "Warehouse" => 0.0,
            "FlatLoad" => 608.6556221
        ),
        "Minneapolis" => Dict(
            "FastFoodRest" => 89.48929949,
            "FullServiceRest" => 277.8184269,
            "Hospital" => 790.9262388,
            "LargeHotel" => 7563.607619,
            "LargeOffice" => 258.6874644,
            "MediumOffice" => 37.28641454,
            "MidriseApartment" => 321.1473562,
            "Outpatient" => 48.51884975,
            "PrimarySchool" => 191.9480118,
            "RetailStore" => 0.0,
            "SecondarySchool" => 491.5554097,
            "SmallHotel" => 433.8738637,
            "SmallOffice" => 11.42620649,
            "StripMall" => 0.0,
            "Supermarket" => 25.53218144,
            "Warehouse" => 0.0,
            "FlatLoad" => 658.8635839
        ),
        "Helena" => Dict(
            "FastFoodRest" => 90.44011877,
            "FullServiceRest" => 280.9757902,
            "Hospital" => 800.0940058,
            "LargeHotel" => 7665.023574,
            "LargeOffice" => 262.1461576,
            "MediumOffice" => 37.68905029,
            "MidriseApartment" => 325.4421541,
            "Outpatient" => 49.09222188,
            "PrimarySchool" => 193.4573283,
            "RetailStore" => 0.0,
            "SecondarySchool" => 494.7393735,
            "SmallHotel" => 438.9398731,
            "SmallOffice" => 11.46564268,
            "StripMall" => 0.0,
            "Supermarket" => 25.72866824,
            "Warehouse" => 0.0,
            "FlatLoad" => 667.2021224
        ),
        "Duluth" => Dict(
            "FastFoodRest" => 98.10641517,
            "FullServiceRest" => 306.4772907,
            "Hospital" => 874.2611723,
            "LargeHotel" => 8484.906093,
            "LargeOffice" => 290.0193773,
            "MediumOffice" => 40.92475821,
            "MidriseApartment" => 360.161261,
            "Outpatient" => 53.53681127,
            "PrimarySchool" => 211.2386551,
            "RetailStore" => 0.0,
            "SecondarySchool" => 543.3733772,
            "SmallHotel" => 479.7414481,
            "SmallOffice" => 11.79316054,
            "StripMall" => 0.0,
            "Supermarket" => 27.3451629,
            "Warehouse" => 0.0,
            "FlatLoad" => 736.3678114
        ),
        "Fairbanks" => Dict(
            "FastFoodRest" => 108.5335945,
            "FullServiceRest" => 341.1572799,
            "Hospital" => 975.1062178,
            "LargeHotel" => 9600.267161,
            "LargeOffice" => 327.8820873,
            "MediumOffice" => 45.32138512,
            "MidriseApartment" => 407.3910855,
            "Outpatient" => 59.6203514,
            "PrimarySchool" => 234.2595741,
            "RetailStore" => 0.0,
            "SecondarySchool" => 604.6838786,
            "SmallHotel" => 535.2525234,
            "SmallOffice" => 12.23744003,
            "StripMall" => 0.0,
            "Supermarket" => 29.53958045,
            "Warehouse" => 0.0,
            "FlatLoad" => 830.07826
        )
    )
    if isempty(city)
        city = find_ashrae_zone_city(latitude, longitude)
    end
    if !(buildingtype in default_buildings)
        error("buildingtype $(buildingtype) not in $(default_buildings).")
    end
    if isnothing(annual_mmbtu)
        annual_mmbtu = dhw_annual_mmbtu[city][buildingtype]
    end
    if length(monthly_mmbtu) == 12
        monthly_mmbtu = convert(Vector{Real}, monthly_mmbtu)
    end
    built_in_load("domestic_hot_water", city, buildingtype, year, annual_mmbtu, monthly_mmbtu)
end

function BuiltInSpaceHeatingLoad(
    city::String,
    buildingtype::String,
    latitude::Real,
    longitude::Real,
    year::Int,
    annual_mmbtu::Union{<:Real, Nothing}=nothing,
    monthly_mmbtu::Vector{<:Real}=Real[],
    )
    spaceheating_annual_mmbtu = Dict(
        "Miami" => Dict(
            "FastFoodRest" => 5.426780867,
            "FullServiceRest" => 12.03181471,
            "Hospital" => 6248.413294,
            "LargeHotel" => 198.0691407,
            "LargeOffice" => 168.9731637,
            "MediumOffice" => 0.036985655,
            "MidriseApartment" => 38.70606161,
            "Outpatient" => 2559.185872,
            "PrimarySchool" => 49.78021153,
            "RetailStore" => 12.12015432,
            "SecondarySchool" => 203.5185485,
            "SmallHotel" => 9.098564901,
            "SmallOffice" => 0.312524873,
            "StripMall" => 20.73216748,
            "Supermarket" => 101.2785324,
            "Warehouse" => 56.0796017,
            "FlatLoad" => 605.2352137
        ),
        "Houston" => Dict(
            "FastFoodRest" => 85.49111065,
            "FullServiceRest" => 199.7942842,
            "Hospital" => 8732.10385,
            "LargeHotel" => 1307.035548,
            "LargeOffice" => 2229.971744,
            "MediumOffice" => 16.25994314,
            "MidriseApartment" => 386.0269973,
            "Outpatient" => 2829.324307,
            "PrimarySchool" => 469.2532935,
            "RetailStore" => 289.0470815,
            "SecondarySchool" => 2011.1678969999998,
            "SmallHotel" => 108.9825885,
            "SmallOffice" => 19.55157672,
            "StripMall" => 292.23235389999996,
            "Supermarket" => 984.7374347000001,
            "Warehouse" => 475.9377273,
            "FlatLoad" => 1277.307359
        ),
        "Phoenix" => Dict(
            "FastFoodRest" => 57.89972381,
            "FullServiceRest" => 147.2569493,
            "Hospital" => 9382.021026,
            "LargeHotel" => 896.790817,
            "LargeOffice" => 1584.061452,
            "MediumOffice" => 1.922551528,
            "MidriseApartment" => 290.9887152,
            "Outpatient" => 3076.340876,
            "PrimarySchool" => 305.573525,
            "RetailStore" => 208.66330580000002,
            "SecondarySchool" => 1400.638544,
            "SmallHotel" => 83.98084516,
            "SmallOffice" => 9.988210938,
            "StripMall" => 230.16060699999997,
            "Supermarket" => 972.3008295,
            "Warehouse" => 362.42249280000004,
            "FlatLoad" => 1188.188154
        ),
        "Atlanta" => Dict(
            "FastFoodRest" => 168.8402371,
            "FullServiceRest" => 379.5865464,
            "Hospital" => 10467.659959999999,
            "LargeHotel" => 2427.6589879999997,
            "LargeOffice" => 3624.593975,
            "MediumOffice" => 49.00635733,
            "MidriseApartment" => 718.9316697,
            "Outpatient" => 3186.250588,
            "PrimarySchool" => 931.7212450999999,
            "RetailStore" => 627.2489826000001,
            "SecondarySchool" => 3968.4936420000004,
            "SmallHotel" => 202.26124219999997,
            "SmallOffice" => 42.74797302,
            "StripMall" => 615.2240506,
            "Supermarket" => 1880.5304489999999,
            "Warehouse" => 930.9449202,
            "FlatLoad" => 1888.856302
        ),
        "LasVegas" => Dict(
            "FastFoodRest" => 100.0877773,
            "FullServiceRest" => 247.21791319999997,
            "Hospital" => 9100.302056,
            "LargeHotel" => 1500.581408,
            "LargeOffice" => 2479.152321,
            "MediumOffice" => 5.220181581,
            "MidriseApartment" => 487.43122850000003,
            "Outpatient" => 2924.8220460000002,
            "PrimarySchool" => 499.6562223,
            "RetailStore" => 386.0185744,
            "SecondarySchool" => 2277.7410649999997,
            "SmallHotel" => 138.4427074,
            "SmallOffice" => 19.16330622,
            "StripMall" => 389.30494280000005,
            "Supermarket" => 1479.302604,
            "Warehouse" => 579.7671637999999,
            "FlatLoad" => 1413.3882199999998
        ),
        "LosAngeles" => Dict(
            "FastFoodRest" => 40.90390152,
            "FullServiceRest" => 97.94277036,
            "Hospital" => 10346.1713,
            "LargeHotel" => 707.848762,
            "LargeOffice" => 1458.148818,
            "MediumOffice" => 0.12342009699999999,
            "MidriseApartment" => 265.2851759,
            "Outpatient" => 3417.120585,
            "PrimarySchool" => 318.73600980000003,
            "RetailStore" => 175.104083,
            "SecondarySchool" => 1198.276619,
            "SmallHotel" => 72.42852638,
            "SmallOffice" => 5.898878347,
            "StripMall" => 193.18730269999998,
            "Supermarket" => 1040.273464,
            "Warehouse" => 323.96697819999997,
            "FlatLoad" => 1228.8385369999999
        ),
        "SanFrancisco" => Dict(
            "FastFoodRest" => 127.22328700000001,
            "FullServiceRest" => 362.48645889999995,
            "Hospital" => 11570.9155,
            "LargeHotel" => 1713.3629670000003,
            "LargeOffice" => 2690.1191,
            "MediumOffice" => 3.8159670660000002,
            "MidriseApartment" => 648.4472797999999,
            "Outpatient" => 3299.0539519999998,
            "PrimarySchool" => 818.2159102,
            "RetailStore" => 569.8081034,
            "SecondarySchool" => 3414.148347,
            "SmallHotel" => 189.0244446,
            "SmallOffice" => 27.53039453,
            "StripMall" => 526.2320428,
            "Supermarket" => 2301.616069,
            "Warehouse" => 675.6758453,
            "FlatLoad" => 1808.604729
        ),
        "Baltimore" => Dict(
            "FastFoodRest" => 305.2671204,
            "FullServiceRest" => 657.1337578,
            "Hospital" => 11253.61694,
            "LargeHotel" => 3731.0254619999996,
            "LargeOffice" => 5109.311943,
            "MediumOffice" => 116.8101842,
            "MidriseApartment" => 1132.964052,
            "Outpatient" => 3285.227941,
            "PrimarySchool" => 1428.239177,
            "RetailStore" => 1068.034778,
            "SecondarySchool" => 6557.634924,
            "SmallHotel" => 346.3683857,
            "SmallOffice" => 63.29818348,
            "StripMall" => 1075.39546,
            "Supermarket" => 2929.182261,
            "Warehouse" => 1568.722061,
            "FlatLoad" => 2539.2645399999997
        ),
        "Albuquerque" => Dict(
            "FastFoodRest" => 199.73581399999998,
            "FullServiceRest" => 398.5712205,
            "Hospital" => 8371.240776999999,
            "LargeHotel" => 2750.8382260000003,
            "LargeOffice" => 3562.0023950000004,
            "MediumOffice" => 47.49307973,
            "MidriseApartment" => 805.0965778,
            "Outpatient" => 2971.868562,
            "PrimarySchool" => 981.4176700999999,
            "RetailStore" => 755.4523907,
            "SecondarySchool" => 4338.227865999999,
            "SmallHotel" => 232.2194443,
            "SmallOffice" => 43.25360481,
            "StripMall" => 760.0982018,
            "Supermarket" => 2302.228741,
            "Warehouse" => 1151.250885,
            "FlatLoad" => 1854.437216
        ),
        "Seattle" => Dict(
            "FastFoodRest" => 255.5992711,
            "FullServiceRest" => 627.5634984000001,
            "Hospital" => 11935.157290000001,
            "LargeHotel" => 3343.683348,
            "LargeOffice" => 5266.970348,
            "MediumOffice" => 28.97979768,
            "MidriseApartment" => 1117.5465470000001,
            "Outpatient" => 3468.128914,
            "PrimarySchool" => 1263.541878,
            "RetailStore" => 952.2758742000001,
            "SecondarySchool" => 6367.850187,
            "SmallHotel" => 310.8087307,
            "SmallOffice" => 49.34878545,
            "StripMall" => 969.1074739000001,
            "Supermarket" => 3004.1844929999997,
            "Warehouse" => 1137.398514,
            "FlatLoad" => 2506.1340600000003
        ),
        "Chicago" => Dict(
            "FastFoodRest" => 441.93439000000006,
            "FullServiceRest" => 888.3312571,
            "Hospital" => 12329.57943,
            "LargeHotel" => 5104.848129,
            "LargeOffice" => 7706.028917000001,
            "MediumOffice" => 216.01411800000002,
            "MidriseApartment" => 1482.040156,
            "Outpatient" => 3506.5381090000005,
            "PrimarySchool" => 2006.0002120000001,
            "RetailStore" => 1472.8704380000001,
            "SecondarySchool" => 8962.172873,
            "SmallHotel" => 479.4653436000001,
            "SmallOffice" => 94.19308949,
            "StripMall" => 1497.556168,
            "Supermarket" => 3696.2112950000005,
            "Warehouse" => 2256.477231,
            "FlatLoad" => 3258.766323
        ),
        "Boulder" => Dict(
            "FastFoodRest" => 306.8980525,
            "FullServiceRest" => 642.8843574,
            "Hospital" => 9169.381845,
            "LargeHotel" => 3975.1080020000004,
            "LargeOffice" => 5027.882454,
            "MediumOffice" => 124.26913059999998,
            "MidriseApartment" => 1098.944993,
            "Outpatient" => 3087.969786,
            "PrimarySchool" => 1356.396807,
            "RetailStore" => 1086.9187570000001,
            "SecondarySchool" => 6268.036872,
            "SmallHotel" => 342.77800099999996,
            "SmallOffice" => 65.95714912,
            "StripMall" => 1093.093638,
            "Supermarket" => 2966.790122,
            "Warehouse" => 1704.8648210000001,
            "FlatLoad" => 2394.8859239999997
        ),
        "Minneapolis" => Dict(
            "FastFoodRest" => 588.8854722,
            "FullServiceRest" => 1121.229499,
            "Hospital" => 13031.2313,
            "LargeHotel" => 6359.946704,
            "LargeOffice" => 10199.279129999999,
            "MediumOffice" => 394.1525556,
            "MidriseApartment" => 1814.148381,
            "Outpatient" => 3661.1462229999997,
            "PrimarySchool" => 2600.964302,
            "RetailStore" => 1869.8106289999998,
            "SecondarySchool" => 11963.323859999999,
            "SmallHotel" => 618.0427338999999,
            "SmallOffice" => 128.12525349999999,
            "StripMall" => 1952.731917,
            "Supermarket" => 4529.776664,
            "Warehouse" => 3231.223746,
            "FlatLoad" => 4004.001148
        ),
        "Helena" => Dict(
            "FastFoodRest" => 468.8276835,
            "FullServiceRest" => 934.8994934,
            "Hospital" => 10760.57411,
            "LargeHotel" => 5554.910785,
            "LargeOffice" => 7373.056709,
            "MediumOffice" => 239.8330306,
            "MidriseApartment" => 1531.102079,
            "Outpatient" => 3390.42972,
            "PrimarySchool" => 2097.777112,
            "RetailStore" => 1494.85988,
            "SecondarySchool" => 9535.484059,
            "SmallHotel" => 499.85992930000003,
            "SmallOffice" => 98.85818175,
            "StripMall" => 1604.0043970000002,
            "Supermarket" => 3948.5338049999996,
            "Warehouse" => 2504.784991,
            "FlatLoad" => 3252.362248
        ),
        "Duluth" => Dict(
            "FastFoodRest" => 738.1353594999999,
            "FullServiceRest" => 1400.36692,
            "Hospital" => 14179.84149,
            "LargeHotel" => 7781.9012760000005,
            "LargeOffice" => 12504.64187,
            "MediumOffice" => 468.2112216,
            "MidriseApartment" => 2204.85149,
            "Outpatient" => 3774.3233130000003,
            "PrimarySchool" => 3160.1200719999997,
            "RetailStore" => 2298.8242920000002,
            "SecondarySchool" => 14468.64346,
            "SmallHotel" => 772.5386662000001,
            "SmallOffice" => 155.8350887,
            "StripMall" => 2411.847491,
            "Supermarket" => 5587.977185,
            "Warehouse" => 3962.122014,
            "FlatLoad" => 4741.886326
        ),
        "Fairbanks" => Dict(
            "FastFoodRest" => 1245.3608279999999,
            "FullServiceRest" => 2209.293209,
            "Hospital" => 20759.042680000002,
            "LargeHotel" => 12298.7791,
            "LargeOffice" => 23214.51532,
            "MediumOffice" => 949.8812392000001,
            "MidriseApartment" => 3398.039504,
            "Outpatient" => 4824.076322999999,
            "PrimarySchool" => 6341.861225,
            "RetailStore" => 3869.670979,
            "SecondarySchool" => 25619.149269999998,
            "SmallHotel" => 1264.41064,
            "SmallOffice" => 297.08593010000004,
            "StripMall" => 3934.89178,
            "Supermarket" => 8515.422039000001,
            "Warehouse" => 6882.6512680000005,
            "FlatLoad" => 7851.508208
        )
    )
    if isempty(city)
        city = find_ashrae_zone_city(latitude, longitude)
    end
    if !(buildingtype in default_buildings)
        error("buildingtype $(buildingtype) not in $(default_buildings).")
    end
    if isnothing(annual_mmbtu)
        annual_mmbtu = spaceheating_annual_mmbtu[city][buildingtype]
    end
    built_in_load("space_heating", city, buildingtype, year, annual_mmbtu, monthly_mmbtu)
end

function BuiltInCoolingLoad(
    city::String,
    buildingtype::String,
    latitude::Float64,
    longitude::Float64,
    year::Int,
    annual_tonhour::Union{<:Real, Nothing}=nothing,
    monthly_tonhour::Vector{<:Real}=Real[],
    existing_chiller_cop::Union{<:Real, Nothing}=nothing
    )
    cooling_annual_kwh = Dict(
        "Albuquerque" => Dict(
            "MidriseApartment" => 34088,
            "FastFoodRest" => 8466,
            "LargeOffice" => 442687,
            "Outpatient" => 357615,
            "PrimarySchool" => 99496,
            "SmallHotel" => 109982,
            "Supermarket" => 40883,
            "SmallOffice" => 8266,
            "StripMall" => 36498,
            "MediumOffice" => 103240,
            "LargeHotel" => 475348,
            "Warehouse" => 5831,
            "SecondarySchool" => 355818,
            "Hospital" => 1382003,
            "RetailStore" => 37955,
            "FullServiceRest" => 25534,
            "FlatLoad" => 220232.0
        ),
        "Helena" => Dict(
            "StripMall" => 13260,
            "RetailStore" => 13952,
            "Outpatient" => 233034,
            "FullServiceRest" => 8420,
            "SmallHotel" => 63978,
            "MediumOffice" => 44528,
            "Hospital" => 1075876,
            "LargeHotel" => 318965,
            "PrimarySchool" => 34468,
            "Supermarket" => 14879,
            "SmallOffice" => 3442,
            "SecondarySchool" => 132887,
            "FastFoodRest" => 2919,
            "MidriseApartment" => 11101,
            "LargeOffice" => 226066,
            "Warehouse" => 1566,
            "FlatLoad" => 137459.0
        ),
        "LosAngeles" => Dict(
            "MidriseApartment" => 22912,
            "Hospital" => 1987427,
            "Outpatient" => 435175,
            "RetailStore" => 22663,
            "SmallHotel" => 116966,
            "StripMall" => 26514,
            "FastFoodRest" => 3452,
            "LargeHotel" => 526169,
            "Supermarket" => 13630,
            "SmallOffice" => 8351,
            "Warehouse" => 673,
            "PrimarySchool" => 126130,
            "FullServiceRest" => 11456,
            "SecondarySchool" => 351337,
            "MediumOffice" => 121320,
            "LargeOffice" => 707617,
            "FlatLoad" => 280112.0
        ),
        "Boulder" => Dict(
            "Warehouse" => 3574,
            "Hospital" => 1223116,
            "PrimarySchool" => 57692,
            "SmallOffice" => 5097,
            "Supermarket" => 24684,
            "LargeHotel" => 384065,
            "LargeOffice" => 316540,
            "MidriseApartment" => 18535,
            "MediumOffice" => 68389,
            "RetailStore" => 23580,
            "FullServiceRest" => 14957,
            "SecondarySchool" => 220472,
            "Outpatient" => 291603,
            "FastFoodRest" => 5270,
            "StripMall" => 23350,
            "SmallHotel" => 84804,
            "FlatLoad" => 172858.0
        ),
        "Chicago" => Dict(
            "Hospital" => 1721142,
            "SmallHotel" => 97554,
            "SecondarySchool" => 336986,
            "SmallOffice" => 6365,
            "Supermarket" => 30785,
            "StripMall" => 32544,
            "FastFoodRest" => 6320,
            "Outpatient" => 328675,
            "LargeOffice" => 533117,
            "MidriseApartment" => 25663,
            "RetailStore" => 31774,
            "LargeHotel" => 469484,
            "MediumOffice" => 85839,
            "Warehouse" => 4370,
            "PrimarySchool" => 82619,
            "FullServiceRest" => 17222,
            "FlatLoad" => 238154.0
        ),
        "Houston" => Dict(
            "PrimarySchool" => 282346,
            "LargeHotel" => 1095441,
            "Hospital" => 2942336,
            "RetailStore" => 115410,
            "StripMall" => 111025,
            "FastFoodRest" => 24525,
            "LargeOffice" => 1471054,
            "Warehouse" => 17665,
            "MidriseApartment" => 97156,
            "SecondarySchool" => 1151985,
            "Outpatient" => 630778,
            "FullServiceRest" => 63610,
            "Supermarket" => 125580,
            "SmallOffice" => 19375,
            "SmallHotel" => 220470,
            "MediumOffice" => 239477,
            "FlatLoad" => 538015.0
        ),
        "Phoenix" => Dict(
            "SecondarySchool" => 1204514,
            "StripMall" => 122803,
            "PrimarySchool" => 302372,
            "LargeOffice" => 1124585,
            "SmallOffice" => 23720,
            "Supermarket" => 145706,
            "Outpatient" => 640137,
            "SmallHotel" => 233842,
            "MediumOffice" => 271151,
            "Hospital" => 2288973,
            "FullServiceRest" => 69402,
            "LargeHotel" => 1024137,
            "MidriseApartment" => 135120,
            "RetailStore" => 123508,
            "Warehouse" => 38561,
            "FastFoodRest" => 29956,
            "FlatLoad" => 486155.0
        ),
        "Fairbanks" => Dict(
            "Warehouse" => 97,
            "LargeOffice" => 147867,
            "Hospital" => 753420,
            "MidriseApartment" => 5035,
            "MediumOffice" => 24183,
            "LargeHotel" => 234082,
            "FullServiceRest" => 2025,
            "SmallOffice" => 1877,
            "Supermarket" => 2500,
            "PrimarySchool" => 15166,
            "RetailStore" => 3289,
            "SecondarySchool" => 53222,
            "Outpatient" => 143322,
            "FastFoodRest" => 656,
            "StripMall" => 3592,
            "SmallHotel" => 48066,
            "FlatLoad" => 89900.0
        ),
        "Seattle" => Dict(
            "StripMall" => 7664,
            "SmallHotel" => 60651,
            "SmallOffice" => 2503,
            "Supermarket" => 4882,
            "Outpatient" => 266075,
            "SecondarySchool" => 89468,
            "MediumOffice" => 38713,
            "MidriseApartment" => 6380,
            "LargeOffice" => 249231,
            "FullServiceRest" => 4144,
            "PrimarySchool" => 28984,
            "Hospital" => 1432385,
            "RetailStore" => 6223,
            "LargeHotel" => 292028,
            "FastFoodRest" => 1245,
            "Warehouse" => 329,
            "FlatLoad" => 155682.0
        ),
        "Duluth" => Dict(
            "Hospital" => 1149799,
            "StripMall" => 8809,
            "MediumOffice" => 37547,
            "RetailStore" => 8489,
            "LargeHotel" => 300364,
            "SecondarySchool" => 99148,
            "PrimarySchool" => 25136,
            "LargeOffice" => 249416,
            "MidriseApartment" => 8181,
            "FullServiceRest" => 4919,
            "Warehouse" => 889,
            "SmallHotel" => 61101,
            "FastFoodRest" => 1709,
            "Supermarket" => 9788,
            "SmallOffice" => 2563,
            "Outpatient" => 198121,
            "FlatLoad" => 135374.0
        ),
        "Minneapolis" => Dict(
            "Warehouse" => 2923,
            "SmallHotel" => 91962,
            "LargeOffice" => 475102,
            "Outpatient" => 289706,
            "Hospital" => 1516644,
            "SmallOffice" => 5452,
            "Supermarket" => 25962,
            "MediumOffice" => 74813,
            "PrimarySchool" => 62017,
            "FullServiceRest" => 14996,
            "LargeHotel" => 445049,
            "SecondarySchool" => 262353,
            "FastFoodRest" => 5429,
            "MidriseApartment" => 22491,
            "RetailStore" => 26859,
            "StripMall" => 26489,
            "FlatLoad" => 209265.0
        ),
        "Baltimore" => Dict(
            "SmallHotel" => 119859,
            "Outpatient" => 406989,
            "RetailStore" => 44703,
            "FullServiceRest" => 25140,
            "SecondarySchool" => 469545,
            "PrimarySchool" => 114110,
            "StripMall" => 43896,
            "MediumOffice" => 124064,
            "MidriseApartment" => 36890,
            "LargeOffice" => 822177,
            "Warehouse" => 7156,
            "LargeHotel" => 594327,
            "Hospital" => 2085465,
            "Supermarket" => 48913,
            "SmallOffice" => 8083,
            "FastFoodRest" => 9287,
            "FlatLoad" => 310038.0
        ),
        "LasVegas" => Dict(
            "RetailStore" => 89484,
            "LargeHotel" => 798740,
            "MediumOffice" => 207085,
            "StripMall" => 84138,
            "PrimarySchool" => 218688,
            "Warehouse" => 28928,
            "SmallOffice" => 16824,
            "Supermarket" => 114543,
            "MidriseApartment" => 92778,
            "FullServiceRest" => 53259,
            "Outpatient" => 513377,
            "FastFoodRest" => 22811,
            "SmallHotel" => 175011,
            "LargeOffice" => 792463,
            "SecondarySchool" => 858081,
            "Hospital" => 1831966,
            "FlatLoad" => 368636.0
        ),
        "Atlanta" => Dict(
            "LargeOffice" => 968973,
            "LargeHotel" => 704148,
            "Supermarket" => 69686,
            "StripMall" => 62530,
            "SmallOffice" => 11259,
            "MediumOffice" => 155865,
            "Hospital" => 2328054,
            "FullServiceRest" => 33832,
            "SmallHotel" => 149237,
            "SecondarySchool" => 608831,
            "Warehouse" => 8214,
            "Outpatient" => 489234,
            "FlatLoad" => 367407.0,
            "FastFoodRest" => 12565,
            "PrimarySchool" => 160138,
            "MidriseApartment" => 50907,
            "RetailStore" => 65044
        ),
        "SanFrancisco" => Dict(
            "SecondarySchool" => 121914,
            "LargeOffice" => 297494,
            "StripMall" => 6626,
            "Supermarket" => 3252,
            "SmallOffice" => 2539,
            "Hospital" => 1427329,
            "PrimarySchool" => 48323,
            "LargeHotel" => 296230,
            "MediumOffice" => 43956,
            "MidriseApartment" => 4484,
            "FullServiceRest" => 3062,
            "Warehouse" => 365,
            "FlatLoad" => 164036.0,
            "SmallHotel" => 71414,
            "RetailStore" => 4515,
            "FastFoodRest" => 838,
            "Outpatient" => 292238
        ),
        "Miami" => Dict(
            "LargeHotel" => 1467216,
            "SmallOffice" => 28235,
            "FastFoodRest" => 36779,
            "FlatLoad" => 700779.0,
            "Supermarket" => 150368,
            "FullServiceRest" => 101567,
            "PrimarySchool" => 434031,
            "Warehouse" => 20108,
            "LargeOffice" => 1878642,
            "SmallHotel" => 318595,
            "MidriseApartment" => 176446,
            "RetailStore" => 163290,
            "Outpatient" => 811571,
            "Hospital" => 3371923,
            "StripMall" => 180645,
            "SecondarySchool" => 1735906,
            "MediumOffice" => 337147
        )
    )
    if isempty(city)
        city = find_ashrae_zone_city(latitude, longitude)
    end
    if !(buildingtype in default_buildings)
        error("buildingtype $(buildingtype) not in $(default_buildings).")
    end
    # Set initial existing_chiller_cop to "cop_unknown_thermal" if not passed in; we will update existing_chiller_cop once the load profile is determined
    if isnothing(existing_chiller_cop)
        existing_chiller_cop = get_existing_chiller_default_cop()
    end
    if isnothing(annual_tonhour)
        annual_kwh = cooling_annual_kwh[city][buildingtype]
    else
        annual_kwh = annual_tonhour * KWH_THERMAL_PER_TONHOUR / existing_chiller_cop
    end
    monthly_kwh = Real[]
    if length(monthly_tonhour) == 12
        monthly_kwh = monthly_tonhour * KWH_THERMAL_PER_TONHOUR / existing_chiller_cop
    end
    built_in_load("cooling", city, buildingtype, year, annual_kwh, monthly_kwh)
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
"""
    BatchNormV2(channels::Integer, λ=identity;
                initβ=zeros32, initγ=ones32,
                affine = true, track_stats = true,
                ϵ=1f-5, momentum=0.1f0)

Same as [`Flux.BatchNorm`](@ref) but with more sensible fallbacks.

# Examples
```julia
m = Chain(Dense(28^2, 64), BatchNormV2(64, relu), Dense(64, 10), BatchNormV2(10), softmax)
```
"""
struct BatchNormV2{F,V,N,W}
    λ::F  # activation function
    β::V  # bias
    γ::V  # scale
    μ::W     # moving mean
    σ²::W    # moving std
    ϵ::N
    momentum::N
    chs::Int # number of channels
    attrs::NormAttributes
end

Flux.hasaffine(gn::BatchNormV2) = Flux.hasaffine(gn.attrs)
Flux._isactive(gn::BatchNormV2) = Flux._isactive(gn.attrs)

Flux.@functor BatchNormV2

Flux.trainable(gn::BatchNormV2) = hasaffine(gn) ? (gn.β, gn.γ) : ()

function BatchNormV2(chs::Int, λ=identity; initβ=zeros32, initγ=ones32, affine=true, track_stats=true, ϵ=1.0f-5,
                     momentum=0.1f0)
    β = affine ? initβ(chs) : nothing
    γ = affine ? initγ(chs) : nothing
    μ = track_stats ? zeros32(chs) : nothing
    σ² = track_stats ? ones32(chs) : nothing

    return BatchNormV2(λ, β, γ, μ, σ², ϵ, momentum, chs, NormAttributes(affine, track_stats, nothing))
end

function batchnorm_fallback(BN::BatchNormV2, x::AbstractArray{T,N}) where {T,N}
    @assert size(x, ndims(x) - 1) == BN.chs
    reduce_dims = [1:(N - 2); N]
    affine_shape = ntuple(i -> i == N - 1 ? size(x, N - 1) : 1, N)
    return norm_forward(BN, x, reduce_dims, affine_shape)
end

(BN::BatchNormV2)(x::AbstractArray) = batchnorm_fallback(BN, x)

function testmode!(m::BatchNormV2, mode=true)
    return (m.attrs.active = (isnothing(mode) || mode == :auto) ? nothing : !mode; m)
end

function Base.show(io::IO, l::BatchNormV2)
    print(io, "BatchNormV2($(l.chs)")
    (l.λ == identity) || print(io, ", $(l.λ)")
    hasaffine(l) || print(io, ", affine=false")
    return print(io, ")")
end

# Use our kernels for now
function (BN::BatchNormV2)(x::Union{CuArray{T,2},CuArray{T,4},CuArray{T,5}}) where {T<:Union{Float32,Float64}}
    (!hasaffine(BN) || !BN.attrs.track_stats) && return batchnorm_fallback(BN, x)
    return BN.λ.(batchnorm(BN.γ, BN.β, x, BN.μ, BN.σ², BN.momentum; alpha=1, beta=0, eps=BN.ϵ,
                           training=Flux._isactive(BN)))
end

ext: jl


text:
# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
# licensing details.

# Depth-first visit / traversal

#################################################
#
#  Depth-first visit
#
#################################################
"""
**Conventions in Breadth First Search and Depth First Search**
VertexColorMap :
- color == 0    => unseen
- color < 0     => examined but not closed
- color > 0     => examined and closed

EdgeColorMap :
- color == 0    => unseen
- color == 1     => examined
"""
struct DepthFirst <: SimpleGraphVisitAlgorithm end

function depth_first_visit_impl!(
    g::G,      # the graph
    stack,                          # an (initialized) stack of vertex
    vcolormap::AVertexMap,    # an (initialized) color-map to indicate status of vertices
    ecolormap::AEdgeMap,      # an (initialized) color-map to indicate status of edges
    visitor::SimpleGraphVisitor) where G<:AGraphOrDiGraph  # the visitor

    while !isempty(stack)
        # u, udsts, tstate = pop!(stack)
        u, udsts, it = pop!(stack)
        found_new_vertex = false

        # while !Base.done(udsts, tstate) && !found_new_vertex
        while it !== nothing && !found_new_vertex
            # v, tstate = Base.next(udsts, tstate)
            v, tstate = it
            ucolor = get(vcolormap, u, 0)
            vcolor = get(vcolormap, v, 0)
            v_edge = Edge(G,u,v) # ordering u,v for undirected graphs
            ecolor = get(ecolormap, v_edge, 0)
            examine_neighbor!(visitor, u, v, ucolor, vcolor, ecolor) #no return here

            ecolormap[v_edge] = 1
            if vcolor == 0
                found_new_vertex = true
                vcolormap[v] = vcolormap[u] - 1 #negative numbers
                discover_vertex!(visitor, v) || return
                # push!(stack, (u, udsts, tstate))
                it = iterate(udsts, tstate)
                if it !== nothing
                    push!(stack, (u, udsts, it))
                end

                open_vertex!(visitor, v)
                vdsts = out_neighbors(g, v)
                if !isempty(vdsts)
                    push!(stack, (v, vdsts, iterate(vdsts)))
                end
            end
        end

        if !found_new_vertex
            close_vertex!(visitor, u)
            vcolormap[u] *= -1 #revert to positive
        end
    end
end

function traverse_graph!(
        g::AGraphOrDiGraph,
        alg::DepthFirst,
        s,
        visitor::SimpleGraphVisitor;
        vcolormap::AVertexMap = VertexMap(g, Int),
        ecolormap::AEdgeMap = ConstEdgeMap(g, 0)
    )

    vcolormap[s] = -1
    discover_vertex!(visitor, s) || return

    sdsts = out_neighbors(g, s)
    # sstate = Base.start(sdsts)
    it = iterate(sdsts)
    stack = [(s, sdsts, it)]

    depth_first_visit_impl!(g, stack, vcolormap, ecolormap, visitor)
end

#################################################
#
#  Useful applications
#
#################################################

# Test whether a graph is cyclic

mutable struct DFSCyclicTestVisitor <: SimpleGraphVisitor
    found_cycle::Bool

    DFSCyclicTestVisitor() = new(false)
end

function examine_neighbor!(
    vis::DFSCyclicTestVisitor,
    u,
    v,
    ucolor,
    vcolor,
    ecolor)

    if vcolor < 0 && vcolor != ucolor+1 # seen and not parent
        vis.found_cycle = true
    end
end

discover_vertex!(vis::DFSCyclicTestVisitor, v) = !vis.found_cycle

"""
    has_cycles(g)

Tests whether a graph contains a simple cycle through depth-first search.
See also [`is_tree`](@ref).
"""
function has_cycles(g::AGraphOrDiGraph)
    cmap = VertexMap(g, zeros(Int, nv(g)))
    visitor = DFSCyclicTestVisitor()
    # em = is_directed(g) ? ConstEdgeMap(g, 0) : EdgeMap(g, Dict{Edge,Int}())
    em = ConstEdgeMap(g, 0)
    for s in vertices(g)
        if cmap[s] == 0
            traverse_graph!(g, DepthFirst(), s, visitor, vcolormap=cmap,
                ecolormap=em)
        end
        visitor.found_cycle && return true
    end
    return false
end
# TODO
# now has_cycles for undir graphs
# is true whenever ne(g) > 0, so this has to be changed
#
# """
#     is_tree(g)

# Check whether `g` is a tree.
# Return `false` whenever [`has_cycles`](@ref) returns `true` and viceversa.
# """
# is_tree(g::AGraphOrDiGraph) = !has_cycles(g)

## This is a faster implementation found on
## http://www.geeksforgeeks.org/detect-cycle-in-a-graph/
##  TODO benchmark carefully and eventually use this
# function _has_cycles(g, v, visited, parent)
#     visited[v] = 1
#     for i in neighbors(g, v)
#         if visited[i] == 0
#             _has_cycles(g, i, visited, v) && return true
#         elseif i != parent
#             return true
#         end
#     end
#     return false
# end
#
#
# function has_cycles(g::AGraph)
#     visited = zeros(Int, nv(g))
#     for i=1:nv(g)
#         if visited[i] == 0
#             _has_cycles(g, i, visited, -1) && return true
#         end
#     end
#
#     return false
# end
#

# Topological sort using DFS

mutable struct TopologicalSortVisitor <: SimpleGraphVisitor
    vertices::Vector{Int}

    function TopologicalSortVisitor(n)
        vs = Vector{Int}()
        sizehint!(vs, n)
        new(vs)
    end
end

function examine_neighbor!(visitor::TopologicalSortVisitor, u, v, ucolor, vcolor, ecolor)
    (vcolor < 0 && ecolor == 0) && error("The input graph contains at least one loop.")
end

function close_vertex!(visitor::TopologicalSortVisitor, v)
    push!(visitor.vertices, v)
end

function topological_sort_by_dfs(g::AGraphOrDiGraph)
    nvg = nv(g)
    cmap = VertexMap(g, zeros(Int, nvg))
    visitor = TopologicalSortVisitor(nvg)

    for s in vertices(g)
        if cmap[s] == 0
            traverse_graph!(g, DepthFirst(), s, visitor, vcolormap=cmap)
        end
    end

    reverse(visitor.vertices)
end

mutable struct TreeDFSVisitor{G <:ADiGraph} <:SimpleGraphVisitor
    tree::G
    predecessor::Vector{Int}
end

TreeDFSVisitor(n, ::Type{G}) where {G} = TreeDFSVisitor(digraph(G(n)), zeros(Int,n))

function examine_neighbor!(visitor::TreeDFSVisitor, u, v, ucolor, vcolor, ecolor)
    if (vcolor == 0)
        visitor.predecessor[v] = u
    end
    return true
end

"""
    dfs_tree(g, s)

Provides a depth-first traversal of the graph `g` starting with source vertex `s`,
and returns a directed acyclic graph of vertices in the order they were discovered.
"""
function dfs_tree(g::G, s) where G
    nvg = nv(g)
    visitor = TreeDFSVisitor(nvg, G)
    traverse_graph!(g, DepthFirst(), s, visitor)
    h = digraph(G(nvg))
    for (v, u) in enumerate(visitor.predecessor)
        if u != 0
            add_edge!(h, u, v)
        end
    end
    return h
end


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

module TestHelpers

using Serialization

include("dimensionful.jl")
export Furlong

function open_fake_pty()
    @static if Sys.iswindows()
        error("Unable to create a fake PTY in Windows")
    end

    O_RDWR = Base.Filesystem.JL_O_RDWR
    O_NOCTTY = Base.Filesystem.JL_O_NOCTTY

    fdm = ccall(:posix_openpt, Cint, (Cint,), O_RDWR|O_NOCTTY)
    fdm == -1 && error("Failed to open PTY master")
    rc = ccall(:grantpt, Cint, (Cint,), fdm)
    rc != 0 && error("grantpt failed")
    rc = ccall(:unlockpt, Cint, (Cint,), fdm)
    rc != 0 && error("unlockpt")

    fds = ccall(:open, Cint, (Ptr{UInt8}, Cint),
        ccall(:ptsname, Ptr{UInt8}, (Cint,), fdm), O_RDWR|O_NOCTTY)

    # slave
    slave = RawFD(fds)
    master = Base.TTY(RawFD(fdm); readable = true)
    slave, master
end

function with_fake_pty(f)
    slave, master = open_fake_pty()
    try
        f(slave, master)
    finally
        ccall(:close,Cint,(Cint,),slave) # XXX: this causes the kernel to throw away all unread data on the pty
        close(master)
    end
end

# OffsetArrays (arrays with indexing that doesn't start at 1)

# This test file is designed to exercise support for generic indexing,
# even though offset arrays aren't implemented in Base.

module OAs

using Base: Indices, IndexCartesian, IndexLinear, tail

export OffsetArray

struct OffsetArray{T,N,AA<:AbstractArray} <: AbstractArray{T,N}
    parent::AA
    offsets::NTuple{N,Int}
end
OffsetVector{T,AA<:AbstractArray} = OffsetArray{T,1,AA}

OffsetArray(A::AbstractArray{T,N}, offsets::NTuple{N,Int}) where {T,N} = OffsetArray{T,N,typeof(A)}(A, offsets)
OffsetArray(A::AbstractArray{T,N}, offsets::Vararg{Int,N}) where {T,N} = OffsetArray(A, offsets)

OffsetArray{T,N}(::UndefInitializer, inds::Indices{N}) where {T,N} =
    OffsetArray{T,N,Array{T,N}}(Array{T,N}(undef, map(length, inds)), map(indsoffset, inds))
OffsetArray{T}(::UndefInitializer, inds::Indices{N}) where {T,N} =
    OffsetArray{T,N}(undef, inds)

Base.IndexStyle(::Type{T}) where {T<:OffsetArray} = Base.IndexStyle(parenttype(T))
parenttype(::Type{OffsetArray{T,N,AA}}) where {T,N,AA} = AA
parenttype(A::OffsetArray) = parenttype(typeof(A))

Base.parent(A::OffsetArray) = A.parent

errmsg(A) = error("size not supported for arrays with indices $(axes(A)); see https://docs.julialang.org/en/latest/devdocs/offset-arrays/")
Base.size(A::OffsetArray) = errmsg(A)
Base.size(A::OffsetArray, d) = errmsg(A)
Base.eachindex(::IndexCartesian, A::OffsetArray) = CartesianIndices(axes(A))
Base.eachindex(::IndexLinear, A::OffsetVector) = axes(A, 1)

# Implementations of indices and indices1. Since bounds-checking is
# performance-critical and relies on indices, these are usually worth
# optimizing thoroughly.
@inline Base.axes(A::OffsetArray, d) = 1 <= d <= length(A.offsets) ? Base.Slice(axes(parent(A))[d] .+ A.offsets[d]) : Base.Slice(1:1)
@inline Base.axes(A::OffsetArray) = _indices(axes(parent(A)), A.offsets)  # would rather use ntuple, but see #15276
@inline _indices(inds, offsets) = (Base.Slice(inds[1] .+ offsets[1]), _indices(tail(inds), tail(offsets))...)
_indices(::Tuple{}, ::Tuple{}) = ()
Base.indices1(A::OffsetArray{T,0}) where {T} = Base.Slice(1:1)  # we only need to specialize this one

const OffsetAxis = Union{Integer, UnitRange, Base.Slice{<:UnitRange}}
function Base.similar(A::OffsetArray, T::Type, dims::Dims)
    B = similar(parent(A), T, dims)
end
function Base.similar(A::AbstractArray, T::Type, inds::Tuple{OffsetAxis,Vararg{OffsetAxis}})
    B = similar(A, T, map(indslength, inds))
    OffsetArray(B, map(indsoffset, inds))
end

Base.similar(::Type{T}, shape::Tuple{OffsetAxis,Vararg{OffsetAxis}}) where {T<:AbstractArray} =
    OffsetArray(T(undef, map(indslength, shape)), map(indsoffset, shape))

Base.reshape(A::AbstractArray, inds::Tuple{OffsetAxis,Vararg{OffsetAxis}}) = OffsetArray(reshape(A, map(indslength, inds)), map(indsoffset, inds))

Base.fill(v, inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {N} =
    fill!(OffsetArray(Array{typeof(v), N}(undef, map(indslength, inds)), map(indsoffset, inds)), v)
Base.zeros(::Type{T}, inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {T, N} =
    fill!(OffsetArray(Array{T, N}(undef, map(indslength, inds)), map(indsoffset, inds)), zero(T))
Base.ones(::Type{T}, inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {T, N} =
    fill!(OffsetArray(Array{T, N}(undef, map(indslength, inds)), map(indsoffset, inds)), one(T))
Base.trues(inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {N} =
    fill!(OffsetArray(BitArray{N}(undef, map(indslength, inds)), map(indsoffset, inds)), true)
Base.falses(inds::NTuple{N, Union{Integer, AbstractUnitRange}}) where {N} =
    fill!(OffsetArray(BitArray{N}(undef, map(indslength, inds)), map(indsoffset, inds)), false)

@inline function Base.getindex(A::OffsetArray{T,N}, I::Vararg{Int,N}) where {T,N}
    checkbounds(A, I...)
    @inbounds ret = parent(A)[offset(A.offsets, I)...]
    ret
end
# Vectors don't support one-based linear indexing; they always use the offsets
@inline function Base.getindex(A::OffsetVector, i::Int)
    checkbounds(A, i)
    @inbounds ret = parent(A)[offset(A.offsets, (i,))[1]]
    ret
end
# But multidimensional arrays allow one-based linear indexing
@inline function Base.getindex(A::OffsetArray, i::Int)
    checkbounds(A, i)
    @inbounds ret = parent(A)[i]
    ret
end
@inline function Base.setindex!(A::OffsetArray{T,N}, val, I::Vararg{Int,N}) where {T,N}
    checkbounds(A, I...)
    @inbounds parent(A)[offset(A.offsets, I)...] = val
    val
end
@inline function Base.setindex!(A::OffsetVector, val, i::Int)
    checkbounds(A, i)
    @inbounds parent(A)[offset(A.offsets, (i,))[1]] = val
    val
end
@inline function Base.setindex!(A::OffsetArray, val, i::Int)
    checkbounds(A, i)
    @inbounds parent(A)[i] = val
    val
end

@inline function Base.deleteat!(A::OffsetArray, i::Int)
    checkbounds(A, i)
    @inbounds deleteat!(parent(A), offset(A.offsets, (i,))[1])
end

@inline function Base.deleteat!(A::OffsetArray{T,N}, I::Vararg{Int, N}) where {T,N}
    checkbounds(A, I...)
    @inbounds deleteat!(parent(A), offset(A.offsets, I)...)
end

@inline function Base.deleteat!(A::OffsetArray, i::UnitRange{Int})
    checkbounds(A, first(i))
    checkbounds(A, last(i))
    first_idx = offset(A.offsets, (first(i),))[1]
    last_idx = offset(A.offsets, (last(i),))[1]
    @inbounds deleteat!(parent(A), first_idx:last_idx)
end

# Computing a shifted index (subtracting the offset)
offset(offsets::NTuple{N,Int}, inds::NTuple{N,Int}) where {N} = _offset((), offsets, inds)
_offset(out, ::Tuple{}, ::Tuple{}) = out
@inline _offset(out, offsets, inds) = _offset((out..., inds[1]-offsets[1]), Base.tail(offsets), Base.tail(inds))

indsoffset(r::AbstractRange) = first(r) - 1
indsoffset(i::Integer) = 0
indslength(r::AbstractRange) = Base._length(r)
indslength(i::Integer) = i

Base.resize!(A::OffsetVector, nl::Integer) = (resize!(A.parent, nl); A)

end

# Mimic a quantity with a physical unit that is not convertable to a real number
struct PhysQuantity{n,T}   # n is like the exponent of the unit
    val::T
end
PhysQuantity{n}(x::T) where {n,T} = PhysQuantity{n,T}(x)
Base.zero(::Type{PhysQuantity{n,T}}) where {n,T} = PhysQuantity{n,T}(zero(T))
Base.zero(x::PhysQuantity) = zero(typeof(x))
Base.:+(x::PhysQuantity{n}, y::PhysQuantity{n}) where n = PhysQuantity{n}(x.val + y.val)
Base.:-(x::PhysQuantity{n}, y::PhysQuantity{n}) where n = PhysQuantity{n}(x.val - y.val)
Base.:*(x::PhysQuantity{n,T}, y::Int) where {n,T} = PhysQuantity{n}(x.val*y)
Base.:/(x::PhysQuantity{n,T}, y::Int) where {n,T} = PhysQuantity{n}(x.val/y)
Base.:*(x::PhysQuantity{n1,S}, y::PhysQuantity{n2,T}) where {n1,n2,S,T} =
    PhysQuantity{n1+n2}(x.val*y.val)
Base.:/(x::PhysQuantity{n1,S}, y::PhysQuantity{n2,T}) where {n1,n2,S,T} =
    PhysQuantity{n1-n2}(x.val/y.val)
Base.convert(::Type{PhysQuantity{0,T}}, x::Int) where T = PhysQuantity{0}(convert(T, x))
Base.convert(::Type{P}, ::Int) where P<:PhysQuantity =
    error("Int is incommensurate with PhysQuantity")

end

max_stars_repo_name: perrutquist/julia

max_issues_repo_name: perrutquist/julia

max_forks_repo_name: perrutquist/julia


text:
export prepend_bootswatch, dbc_themes

function prepend_bootswatch(suffix)
  BOOTSWATCH_BASE = "https://stackpath.bootstrapcdn.com/bootswatch/4.5.2/"
  return string(BOOTSWATCH_BASE, suffix)
end

dbc_themes = (
  BOOTSTRAP = "https://stackpath.bootstrapcdn.com/bootstrap/4.5.2/css/bootstrap.min.css",
  GRID = "https://stackpath.bootstrapcdn.com/bootstrap/4.5.2/css/bootstrap-grid.min.css",
  CERULEAN = prepend_bootswatch("cerulean/bootstrap.min.css"),
  COSMO = prepend_bootswatch("cosmo/bootstrap.min.css"),
  CYBORG = prepend_bootswatch("cyborg/bootstrap.min.css"),
  DARKLY = prepend_bootswatch("darkly/bootstrap.min.css"),
  FLATLY = prepend_bootswatch("flatly/bootstrap.min.css"),
  JOURNAL = prepend_bootswatch("journal/bootstrap.min.css"),
  LITERA = prepend_bootswatch("litera/bootstrap.min.css"),
  LUMEN = prepend_bootswatch("lumen/bootstrap.min.css"),
  LUX = prepend_bootswatch("lux/bootstrap.min.css"),
  MATERIA = prepend_bootswatch("materia/bootstrap.min.css"),
  MINTY = prepend_bootswatch("minty/bootstrap.min.css"),
  PULSE = prepend_bootswatch("pulse/bootstrap.min.css"),
  SANDSTONE = prepend_bootswatch("sandstone/bootstrap.min.css"),
  SIMPLEX = prepend_bootswatch("simplex/bootstrap.min.css"),
  SKETCHY = prepend_bootswatch("sketchy/bootstrap.min.css"),
  SLATE = prepend_bootswatch("slate/bootstrap.min.css"),
  SOLAR = prepend_bootswatch("solar/bootstrap.min.css"),
  SPACELAB = prepend_bootswatch("spacelab/bootstrap.min.css"),
  SUPERHERO = prepend_bootswatch("superhero/bootstrap.min.css"),
  UNITED = prepend_bootswatch("united/bootstrap.min.css"),
  YETI = prepend_bootswatch("yeti/bootstrap.min.css")
)

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
# caeser cipher
const n = length('A':'Z')
# function encrypt(ch::Char, shift::Integer)
#     if ch in 'A':'Z'
#         'A' + mod((ch - 'A') + shift, n)
#     elseif ch in 'a':'z'
#         'a' + mod((ch - 'a') + shift, n)
#     else
#         ch
#     end
# end

# function decrypt(ch::Char, shift::Integer)
#     if ch in 'A':'Z'
#         'A' + mod((ch - 'A') - shift, n)
#     elseif ch in 'a':'z'
#         'a' + mod((ch - 'a') - shift, n)
#     else
#         ch
#     end
# end

# more elegant solution based on discussion here:
# https://discourse.julialang.org/t/a-switch-type-function/13491/15
function encryptc(ch::Char, shift::Integer)
    ch in 'A':'Z' ? 'A' + mod((ch - 'A') + shift, n) :
    ch in 'a':'z' ? 'a' + mod((ch - 'a') + shift, n) :
    #= otherwise =# ch
end

function decryptc(ch::Char, shift::Integer)
    ch in 'A':'Z' ? 'A' + mod((ch - 'A') - shift, n) :
    ch in 'a':'z' ? 'a' + mod((ch - 'a') - shift, n) :
    #= otherwise =# ch
end

# Example
# julia> message = "This is the end of the world!"
# "This is the end of the world"
#
# julia> shift = 3
# 3
#
# julia> emessage = map(message) do char
#     encryptc(char, shift)
# end
# "Wklv lv wkh hqg ri wkh zruog!"

# julia> map(emessage) do char
#     decryptc(char, shift)
# end
# "This is the end of the world!"

# substition cipher
using Random
# Engheim's definition is wrong - he uses 'char' for 'ch'
function encrypts(ch::Char, updict::Dict{Char, Char}, lowdict::Dict{Char, Char})
    ch in 'A':'Z' ? get(updict, ch, ch) :
    ch in 'a':'z' ? get(lowdict, ch, ch) :
    #= otherwise =# ch
end

# Also don't think decrypt is confusing as Engheim claims
function decrypts(ch::Char, rupdict::Dict{Char, Char}, rlowdict::Dict{Char, Char})
    ch in 'A':'Z' ? get(rupdict, ch, ch) :
    ch in 'a':'z' ? get(rlowdict, ch, ch) :
    #= otherwise =# ch
end

# example
# upshuf=shuffle('A':'Z')
# lowshuf=shuffle('a':'z')
# julia> updict = Dict(zip('A':'Z', upshuf))
# Dict{Char, Char} with 26 entries:
#   'E' => 'D'
#   'Z' => 'C'
#   'X' => 'G'
#   'B' => 'U'
#   'C' => 'T'
#   'D' => 'Y'
#   'A' => 'S'
#   'R' => 'F'
#   'G' => 'A'
#   'F' => 'P'
#   'N' => 'Q'
#   'M' => 'M'
#   'K' => 'V'
#   ⋮   => ⋮
#
#   julia> lowdict = Dict(zip('a':'z', lowshuf))
#   Dict{Char, Char} with 26 entries:
#     'n' => 'x'
#     'f' => 'b'
#     'w' => 'y'
#     'd' => 'a'
#     'e' => 'f'
#     'o' => 'l'
#     'h' => 'c'
#     'j' => 'm'
#     'i' => 't'
#     'k' => 'z'
#     'r' => 'v'
#     's' => 'q'
#     't' => 'w'
#     ⋮   => ⋮
#
#   julia> encrypts('a', updict, lowdict)
#   'i': ASCII/Unicode U+0069 (category Ll: Letter, lowercase)
#
#   julia> encrypts('?', updict, lowdict)
#   '?': ASCII/Unicode U+003F (category Po: Punctuation, other)
#
#   julia> encrypts('Z', updict, lowdict)
#   'C': ASCII/Unicode U+0043 (category Lu: Letter, uppercase)
#
# The idea for how to reverse the dictionary comes from here:
# https://discourse.julialang.org/t/is-something-like-reversed-dict-findall-x-x-house-cc-is-to-slow/16443/2
# julia> rupdict = Dict(values(updict) .=> keys(updict))
# Dict{Char, Char} with 26 entries:
#   'E' => 'P'
#   'Z' => 'O'
#   'X' => 'Q'
#   'C' => 'Z'
#   'B' => 'Y'
#   'D' => 'E'
#   'A' => 'G'
#   'R' => 'T'
#   'G' => 'X'
#   'F' => 'R'
#   'Q' => 'N'
#   'P' => 'F'
#   'M' => 'M'
#   ⋮   => ⋮
#
# julia> rlowdict = Dict(values(lowdict) .=> keys(lowdict))
# Dict{Char, Char} with 26 entries:
#   'n' => 'g'
#   'f' => 'e'
#   'w' => 't'
#   'd' => 'q'
#   'e' => 'v'
#   'o' => 'm'
#   'h' => 'z'
#   'y' => 'w'
#   'i' => 'a'
#   't' => 'i'
#   's' => 'l'
#   'r' => 'u'
#   'j' => 'x'
#   ⋮   => ⋮
# julia> message = "This is the end of the world!"
# "This is the end of the world!"

# julia> emessage = map(message) do char
#            encrypts(char,updict,lowdict)
#         end
# "Rctq tq wcf fxa lb wcf ylvsa!"

# julia> map(emessage) do char
#            decrypts(char,rupdict,rlowdict)
#         end
# "This is the end of the world!"

max_stars_repo_name: FourMInfo/JuliaforBeginners

max_issues_repo_name: FourMInfo/JuliaforBeginners

max_forks_repo_name: FourMInfo/JuliaforBeginners


text:
scribename(fname::Symbol) = Symbol(string("##", fname, "__scribe"))

rawfuncname(fname::Symbol) = Symbol(string("##", fname, "__raw"))

macro scribename(fname::Symbol)
    fname = Expr(:quote, scribename(fname))
    esc(:($fname))
end

macro scribeof(func::Symbol)
    sname = scribename(func)
    esc(:($sname))
end

macro rawfunc(func::Symbol)
    fname = rawfuncname(func)
    esc(:($fname))
end

macro scribeofrawfunc(func::Symbol)
    sname = scribename(rawfuncname(func))
    esc(:($sname))
end

macro rawcall(fcall::Expr)
    @capture(fcall, fname_(args__))
    fname = rawfuncname(fname)
    esc(:($fname($(args...))))
end


text:
# Configurations and forces
mutable struct NewtonIndicesConfigurationForce{nq,nu,nc,nb,n1,n2,n3} <: NewtonIndices
    nd::Int                                   # implicit dynamics constraint
    nr::Int                                   # size of a one-time-step block
    iq::SizedArray{Tuple{nq},Int,1,1}         # configuration indices
    iu::SizedArray{Tuple{nu},Int,1,1}         # control indices
    iγ::SizedArray{Tuple{nc},Int,1,1}         # impact indices
    ib::SizedArray{Tuple{nb},Int,1,1}         # linear friction indices
    iν::SizedArray{Tuple{n1},Int,1,1}         # implicit dynamics lagrange multiplier
    iz::SizedArray{Tuple{n2},Int,1,1}         # IP solver solution [q2, γ1, b1]
    iθ::SizedArray{Tuple{n3},Int,1,1}         # IP solver data [q0, q1, u1]
    Iq::Vector{SizedArray{Tuple{nq},Int,1,1}} # configuration indices
    Iu::Vector{SizedArray{Tuple{nu},Int,1,1}} # control indices
    Iγ::Vector{SizedArray{Tuple{nc},Int,1,1}} # impact indices
    Ib::Vector{SizedArray{Tuple{nb},Int,1,1}} # linear friction indices
    Iν::Vector{SizedArray{Tuple{n1},Int,1,1}} # implicit dynamics lagrange multiplier
    Iz::Vector{SizedArray{Tuple{n2},Int,1,1}} # IP solver solution [q2, γ1, b1]
    Iθ::Vector{SizedArray{Tuple{n3},Int,1,1}} # IP solver data [q0, q1, u1]
end

function NewtonIndicesConfigurationForce(model::ContactModel, env::Environment, H::Int)
    dim = model.dim
    nq = dim.q # configuration
    nu = dim.u # control
    nw = dim.w # disturbance
    nc = dim.c # contact
    nb = nc * friction_dim(env) # linear friction
    nd = nq + nc + nb # implicit dynamics constraint
    nr = nq + nu + nc + nb + nd # size of a one-time-step block

    off = 0
    iq = SizedVector{nq}(off .+ (1:nq)); off += nq # index of the configuration q2
    iu = SizedVector{nu}(off .+ (1:nu)); off += nu # index of the control u1
    iγ = SizedVector{nc}(off .+ (1:nc)); off += nc # index of the impact γ1
    ib = SizedVector{nb}(off .+ (1:nb)); off += nb # index of the linear friction b1
    iν = SizedVector{nd}(off .+ (1:nd)); off += nd # index of the dynamics lagrange multiplier ν1
    iz = vcat(iq, iγ, ib) # index of the IP solver solution [q2, γ1, b1]
    iθ = vcat(iq .- 2nr, iq .- nr, iu) # index of the IP solver data [q0, q1, u1]

    Iq = [(t - 1) * nr .+ iq for t = 1:H]
    Iu = [(t - 1) * nr .+ iu for t = 1:H]
    Iγ = [(t - 1) * nr .+ iγ for t = 1:H]
    Ib = [(t - 1) * nr .+ ib for t = 1:H]
    Iν = [(t - 1) * nr .+ iν for t = 1:H]
    Iz = [(t - 1) * nr .+ iz for t = 1:H]
    Iθ = [(t - 1) * nr .+ iθ for t = 1:H]

    return NewtonIndicesConfigurationForce{nq,nu,nc,nb,nd,nd,2nq+nu}(
        nd, nr,
        iq, iu, iγ, ib, iν, iz, iθ,
        Iq, Iu, Iγ, Ib, Iν, Iz, Iθ)
end

mutable struct NewtonIndicesConfiguration{nq,nu,n1,n2,n3} <: NewtonIndices
    nd::Int                                   # implicit dynamics constraint
    nr::Int                                   # size of a one-time-step block
    iq::SizedArray{Tuple{nq},Int,1,1}         # configuration indices
    iu::SizedArray{Tuple{nu},Int,1,1}         # control indices
    iν::SizedArray{Tuple{n1},Int,1,1}         # implicit dynamics lagrange multiplier
    iz::SizedArray{Tuple{n2},Int,1,1}         # IP solver solution [q2]
    iθ::SizedArray{Tuple{n3},Int,1,1}         # IP solver data [q0, q1, u1]
    Iq::Vector{SizedArray{Tuple{nq},Int,1,1}} # configuration indices
    Iu::Vector{SizedArray{Tuple{nu},Int,1,1}} # control indices
    Iν::Vector{SizedArray{Tuple{n1},Int,1,1}} # implicit dynamics lagrange multiplier
    Iz::Vector{SizedArray{Tuple{n2},Int,1,1}} # IP solver solution [q2]
    Iθ::Vector{SizedArray{Tuple{n3},Int,1,1}} # IP solver data [q0, q1, u1]
end

function NewtonIndicesConfiguration(model::ContactModel, env::Environment, H::Int)
    dim = model.dim
    nq = dim.q # configuration
    nu = dim.u # control
    nw = dim.w # disturbance
    nc = dim.c # contact
    nd = nq # implicit dynamics constraint
    nr = nq + nu + nd # size of a one-time-step block

    off = 0
    iq = SizedVector{nq}(off .+ (1:nq)); off += nq # index of the configuration q2
    iu = SizedVector{nu}(off .+ (1:nu)); off += nu # index of the control u1
    iν = SizedVector{nd}(off .+ (1:nd)); off += nd # index of the dynamics lagrange multiplier ν1
    iz = copy(iq)                                  # index of the IP solver solution [q2]
    iθ = vcat(iq .- 2nr, iq .- nr, iu) # index of the IP solver data [q0, q1, u1]

    Iq = [(t - 1) * nr .+ iq for t = 1:H]
    Iu = [(t - 1) * nr .+ iu for t = 1:H]
    Iν = [(t - 1) * nr .+ iν for t = 1:H]
    Iz = [(t - 1) * nr .+ iz for t = 1:H]
    Iθ = [(t - 1) * nr .+ iθ for t = 1:H]

    return NewtonIndicesConfiguration{nq,nu,nd,nd,2nq+nu}(
        nd, nr,
        iq, iu, iν, iz, iθ,
        Iq, Iu, Iν, Iz, Iθ)
end

function NewtonIndices(model::ContactModel, env::Environment, H::Int;
    mode = :configurationforce)

    if mode == :configurationforce
        return NewtonIndicesConfigurationForce(model, env, H)
    elseif mode == :configuration
        NewtonIndicesConfiguration(model, env, H)
    else
        @error "mode not implemented"
    end
end


text:
using OffsetArrays, LinearAlgebra
function mydotavx(a, b)
    s = zero(eltype(a))
    @avxt for i ∈ eachindex(a,b)
        s += a[i]*b[i]
    end
    s
end
function AmulB!(C,A,B)
    @avxt for n in indices((B,C),2), m in indices((A,C),1)
        Cmn = zero(eltype(C))
        for k in indices((A,B),(2,1))
            Cmn += A[m,k] * B[k,n]
        end
        C[m,n] = Cmn
    end
    C
end
function dot3(x::AbstractVector{Complex{T}}, A::AbstractMatrix{Complex{T}}, y::AbstractVector{Complex{T}}) where {T}
    xr = reinterpret(reshape, T, x);
    yr = reinterpret(reshape, T, y);
    Ar = reinterpret(reshape, T, A);
    sre = zero(T)
    sim = zero(T)
    @avxt for n in axes(Ar,3)
        tre = zero(T)
        tim = zero(T)
        for m in axes(Ar,2)
            tre += xr[1,m] * Ar[1,m,n] + xr[2,m] * Ar[2,m,n]
            tim += xr[1,m] * Ar[2,m,n] - xr[2,m] * Ar[1,m,n]
        end
        sre += tre * yr[1,n] - tim * yr[2,n]
        sim += tre * yr[2,n] + tim * yr[1,n]
    end
    sre + im*sim
end
function conv!(out, A, kern)
    @avxt for n ∈ axes(out,2), m ∈ axes(out,1)
        tmp = zero(eltype(out))
        for nᵢ ∈ -1:1, mᵢ ∈ -1:1
            tmp += A[mᵢ+m, nᵢ+n]*kern[mᵢ, nᵢ]
        end
        out[m,n] = tmp
    end
    out
end
function conv_baseline!(out, A, kern)
    for n ∈ axes(out,2), m ∈ axes(out,1)
        tmp = zero(eltype(out))
        for nᵢ ∈ -1:1, mᵢ ∈ -1:1
            tmp += A[mᵢ+m, nᵢ+n]*kern[mᵢ, nᵢ]
        end
        out[m,n] = tmp
    end
    out
end

struct DenseConvDims{N,K,C_in,C_out} end

function kernaxes(::DenseConvDims{2,K,C_in, C_out}) where {K,C_in, C_out}
    K₁ =  LoopVectorization.StaticInt(1):LoopVectorization.StaticInt(K[1])
    K₂ =  LoopVectorization.StaticInt(1):LoopVectorization.StaticInt(K[2])
    Cᵢₙ =  LoopVectorization.StaticInt(1):LoopVectorization.StaticInt(C_in)
    Cₒᵤₜ = LoopVectorization.StaticInt(1):LoopVectorization.StaticInt(C_out)
    (K₁, K₂, Cᵢₙ, Cₒᵤₜ)
end

function convlayer!(
    out::AbstractArray{<:Any,4}, img, kern,
    dcd::DenseConvDims{2, <:Any, <:Any, <:Any}
)
    (K₁, K₂, Cᵢₙ, Cₒᵤₜ) = kernaxes(dcd)
    @avxt for j₁ ∈ axes(out,1), j₂ ∈ axes(out,2), d ∈ axes(out,4), o ∈ Cₒᵤₜ
        s = zero(eltype(out))
        for k₁ ∈ K₁, k₂ ∈ K₂, i ∈ Cᵢₙ
            s += img[j₁ + k₁ - 1, j₂ + k₂ - 1, i, d] * kern[k₁, k₂, i, o]
        end
        out[j₁, j₂, o, d] = s
    end
    out
end
function convlayer_direct!(
    out::AbstractArray{<:Any,4}, img, kern,
    dcd::DenseConvDims{2, <:Any, <:Any, <:Any}
)
    (K₁, K₂, Cᵢₙ, Cₒᵤₜ) = kernaxes(dcd)
    @inbounds @fastmath for j₁ ∈ axes(out,1), j₂ ∈ axes(out,2), d ∈ axes(out,4), o ∈ Cₒᵤₜ
        s = zero(eltype(out))
        for k₁ ∈ K₁, k₂ ∈ K₂, i ∈ Cᵢₙ
            s += img[j₁ + k₁ - 1, j₂ + k₂ - 1, i, d] * kern[k₁, k₂, i, o]
        end
        out[j₁, j₂, o, d] = s
    end
    out
end

@testset "Threading" begin
    dcd = DenseConvDims{2,(5,5),3,6}()
    kern4 = rand(Float32, 5, 5, 3, 6);
    for M ∈ 17:50:267
        img = rand(Float32, M, M, 3, 100);
        outimage1 = Array{Float32}(undef, size(img,1)+1-size(kern4,1), size(img,2)+1-size(kern4,2), size(kern4,4), size(img,4));
        outimage2 = similar(outimage1);

        convlayer!(outimage1, img, kern4, dcd);
        convlayer_direct!(outimage2, img, kern4, dcd);
        @test outimage1 ≈ outimage2
    end

    for M ∈ 17:399
        # @show M
        K = M; N = M;
        A = rand(M,K); B = rand(K,N);
        @test dot(A,B) ≈ mydotavx(A,B)

        C1 = A * B; C0 = similar(C1);
        @test AmulB!(C0, A, B) ≈ C1

        x = randn(Complex{Float64}, 3M-1);
        W = randn(Complex{Float64}, 3M-1, 3M+1);
        y = randn(Complex{Float64}, 3M+1);
        @test dot(x,W,y) ≈ dot3(x,W,y)

        kern = OffsetArray(randn(3,3),-2,-2)
        out1 = OffsetArray(randn(size(A) .- 2), 1, 1)
        out2 = similar(out1);
        @test conv!(out1, A, kern) ≈ conv_baseline!(out2, A, kern)
    end
end


text:
struct GCNConv{T,F}
    weight::AbstractMatrix{T}
    bias::AbstractMatrix{T}
    norm::AbstractMatrix{T}
    σ::F
end

function GCNConv(adj::AbstractMatrix, ch::Pair{<:Integer,<:Integer}, σ = identity;
                 init = glorot_uniform, T::DataType=Float32, bias::Bool=true)
    N = size(adj, 1)
    b = bias ? param(init(ch[2], N)) : zeros(T, ch[2], N)
    GCNConv(param(init(ch[2], ch[1])), b, normalized_laplacian(adj+I, T), σ)
end

@treelike GCNConv

(g::GCNConv)(X::AbstractMatrix) = g.σ.(g.weight * X * g.norm + g.bias)

struct ChebConv{T}
    weight::AbstractArray{T,3}
    bias::AbstractMatrix{T}
    L̃::AbstractMatrix{T}
    k::Integer
    in_channel::Integer
    out_channel::Integer
end

function ChebConv(adj::AbstractMatrix, ch::Pair{<:Integer,<:Integer}, k::Integer;
                  init = glorot_uniform, T::DataType=Float32, bias::Bool=true)
    N = size(adj, 1)
    b = bias ? param(init(ch[2], N)) : zeros(T, ch[2], N)
    L̃ = T(2. / eigmax(adj)) * normalized_laplacian(adj, T) - I
    ChebConv(param(init(ch[2], ch[1], k)), b, L̃, k, ch[1], ch[2])
end

@treelike ChebConv

function (c::ChebConv)(X::AbstractMatrix{T}) where {T<:Real}
    fin = c.in_channel
    @assert size(X, 1) == fin "Input feature size must match input channel size."
    N = size(c.L̃, 1)
    @assert size(X, 2) == N "Input vertex number must match Laplacian matrix size."
    fout = c.out_channel

    Z = Array{T}(undef, fin, N, c.k)
    Z[:,:,1] = X
    Z[:,:,2] = X * c.L̃
    for k = 3:c.k
        Z[:,:,k] = 2*view(Z, :, :, k-1)*c.L̃ - view(Z, :, :, k-2)
    end

    Y = view(c.weight, :, :, 1) * view(Z, :, :, 1)
    for k = 2:c.k
        Y += view(c.weight, :, :, k) * view(Z, :, :, k)
    end
    Y += c.bias
    return Y
end

struct GraphConv{V,T} <: MessagePassing
    adjlist::V
    weight1::AbstractMatrix{T}
    weight2::AbstractMatrix{T}
    bias::AbstractMatrix{T}
    aggr::Symbol
end

function GraphConv(el::AbstractVector{<:AbstractVector{<:Integer}},
                   ch::Pair{<:Integer,<:Integer}, aggr=:add;
                   init = glorot_uniform, bias::Bool=true)
    N = size(el, 1)
    b = bias ? param(init(ch[2], N)) : zeros(T, ch[2], N)
    GraphConv(el, param(init(ch[2], ch[1])), param(init(ch[2], ch[1])), b, aggr)
end

function GraphConv(adj::AbstractMatrix, ch::Pair{<:Integer,<:Integer}, aggr=:add;
                   init = glorot_uniform, bias::Bool=true)
    N = size(adj, 1)
    b = bias ? param(init(ch[2], N)) : zeros(T, ch[2], N)
    GraphConv(neighbors(adj), param(init(ch[2], ch[1])), param(init(ch[2], ch[1])), b, aggr)
end

@treelike GraphConv

message(g::GraphConv; x_i=zeros(0), x_j=zeros(0)) = g.weight2 * x_j
update(g::GraphConv; X=zeros(0), M=zeros(0)) = g.weight1*X + M + g.bias
(g::GraphConv)(X::AbstractMatrix) = propagate(g, X=X, aggr=:add)

struct GATConv{V,T} <: MessagePassing
    adjlist::V
    weight::AbstractMatrix{T}
    bias::AbstractMatrix{T}
    a::AbstractArray
    negative_slope::Real
end

function GATConv(adj::AbstractMatrix, ch::Pair{<:Integer,<:Integer}; heads=1,
                 concat=true, negative_slope=0.2, init=glorot_uniform, bias::Bool=true)
    N = size(adj, 1)
    b = bias ? param(init(ch[2], N)) : zeros(T, ch[2], N)
    GATConv(neighbors(adj), param(init(ch[2], ch[1])), b, param(init(2 * ch[2])), negative_slope)
end

@treelike GATConv

function message(g::GATConv; x_i=zeros(0), x_j=zeros(0))
    n = size(x_j, 2)
    α = leakyrelu.(g.a' * vcat(repeat(x_i, outer=(1,n)), x_j), g.negative_slope)
    α = asoftmax(α)
    α .* x_j
end
update(g::GATConv; X=zeros(0), M=zeros(0)) = M + g.bias
(g::GATConv)(X::AbstractMatrix) = propagate(g, X=g.weight*X, aggr=:add)

function asoftmax(xs)
    xs = [exp.(x) for x in xs]
    s = sum(xs)
    return [x ./ s for x in xs]
end

struct GatedGraphConv{V,T,R} <: MessagePassing
    adjlist::V
    weight::AbstractArray{T}
    gru::R
    out_ch::Integer
    num_layers::Integer
    aggr::Symbol
end

function GatedGraphConv(adj::AbstractMatrix, out_ch::Integer, num_layers::Integer;
                        aggr=:add, init=glorot_uniform)
    N = size(adj, 1)
    w = param(init(out_ch, out_ch, num_layers))
    gru = GRUCell(out_ch, out_ch)
    GatedGraphConv(neighbors(adj), w, gru, out_ch, num_layers, aggr)
end

@treelike GatedGraphConv

message(g::GatedGraphConv; x_i=zeros(0), x_j=zeros(0)) = x_j
update(g::GatedGraphConv; X=zeros(0), M=zeros(0)) = M
function (g::GatedGraphConv)(X::AbstractMatrix{T}) where {T<:Real}
    H = X
    m, n = size(H)
    @assert (m <= g.out_ch) "number of input features must less or equals to output features."
    (m < g.out_ch) && (H = vcat(H, zeros(T, g.out_ch - m, n)))

    for i = 1:g.num_layers
        M = view(g.weight, :, :, i) * H
        M = propagate(g, X=M, aggr=g.aggr)
        H, _ = g.gru(H, M)
    end
    H
end

struct EdgeConv{V} <: MessagePassing
    adjlist::V
    nn
    aggr::Symbol
end

function EdgeConv(adj::AbstractMatrix, nn; aggr::Symbol=:max)
    EdgeConv(neighbors(adj), nn, aggr)
end

@treelike EdgeConv

function message(e::EdgeConv; x_i=zeros(0), x_j=zeros(0))
    n = size(x_j, 2)
    e.nn(vcat(repeat(x_i, outer=(1,n)), x_j .- x_i))
end
update(e::EdgeConv; X=zeros(0), M=zeros(0)) = M
(e::EdgeConv)(X::AbstractMatrix) = propagate(e, X=X, aggr=e.aggr)