def disagg_prec_cascade(precip_daily, 
                        cascade_options, hourly=True,level=9,
                        shift=0,
                        test=False):
    """Precipitation disaggregation with cascade model (Olsson, 1998)

    Parameters
    ----------
    precip_daily : pd.Series
        daily data
    hourly: Boolean (for an hourly resolution disaggregation)
        if False, then returns 5-min disaggregated precipitation 
        (disaggregation level depending on the "level" variable)
    cascade_options : cascade object
        including statistical parameters for the cascade model
    shift : int
        shifts the precipitation data by shift steps (eg +7 for 7:00 to
        6:00)
    test : bool
        test mode, returns time series of each cascade level
    """

    if len(precip_daily) < 2:
        raise ValueError('Input data must have at least two elements.')

    # set missing values to zero:
    precip_daily = precip_daily.copy()
    missing_days = precip_daily.index[precip_daily.isnull()]
    precip_daily[missing_days] = 0

    if hourly:
        si = 5  # index of first level
    else:
        si = level
    # statistics for branching into two bins
    wxxcum = np.zeros((7, 2, 4))

    if isinstance(cascade_options, melodist.cascade.CascadeStatistics):
        # this is the standard case considering one data set for all levels
        # get cumulative probabilities for branching
        overwrite_stats = False
        for k in range(0, 7):
            wxxcum[k, :, :] = cascade_options.wxx[k, :, :]
            if k > 0:
                wxxcum[k, :, :] = wxxcum[k-1, :, :] + wxxcum[k, :, :]
    elif isinstance(cascade_options, list):
        if len(cascade_options) == si:#5
            overwrite_stats = True
            list_casc = cascade_options
        else:
            raise ValueError('Cascade statistics list must have %s elements!' % si)
    else:
        raise TypeError('cascade_options has invalid type')

    # arrays for each level
    n = len(precip_daily)
    vdn1 = np.zeros(n*2)
    vdn2 = np.zeros(n*4)
    vdn3 = np.zeros(n*8)
    vdn4 = np.zeros(n*16)
    vdn5 = np.zeros(n*32)
    if not hourly:
        vdn6 = np.zeros(n*64)
        vdn7 = np.zeros(n*128)
        vdn8 = np.zeros(n*256)
        vdn9 = np.zeros(n*512)
    if level == 10 or level == 11:
        vdn10 = np.zeros(n*1024)
        if level == 11:
            vdn11 = np.zeros(n*2048)
            
    # class boundaries for histograms
    wclassbounds = np.array([0.0,
                             0.1429,
                             0.2857,
                             0.4286,
                             0.5714,
                             0.7143,
                             0.8571,
                             1.0])

    # disaggregation for each level
    for l in range(1, si+1):
        if l == 1:
            vdn_in = precip_daily
            vdn_out = vdn1
        elif l == 2:
            vdn_in = vdn_out
            vdn_out = vdn2
        elif l == 3:
            vdn_in = vdn_out
            vdn_out = vdn3
        elif l == 4:
            vdn_in = vdn_out
            vdn_out = vdn4
        elif l == 5:
            vdn_in = vdn_out
            vdn_out = vdn5
       
        elif l == 6:
            vdn_in = vdn_out
            vdn_out = vdn6
        elif l == 7:
            vdn_in = vdn_out
            vdn_out = vdn7
        elif l == 8:
            vdn_in = vdn_out
            vdn_out = vdn8
        elif l == 9:
            vdn_in = vdn_out
            vdn_out = vdn9
        elif l == 10:
            vdn_in = vdn_out
            vdn_out = vdn10
        elif l == 11:
            vdn_in = vdn_out
            vdn_out = vdn11
                
        si -= 1
        if overwrite_stats:
            cascade_options = list_casc[si]
            for k in range(0, 7):
                wxxcum[k, :, :] = cascade_options.wxx[k, :, :]
                if k > 0:
                    wxxcum[k, :, :] = wxxcum[k-1, :, :] + wxxcum[k, :, :]
            meanvol = cascade_options.threshold[0]
        else:
            meanvol = cascade_options.threshold[si]

        # evaluate mean rainfall intensity for wet boxes
        # these values should be determined during the aggregation phase!!!!!
        # mean volume threshold
        # meanvol = np.mean(vdn_in[vdn_in>0.])

        # use values derived parameter by parameter estimation instead
        # see above

        j = 0
        for i in range(0, len(vdn_in)):
            # it's raining now?
            if vdn_in[i] > 0:
                # determine type of box
                if i == 0:  # only starting or isolated
                    if vdn_in[i+1] > 0:
                        vbtype = cascade.BoxTypes.starting
                    else:
                        vbtype = cascade.BoxTypes.isolated
                elif i == len(vdn_in)-1:  # only ending or isolated
                    if vdn_in[i-1] > 0:
                        vbtype = cascade.BoxTypes.ending
                    else:
                        vbtype = cascade.BoxTypes.isolated
                else:  # neither at at the end nor at the beginning
                    if vdn_in[i-1] == 0 and vdn_in[i+1] == 0:
                        vbtype = cascade.BoxTypes.isolated

                    if vdn_in[i-1] == 0 and vdn_in[i+1] > 0:
                        vbtype = cascade.BoxTypes.starting

                    if vdn_in[i-1] > 0 and vdn_in[i+1] > 0:
                        vbtype = cascade.BoxTypes.enclosed

                    if vdn_in[i-1] > 0 and vdn_in[i+1] == 0:
                        vbtype = cascade.BoxTypes.ending

                # above or below mean?
                if vdn_in[i] > meanvol:
                    belowabove = 1  # above mean
                else:
                    belowabove = 0  # below mean

                #
                p = np.zeros((3, 1))
                p[0] = cascade_options.p01[belowabove, vbtype-1]  # index changed!
                p[1] = cascade_options.p10[belowabove, vbtype-1]
                p[2] = cascade_options.pxx[belowabove, vbtype-1]

                # draw a random number to determine the braching type
                rndp = np.random.random()

                if rndp <= p[0]:
                    # first box 0, second box: 1 P(0/1)
                    vdn_out[j] = 0.0
                    j = j + 1
                    vdn_out[j] = vdn_in[i]
                    j = j + 1
                elif rndp > p[0] and rndp <= p[0] + p[1]:
                    # first box 1, second box: 0 P(1/0)
                    vdn_out[j] = vdn_in[i]
                    j = j + 1
                    vdn_out[j] = 0.0
                    j = j + 1
                else:
                    # both boxes wet
                    # we need a new random number
                    rndw = np.random.random()

                    # guess w1:
                    for k in range(0, 7):
                        if rndw <= wxxcum[k, belowabove, vbtype-1]:
                            w1 = wclassbounds[k+1] - 1./14.  # class center
                            break

                    vdn_out[j] = w1 * vdn_in[i]
                    j = j + 1
                    vdn_out[j] = (1. - w1) * vdn_in[i]
                    j = j + 1

                    # check results (in the previous version this error has never been observed)
                    if w1 < 0 or w1 > 1:
                        print('error')
                        return
            else:
                # add two dry boxes
                vdn_out[j] = 0.0
                j = j + 1
                vdn_out[j] = 0.0
                j = j + 1
    
    if hourly:
        # uniformly disaggregate 0.75 h values to 0.25 h values
        vdn_025 = np.zeros(len(vdn_out)*3)
        j = 0
        for i in range(0, len(vdn_out)):
            for m in range(0, 3):
                vdn_025[j+m] = vdn_out[i] / 3.
            j = j + 3
    
        # aggregate to hourly time steps
        vdn_025cs = np.cumsum(vdn_025)
        vdn = np.zeros(int(len(vdn_025)/4))
        for i in range(0, len(vdn)+1):
            # for first hour take 4th item
            if i == 0:
                vdn[i] = vdn_025cs[3]
            elif i == 1:
                pass
            else:
                # >1 (starting with 2-1 = 1 item)
                vdn[i-1] = vdn_025cs[(i*4)-1] - vdn_025cs[(i*4)-5]
    
        disagg_precip = pd.Series(index=melodist.util.hourly_index(precip_daily.index), data=vdn)
    
    else:
        precip_sn = pd.Series(index= sub_level_index(precip_daily.index, level=level, fill_gaps=False), data=vdn_out)
        disagg_precip = precip_sn.resample('5min').sum()
    
    # set missing days to nan again:
    for date in missing_days:
         disagg_precip[ disagg_precip.index.date == date.date()] = np.nan

    # shifts the data by shift steps (fills with nan/cuts edge data )
    if shift != 0:
         disagg_precip =  disagg_precip.shift(shift) #? freq='1U')

    # return time series
    if test:
        if hourly:
            return vdn1, vdn2, vdn3, vdn4, vdn5, vdn_025, disagg_precip
        else:
            if level == 9:
                return vdn1, vdn2, vdn3, vdn4, vdn5, vdn6, vdn7, vdn8, vdn9, precip_sn, disagg_precip
            elif level == 10:
                return vdn1, vdn2, vdn3, vdn4, vdn5, vdn6, vdn7, vdn8, vdn9, vdn10, precip_sn, disagg_precip
            else:
                return vdn1, vdn2, vdn3, vdn4, vdn5, vdn6, vdn7, vdn8, vdn9, vdn10, vdn11, precip_sn, disagg_precip
    else:
        return  disagg_precip
