def my_cumulative(s, d, r, b):
    """
    Decompositon of cumulative.
    
    Inspired by the MiniZinc implementation.
    
    The MiniZinc decomposition is discussed in the paper:
    A. Schutt, T. Feydy, P.J. Stuckey, and M. G. Wallace.
    'Why cumulative decomposition is not as bad as it sounds.'

    Parameters:

    s: start_times    assumption: array of varint
    d: durations      assumption: array of int
    r: resources      assumption: array of int
    b: resource limit assumption: varint or int
    """
    
    constraints = []
    max_d = max(d)
    tasks = [i for i in range(len(s)) if r[i] > 0 and d[i] > 0]
    times_min = min([s[i].lb for i in tasks])
    times_max = max([s[i].ub + max_d for i in tasks])
    for t in range(times_min, times_max + 1):
        constraints += [ b >= sum([((s[i] <= t) & (t < s[i] + d[i])) * r[i] for i in tasks])]
            
    # Somewhat experimental:
    # This constraint is needed to contrain the upper limit of b.
    if not isinstance(b, int):
        constraints += [b <= sum(r)]
        
    return constraints