"""
cancel out the basic system load and replace it with a convex combination of bids

note: the demand_module (or some subsidiary module) may store calibration data
at the module level (not in the model), so this module should only be used with one 
model at a time. An alternative approach would be to receive a calibration_data 
object back from demand_module.calibrate(), then add that to the model and pass
it back to the bid function when needed.

note: we also take advantage of this assumption and store a reference to the
current demand_module in this module (rather than storing it in the model itself)
"""

import os, sys
from pprint import pprint
from pyomo.environ import *
import switch_mod.utilities as utilities
demand_module = None    # will be set via command-line options

import util
from util import get

def define_arguments(argparser):
    argparser.add_argument("--dr_flat_pricing", action='store_true', default=False,
        help="Charge a constant (average) price for electricity, rather than varying hour by hour")
    argparser.add_argument("--dr_total_cost_pricing", action='store_true', default=False,
        help="Include both marginal and non-marginal(fixed) costs when setting prices")
    argparser.add_argument("--dr_elasticity_scenario", type=int, default=3,
        help="Choose a scenario of customer elasticity (1-3), defined in the demand_module")
    argparser.add_argument("--dr_demand_module", default=None,
        help="Name of module to use for demand-response bids. This should also be "
        "specified in the modules list, and should provide calibrate() and bid() functions. "
        "Pre-written options include constant_elasticity_demand_system or r_demand_system. "
        "Specify one of these in the modules list and use --help again to see module-specific options.")