def disagg_prec_cascade ( precip_daily , cascade_options , hourly = True , level = 9 , shift = 0 , test = False ) : if len ( precip_daily ) < 2 : raise ValueError ( 'Input data must have at least two elements.' ) precip_daily = precip_daily . copy ( ) missing_days = precip_daily . index [ precip_daily . isnull ( ) ] precip_daily [ missing_days ] = 0 if hourly : si = 5 else : si = level wxxcum = np . zeros ( ( 7 , 2 , 4 ) ) if isinstance ( cascade_options , melodist . cascade . CascadeStatistics ) : overwrite_stats = False for k in range ( 0 , 7 ) : wxxcum [ k , : , : ] = cascade_options . wxx [ k , : , : ] if k > 0 : wxxcum [ k , : , : ] = wxxcum [ k - 1 , : , : ] + wxxcum [ k , : , : ] elif isinstance ( cascade_options , list ) : if len ( cascade_options ) == si : overwrite_stats = True list_casc = cascade_options else : raise ValueError ( 'Cascade statistics list must have %s elements!' % si ) else : raise TypeError ( 'cascade_options has invalid type' ) n = len ( precip_daily ) vdn1 = np . zeros ( n * 2 ) vdn2 = np . zeros ( n * 4 ) vdn3 = np . zeros ( n * 8 ) vdn4 = np . zeros ( n * 16 ) vdn5 = np . zeros ( n * 32 ) if not hourly : vdn6 = np . zeros ( n * 64 ) vdn7 = np . zeros ( n * 128 ) vdn8 = np . zeros ( n * 256 ) vdn9 = np . zeros ( n * 512 ) if level == 10 or level == 11 : vdn10 = np . zeros ( n * 1024 ) if level == 11 : vdn11 = np . zeros ( n * 2048 ) wclassbounds = np . array ( [ 0.0 , 0.1429 , 0.2857 , 0.4286 , 0.5714 , 0.7143 , 0.8571 , 1.0 ] ) for l in range ( 1 , si + 1 ) : if l == 1 : vdn_in = precip_daily vdn_out = vdn1 elif l == 2 : vdn_in = vdn_out vdn_out = vdn2 elif l == 3 : vdn_in = vdn_out vdn_out = vdn3 elif l == 4 : vdn_in = vdn_out vdn_out = vdn4 elif l == 5 : vdn_in = vdn_out vdn_out = vdn5 elif l == 6 : vdn_in = vdn_out vdn_out = vdn6 elif l == 7 : vdn_in = vdn_out vdn_out = vdn7 elif l == 8 : vdn_in = vdn_out vdn_out = vdn8 elif l == 9 : vdn_in = vdn_out vdn_out = vdn9 elif l == 10 : vdn_in = vdn_out vdn_out = vdn10 elif l == 11 : vdn_in = vdn_out vdn_out = vdn11 si -= 1 if overwrite_stats : cascade_options = list_casc [ si ] for k in range ( 0 , 7 ) : wxxcum [ k , : , : ] = cascade_options . wxx [ k , : , : ] if k > 0 : wxxcum [ k , : , : ] = wxxcum [ k - 1 , : , : ] + wxxcum [ k , : , : ] meanvol = cascade_options . threshold [ 0 ] else : meanvol = cascade_options . threshold [ si ] j = 0 for i in range ( 0 , len ( vdn_in ) ) : if vdn_in [ i ] > 0 : if i == 0 : if vdn_in [ i + 1 ] > 0 : vbtype = cascade . BoxTypes . starting else : vbtype = cascade . BoxTypes . isolated elif i == len ( vdn_in ) - 1 : if vdn_in [ i - 1 ] > 0 : vbtype = cascade . BoxTypes . ending else : vbtype = cascade . BoxTypes . isolated else : if vdn_in [ i - 1 ] == 0 and vdn_in [ i + 1 ] == 0 : vbtype = cascade . BoxTypes . isolated if vdn_in [ i - 1 ] == 0 and vdn_in [ i + 1 ] > 0 : vbtype = cascade . BoxTypes . starting if vdn_in [ i - 1 ] > 0 and vdn_in [ i + 1 ] > 0 : vbtype = cascade . BoxTypes . enclosed if vdn_in [ i - 1 ] > 0 and vdn_in [ i + 1 ] == 0 : vbtype = cascade . BoxTypes . ending if vdn_in [ i ] > meanvol : belowabove = 1 else : belowabove = 0 p = np . zeros ( ( 3 , 1 ) ) p [ 0 ] = cascade_options . p01 [ belowabove , vbtype - 1 ] p [ 1 ] = cascade_options . p10 [ belowabove , vbtype - 1 ] p [ 2 ] = cascade_options . pxx [ belowabove , vbtype - 1 ] rndp = np . random . random ( ) if rndp <= p [ 0 ] : vdn_out [ j ] = 0.0 j = j + 1 vdn_out [ j ] = vdn_in [ i ] j = j + 1 elif rndp > p [ 0 ] and rndp <= p [ 0 ] + p [ 1 ] : vdn_out [ j ] = vdn_in [ i ] j = j + 1 vdn_out [ j ] = 0.0 j = j + 1 else : rndw = np . random . random ( ) for k in range ( 0 , 7 ) : if rndw <= wxxcum [ k , belowabove , vbtype - 1 ] : w1 = wclassbounds [ k + 1 ] - 1. / 14. break vdn_out [ j ] = w1 * vdn_in [ i ] j = j + 1 vdn_out [ j ] = ( 1. - w1 ) * vdn_in [ i ] j = j + 1 if w1 < 0 or w1 > 1 : print ( 'error' ) return else : vdn_out [ j ] = 0.0 j = j + 1 vdn_out [ j ] = 0.0 j = j + 1 if hourly : vdn_025 = np . zeros ( len ( vdn_out ) * 3 ) j = 0 for i in range ( 0 , len ( vdn_out ) ) : for m in range ( 0 , 3 ) : vdn_025 [ j + m ] = vdn_out [ i ] / 3. j = j + 3 vdn_025cs = np . cumsum ( vdn_025 ) vdn = np . zeros ( int ( len ( vdn_025 ) / 4 ) ) for i in range ( 0 , len ( vdn ) + 1 ) : if i == 0 : vdn [ i ] = vdn_025cs [ 3 ] elif i == 1 : pass else : vdn [ i - 1 ] = vdn_025cs [ ( i * 4 ) - 1 ] - vdn_025cs [ ( i * 4 ) - 5 ] disagg_precip = pd . Series ( index = melodist . util . hourly_index ( precip_daily . index ) , data = vdn ) else : precip_sn = pd . Series ( index = sub_level_index ( precip_daily . index , level = level , fill_gaps = False ) , data = vdn_out ) disagg_precip = precip_sn . resample ( '5min' ) . sum ( ) for date in missing_days : disagg_precip [ disagg_precip . index . date == date . date ( ) ] = np . nan if shift != 0 : disagg_precip = disagg_precip . shift ( shift ) if test : if hourly : return vdn1 , vdn2 , vdn3 , vdn4 , vdn5 , vdn_025 , disagg_precip else : if level == 9 : return vdn1 , vdn2 , vdn3 , vdn4 , vdn5 , vdn6 , vdn7 , vdn8 , vdn9 , precip_sn , disagg_precip elif level == 10 : return vdn1 , vdn2 , vdn3 , vdn4 , vdn5 , vdn6 , vdn7 , vdn8 , vdn9 , vdn10 , precip_sn , disagg_precip else : return vdn1 , vdn2 , vdn3 , vdn4 , vdn5 , vdn6 , vdn7 , vdn8 , vdn9 , vdn10 , vdn11 , precip_sn , disagg_precip else : return disagg_precip
Precipitation disaggregation with cascade model (Olsson, 1998)

    Parameters
    ----------
    precip_daily : pd.Series
        daily data
    hourly: Boolean (for an hourly resolution disaggregation)
        if False, then returns 5-min disaggregated precipitation 
        (disaggregation level depending on the "level" variable)
    cascade_options : cascade object
        including statistical parameters for the cascade model
    shift : int
        shifts the precipitation data by shift steps (eg +7 for 7:00 to
        6:00)
    test : bool
        test mode, returns time series of each cascade level