rule_fcst_par_conc_futrs_notdone.py

Go to the documentation of this file.

# -*- coding: utf-8 -*-
"""
Created on Mon May 18 19:48:52 2020
 
@author: cd
"""
 
 
import abc
import numpy as np
from efo.rule import RuleBase, RuleRlsSpecified
from efo.junction import ReservoirJunction
# import multiprocessing as mp
# from multiprocessing import Pool as p
import concurrent.futures
# from collections import deque
 
# TODO: Numpy array variable should be passed as reference or pointer
# TODO: Function arguments should be labeled "ref" or "ptr" at beginning for clarification
# TODO: Function argments that are array copies should be labeled "cpy"
 
 
class RuleFcstBase(RuleBase):
    def __init__(self, name, timeFcst):
        # Call super class constructor
        super().__init__(name, timeFcst)
        # self.T = timeFcst
        self.releaserelease = np.full(self.T.Tcont.nSteps, np.nan)
        self.releaseSchedreleaseSched = np.empty(self.T.nSteps)
 
    def calc_release(self, rlsProposed, rlsUnCtrl, stor, rlsPrev=None, qIn=None):
        rlsProposed = super().calc_release(
            rlsProposed, rlsPrev=rlsPrev, rlsUnCtrl=rlsUnCtrl, stor=stor, qIn=qIn)
        # TODO: Check to see if the current date is a forecast date
        self.releaseSchedreleaseSched = self.calc_release_fcstcalc_release_fcst(rlsProposed, rlsUnCtrl, stor, qIn)
        return self.releaseSchedreleaseSched[1]
 
    @abc.abstractmethod
    def calc_release_fcst(self, rlsProposed, rlsUnCtrl, stor, qIn):
        # TODO: This should return a tuple that has the release for the current time step
        # but also the release schedule and possilby the max and min?
        # TODO: In the world of OOP this should be a protected method because we only
        # want to call this from other forecast rules, but in python not sure.
        pass
 
    @classmethod
    def __subclasshook__(cls, C):
        if cls is RuleFcstBase:
            attrs = set(dir(C))
            if set(cls.__abstractmethods__) <= attrs:
                return True
        return NotImplemented
 
 