def define_components(m):

    ###################
    # Choose the right demand module.
    # NOTE: we assume only one model will be run at a time, so it's safe to store
    # the setting in this module instead of in the model.
    ##################

    global demand_module
    if m.options.dr_demand_module is None:
        raise RuntimeError(
            "No demand module was specified for the demand_response system; unable to continue. "
            "Please use --dr_demand_module <module_name> in options.txt, scenarios.txt or on "
            "the command line. "
            "You should also add this module to the list of modules to load "
            " via modules.txt or --include_module <module_name>."
        )
    if m.options.dr_demand_module not in sys.modules:
        raise RuntimeError(
            "Demand module {mod} cannot be used because it has not been loaded. "
            "Please add this module to the the modules list (usually modules.txt) "
            "or specify --include_module {mod} in options.txt or on the command line."
            "".format(mod=m.options.dr_demand_module)
        )
    demand_module = sys.modules[m.options.dr_demand_module]
    
    # Make sure the model has a dual suffix
    if not hasattr(m, "dual"):
        m.dual = Suffix(direction=Suffix.IMPORT)
    
    ###################
    # Unserved load, with a penalty.
    # to ensure the model is always feasible, no matter what demand bids we get
    ##################

    # cost per MWh for unserved load (high)
    m.dr_unserved_load_penalty_per_mwh = Param(initialize=10000)
    # amount of unserved load during each timepoint
    m.DRUnservedLoad = Var(m.LOAD_ZONES, m.TIMEPOINTS, within=NonNegativeReals)
    # total cost for unserved load
    m.DR_Unserved_Load_Penalty = Expression(m.TIMEPOINTS, rule=lambda m, tp:
        sum(m.DRUnservedLoad[lz, tp] * m.dr_unserved_load_penalty_per_mwh for lz in m.LOAD_ZONES)
    )
    # add the unserved load to the model's energy balance
    m.LZ_Energy_Components_Produce.append('DRUnservedLoad')
    # add the unserved load penalty to the model's objective function
    m.cost_components_tp.append('DR_Unserved_Load_Penalty')

    ###################
    # Price Responsive Demand bids
    ##################
    
    # list of all bids that have been received from the demand system
    m.DR_BID_LIST = Set(initialize = [], ordered=True)
    # we need an explicit indexing set for everything that depends on DR_BID_LIST
    # so we can reconstruct it (and them) each time we add an element to DR_BID_LIST
    # (not needed, and actually doesn't work -- reconstruct() fails for sets)
    # m.DR_BIDS_LZ_TP = Set(initialize = lambda m: m.DR_BID_LIST * m.LOAD_ZONES * m.TIMEPOINTS)
    # m.DR_BIDS_LZ_TS = Set(initialize = lambda m: m.DR_BID_LIST * m.LOAD_ZONES * m.TIMESERIES)
    
    # data for the individual bids; each load_zone gets one bid for each timeseries,
    # and each bid covers all the timepoints in that timeseries. So we just record 
    # the bid for each timepoint for each load_zone.
    m.dr_bid = Param(m.DR_BID_LIST, m.LOAD_ZONES, m.TIMEPOINTS, mutable=True)

    # price used to get this bid (only kept for reference)
    m.dr_price = Param(m.DR_BID_LIST, m.LOAD_ZONES, m.TIMEPOINTS, mutable=True)

    # the private benefit of serving each bid
    m.dr_bid_benefit = Param(m.DR_BID_LIST, m.LOAD_ZONES, m.TIMESERIES, mutable=True)

    # weights to assign to the bids for each timeseries when constructing an optimal demand profile
    m.DRBidWeight = Var(m.DR_BID_LIST, m.LOAD_ZONES, m.TIMESERIES, within=NonNegativeReals)
    
    # def DR_Convex_Bid_Weight_rule(m, lz, ts):
    #     if len(m.DR_BID_LIST) == 0:
    #         print "no items in m.DR_BID_LIST, skipping DR_Convex_Bid_Weight constraint"
    #         return Constraint.Skip
    #     else:
    #         print "constructing DR_Convex_Bid_Weight constraint"
    #         return (sum(m.DRBidWeight[b, lz, ts] for b in m.DR_BID_LIST) == 1)
    # 
    # choose a convex combination of bids for each zone and timeseries
    m.DR_Convex_Bid_Weight = Constraint(m.LOAD_ZONES, m.TIMESERIES, rule=lambda m, lz, ts: 
        Constraint.Skip if len(m.DR_BID_LIST) == 0 
            else (sum(m.DRBidWeight[b, lz, ts] for b in m.DR_BID_LIST) == 1)
    )
    
    # Since we don't have differentiated prices for each zone, we have to use the same
    # weights for all zones. (Otherwise the model will try to micromanage load in each
    # zone, but that won't be reflected in the prices we report.)
    # Note: LOAD_ZONES is not an ordered set, so we have to use a trick to get a single
    # arbitrary one to refer to (next(iter(m.LOAD_ZONES)) would also work).
    m.DR_Load_Zone_Shared_Bid_Weight = Constraint(
        m.DR_BID_LIST, m.LOAD_ZONES, m.TIMESERIES, rule=lambda m, b, lz, ts: 
            m.DRBidWeight[b, lz, ts] == m.DRBidWeight[b, list(m.LOAD_ZONES)[0], ts]
    )

    # For flat-price models, we have to use the same weight for all timeseries within the
    # same year (period), because there is only one price for the whole period, so it can't
    # induce different adjustments in individual timeseries.
    if m.options.dr_flat_pricing:
        m.DR_Flat_Bid_Weight = Constraint(
            m.DR_BID_LIST, m.LOAD_ZONES, m.TIMESERIES, rule=lambda m, b, lz, ts: 
                m.DRBidWeight[b, lz, ts] 
                == m.DRBidWeight[b, lz, m.tp_ts[m.PERIOD_TPS[m.ts_period[ts]].first()]]
        )
                
    
    # Optimal level of demand, calculated from available bids (negative, indicating consumption)
    m.FlexibleDemand = Expression(m.LOAD_ZONES, m.TIMEPOINTS, 
        rule=lambda m, lz, tp:
            sum(m.DRBidWeight[b, lz, m.tp_ts[tp]] * m.dr_bid[b, lz, tp] for b in m.DR_BID_LIST)
    )

    # # FlexibleDemand reported as an adjustment (negative equals more demand)
    # # We have to do it this way because there's no way to remove the lz_demand_mw from the model
    # # without changing the core code.
    # m.DemandPriceResponse = Expression(m.LOAD_ZONES, m.TIMEPOINTS, 
    #     rule=lambda m, lz, tp: m.lz_demand_mw[lz, tp] - m.FlexibleDemand[lz, tp]
    # )
    
    # private benefit of the electricity consumption 
    # (i.e., willingness to pay for the current electricity supply)
    # reported as negative cost, i.e., positive benefit
    # also divide by number of timepoints in the timeseries
    # to convert from a cost per timeseries to a cost per timepoint.
    m.DR_Welfare_Cost = Expression(m.TIMEPOINTS, rule=lambda m, tp:
        (-1.0) 
        * sum(m.DRBidWeight[b, lz, m.tp_ts[tp]] * m.dr_bid_benefit[b, lz, m.tp_ts[tp]] 
            for b in m.DR_BID_LIST for lz in m.LOAD_ZONES) 
        * m.tp_duration_hrs[tp] / m.ts_num_tps[m.tp_ts[tp]]
    )

    # add the private benefit to the model's objective function
    m.cost_components_tp.append('DR_Welfare_Cost')

    # annual costs, recovered via baseline prices
    # but not included in switch's calculation of costs
    m.other_costs = Param(m.PERIODS, mutable=True, default=0.0)
    m.cost_components_annual.append('other_costs')
    
    # variable to store the baseline data
    m.base_data = None

