プログラムのソースコード Source Code of Programs


多状態スキーレンタル問題 Multislope Ski-Rental Problem


# The best possible strategy for the multislope ski-rental problem,
# coded by Ryota Uchida and Hiroshi Fujiwara for Python 3.x.
#
# The algorithm is due to Section 5 of the FOCS paper
# [Augustine&Irani&Swamy2004]. The runtime is O(k^2 log k log
# 1/epsilon) for a (k+1)-slope ski-rental instance.  The notation is
# due to the STACS paper [Lotker&Patt-Shamir&Rawitz2008].
#
# Example:
# 
# best_possible_skirental_strategy([1.0, 0.6, 0.0],
#                                  [[0.0, 0.4, 1.0],
#                                   [0.0, 0.0, 1.0]])
# 
# best_possible_skirental_strategy([1.0, 0.6, 0.3, 0.0],
#                                  [[0.0, 0.2, 0.5, 1.0],
#                                   [0.0, 0.0, 0.3, 0.8],
#                                   [0.0, 0.0, 0.0, 0.5]])
#
# (The latter is the example in Section 2.4 in the JOCO paper
# [Fujiwara&Kitano&Fujito2016].)
 
def a(r, b, j):
    """
    For multislope ski-rental instance (r, b), the function returns
    the root of equation OFF_{j-1}(t) = OFF_{j}(t), where OFF_{j}(t) =
    r_{j}*t + b_{0,j}. Namely, the returned value is the time at which
    OFF_j is optimal for the first time.
    In [Augustine&Irani&Swamy2004], a(r, b, j) is denoted by b_j.
    """
    return (b[0][j] - b[0][j-1]) / (r[j-1] - r[j])
 
def OPT(r, b, t):
    """
    For multislope ski-rental instance (r, b),
    the function returns min_i (r[i]*t + b[0][i]).
    """
 
    k = len(r) - 1
 
    for i in range(k):
        # if a(r, b, i) <= t < a(r, b, i+1)
        if t < a(r, b, i+1):
            return r[i]*t + b[0][i]
 
    # if a(r, b, k) <= t
    return r[k]*t + b[0][k]   
 
def opt_state_search(r, b, p, i, j, t_i, left_most_state, right_most_state):
    """
    For multislope ski-rental instance (r, b), suppose that a strategy makes a
    p-eager transition to state i at time t_i, and next makes a
    p-eager transition to state j.
 
    (A) If there is l between left_most_state and right_most_state
    such that the total cost so far = p OFF_{l}(t) holds, the function
    returns this l, where OFF_{j}(t) = r_{j}*t + b_{0,j}.
 
    (B) Otherwise (i.e., if a p-eager transition to state j is
    impossible), the function returns None.
 
    Assumption:
    For each i = 0, 1, ..., j-1, E[i][0] and E[i][1] are already
    calculated, where
    E[i][0] = the earliest time of eager p-transition to state i, and
    E[i][1] = the list of states that the strategy passed through.
    For example,
    E = [(0, [0]), (0.2, [0, 1]), (0.8, [0, 2]), (None, None), (1.3, [0, 1, 3])].
 
    Note:
    The algorithm is based on bisection search with recursive function
    call, due to Lemma 5 in [Augustine&Irani&Swamy2004].  The
    bisection search takes O(log k) time.
    """
 
    middle_state = (left_most_state + right_most_state) // 2
 
    # breadth of bisection search = 1
    if left_most_state == middle_state:
        l = p*OPT(r,b,t_i) + r[i]*(a(r,b,middle_state+1) - t_i) + b[i][j]
 
        # earliest p-eager transition is possible
        if l < p*OPT(r,b,a(r,b,middle_state+1)):
            return middle_state
        # earliest p-eager transition is impossible
        else: 
            return None 
 
    # breadth of bisection search > 1
    l = p*OPT(r,b,t_i) + r[i]*(a(r,b,middle_state) - t_i) + b[i][j]
    if l < p*OPT(r,b,a(r,b,middle_state)) or \
            (l >= p*OPT(r,b,a(r,b,middle_state)) and p*r[middle_state] <= r[i]):
        return opt_state_search(r, b, p, i, j, t_i, left_most_state, middle_state)
    else:
        return opt_state_search(r, b, p, i, j, t_i, middle_state, right_most_state)
 
 