class RuleFcstEfo(RuleFcstBase):
    def __init__(self, name, timeFcst, efoResrvr, storTh, riskTol, constants, *,
                  efoNetwork=None, ruleRlsSpecified=None, efoRlsSched=None):
        # Call super class constructor
        super().__init__(name, timeFcst)
        # Set sub class properties
        self.efoNetworkefoNetwork = efoNetwork
        self.efoResrvrefoResrvr = efoResrvr
        self.ruleRlsSpecifiedruleRlsSpecified = ruleRlsSpecified
        self.efoRlsSchedefoRlsSched = efoRlsSched
        self.nMembersnMembers = len(efoResrvr)
        # TODO: This shoud be a storage guide curve rule for flexibility
        self.storThstorTh = float(storTh)
        self.riskTolriskTol = riskTol
        self.constantsconstants = constants
        self.concTaskerconcTasker = ConcurrentTasker()
 
    def calc_release_fcst(self, rlsProposed, rlsUnCtrl, stor, qIn=None):
        # Set initial conditions of the ensemble networks to the current observed
        # TODO: Need to create efo network that has the primary network and the ensembles
        # TODO: Should this have max and min release properties for easier access for results?
        if self.efoNetworkefoNetwork is not None: 
            # TODO: Here it should just return a copy with the data for today's forecast
            # TODO: This will avoid pickling all of the data in parallel processing
            # TODO: May need to add a copy method to allow this
            for net in self.efoNetworkefoNetwork: net.set_init_cond()
        # Initialize variables
        t = self.T.get_fcst_time()
        # self.storFcst = np.empty(t.nSteps, self.nMembers, t.nSteps - 1)
        # self.storFcst[:] = np.nan
        self.storFcststorFcst = np.full((t.nSteps, self.nMembersnMembers, t.nSteps-1), np.nan)
        # self.storFcst[0, :, :] = [curRes.stor[0] for curRes in self.efoResrvr]
        self.storFcststorFcst[0, :, :] = self.efoResrvrefoResrvr[0].stor[0]
        # self.rlsFcst = np.full((t.nSteps, self.nMembers, t.nSteps - 1), np.nan)
        self.rlsUnCtrlFcstrlsUnCtrlFcst = np.full((t.nSteps, self.nMembersnMembers, t.nSteps-1), np.nan)
        self.rlsUnCtrlFcstrlsUnCtrlFcst[0, :, :] = self.efoResrvrefoResrvr[0].rlsUnCtrl[0]
        self.rlsFcstrlsFcst = np.zeros((t.nSteps, self.nMembersnMembers, t.nSteps-1))
        for i in range(t.nSteps-1):
            self.rlsFcstrlsFcst[i+2:, 0:self.nMembersnMembers, i] = np.nan
        # self.rlsFcst[0, :, :] = [curRes.rlsCtrl[0] for curRes in self.efoResrvr]
        # TODO: Maybe use the passed rlsPrev here?
        self.rlsFcstrlsFcst[0, :, :] = self.efoResrvrefoResrvr[0].rlsCtrl[0]
        self.rlsMaxrlsMax = np.full((t.nSteps, self.nMembersnMembers, t.nSteps-1), np.nan)
        self.rlsNoConstraintrlsNoConstraint = np.copy(self.rlsFcstrlsFcst)
        self.prExcThprExcTh = np.full(t.nSteps, np.nan)
        self.prExcThPreRlsprExcThPreRls = np.copy(self.prExcThprExcTh)
        self.iRiskyMbrsiRiskyMbrs = np.full((t.nSteps, self.nMembersnMembers), False)
        riskySelect = np.arange(0, self.nMembersnMembers, dtype=int)
        riskyMbr = np.copy(riskySelect)
        iRisky = np.full(self.nMembersnMembers, False)
        iRiskyChk = np.full(self.nMembersnMembers, True)
        volAbvThresh = np.zeros([t.nSteps, self.nMembersnMembers])
        vol2Rls = np.zeros([t.nSteps, self.nMembersnMembers])
        rlsFcst = np.full(t.nSteps, np.nan)
        rlsFcst[0] = self.efoResrvrefoResrvr[0].rlsCtrl[0]
        rlsMax = np.full(t.nSteps, np.nan)
        rlsToday = np.empty((t.nSteps-1, self.nMembersnMembers))
        rlsSched = np.zeros(t.nSteps)
        for tsHoriz in t.steps[1:]:
            self.T.step = tsHoriz
            if tsHoriz > 1:
                self.storFcststorFcst[:, :, tsHoriz-1] = self.storFcststorFcst[:, :, tsHoriz-2]
                self.rlsUnCtrlFcstrlsUnCtrlFcst[:, :, tsHoriz-1] = self.rlsUnCtrlFcstrlsUnCtrlFcst[:, :, tsHoriz-2]
            future2ensMbr = dict()
            with concurrent.futures.ProcessPoolExecutor(max_workers=6) as executor:
                for j in range(self.nMembersnMembers):
                    # TODO: Maybe shouldn't set the override to 0 here?
                    if issubclass(type(self.ruleRlsSpecifiedruleRlsSpecified[riskySelect[j]]), RuleRlsSpecified):
                        self.ruleRlsSpecifiedruleRlsSpecified[riskySelect[j]].set_release(tsHoriz, 0)
                    # self.futures.append(executor.submit(self.efoResrvr[j].calc_qout()))
                    # futures.append(executor.submit(self.efoResrvr[j].calc_qout()))
                    future2ensMbr.update( {executor.submit(self.efoResrvrefoResrvr[j].calc_qout()): j} )
                    # self.efoResrvr[j].calc_qout()
                    # TODO: In a network this should just process up to the efoReseroir...
                    # TODO: and not points downstream to save compute time
                for curFuture in concurrent.futures.as_completed(future2ensMbr):
                    curMbr = future2ensMbr[curFuture]
                    self.storFcststorFcst[tsHoriz, riskySelect[curMbr], tsHoriz-1] = self.efoResrvrefoResrvr[curMbr].stor[tsHoriz]
                    self.rlsUnCtrlFcstrlsUnCtrlFcst[tsHoriz, riskySelect[curMbr], tsHoriz-1] = \
                        self.efoResrvrefoResrvr[curMbr].rlsUnCtrl[tsHoriz]
                # while futures:
                #     future = futures.popleft()
                #     if future.exception():
                #         continue
                #     elif future.done():
                #         self.storFcst[tsHoriz, riskySelect[j], tsHoriz-1] = self.efoResrvr[j].stor[tsHoriz]
                #         self.rlsUnCtrlFcst[tsHoriz, riskySelect[j], tsHoriz-1] = \
                #             self.efoResrvr[j].rlsUnCtrl[tsHoriz]
            risky = True
            while risky:
                storPlusUnCtrl = self.storFcststorFcst[tsHoriz, :, tsHoriz-1]\
                    + self.rlsUnCtrlFcstrlsUnCtrlFcst[tsHoriz, :, tsHoriz-1]*self.constantsconstants.cfs2af
                iRisky = (storPlusUnCtrl > self.storThstorTh) & iRiskyChk
                self.prExcThprExcTh[tsHoriz] = np.sum(iRisky)/self.nMembersnMembers
                if self.prExcThprExcTh[tsHoriz] > self.riskTolriskTol[tsHoriz]:
                    self.prExcThPreRlsprExcThPreRls[tsHoriz] = self.prExcThprExcTh[tsHoriz]
                    iRiskyChk = np.copy(iRisky)
                    # riskyMbr = np.where(self.storFcst[tsHoriz, :, tsHoriz-1] > self.storTh)[0]
                    riskyMbr = np.where(iRisky)[0]
                    nMbrs2Reduce = int(len(riskyMbr) - np.floor(self.riskTolriskTol[tsHoriz]*self.nMembersnMembers))
                    volAbvThresh[tsHoriz, riskyMbr] = \
                        self.storFcststorFcst[tsHoriz, riskyMbr, tsHoriz-1]\
                            + sum(self.rlsUnCtrlFcstrlsUnCtrlFcst[1:tsHoriz+1, riskyMbr, tsHoriz-1])*self.constantsconstants.cfs2af\
                                - self.storThstorTh
                    mbrSorted = np.argsort(volAbvThresh[tsHoriz, riskyMbr])
                    riskySelect = riskyMbr[mbrSorted[0:nMbrs2Reduce]]
                    self.iRiskyMbrsiRiskyMbrs[tsHoriz, riskySelect] = True
                    # TODO: Should add a while loop so you can account for evaportion iterations
                    # TODO: This percent increase (1.01) should be a property
                    # vol2Rls[tsHoriz, riskySelect] = \
                    #     1.01*(self.storFcst[tsHoriz, riskySelect, tsHoriz-1] - self.storTh)
                    vol2Rls[tsHoriz, riskySelect] = 1.01*volAbvThresh[tsHoriz, riskySelect]
                    tasks = []
                    for j in range(len(riskySelect)):
                        rlsFcst[1:tsHoriz+1] = np.sum(
                                vol2Rls[:, riskySelect[j]])*self.constantsconstants.af2cfs/tsHoriz
                        self.rlsNoConstraintrlsNoConstraint[1:tsHoriz+1, riskySelect[j], tsHoriz-1] = np.copy(rlsFcst[1:tsHoriz+1])
                        rlsFcst[1:tsHoriz+1], rlsMax[1:tsHoriz+1] = self._calc_release_sched_calc_release_sched(
                            rlsFcst[1:tsHoriz+1], tsHoriz, riskySelect[j])
                        self.rlsMaxrlsMax[1:tsHoriz+1, riskySelect[j], tsHoriz-1] = np.copy(rlsMax[1:tsHoriz+1])
                        if np.sum(rlsFcst[1:tsHoriz+1])*self.constantsconstants.cfs2af < np.sum(vol2Rls[:, riskySelect[j]]):
                            iRiskyChk[riskySelect[j]] = False
                            vol2Rls[tsHoriz, riskySelect[j]] = \
                                vol2Rls[tsHoriz, riskySelect[j]] - \
                                (np.sum(vol2Rls[:, riskySelect[j]]) - np.sum(rlsFcst[1:tsHoriz+1])*self.constantsconstants.cfs2af)
                        self.rlsFcstrlsFcst[1:tsHoriz+1, riskySelect[j], tsHoriz-1] = np.copy(rlsFcst[1:tsHoriz+1])
                    risky = True
                else:
                    if np.isnan(self.prExcThPreRlsprExcThPreRls[tsHoriz]):
                        self.prExcThPreRlsprExcThPreRls[tsHoriz] = self.prExcThprExcTh[tsHoriz]
                    # TODO: Create a function that initializes these variables, would have to return...
                    # TODO: a tuple of multiple variables. Not sure this would save any code redundancy
                    riskySelect = np.arange(0, self.nMembersnMembers, dtype=int)
                    riskyMbr = np.copy(riskySelect)
                    iRiskyChk = np.full(self.nMembersnMembers, True)
                    risky = False
        rlsToday[:, :] = self.rlsFcstrlsFcst[1, :, :].T
        if np.any(rlsToday[:] > 0):
            # TODO: You could have a count of the number of ens mbrs and ts that have the same day 1 release
            # TODO: This could be used to provide multiple release schedules
            tsMax, mbrMax  = np.where(rlsToday==np.max(rlsToday))
            # TODO: Does this need to be set since you return the release schedule and it is reset by super class?
            self.releaserelease[self.T.Tcont.step] = self.rlsFcstrlsFcst[1, mbrMax[0], tsMax[0]]
            rlsSched = self.rlsFcstrlsFcst[:, mbrMax[0], tsMax[0]]
        else:
            # TODO: Does this need to be set since you return the release schedule and it is reset by super class?
            self.releaserelease[self.T.Tcont.step] = 0
        return rlsSched
    
    def _calc_release_sched(self, rlsFcst, tsHoriz, ensMbr):
        pctConvg = 0.01
        rlsTot = np.nansum(rlsFcst)
        maxItr = 100
        curItr = 0
        # Create reference variables to ensemble reservoir objects
        rlsMax = self.efoResrvrefoResrvr[ensMbr].rlsMax[1:tsHoriz+1]
        rlsApplied = self.efoResrvrefoResrvr[ensMbr].rlsCtrl[1:tsHoriz+1]
        rlsMaxPrev = np.zeros(len(rlsMax))
        rlsFcstInit = np.copy(rlsFcst)
        # idxZero = 1
        while True:
            # Recaculate storage
            for i in range(1, tsHoriz+1):
                if issubclass(type(self.efoResrvrefoResrvr[ensMbr]), ReservoirJunction):
                    self.T.step = i
                    if issubclass(type(self.ruleRlsSpecifiedruleRlsSpecified[ensMbr]), RuleRlsSpecified):
                        self.ruleRlsSpecifiedruleRlsSpecified[ensMbr].set_release(i, np.copy(rlsFcst[i-1]))
                    # TODO: In a network this should just process up to the efoReseroir...
                    # TODO: and not points downstream to save compute time
                    self.efoNetworkefoNetwork[ensMbr].process_fcst_junctions()
                    self.storFcststorFcst[i, ensMbr, tsHoriz-1] = np.copy(self.efoResrvrefoResrvr[ensMbr].stor[i])
            # TODO: If you are maintaining forecasted storage levels then all future timesteps...
            # TODO: should be at the max release schedule
            if rlsTot - np.sum(rlsApplied) < pctConvg*rlsTot\
                or np.all(np.absolute(rlsMaxPrev - rlsMax) < 0.01) or curItr > maxItr:
                break
            # zChk = np.array((rlsFcst, self.efoResrvr[ensMbr].rlsCtrl, self.storFcst[:, ensMbr, tsHoriz-1], rlsMax)).T
            rlsFcst = self._get_adjusted_release_eqldist_get_adjusted_release_eqldist(rlsFcstInit, rlsMax)
            rlsMaxPrev = np.copy(rlsMax)
            curItr += 1
        return rlsFcst, np.copy(rlsMax)
    
    def _get_adjusted_release_eqldist(self, rlsFcstInit, rlsMax):
        rlsFcst = np.copy(rlsFcstInit)
        rlsFcstCheck = np.zeros(len(rlsFcst))
        iVal = ~np.isnan(rlsFcst) & ~np.isnan(rlsMax)
        iMaxed = np.full(len(rlsFcst), False)
        iMaxed[0] = True
        # iMaxed = np.copy(iZeroStor)
        while any(iMaxed[1:]==False) and any(rlsFcst[iVal] - rlsFcstCheck[iVal] != 0):
            rlsFcstCheck = np.copy(rlsFcst)
            iGrtrMax = rlsFcst > rlsMax
            iMaxed = iMaxed | iGrtrMax
            rlsDiffVol = np.sum(rlsFcst[iGrtrMax] - rlsMax[iGrtrMax])
            rlsFcst[iGrtrMax] = rlsMax[iGrtrMax]
            rlsAdd = rlsDiffVol/max(1,np.sum(~iMaxed & iVal))
            rlsFcst[~iMaxed] += rlsAdd
        return rlsFcst
    
    def _get_adjusted_release_rollback(self, rlsFcstInit, tsHoriz, rlsMax):
        rlsFcst = np.copy(rlsFcstInit)
        rlsFcstCheck = np.copy(rlsFcst)
        carryOver = 0
        for i in range(tsHoriz, 0, -1):
            rlsFcst[i] += carryOver
            rlsFcstCheck[i] = rlsFcst[i]
            if rlsFcst[i] > rlsMax[i]:
                rlsFcst[i] = rlsMax[i]
            # Readjust carryover
            carryOver += (rlsFcstCheck[i] -  rlsFcst[i])/(i-1) if i > 1 else 0
        return rlsFcst
 
 