def post_iterate(m):
    print "\n\n======================================================="
    print "Solved model"
    print "======================================================="
    print "Total cost: ${v:,.0f}".format(v=value(m.SystemCost))
    print "marginal costs (first day):"
    print [
        electricity_marginal_cost(m, lz, tp) 
            for lz in m.LOAD_ZONES
                for tp in m.TS_TPS[m.TIMESERIES[1]]
    ]
    print "marginal costs (second day):"
    print [
        electricity_marginal_cost(m, lz, tp) 
            for lz in m.LOAD_ZONES
                for tp in m.TS_TPS[m.TIMESERIES[2]]
    ]

    # if m.iteration_number % 5 == 0:
    #     # save time by only writing results every 5 iterations
    # write_results(m)

    # Retrieve SystemCost before calling update_demand()
    # because that breaks SystemCost until the next solve
    old_SystemCost = getattr(m, "last_SystemCost", None)
    new_SystemCost = value(m.SystemCost)
    m.last_SystemCost = new_SystemCost
    if m.iteration_number > 0:
        # store cost of current solution before it gets altered by update_demand()
        current_cost = value(sum(
            (
                sum(
                    electricity_marginal_cost(m, lz, tp) * electricity_demand(m, lz, tp)
                        for lz in m.LOAD_ZONES 
                ) + m.DR_Welfare_Cost[tp]
            ) * m.bring_timepoint_costs_to_base_year[tp]
                for ts in m.TIMESERIES
                    for tp in m.TS_TPS[ts]
        ))
    
    update_demand(m)

    if m.iteration_number > 0:
        # get an estimate of best possible net cost of serving load
        # (if we could completely serve the last bid at the prices we quoted,
        # that would be an optimum; the actual cost may be higher but never lower)
        b = m.DR_BID_LIST.last()
        best_cost = value(sum(
            sum(
                electricity_marginal_cost(m, lz, tp) * m.dr_bid[b, lz, tp] 
                - m.dr_bid_benefit[b, lz, ts] * m.tp_duration_hrs[tp] / m.ts_num_tps[ts]
                    for lz in m.LOAD_ZONES 
            ) * m.bring_timepoint_costs_to_base_year[tp]
                for ts in m.TIMESERIES
                    for tp in m.TS_TPS[ts]
        ))
        print "last_SystemCost={}, SystemCost={}, ratio={}".format(
            old_SystemCost, new_SystemCost, new_SystemCost/old_SystemCost)
        print "lower bound={}, current cost={}, ratio={}".format(
            best_cost, current_cost, current_cost/best_cost)
        print "discount factors: " + " ".join([
            "{}={}".format(p, m.bring_timepoint_costs_to_base_year[m.PERIOD_TPS[p].first()])
                for p in m.PERIODS
        ])

    # Check for convergence (no progress during the last iteration)
    converged = (m.iteration_number > 0 and new_SystemCost == old_SystemCost)
        
    return converged