def earliest_transition_time(r, b, p, E, j):
    """
    For multislope ski-rental instance (r, b):
 
    (A) If an eager p-transition transition time to state j is
    possible, the function return a tuple:
    (the earliest time, the used states (including state j)),
    such as (1.3, [0, 1, 3]).
 
    (B) Otherwise, the function returns a tuple: (None, None).
 
    Assumption:
    For each i = 0, 1, ..., j-1, E[i][0] and E[i][1] are already
    calculated, where
    E[i][0] = the earliest time of eager p-transition to state i, and
    E[i][1] = the list of states that the strategy passed through.
    For example,
    E = [(0, [0]), (0.2, [0, 1]), (0.8, [0, 2]), (None, None), (1.3, [0, 1, 3])].
 
    Note:
    The function stands for a recursive relation for dynamic
    programming; one can calculate all E's recursively by assigning
    E[j] = earliest_transition_time(r, b, p, E, j).  The function
    calls function opt_state_search O(k^2) times.  Therefore, the
    runtime is O(k^2 log k). The procedure is referred to as procedure
    EXISTS in the mention after Lemma 6 in
    [Augustine&Irani&Swamy2004], where p is denoted by rho.
    """
 
    k = len(r) - 1
 
    earliest = float("inf")
    origin_state = None
 
    # seek best origin state finally to state j
    # candidates are states 0, 1, ..., j-1
    for i in range (0,j):
        if E[i][0] is not None:
            opt_state = opt_state_search(r, b, p, i, j, E[i][0], 0, k)
 
            # p-transition from i to j is possible
            if opt_state is not None:
                # calculate transition time
                t = (p*OPT(r,b,E[i][0])-r[i]*E[i][0]+b[i][j]-p*b[0][opt_state]) / \
                    (p*r[opt_state]-r[i])
 
                # t is earliest found so far
                if t < earliest:
                    earliest = t
                    origin_state = i
 
    # p-transition is possible for some i
    if earliest < float("inf"):
        # construct a new list of used states
        used_states = list(E[origin_state][1])
        used_states.append(j)
 
        # for example, (1.3, [0, 1, 3])
        return (earliest, used_states)
 
    # p-transition is impossible for any i
    else:
        return (None, None)
 
def best_possible_skirental_strategy(
    r, b, epsilon=1e-14, upper_bound_on_CR=4.0, lower_bound_on_CR=1.0):
    """
    The function returns, for multislope ski-rental instance (r, b),
    a tuple:
    (the competitive ratio, the strategy, the states that the strategy uses).
    For the instance, the returned strategy is optimal with precision epsilon.
    The search for the competitive ratio is done between 
    lower_bound_on_CR and upper_bound_on_CR.
 
    Assumption:
    A multislope ski-rental instance is given as follows:
    For i = 0, 1, ..., k,
    r[i] = per-time cost for staying state i.
    For 0 <= i < j <= k,
    b[i][j] = transition cost from state i to state j.
    r is a (k+1)-element list, and
    b is a k-element list of (k+1)-element lists.
    In [Augustine&Irani&Swamy2004],
    r[i] and b[i][j] are respectively denoted by kappa_i and d_i,j.
 
    An example of a 3-slope ski-rental instance:
    r = [1.0, 0.6, 0.0], and
    b = [[0.0, 0.4, 1.0], [0.0, 0.0, 1.0]].
 
    Note:
    The algorithm is based on bisection search mentioned after Lemma 6
    in [Augustine&Irani&Swamy2004].  The function calls function
    earliest_transition_time O(log 1/epsilon) times.  Therefore, the
    runtime is O(k^2 log k log 1/epsilon).
    """
 
    k = len(r) - 1
 
    previous_target_CR = lower_bound_on_CR
    while True: 
        target_CR = (upper_bound_on_CR + lower_bound_on_CR) / 2.0
 
        # E is a table for dynamic programming:
        E = [(0.0, [0])]
 
        # calculate E[1], E[2], ..., E[k] by dynamic programming
        for i in range (1, k+1):
            # obtain E[i] from E[1], E[2], ..., E[i-1]
            E.append(earliest_transition_time(r, b, target_CR, E, i))
 
        # (A) there is a feasible strategy
        if E[k][0] is not None:
            # (A-1) within desired precision: bisection search completed
            if abs(target_CR - previous_target_CR) < epsilon:
                # x: best possible strategy
                x = []
 
                # trick just for setting x_0 to be zero
                current_state = -1
 
                # for over state used
                for state in E[k][1]:
                    # if the strategy skips a state,
                    # multiple entries have the same value.
                    while current_state < state:
                        x.append(E[state][0])
                        current_state = current_state + 1
 
                return (target_CR, x, E[k][1])
 
            # (A-2) out of desired precision
            else:
                previous_target_CR = target_CR
                upper_bound_on_CR = target_CR
 
        # (B) there is no feasible strategy
        else:
            previous_target_CR = target_CR
            lower_bound_on_CR = target_CR