class ConcurrentTasker():
    def __init__(self):
        pass
    
    def _calc_release_sched(self, rlsFcst, tsHoriz, ensMbr):
        pctConvg = 0.01
        rlsTot = np.nansum(rlsFcst)
        maxItr = 100
        curItr = 0
        # Create reference variables to ensemble reservoir objects
        rlsMax = self.efoResrvr[ensMbr].rlsMax[1:tsHoriz+1]
        rlsApplied = self.efoResrvr[ensMbr].rlsCtrl[1:tsHoriz+1]
        rlsMaxPrev = np.zeros(len(rlsMax))
        rlsFcstInit = np.copy(rlsFcst)
        # idxZero = 1
        while True:
            # Recaculate storage
            for i in range(1, tsHoriz+1):
                if issubclass(type(self.efoResrvr[ensMbr]), ReservoirJunction):
                    self.T.step = i
                    if issubclass(type(self.ruleRlsSpecified[ensMbr]), RuleRlsSpecified):
                        self.ruleRlsSpecified[ensMbr].set_release(i, np.copy(rlsFcst[i-1]))
                    # TODO: In a network this should just process up to the efoReseroir...
                    # TODO: and not points downstream to save compute time
                    self.efoNetwork[ensMbr].process_fcst_junctions()
                    self.storFcst[i, ensMbr, tsHoriz-1] = np.copy(self.efoResrvr[ensMbr].stor[i])
            # TODO: If you are maintaining forecasted storage levels then all future timesteps...
            # TODO: should be at the max release schedule
            if rlsTot - np.sum(rlsApplied) < pctConvg*rlsTot\
                or np.all(np.absolute(rlsMaxPrev - rlsMax) < 0.01) or curItr > maxItr:
                break
            # zChk = np.array((rlsFcst, self.efoResrvr[ensMbr].rlsCtrl, self.storFcst[:, ensMbr, tsHoriz-1], rlsMax)).T
            rlsFcst = self._get_adjusted_release_eqldist(rlsFcstInit, rlsMax)
            rlsMaxPrev = np.copy(rlsMax)
            curItr += 1
        return rlsFcst, np.copy(rlsMax)
    
    def _get_adjusted_release_eqldist(self, rlsFcstInit, rlsMax):
        rlsFcst = np.copy(rlsFcstInit)
        rlsFcstCheck = np.zeros(len(rlsFcst))
        iVal = ~np.isnan(rlsFcst) & ~np.isnan(rlsMax)
        iMaxed = np.full(len(rlsFcst), False)
        iMaxed[0] = True
        # iMaxed = np.copy(iZeroStor)
        while any(iMaxed[1:]==False) and any(rlsFcst[iVal] - rlsFcstCheck[iVal] != 0):
            rlsFcstCheck = np.copy(rlsFcst)
            iGrtrMax = rlsFcst > rlsMax
            iMaxed = iMaxed | iGrtrMax
            rlsDiffVol = np.sum(rlsFcst[iGrtrMax] - rlsMax[iGrtrMax])
            rlsFcst[iGrtrMax] = rlsMax[iGrtrMax]
            rlsAdd = rlsDiffVol/max(1,np.sum(~iMaxed & iVal))
            rlsFcst[~iMaxed] += rlsAdd
        return rlsFcst
 
 