def update_demand(m):
    """
    This should be called after solving the model, in order to calculate new bids
    to include in future runs. The first time through, it also uses the fixed demand
    and marginal costs to calibrate the demand system, and then replaces the fixed
    demand with the flexible demand system.
    """
    first_run = (m.base_data is None)
    outputs_dir = m.options.outputs_dir
    tag = m.options.scenario_name

    print "attaching new demand bid to model"
    if first_run:
        calibrate_model(m)

        util.create_table(
            output_file=os.path.join(outputs_dir, "bid_weights_{t}.tsv".format(t=tag)), 
            headings=("iteration", "load_zone", "timeseries", "bid_num", "weight")
        )
    else:   # not first run
        # print "m.DRBidWeight (first day):"
        # print [(b, lz, ts, value(m.DRBidWeight[b, lz, ts])) 
        #     for b in m.DR_BID_LIST
        #     for lz in m.LOAD_ZONES
        #     for ts in m.TIMESERIES]
        print "m.DRBidWeight:"
        pprint([(lz, ts, [(b, value(m.DRBidWeight[b, lz, ts])) for b in m.DR_BID_LIST])
            for lz in m.LOAD_ZONES
            for ts in m.TIMESERIES])
        #print "DR_Convex_Bid_Weight:"
        #m.DR_Convex_Bid_Weight.pprint()

        # store the current bid weights for future reference
        # This should be done before adding the new bid.
        util.append_table(m, m.LOAD_ZONES, m.TIMESERIES, m.DR_BID_LIST, 
            output_file=os.path.join(outputs_dir, "bid_weights_{t}.tsv".format(t=tag)), 
            values=lambda m, lz, ts, b: (len(m.DR_BID_LIST), lz, ts, b, m.DRBidWeight[b, lz, ts])
        )

    # get new bids from the demand system at the current prices
    bids = get_bids(m)
    
    print "adding bids to model"
    # print "first day (lz, ts, prices, demand, wtp) ="
    # pprint(bids[0])
    # add the new bids to the model
    add_bids(m, bids)
    print "m.dr_bid_benefit (first day):"
    pprint([(b, lz, ts, value(m.dr_bid_benefit[b, lz, ts])) 
        for b in m.DR_BID_LIST
        for lz in m.LOAD_ZONES
        for ts in [m.TIMESERIES.first()]])
    
    # print "m.dr_bid (first day):"
    # print [(b, lz, ts, value(m.dr_bid[b, lz, ts]))
    #     for b in m.DR_BID_LIST
    #     for lz in m.LOAD_ZONES 
    #     for ts in m.TS_TPS[m.TIMESERIES.first()]]
    
    if first_run:
        # replace lz_demand_mw with FlexibleDemand in the energy balance constraint
        # note: it is easiest to do this after retrieving the bids because this
        # destroys the dual values which are needed for calculating the bids
        # note: the first two lines are simpler than the method I use, but my approach
        # preserves the ordering of the list, which is nice for reporting.
        # m.LZ_Energy_Components_Consume.remove('lz_demand_mw')
        # m.LZ_Energy_Components_Consume.append('FlexibleDemand')
        ecc = m.LZ_Energy_Components_Consume
        ecc[ecc.index('lz_demand_mw')] = 'FlexibleDemand'
        reconstruct_energy_balance(m)

