<jupyter_start><jupyter_text>InfoThis notebook includes a solver modified to simulate forward different resource arrival scenarios. The solver contains a few small differences, subtracting the number of allocated copies still pending at a location.<jupyter_code>import numpy as np
import pickle
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
from gurobi import *<jupyter_output><empty_output><jupyter_text>General utility functions<jupyter_code>def print_tableau(matches):
    for i in range(matches.shape[2]):
        print("t=",i)
        print(matches[:,:,i])<jupyter_output><empty_output><jupyter_text>Primal-Dual Solver Implementation  Solver utility functions<jupyter_code>def in_constraint(v):
    if v[1]:
        return True
    else:
        return False

def make_alpha_mapping(I,J,T,alphas,valid_matches):
    '''Creates an index into the alpha array for each position in valid matches'''
    constraints_d, _ = dual_constraint_matrix(valid_matches,pairing_weights,I,J,T,k)
    constraints_d = constraints_d[:,:alphas.size]
    alpha_map = np.zeros((*valid_matches.shape,constraints_d.shape[1]),dtype=np.bool)

    cix=0
    for i in range(I):
            for j in range(J):
                for t in range(T):
                    if valid_matches[i][j][t]:
                        alpha_map[i,j,t,:] = constraints_d[cix,:]
                        cix += 1
    return alpha_map

def sum_alpha_it(alphas,i,t,k):
    '''Return the alpha terms summed at a given i,t'''
    alpha_it = alphas.reshape(I,T)
    startit = max(t-k[i]+1,0)
    return np.sum(alpha_it[i,startit:t+1])<jupyter_output><empty_output><jupyter_text>Convert the problem formulation to matricies of the form `Ax = b`To make the dual constraint matrix: - Create a constraint map to see which alphas/betas apply at a given location in the primal. Each valid location will correspond to a constraint in the dual, and the variables == 1 will be those in the `cmap[i][j][t]`<jupyter_code>def primal_constraint_matrix(valid_matches,I,J,T,k):
    '''
    Primal constraint matrix: Number constraints * number of IJT positions
             - Each row corresponds to an IJT position in the grid
             - Each column corresponds to a dual variable 
             
    This matrix is formed by constructing all primal constraints, then boolean indexing each IJT position
    '''

    constraints = np.zeros((T*I+J,valid_matches.size),dtype=np.float128)
    cix = 0
    
    #constraints limiting to one resource allocation in the time interval
    for i in range(I):
        for t in range(T):
            constraint = np.zeros((I,J,T), np.int)
            valid_mask = constraint.copy()
            endix = min(t+k[i],T)
            valid_mask[i,:,t:endix] = 1 
            constraint[(valid_mask == 1)] = 1
            constraints[cix,:] = constraint.reshape((1, constraint.shape[0] * constraint.shape[1] * constraint.shape[2]))
            cix += 1

    #constraints limiting each agent to only match once            
    for j in range(J):
        constraint = np.zeros((I,J,T), np.int)
        valid_mask = constraint.copy()
        valid_mask[1:,j,:] = 1

        constraint[(valid_matches == 1) & (valid_mask ==1)] = 1
        constraints[cix+j,:] = constraint.reshape((1, constraint.shape[0] * constraint.shape[1] * constraint.shape[2]))
    
    return constraints