class RuleFcstPerfectEfoRisk(RuleFcstBase):
    def __init__(self, name, rulePfo, ruleEfoNoRls):
         # Call super class constructor
        super().__init__(name, rulePfo.T)
        self.rulePforulePfo = rulePfo
        self.ruleEfoNoRlsruleEfoNoRls = ruleEfoNoRls
        self.riskEforiskEfo = np.full((self.T.Tcont.nSteps, self.T.nSteps), 0.)
        
    def calc_release_fcst(self, rlsProposed, rlsUnCtrl, stor, qIn=None):
        zeroRls = self.ruleEfoNoRlsruleEfoNoRls.calc_release_fcst(rlsProposed, rlsUnCtrl, stor, qIn)
        self.riskEforiskEfo[self.T.Tcont.step,1:]  = self.ruleEfoNoRlsruleEfoNoRls.prExcTh[1:]
        rlsSched = self.rulePforulePfo.calc_release_fcst(rlsProposed, rlsUnCtrl, stor, qIn)
        # TODO: Does this need to be set since you return the release schedule and it is reset by super class?
        self.releaserelease[self.T.Tcont.step] = rlsSched[1]
        return rlsSched
 
 
# NOT DONE
# TODO: This should keep track of the mapping between the fcst reservoir and the network.
class RuleFcstRlsSched(RuleFcstBase):
    def __init__(self, name, timeFcst, fcstNetwork):
        # Call super class constructor
        super().__init__(name, timeFcst)
        self.fcstNetworkfcstNetwork = fcstNetwork
 
    def calc_release_fcst(self, rlsProposed, rlsUnCtrl, stor):
        pass