def sum_product(vector1, vector2):
    return sum(v1*v2 for (v1, v2) in zip(vector1, vector2))

def total_direct_costs_per_year(m, period):
    """Return undiscounted total cost per year, during each period, as calculated by SWITCH,
    including everything except DR_Welfare_Cost.

    note: during the first iteration, this doesn't include "other_costs", but in later 
    iterations it does.

    This code comes from financials.calc_sys_costs_per_period(), excluding discounting
    and upscaling to the period.
    
    NOTE: ideally this would give costs by zone and period, to allow calibration for different
    utilities within a large study. But the cost components don't distinguish that way.
    (They probably should, since that would allow us to discuss average electricity costs
    in each zone.)
    """
    return value(
        sum(getattr(m, annual_cost)[period] for annual_cost in m.cost_components_annual)
        + sum(
            getattr(m, tp_cost)[t] * m.tp_weight_in_year[t]
            for t in m.PERIOD_TPS[period]
                for tp_cost in m.cost_components_tp
                    if tp_cost != "DR_Welfare_Cost"
        )
    )    

def electricity_marginal_cost(m, lz, tp):
    """Return marginal cost of production per MWh in load_zone lz during timepoint tp."""
    return m.dual[m.Energy_Balance[lz, tp]]/m.bring_timepoint_costs_to_base_year[tp]

def electricity_demand(m, lz, tp):
    """Return total electricity consumption by customers in load_zone lz during timepoint tp."""
    return value(sum(
        getattr(m, component)[lz, tp]
            for component in ('lz_demand_mw', 'FlexibleDemand')
                if component in m.LZ_Energy_Components_Consume
    ))
    
def make_prices(m):
    """Calculate hourly prices for customers, based on the current model configuration.
    These may be any combination of marginal vs. total-cost and flat vs. dynamic.
    """

    if m.options.dr_total_cost_pricing:
        # rescale (long-run) marginal costs to recover all costs 
        # (sunk and other, in addition to marginal)
        # calculate the ratio between potential revenue 
        # at marginal-cost pricing and total costs for each period
        mc_annual_revenue = {
            (lz, p): 
            sum(
                electricity_demand(m, lz, tp) 
                * electricity_marginal_cost(m, lz, tp) 
                * m.tp_weight_in_year[tp]
                for tp in m.PERIOD_TPS[p]
            )
            for lz in m.LOAD_ZONES for p in m.PERIODS
        }
        # note: it would be nice to do this on a zonal basis, but production costs
        # are only available model-wide.
        price_scalar = {
            p: total_direct_costs_per_year(m, p) 
                / sum(mc_annual_revenue[lz, p] for lz in m.LOAD_ZONES) 
            for p in m.PERIODS
        }
    else:
        # use marginal costs directly as prices
        price_scalar = {p: 1.0 for p in m.PERIODS}
        
    # calculate hourly prices
    hourly_prices = {
        (lz, tp): price_scalar[m.tp_period[tp]] * electricity_marginal_cost(m, lz, tp)
            for lz in m.LOAD_ZONES for tp in m.TIMEPOINTS
    }
    
    if m.options.dr_flat_pricing:
        # use flat prices each year
        # calculate annual average prices (total revenue / total kWh)
        average_prices = {
            (lz, p): 
            sum(
                hourly_prices[lz, tp] 
                * electricity_demand(m, lz, tp) 
                * m.tp_weight_in_year[tp] 
                for tp in m.PERIOD_TPS[p]
            ) 
            / 
            sum(
                electricity_demand(m, lz, tp) 
                * m.tp_weight_in_year[tp] 
                for tp in m.PERIOD_TPS[p]
            )
            for lz in m.LOAD_ZONES for p in m.PERIODS
        }
        prices = {
            (lz, tp): average_prices[lz, m.tp_period[tp]]
            for lz in m.LOAD_ZONES for tp in m.TIMEPOINTS
        }