def getItemReadVariables(events, item_reads, surface):
    """Get ItemRead values.
    
    Return all variables used in ItemRead optimizations.

    Args:
        events: List of Event dicts.
        item_reads: List of ItemRead dicts.
        surface: Surface for which to return variables.
    
    Returns:
        Tuple of variables consisting of the following:
            surface_events: Same as events but only for this surface.
            surface_item_reads: Same as item_reads but only for this surface.
            W: R^N weight vector where N is the number of surface ItemReads.
            A: R^{TxN} arrival matrix where A[t][i] is whether ItemRead i 
                arrived at time t.
            R: R^{TxN} removal matrix where R[t][i] is whether ItemRead i
                was removed at time t.
            S: R^T scale weight vector.
    """
    # Get all events for this surface and for events 
    # not associated with any surface (no weight data).
    surface_events = [v for v in events if v["Event.surface"]==surface or v["Event.surface"]==None]
    # Get all item reads that arrived at and/or were removed from this surface.
    surface_item_reads = [v for v in item_reads if v["ArrivalEvent.surface"]==surface or \
                          ("RemovalEvent.surface" in v and v["RemovalEvent.surface"]==surface)]
    T, N = len(surface_events), len(surface_item_reads)
    times, S = zip(*sorted([(str(v["Event.time"]), v["Event.event_force"]) for v in surface_events]))
    # Arrival and removal matrices
    A, R = [[[0 for i in range(N)] for t in range(T)] for v in range(2)]
    for (i, item) in enumerate(surface_item_reads):
        t_arrival = times.index(item["ArrivalEvent.time"])
        A[t_arrival][i] = 1
        if "RemovalEvent.time" in item:
            t_removal = times.index(item["RemovalEvent.time"])
            R[t_removal][i] = 1
    # Get optimization variables.
    W = [v["ItemRead.weight"] for v in surface_item_reads]
    return surface_events, surface_item_reads, W, A, R, S