def dual_constraint_matrix(valid_matches,pairing_weights,I,J,T,k):
    '''
    Dual constraint matrix: Number IJT positions * number dual variables
             - Each row corresponds to an IJT position in the grid
             - Each column corresponds to a dual variable 
             
    This matrix is formed by constructing all primal constraints, then boolean indexing each IJT position
    '''
    num_positions = I*J*T
    num_primal_constraints = I*T+J
    dual_constraint_matrix = np.zeros((num_positions, num_primal_constraints))
    
    inequalities = np.zeros(num_positions)
    constraint_map = np.zeros((I,J,T,num_primal_constraints), np.int) 
    
    cix = 0

    #constraints limiting to one resource allocation in the time interval
    for i in range(I):
        for t in range(T):
            constraint = np.zeros((I,J,T), np.int)
            valid_mask = constraint.copy()

            endix = min(t+k[i],T)
            valid_mask[i,:,t:endix] = 1 
            constraint[(valid_mask == 1)] = 1
            constraint_map[:,:,:,cix] = constraint.copy()
            cix += 1

    #constraints limiting each agent to only match once            
    for j in range(J):
        constraint = np.zeros((I,J,T), np.int)
        valid_mask = constraint.copy()
        valid_mask[1:,j,:] = 1
        constraint[valid_mask ==1] = 1
        constraint_map[:,:,:,cix] = constraint.copy()
        cix += 1
    
    cix = 0
    for i in range(I):
        for j in range(J):
            for t in range(T):
                dual_constraint_matrix[cix,:] = constraint_map[i,j,t,:] 
                inequalities[cix] = pairing_weights[i,j,t]
                cix += 1
    
    return dual_constraint_matrix, inequalities<jupyter_output><empty_output><jupyter_text>Solving the primal and dual via the Gurobi library<jupyter_code>def primal_solutions(valid_matches, pairing_weights, I, J, T, k, c, pending_allocs=None):
    '''
        Sets up and solves the primal linear program using the gurobi library API. 
        Defers most business logic for for primal_constraint_matrix
        
        pending_allocs: if specified, accounts for previously made allocations in primal inequality, 
            which is necessary in the exaustive search simulator
    '''
   
    m = Model("dynamicmatch_primal")
    m.modelSense = GRB.MAXIMIZE
    m.setParam( 'OutputFlag', False )
    m.setParam( 'NumericFocus', 3)
    
    weights = pairing_weights.reshape(pairing_weights.shape[0] * pairing_weights.shape[1] * pairing_weights.shape[2])
    constraints = primal_constraint_matrix(valid_matches,I,J,T,k)
    
    keys = range(constraints.shape[1])
    variables = m.addVars(keys,
                    vtype=GRB.CONTINUOUS,
                     obj=weights,
                     name="primal",
                     lb=0)

    for cix, constraint in enumerate(constraints):
        ix = cix // T
        tx = cix % T
        
        used_copies = pending_allocs[ix,tx] if pending_allocs is not None and cix < T * I else 0
        equality = c[cix // T] - used_copies if cix < T * I else 1
        m.addConstr(sum(variables[o]*c for o,c in filter(in_constraint, zip(variables,constraint))) <= equality)

    m.optimize()
    allocations = np.array([variables[var].X for var in variables], dtype=np.float128).reshape(pairing_weights.shape)

    return m.objVal, allocations


def dual_solutions(valid_matches, pairing_weights, I, J, T, k, c, pending_allocs=False):
    '''
        Sets up and solves the dual linear program using the gurobi library API. 
        Defers most business logic for for dual_constraint_matrix

        pending_allocs: if specified, accounts for previously made allocations in dual coefficients, 
                which is necessary in the exaustive search simulator
    '''

    md = Model("dynamicmatch_dual")
    md.modelSense = GRB.MINIMIZE
    md.setParam( 'OutputFlag', False )
    md.setParam( 'NumericFocus', 3)

    constraints_d, inequalities = dual_constraint_matrix(valid_matches,pairing_weights,I,J,T,k)
    variables_d = np.ones(constraints_d.shape[1])
    
    for cix in range(constraints_d.shape[1]):
        ix = cix // T
        tx = cix % T
        used_copies = pending_allocs[ix,tx] if pending_allocs is not None and cix < T * I else 0
        
        variables_d[ix] = c[ix // T] - used_copies if ix < T * I else 1
    
    keys = range(constraints_d.shape[1])
    variables = md.addVars(keys,
                    vtype=GRB.CONTINUOUS,
                    obj=variables_d,
                    name="dual",
                    lb=0)

    for cix, constraint in enumerate(constraints_d):
        constr = inequalities[cix] + sum( variables[o]*-1*c for o,c in filter(in_constraint, zip(variables,constraint))) <= 0
        md.addConstr(constr)
        
    md.optimize()
    duals = np.array([variables[var].X for var in variables],dtype=np.longdouble)
    betas = duals[duals.size - J:]
    alphas = duals[:duals.size - J]
    
    return md.objVal, alphas, betas

def create_simulation_scenario(I,J,T, scenario=None, arrivals=''):
    "Defines tableau for 2 x 2 model according to arrivals and agent scenario distribution"
    if scenario:
        assert I == scenario['k'].size and I == scenario['c'].size and I == len(scenario['resource_configs'])+1, 'Scenario provided is not valid'
        assert arrivals != '', 'Arrival model not specified'
    
    valid_matches = np.zeros((I,J,T), np.int)
    weight_noise = np.random.random(valid_matches.shape) 

    #Only consider valid indicies
    for t in range(T):
        valid_matches[:,max(0,t-d):min(t+1,J),t] = 1

    if scenario == None:
        pairing_weights = np.random.random(valid_matches.shape) 

    #Construct matching matrix based on scenario configuration
    else:
        pairing_weights = np.zeros(valid_matches.shape)
        
        for resource_ix, config in enumerate(scenario['resource_configs']):
            
            if arrivals == "Fixed":
                agentAs = np.arange(J)[::2]
                agentBs = np.arange(J)[1::2]

            elif arrivals == "Stochastic":
                pA = scenario['resource_configs'][resource_ix]['probability-A']
                agentAs = np.arange(J)[np.random.choice([1, 0], size=J, p=[pA, 1-pA]) == 1]
                agentBs = np.setdiff1d(np.arange(J), agentAs)

            elif arrivals == "Adversarial":
                agentAs = np.arange(J)[0:int(J/2)]
                agentBs = np.setdiff1d(np.arange(J), agentAs)
                
            pairing_weights[resource_ix+1,agentAs,:] = scenario['resource_configs'][resource_ix]['agent-A']
            pairing_weights[resource_ix+1,agentBs,:] = scenario['resource_configs'][resource_ix]['agent-B']

    #Random noise to break ties
    for t in range(1,pairing_weights.shape[2]):
        pairing_weights[:,:,t] = pairing_weights[:,:,t] + .0001 * weight_noise[:,:,t]  
    
    #Modify pairing weights based on problem formulation 
    pairing_weights[0,:,:] =  0
    pairing_weights[valid_matches == 0] = -1
    pairing_weights[1:,:,:] = np.round(pairing_weights[1:,:,:],decimals=8)

    return valid_matches, pairing_weights<jupyter_output><empty_output><jupyter_text>init completed_allocsfor t in range(t)[0::sim_interval]:      for n in range(allocs_per_sim_interval):        pending_resources = get_pending_resource_counts(completed_allocs)                for sim in nsims:            generate random simulatin for rest of arrivals based on model            if primal:                sim_allocs[i] = primal_solutions_exhaustive()            elif dual:                 objd, alphas, betas = dual_solutions(valid_matches,pairing_weights, I, J, T, k, c)                objo, sim_allocs[i] = online_matching(I,J,T,k,c,alphas,betas,allocs,valid_matches,pairing_weights)        aix,ajx,atx = argmax(sim_allocs)        if atx < t + sim_interval            completed_allocs[aix, ajx, atx] = 1<jupyter_code>def exhaustive_simulate_forward(valid_matches_org, pairing_weights_org, I, J, T, scenario,
                                   sim_interval=1, nsims=10, allocs_per_interval=1, method='primal'):
     
    c , k = scenario['c'], scenario['k']
    completed_allocs = np.zeros(valid_matches_org.shape)
    
    for tx in range(T)[0::sim_interval]:  
        for n in range(allocs_per_interval):

            sim_allocs = np.zeros((nsims,I,J,T-tx))
            pending_resources = get_resource_uses_from_alloc(I,J,T,c,k,completed_allocs)
            
            for sim in range(nsims):    
                valid_matches_t, pairing_weights_t = simulate_future_arrivals(I,J,T,tx,valid_matches_org, 
                                                                    pairing_weights_org, completed_allocs, scenario=scenario)           
                
                if method == 'primal':
                    _, sim_allocs[sim,:,:,:] = primal_solutions(valid_matches_t, pairing_weights_t, I, J, 
                                                    T-tx, k, c, pending_allocs=pending_resources[:,tx:])
            
            aix, ajx, atx = np.unravel_index(np.argmax(sim_allocs.sum(axis=0), axis=None), (I,J,T-tx))
            
            if atx < sim_interval and aix > 0:
                completed_allocs[aix, ajx, atx+tx] = 1
    
    obj = 0
    for i in range(I):
        for j in range(J):
            for t in range(T):
                obj += pairing_weights_org[i,j,t] * completed_allocs[i,j,t]
    
    return obj, completed_allocs
            
def get_resource_uses_from_alloc(I,J,T,c,k,allocs):

    resources_used = np.zeros((I,T),dtype=np.float)

    for t in range(T):
        for i in range(I):
            resources_used[i,t:t+k[i]] += allocs[i,:,t].sum()    
    return resources_used


def simulate_future_arrivals(I,J,T,tx,valid_matches_base, pairing_weights_base, allocs, scenario=None):
    '''
    Defines tableau for 2 x 2 model according to stochastic arrivals and agent scenario distribution
    tx: the index of the current arrival
    For a given t: 
        - Only agents ariving after time t have weights simulated
        - Previously-matched agents have weights set to zero 
    '''
    
    arrived_agents = np.arange(J)[np.array(np.arange(J) <= tx, dtype=np.bool)]
    matched_agents = np.arange(J)[np.array(allocs.sum(axis=2).sum(axis=0) == 1, dtype=np.bool)]
    
    _, pairing_weights_sim = create_simulation_scenario(I,J,T, scenario=scenario, arrivals='Stochastic')
    
    #Agents who already arrived remain the same
    pairing_weights_sim[:,arrived_agents, :] = pairing_weights_base[:,arrived_agents, :]
    
    #Previously matched agents have zero utility
    pairing_weights_sim[:,matched_agents,:] = 0
    
    #Mantain problem formulation restrictions
    pairing_weights_sim[0,:,:] =  0
    pairing_weights_sim[valid_matches_base == 0] = -1

    return valid_matches_base[:,:,tx:], pairing_weights_sim[:,:,tx:]<jupyter_output><empty_output><jupyter_text>Creating a scenario and calling the solver<jupyter_code>J=30
I=3
d=2
T = J+d 

scenarios = {}
scenarios['1'] = {}
scenarios['1']['k'] = np.array([1,2,2])
scenarios['1']['c'] =  np.array([J,1, 1])
scenarios['1']['resource_configs'] = []
scenarios['1']['resource_configs'].append({'agent-A': .05,'agent-B':.05,'probability-A': 1})
scenarios['1']['resource_configs'].append({'agent-A': .1,'agent-B':.9,'probability-A': .5})
valid_matches, pairing_weights = create_simulation_scenario(I,J,T,scenario=scenarios['1'],arrivals='Stochastic')

objpe, allocs_exhaustive = exhaustive_simulate_forward(valid_matches,pairing_weights,I,J,T,scenarios['1'],
                                   sim_interval=30, nsims=10, allocs_per_interval=1, method='primal')

objp, allocs = primal_solutions(valid_matches, pairing_weights, 
                                    I, J, T, scenarios['1']['k'], scenarios['1']['c'])

print('primal', objp)
print('simulate forward', objpe)<jupyter_output>primal 11.10177095
simulate forward 0.05005791
StarcoderdataJupyterStructuredCleanDedup