From 2cd68b4cc2a6246ec85151e12dda408b2e1a84eb Mon Sep 17 00:00:00 2001 From: ishikaghosh2201 <112980412+ishikaghosh2201@users.noreply.github.com> Date: Thu, 4 Sep 2025 16:16:48 -0400 Subject: [PATCH 1/7] ILP optimization without the distance matrix --- .../cereeberus/distance/ilp_without_D.py | 404 ++++++++++++++++++ 1 file changed, 404 insertions(+) create mode 100644 cereeberus/cereeberus/distance/ilp_without_D.py diff --git a/cereeberus/cereeberus/distance/ilp_without_D.py b/cereeberus/cereeberus/distance/ilp_without_D.py new file mode 100644 index 0000000..3cd1e51 --- /dev/null +++ b/cereeberus/cereeberus/distance/ilp_without_D.py @@ -0,0 +1,404 @@ +# from cereeberus import ReebGraph, MapperGraph, Interleave +from .labeled_blocks import LabeledBlockMatrix as LBM +from .labeled_blocks import LabeledMatrix as LM +# import cereeberus.data.ex_mappergraphs as ex_mg + +import matplotlib.pyplot as plt +import numpy as np + +import pulp #for ILP optimization + + + +# function to build the phi and psi matrices after the ILP optimization +def build_map_matrices(myAssgn, map_results): + """ + Function to build the map matrices after the ILP optimization. The end result is a dictionary with keys + - 'phi_0_V', 'phi_0_E' + - 'phi_n_V', 'phi_n_E' + - 'psi_0_V', 'psi_0_E' + - 'psi_n_V', 'psi_n_E'. + Each key corresponds to a LabeledBlockMatrix containing the respective map matrices for the given thickening. + + Parameters: + myAssgn (Assignment): the Assignment object + map_results (dict): the dictionary containing the results of the ILP optimization from the solve_ilp function + + Returns: + dict: the dictionary containing the final map matrices + """ + + # get the function values + func_vals = myAssgn.all_func_vals() + + # create a dictionary to store the final matrices + final_LBMs = {} + + for thickening in ['0', 'n']: + for obj_type in ['V', 'E']: + for map_type in ['phi', 'psi']: + + # Gets the full LBM for the relevant interleaving map + M = myAssgn.get_interleaving_map(maptype=map_type, key=thickening, obj_type=obj_type) + + for i in M.get_all_block_indices(): + block_shape = M[i].array.shape + + # Reset to zero + M[i].array = np.zeros(block_shape) + + # Set the values according to the ILP optimization results + for a in range(block_shape[0]): + for b in range(block_shape[1]): + # Get the relevant value from the ILP optimization results + if a<0 or b<0: + print(f"Warning: a={a}, b={b} for block {i} in {map_type}_{thickening}_{obj_type}.") + try: + M[i].array[a,b] = pulp.value(map_results[map_type+'_vars'][i][thickening][obj_type][(a,b)]) + except KeyError: + # TODO: This makes me nervous, but I think the issue is just that there are empty rows/columns that we want to skip + # print(f"KeyError: {map_type+'_vars'} for block {i} in {map_type}_{thickening}_{obj_type}.") + continue + + + output_key = f"{map_type}_{thickening}_{obj_type}" + final_LBMs[output_key] = M + + return final_LBMs + +##------------------- ILP Optimization -------------------## + +def solve_ilp(myAssgn, pulp_solver = None, verbose=False): + + """ + Function to solve the ILP optimization problem for interleaving maps. The function creates a linear programming problem using the PuLP library and solves it to find the optimal interleaving maps. + + Parameters: + myAssgn (Assignment): the Assignment object containing the interleaving maps and other relevant data + pulp_solver (pulp.LpSolver): the solver to use for the ILP optimization. If None, the default solver is used. + verbose (bool): whether to print the optimization status and results + + Returns: + tuple: a tuple containing the final interleaving maps (as a dictionary of LabeledBlockMatrices) and the optimized loss value + """ + # function values + func_vals = myAssgn.all_func_vals() + + # create the ILP problem + prob = pulp.LpProblem("Interleave_Optimization_Problem", pulp.LpMinimize) + + # create empty dictionaries to store the decision variables + phi_vars = {block : {thickening:{obj_type:{} for obj_type in ['V', 'E']} for thickening in ['0', 'n']} for block in func_vals} + psi_vars = {block : {thickening:{obj_type:{} for obj_type in ['V', 'E']} for thickening in ['0', 'n']} for block in func_vals} + + + # create the decision variables (NOTE: these are all the decision variables for all the diagrams) + for thickening in ['0', 'n']: + for obj_type in ['V', 'E']: + for block in (func_vals[:-1] if obj_type == 'E' else func_vals): + + # create lp variables for phi + n_rows = myAssgn.phi(thickening, obj_type)[block].get_array().shape[0] + n_cols = myAssgn.phi(thickening, obj_type)[block].get_array().shape[1] + + phi_vars[block][thickening][obj_type] = pulp.LpVariable.dicts('phi_'+thickening+obj_type+'_'+str(block), ((a, b) for a in range(n_rows) for b in range(n_cols)), cat='Binary') + + # # set the initial values + # for a in range(n_rows): + # for b in range(n_cols): + # prob += phi_vars[block][thickening][obj_type][(a,b)] == myAssgn.phi(thickening, obj_type)[block].get_array()[a][b] + + + # create lp variables for psi + n_rows = myAssgn.psi(thickening, obj_type)[block].get_array().shape[0] + n_cols = myAssgn.psi(thickening, obj_type)[block].get_array().shape[1] + + psi_vars[block][thickening][obj_type] = pulp.LpVariable.dicts('psi_'+thickening+obj_type+'_'+str(block), ((a, b) for a in range(n_rows) for b in range(n_cols)), cat='Binary') + + # # set the initial values + # for a in range(n_rows): + # for b in range(n_cols): + # prob += psi_vars[block][thickening][obj_type][(a,b)] == myAssgn.psi(thickening, obj_type)[block].get_array()[a][b] + + # create the other decision variables + z_vars = {block : {starting_map: {obj_type: {} for obj_type in ['V', 'E']} for starting_map in ['F', 'G']} for block in func_vals} + + map_product_vars = {block : {starting_map: {obj_type: {} for obj_type in ['V', 'E']} for starting_map in ['F', 'G']} for block in func_vals} + for obj_type in ['V', 'E']: + for starting_map in ['F', 'G']: + for block in (func_vals[:-1] if obj_type == 'E' else func_vals): + + if starting_map == 'F': + n_rows_1 = myAssgn.psi('n', obj_type)[block].get_array().shape[0] + n_cols_1 = myAssgn.phi('0', obj_type)[block].get_array().shape[1] + + # set the map product variables + map_product_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts(starting_map+'_'+obj_type+'_'+str(block), ((a, b) for a in range(n_rows_1) for b in range(n_cols_1)), cat='Integer') + + n_rowcol_1 = myAssgn.psi('n', obj_type)[block].get_array().shape[1] + + # set the z variables + z_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts('z_'+starting_map+'_'+obj_type+'_'+str(block), ((a, b, c) for a in range(n_rows_1) for b in range(n_rowcol_1) for c in range(n_cols_1)), cat='Binary') + + else: + n_rows_1 = myAssgn.phi('n', obj_type)[block].get_array().shape[0] + n_cols_1 = myAssgn.psi('0', obj_type)[block].get_array().shape[1] + + # set the map product variables + map_product_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts(starting_map+'_'+obj_type+'_'+str(block), ((a, b) for a in range(n_rows_1) for b in range(n_cols_1)), cat='Integer') + + n_rowcol_2 = myAssgn.phi('n', obj_type)[block].get_array().shape[1] + + z_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts('z_'+starting_map+'_'+obj_type+'_'+str(block), ((a, b, c) for a in range(n_rows_1) for b in range(n_rowcol_2) for c in range(n_cols_1)), cat='Binary') + + # decison variables for the triangles (to make the ceiling(expression/2 work) + aux_vars = {block : {starting_map: {obj_type: {} for obj_type in ['V', 'E']} for starting_map in ['F', 'G']} for block in func_vals} + + for obj_type in ['V', 'E']: + for starting_map in ['F', 'G']: + # for block in func_vals: + for block in (func_vals[:-1] if obj_type == 'E' else func_vals): + if starting_map == 'F': + shape_m_tri = myAssgn.D('F', '2n', obj_type)[block].get_array().shape[0] + else: + shape_m_tri = myAssgn.D('G', '2n', obj_type)[block].get_array().shape[0] + + aux_vars[block][starting_map][obj_type] = pulp.LpVariable('aux_'+starting_map+'_'+obj_type+'_'+str(block), cat='Integer') + + + + + + # create the constraints + for block in func_vals: + for starting_map in ['F', 'G']: + # set the other map based on starting map + if starting_map == 'F': + other_map = 'G' + else: + other_map = 'F' + + for up_or_down in ['up', 'down']: # deals with 1 (up, down) and 2 (up, down) + + if block == func_vals[-1]: # skip the last block for this type of diagrams + continue + + #set the matrices + if up_or_down == 'up': #NOTE: the change in block indices + bou_n = myAssgn.B_up(other_map, 'n')[block].get_array() + bou_0 = myAssgn.B_up(starting_map, '0')[block].get_array() + if starting_map == 'F': + map_V = myAssgn.phi('0', 'V')[block+1].get_array() + map_E = myAssgn.phi('0', 'E')[block].get_array() + map_V_vars = phi_vars[block+1]['0']['V'] + map_E_vars = phi_vars[block]['0']['E'] + else: + map_V = myAssgn.psi('0', 'V')[block+1].get_array() + map_E = myAssgn.psi('0', 'E')[block].get_array() + map_V_vars = psi_vars[block+1]['0']['V'] + map_E_vars = psi_vars[block]['0']['E'] + else: + bou_n = myAssgn.B_down(other_map, 'n')[block].get_array() + bou_0 = myAssgn.B_down(starting_map, '0')[block].get_array() + if starting_map == 'F': + map_V = myAssgn.phi('0', 'V')[block].get_array() + map_E = myAssgn.phi('0', 'E')[block].get_array() + map_V_vars = phi_vars[block]['0']['V'] + map_E_vars = phi_vars[block]['0']['E'] + else: + map_V = myAssgn.psi('0', 'V')[block].get_array() + map_E = myAssgn.psi('0', 'E')[block].get_array() + map_V_vars = psi_vars[block]['0']['V'] + map_E_vars = psi_vars[block]['0']['E'] + + # set the dimensions + shape_m_mix = map_V.shape[0] + shape_n_mix = map_V.shape[1] + shape_o_mix = bou_n.shape[1] + shape_p_mix = map_E.shape[1] + + # constraint 1: loss is bigger than the absolute value of each matrix elements + + # changed below + for i in range(shape_m_mix): + for k in range(shape_p_mix): + # inner difference + first_term = pulp.lpSum([map_V_vars[i,j] * bou_0[j][k] for j in range(shape_n_mix)]) + second_term = pulp.lpSum([bou_n[i][l] * map_E_vars[l,k] for l in range(shape_o_mix)]) + + # total expression + mixed_expression = (first_term - second_term) + + prob += mixed_expression == 0 + + + + + # constraint 2: each column sums to 1 + for j in range(shape_n_mix): + prob += pulp.lpSum(map_V_vars[h,j] for h in range(shape_m_mix)) == 1 + + for k in range(shape_p_mix): + prob += pulp.lpSum(map_E_vars[l, k] for l in range(shape_o_mix)) == 1 + + + + for obj_type in ['V', 'E']: # deals with 3, 4, 5, 6, 7, 8, 9, 10 + if obj_type == 'E' and block == func_vals[-1]: # skip the last block for this type of diagrams. This is because we don't have an edge with the highest function value + continue + + # multiply inclusion matrices. Needed for the triangles + i_n_i_0 = myAssgn.I(starting_map, 'n', obj_type)[block].get_array() @ myAssgn.I(starting_map, '0', obj_type)[block].get_array() + + # write inclusion matrices for easier reference. Needed for the parallelograms + inc_0_para = myAssgn.I(starting_map, '0', obj_type)[block].get_array() + inc_n_para = myAssgn.I(other_map, 'n', obj_type)[block].get_array() + + + # set map matrices for easier reference + if starting_map == 'F': + map_0_para_vars = phi_vars[block]['0'][obj_type] + map_n_para_vars = phi_vars[block]['n'][obj_type] + + if starting_map == 'G': + map_0_para_vars = psi_vars[block]['0'][obj_type] + map_n_para_vars = psi_vars[block]['n'][obj_type] + + # set the dimensions + shape_m_tri = i_n_i_0.shape[0] # for triangles + + shape_o_para = myAssgn.I(other_map, 'n', obj_type)[block].get_array().shape[1] # for parallelograms + + shape_n_para = myAssgn.I(starting_map, '0', obj_type)[block].get_array().shape[0] # for parallelograms + + if starting_map == 'F': + shape_n_tri = myAssgn.psi('n', obj_type)[block].get_array().shape[1] # for triangles + shape_o_tri = myAssgn.phi('0', obj_type)[block].get_array().shape[1] # for triangles + + # changed shape_n_para = myAssgn.phi('n', obj_type)[block].get_array().shape[1] # for parallelograms + + shape_m_para = myAssgn.phi('n', obj_type)[block].get_array().shape[0] # for parallelograms + shape_p_para = myAssgn.phi('0', obj_type)[block].get_array().shape[1] # for parallelograms + else: + shape_n_tri = myAssgn.phi('n', obj_type)[block].get_array().shape[1] # for triangles + shape_o_tri = myAssgn.psi('0', obj_type)[block].get_array().shape[1] # for triangles + + # changed shape_n_para = myAssgn.psi('n', obj_type)[block].get_array().shape[1] # for parallelograms + + shape_m_para = myAssgn.psi('n', obj_type)[block].get_array().shape[0] # for parallelograms + shape_p_para = myAssgn.psi('0', obj_type)[block].get_array().shape[1] # for parallelograms + + + + + + # constraint 1: loss is bigger than the absolute value of each matrix elements + + # for triangles + + # changed below + for h in range(shape_m_tri): + + tri_expression = pulp.lpSum( (i_n_i_0[h,k] - map_product_vars[block][starting_map][obj_type][h,k]) for k in range(shape_o_tri)) + + + # prob += aux_vars[block][starting_map][obj_type] * 2 >= tri_expression # ceiling of half of the expression + + + prob += tri_expression == 0 + + + + + # for parallelograms + + # changed below + for i in range(shape_m_para): + for k in range(shape_p_para): + # inner difference + first_term = pulp.lpSum([map_n_para_vars[i,j] * inc_0_para[j][k] for j in range(shape_n_para)]) + second_term = pulp.lpSum([inc_n_para[i][l] * map_0_para_vars[l,k] for l in range(shape_o_para)]) + + + # total expression + para_expression = first_term - second_term + + prob += para_expression == 0 + + + # constraint 2: map_multiplication and z relation. This is for triangles + for i in range(shape_m_tri): + for k in range(shape_o_tri): + prob += map_product_vars[block][starting_map][obj_type][i,k] == pulp.lpSum(z_vars[block][starting_map][obj_type][i,j,k] for j in range(shape_n_tri)) + + # constraint 3: z is less than either of the map values and greater than sum of the map values minus 1. This is for triangles + for i in range(shape_m_tri): + for j in range(shape_n_tri): + for k in range(shape_o_tri): + if starting_map == 'F': + prob += z_vars[block][starting_map][obj_type][i,j,k] <= psi_vars[block]['n'][obj_type][i,j] + prob += z_vars[block][starting_map][obj_type][i,j,k] <= phi_vars[block]['0'][obj_type][j,k] + prob += z_vars[block][starting_map][obj_type][i,j,k] >= psi_vars[block]['n'][obj_type][i,j] + phi_vars[block]['0'][obj_type][j,k] - 1 + else: + prob += z_vars[block][starting_map][obj_type][i,j,k] <= phi_vars[block]['n'][obj_type][i,j] + prob += z_vars[block][starting_map][obj_type][i,j,k] <= psi_vars[block]['0'][obj_type][j,k] + prob += z_vars[block][starting_map][obj_type][i,j,k] >= phi_vars[block]['n'][obj_type][i,j] + psi_vars[block]['0'][obj_type][j,k] - 1 + + # constraint 4: each column sums to 1 + if starting_map == 'F': + # for triangles + for j in range(shape_n_tri): + prob += pulp.lpSum(psi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_tri)) == 1 + for k in range(shape_o_tri): + prob += pulp.lpSum(phi_vars[block]['0'][obj_type][j,k] for j in range(shape_n_tri)) == 1 + + # for parallelograms + for j in range(shape_n_para): + prob += pulp.lpSum(phi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_para)) == 1 + + for k in range(shape_p_para): + prob += pulp.lpSum(phi_vars[block]['0'][obj_type][j,k] for j in range(shape_o_para)) == 1 + + else: + # for triangles + for j in range(shape_n_tri): + prob += pulp.lpSum(phi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_tri)) == 1 + + for k in range(shape_o_tri): + prob += pulp.lpSum(psi_vars[block]['0'][obj_type][j,k] for j in range(shape_n_tri)) == 1 + + # for parallelograms + for j in range(shape_n_para): + prob += pulp.lpSum(psi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_para)) == 1 + + for k in range(shape_p_para): + prob += pulp.lpSum(psi_vars[block]['0'][obj_type][j,k] for j in range(shape_o_para)) == 1 + + + # Set the objective function + prob += 0 + + # solve the problem + if pulp_solver == 'GUROBI': + prob.solve(pulp.GUROBI_CMD(msg=0)) + else: + prob.solve(pulp.PULP_CBC_CMD(msg=0)) + + + # prob.solve(pulp.GUROBI_CMD(msg=0)) + if prob.status != 1: + return prob.status + + # create a dictionary to store the results + map_results = {'phi_vars': phi_vars, 'psi_vars': psi_vars} + + # make the results a LabeledBlockMatrix + final_maps = build_map_matrices(myAssgn, map_results) + + if verbose: + print("Status:", pulp.LpStatus[prob.status]) + # prob.writeLP("model.lp") # Write the model in LP format + + + # return results + return final_maps, prob.status \ No newline at end of file From 109c28a281f6dda94eb6969751f5d7eb17bc41a7 Mon Sep 17 00:00:00 2001 From: ishikaghosh2201 <112980412+ishikaghosh2201@users.noreply.github.com> Date: Fri, 5 Sep 2025 13:47:47 -0400 Subject: [PATCH 2/7] Updating the fit function (computation without dist matrix) --- cereeberus/cereeberus/distance/ilp.py | 71 +++++++------- .../distance/{ilp_without_D.py => ilp_old.py} | 70 ++++++------- cereeberus/cereeberus/distance/interleave.py | 97 ++++++++++++++++++- 3 files changed, 163 insertions(+), 75 deletions(-) rename cereeberus/cereeberus/distance/{ilp_without_D.py => ilp_old.py} (88%) diff --git a/cereeberus/cereeberus/distance/ilp.py b/cereeberus/cereeberus/distance/ilp.py index f7e805f..400193e 100644 --- a/cereeberus/cereeberus/distance/ilp.py +++ b/cereeberus/cereeberus/distance/ilp.py @@ -166,9 +166,7 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): aux_vars[block][starting_map][obj_type] = pulp.LpVariable('aux_'+starting_map+'_'+obj_type+'_'+str(block), cat='Integer') - - # create the minmax variable - minmax_var = pulp.LpVariable('minmax_var', cat='Integer') + # create the constraints @@ -187,7 +185,6 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): #set the matrices if up_or_down == 'up': #NOTE: the change in block indices - dist_n_other = myAssgn.D(other_map, 'n', 'V')[block+1].get_array() bou_n = myAssgn.B_up(other_map, 'n')[block].get_array() bou_0 = myAssgn.B_up(starting_map, '0')[block].get_array() if starting_map == 'F': @@ -201,7 +198,6 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): map_V_vars = psi_vars[block+1]['0']['V'] map_E_vars = psi_vars[block]['0']['E'] else: - dist_n_other = myAssgn.D(other_map, 'n', 'V')[block].get_array() bou_n = myAssgn.B_down(other_map, 'n')[block].get_array() bou_0 = myAssgn.B_down(starting_map, '0')[block].get_array() if starting_map == 'F': @@ -216,13 +212,14 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): map_E_vars = psi_vars[block]['0']['E'] # set the dimensions - shape_m_mix = dist_n_other.shape[0] + shape_m_mix = map_V.shape[0] shape_n_mix = map_V.shape[1] shape_o_mix = bou_n.shape[1] shape_p_mix = map_E.shape[1] # constraint 1: loss is bigger than the absolute value of each matrix elements + # changed below for i in range(shape_m_mix): for k in range(shape_p_mix): # inner difference @@ -230,9 +227,12 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): second_term = pulp.lpSum([bou_n[i][l] * map_E_vars[l,k] for l in range(shape_o_mix)]) # total expression - mixed_expression = pulp.lpSum(dist_n_other[i][h] * (first_term - second_term) for h in range(shape_m_mix)) - - prob += minmax_var >= mixed_expression + mixed_expression = (first_term - second_term) + + prob += mixed_expression == 0 + + + # constraint 2: each column sums to 1 for j in range(shape_n_mix): @@ -254,11 +254,6 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): inc_0_para = myAssgn.I(starting_map, '0', obj_type)[block].get_array() inc_n_para = myAssgn.I(other_map, 'n', obj_type)[block].get_array() - # write dist matrix for easier reference. - # for triangles - dist_2n_starting = myAssgn.D(starting_map, '2n', obj_type)[block].get_array() - # for parallelograms - dist_2n_other = myAssgn.D(other_map, '2n', obj_type)[block].get_array() # set map matrices for easier reference if starting_map == 'F': @@ -270,24 +265,27 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): map_n_para_vars = psi_vars[block]['n'][obj_type] # set the dimensions - shape_m_tri = dist_2n_starting.shape[0] # for triangles - shape_m_para = dist_2n_other.shape[0] # for parallelograms + shape_m_tri = i_n_i_0.shape[0] # for triangles shape_o_para = myAssgn.I(other_map, 'n', obj_type)[block].get_array().shape[1] # for parallelograms + shape_n_para = myAssgn.I(starting_map, '0', obj_type)[block].get_array().shape[0] # for parallelograms if starting_map == 'F': shape_n_tri = myAssgn.psi('n', obj_type)[block].get_array().shape[1] # for triangles shape_o_tri = myAssgn.phi('0', obj_type)[block].get_array().shape[1] # for triangles - shape_n_para = myAssgn.phi('n', obj_type)[block].get_array().shape[1] # for parallelograms + # changed shape_n_para = myAssgn.phi('n', obj_type)[block].get_array().shape[1] # for parallelograms + shape_m_para = myAssgn.phi('n', obj_type)[block].get_array().shape[0] # for parallelograms shape_p_para = myAssgn.phi('0', obj_type)[block].get_array().shape[1] # for parallelograms else: shape_n_tri = myAssgn.phi('n', obj_type)[block].get_array().shape[1] # for triangles shape_o_tri = myAssgn.psi('0', obj_type)[block].get_array().shape[1] # for triangles - shape_n_para = myAssgn.psi('n', obj_type)[block].get_array().shape[1] # for parallelograms + # changed shape_n_para = myAssgn.psi('n', obj_type)[block].get_array().shape[1] # for parallelograms + + shape_m_para = myAssgn.psi('n', obj_type)[block].get_array().shape[0] # for parallelograms shape_p_para = myAssgn.psi('0', obj_type)[block].get_array().shape[1] # for parallelograms @@ -297,17 +295,24 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): # constraint 1: loss is bigger than the absolute value of each matrix elements # for triangles + + # changed below for h in range(shape_m_tri): - tri_expression = pulp.lpSum(dist_2n_starting[i,h] * (i_n_i_0[h,k] - map_product_vars[block][starting_map][obj_type][h,k]) for i in range(shape_m_tri) for k in range(shape_o_tri)) + tri_expression = pulp.lpSum( (i_n_i_0[h,k] - map_product_vars[block][starting_map][obj_type][h,k]) for k in range(shape_o_tri)) - prob += aux_vars[block][starting_map][obj_type] * 2 >= tri_expression # ceiling of half of the expression - prob += minmax_var >= aux_vars[block][starting_map][obj_type] + # prob += aux_vars[block][starting_map][obj_type] * 2 >= tri_expression # ceiling of half of the expression + + + prob += tri_expression == 0 + # for parallelograms + + # changed below for i in range(shape_m_para): for k in range(shape_p_para): # inner difference @@ -316,9 +321,9 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): # total expression - para_expression = pulp.lpSum(dist_2n_other[i][h] * (first_term - second_term) for h in range(shape_m_para)) + para_expression = first_term - second_term - prob += minmax_var >= para_expression + prob += para_expression == 0 # constraint 2: map_multiplication and z relation. This is for triangles @@ -371,7 +376,7 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): # Set the objective function - prob += minmax_var + prob += 0 # solve the problem if pulp_solver == 'GUROBI': @@ -379,11 +384,12 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): else: prob.solve(pulp.PULP_CBC_CMD(msg=0)) - + status_str = pulp.LpStatus[prob.status] + # prob.solve(pulp.GUROBI_CMD(msg=0)) if prob.status != 1: - raise ValueError("The ILP optimization did not converge. Please check the input data and try again.") - + return None, status_str + # create a dictionary to store the results map_results = {'phi_vars': phi_vars, 'psi_vars': psi_vars} @@ -391,14 +397,9 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): final_maps = build_map_matrices(myAssgn, map_results) if verbose: - print(f"The optimized loss is: {pulp.value(minmax_var)}") print("Status:", pulp.LpStatus[prob.status]) - prob.writeLP("model.lp") # Write the model in LP format + # prob.writeLP("model.lp") # Write the model in LP format - # if get_thickened_maps: - # final_maps = build_map_matrices(myAssgn, map_results, thickening = 'n') - - # return final_maps, pulp.value(minmax_var) - + # return results - return final_maps, pulp.value(minmax_var) + return final_maps, status_str \ No newline at end of file diff --git a/cereeberus/cereeberus/distance/ilp_without_D.py b/cereeberus/cereeberus/distance/ilp_old.py similarity index 88% rename from cereeberus/cereeberus/distance/ilp_without_D.py rename to cereeberus/cereeberus/distance/ilp_old.py index 3cd1e51..f7e805f 100644 --- a/cereeberus/cereeberus/distance/ilp_without_D.py +++ b/cereeberus/cereeberus/distance/ilp_old.py @@ -3,7 +3,7 @@ from .labeled_blocks import LabeledMatrix as LM # import cereeberus.data.ex_mappergraphs as ex_mg -import matplotlib.pyplot as plt +import matplotlib.pyplot as plt import numpy as np import pulp #for ILP optimization @@ -166,7 +166,9 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): aux_vars[block][starting_map][obj_type] = pulp.LpVariable('aux_'+starting_map+'_'+obj_type+'_'+str(block), cat='Integer') - + + # create the minmax variable + minmax_var = pulp.LpVariable('minmax_var', cat='Integer') # create the constraints @@ -185,6 +187,7 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): #set the matrices if up_or_down == 'up': #NOTE: the change in block indices + dist_n_other = myAssgn.D(other_map, 'n', 'V')[block+1].get_array() bou_n = myAssgn.B_up(other_map, 'n')[block].get_array() bou_0 = myAssgn.B_up(starting_map, '0')[block].get_array() if starting_map == 'F': @@ -198,6 +201,7 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): map_V_vars = psi_vars[block+1]['0']['V'] map_E_vars = psi_vars[block]['0']['E'] else: + dist_n_other = myAssgn.D(other_map, 'n', 'V')[block].get_array() bou_n = myAssgn.B_down(other_map, 'n')[block].get_array() bou_0 = myAssgn.B_down(starting_map, '0')[block].get_array() if starting_map == 'F': @@ -212,14 +216,13 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): map_E_vars = psi_vars[block]['0']['E'] # set the dimensions - shape_m_mix = map_V.shape[0] + shape_m_mix = dist_n_other.shape[0] shape_n_mix = map_V.shape[1] shape_o_mix = bou_n.shape[1] shape_p_mix = map_E.shape[1] # constraint 1: loss is bigger than the absolute value of each matrix elements - # changed below for i in range(shape_m_mix): for k in range(shape_p_mix): # inner difference @@ -227,12 +230,9 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): second_term = pulp.lpSum([bou_n[i][l] * map_E_vars[l,k] for l in range(shape_o_mix)]) # total expression - mixed_expression = (first_term - second_term) - - prob += mixed_expression == 0 - - - + mixed_expression = pulp.lpSum(dist_n_other[i][h] * (first_term - second_term) for h in range(shape_m_mix)) + + prob += minmax_var >= mixed_expression # constraint 2: each column sums to 1 for j in range(shape_n_mix): @@ -254,6 +254,11 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): inc_0_para = myAssgn.I(starting_map, '0', obj_type)[block].get_array() inc_n_para = myAssgn.I(other_map, 'n', obj_type)[block].get_array() + # write dist matrix for easier reference. + # for triangles + dist_2n_starting = myAssgn.D(starting_map, '2n', obj_type)[block].get_array() + # for parallelograms + dist_2n_other = myAssgn.D(other_map, '2n', obj_type)[block].get_array() # set map matrices for easier reference if starting_map == 'F': @@ -265,27 +270,24 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): map_n_para_vars = psi_vars[block]['n'][obj_type] # set the dimensions - shape_m_tri = i_n_i_0.shape[0] # for triangles + shape_m_tri = dist_2n_starting.shape[0] # for triangles + shape_m_para = dist_2n_other.shape[0] # for parallelograms shape_o_para = myAssgn.I(other_map, 'n', obj_type)[block].get_array().shape[1] # for parallelograms - shape_n_para = myAssgn.I(starting_map, '0', obj_type)[block].get_array().shape[0] # for parallelograms if starting_map == 'F': shape_n_tri = myAssgn.psi('n', obj_type)[block].get_array().shape[1] # for triangles shape_o_tri = myAssgn.phi('0', obj_type)[block].get_array().shape[1] # for triangles - # changed shape_n_para = myAssgn.phi('n', obj_type)[block].get_array().shape[1] # for parallelograms + shape_n_para = myAssgn.phi('n', obj_type)[block].get_array().shape[1] # for parallelograms - shape_m_para = myAssgn.phi('n', obj_type)[block].get_array().shape[0] # for parallelograms shape_p_para = myAssgn.phi('0', obj_type)[block].get_array().shape[1] # for parallelograms else: shape_n_tri = myAssgn.phi('n', obj_type)[block].get_array().shape[1] # for triangles shape_o_tri = myAssgn.psi('0', obj_type)[block].get_array().shape[1] # for triangles - # changed shape_n_para = myAssgn.psi('n', obj_type)[block].get_array().shape[1] # for parallelograms - - shape_m_para = myAssgn.psi('n', obj_type)[block].get_array().shape[0] # for parallelograms + shape_n_para = myAssgn.psi('n', obj_type)[block].get_array().shape[1] # for parallelograms shape_p_para = myAssgn.psi('0', obj_type)[block].get_array().shape[1] # for parallelograms @@ -295,24 +297,17 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): # constraint 1: loss is bigger than the absolute value of each matrix elements # for triangles - - # changed below for h in range(shape_m_tri): - tri_expression = pulp.lpSum( (i_n_i_0[h,k] - map_product_vars[block][starting_map][obj_type][h,k]) for k in range(shape_o_tri)) - - - # prob += aux_vars[block][starting_map][obj_type] * 2 >= tri_expression # ceiling of half of the expression + tri_expression = pulp.lpSum(dist_2n_starting[i,h] * (i_n_i_0[h,k] - map_product_vars[block][starting_map][obj_type][h,k]) for i in range(shape_m_tri) for k in range(shape_o_tri)) + prob += aux_vars[block][starting_map][obj_type] * 2 >= tri_expression # ceiling of half of the expression - prob += tri_expression == 0 - + prob += minmax_var >= aux_vars[block][starting_map][obj_type] # for parallelograms - - # changed below for i in range(shape_m_para): for k in range(shape_p_para): # inner difference @@ -321,9 +316,9 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): # total expression - para_expression = first_term - second_term + para_expression = pulp.lpSum(dist_2n_other[i][h] * (first_term - second_term) for h in range(shape_m_para)) - prob += para_expression == 0 + prob += minmax_var >= para_expression # constraint 2: map_multiplication and z relation. This is for triangles @@ -376,7 +371,7 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): # Set the objective function - prob += 0 + prob += minmax_var # solve the problem if pulp_solver == 'GUROBI': @@ -387,8 +382,8 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): # prob.solve(pulp.GUROBI_CMD(msg=0)) if prob.status != 1: - return prob.status - + raise ValueError("The ILP optimization did not converge. Please check the input data and try again.") + # create a dictionary to store the results map_results = {'phi_vars': phi_vars, 'psi_vars': psi_vars} @@ -396,9 +391,14 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): final_maps = build_map_matrices(myAssgn, map_results) if verbose: + print(f"The optimized loss is: {pulp.value(minmax_var)}") print("Status:", pulp.LpStatus[prob.status]) - # prob.writeLP("model.lp") # Write the model in LP format + prob.writeLP("model.lp") # Write the model in LP format - + # if get_thickened_maps: + # final_maps = build_map_matrices(myAssgn, map_results, thickening = 'n') + + # return final_maps, pulp.value(minmax_var) + # return results - return final_maps, prob.status \ No newline at end of file + return final_maps, pulp.value(minmax_var) diff --git a/cereeberus/cereeberus/distance/interleave.py b/cereeberus/cereeberus/distance/interleave.py index fc4c4bc..326b998 100644 --- a/cereeberus/cereeberus/distance/interleave.py +++ b/cereeberus/cereeberus/distance/interleave.py @@ -147,6 +147,89 @@ def old_fit(self, pulp_solver = None, verbose = False, max_n_for_error = 100, return self.n def fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): + """ + Finds the smallest feasible n using exponential + binary search and stores the optimal assignment in self.assignment and the n in self.n. + + Parameters: + pulp_solver (pulp.LpSolver): + solver to use for ILP. If None, the default solver is used. + verbose (bool, optional): + If True, print the progress of the optimization. Defaults to False. + max_n_for_error (int, optional): + The maximum value of `n` to search for. If the interleaving distance is not found by this value, a ValueError is raised. Defaults to 100. ####NOTE: this can be replaced by the bounding box. + + Returns: + int: smallest feasible n + """ + + # Step 0: Check n = 0 + myAssgn = Assignment(self.F, self.G, n=0) + prob_status = myAssgn.optimize(pulp_solver=pulp_solver) + + if verbose: + print(f"\n-\nTrying n = 0...") + print(f"n = 0, status = {prob_status}") + + if prob_status == 1: + self.n = 0 + self.assignment = myAssgn + return self.n + + # Step 1: Exponential search for first feasible n + low, high = 0, 1 # last infeasible, first candidate for feasible + found_feasible_n = False + + while high <= max_n_for_error: + myAssgn = Assignment(self.F, self.G, n=high) + prob_status = myAssgn.optimize(pulp_solver=pulp_solver) + + if verbose: + print(f"\n-\nTrying n = {high}...") + print(f"n = {high}, status = {prob_status}") + + if prob_status == 1: + found_feasible_n = True + break + low, high = high, high * 2 # double n + else: + if not found_feasible_n: + raise ValueError(f"Interleaving distance not found for n <= {max_n_for_error}.") + + + high = min(high, max_n_for_error) # Clamp to max allowed + + # Step 2: Binary search for smallest feasible n in [low+1, high] + best_n = high + + while low <= high: + mid = (low + high) // 2 + myAssgn = Assignment(self.F, self.G, n=mid) + prob_status = myAssgn.optimize(pulp_solver=pulp_solver) + + if verbose: + print(f"\n-\nTrying n = {mid}...") + print(f"n = {mid}, status = {prob_status}") + + if prob_status == 1: + best_n = mid + high = mid - 1 # try smaller n + + else: + low = mid+1 # update to last infeasible n + + # Step 3: Final validation + self.n = best_n + self.assignment = Assignment(self.F, self.G, n=self.n) + prob_status = self.assignment.optimize(pulp_solver=pulp_solver) + + if prob_status != 1: + raise ValueError(f"Final fit object is not an interleaving. Status = {prob_status} for n = {self.n}.") + + return self.n + + + + def fit_with_D(self, pulp_solver = None, verbose= False, max_n_for_error = 100): """ Compute the interleaving distance between the two Mapper graphs. @@ -1699,11 +1782,15 @@ def optimize(self, pulp_solver = None): Parameters: pulp_solver (pulp.LpSolver): the solver to use for the ILP optimization. If None, the default solver is used. Returns: - float : - The loss value found by the ILP solver. + int or None: + Returns 1 if an optimal solution was found and None otherwise. + """ - - map_dict, loss_val = solve_ilp(self, pulp_solver = pulp_solver) + + map_dict, prob_status = solve_ilp(self, pulp_solver = pulp_solver) + + if prob_status != 'Optimal': + return None self.phi_['0'] = {'V': map_dict['phi_0_V'], 'E': map_dict['phi_0_E']} self.phi_['n'] = {'V': map_dict['phi_n_V'], 'E': map_dict['phi_n_E']} @@ -1711,7 +1798,7 @@ def optimize(self, pulp_solver = None): self.psi_['n'] = {'V': map_dict['psi_n_V'], 'E': map_dict['psi_n_E']} - return loss_val + return 1 From f3802b4bf9df9a05a2269b8b6cd65f2845976c30 Mon Sep 17 00:00:00 2001 From: ishikaghosh2201 <112980412+ishikaghosh2201@users.noreply.github.com> Date: Mon, 22 Sep 2025 13:42:49 -0400 Subject: [PATCH 3/7] update to interleave file --- cereeberus/cereeberus/distance/interleave.py | 386 ++++++++++++------- 1 file changed, 250 insertions(+), 136 deletions(-) diff --git a/cereeberus/cereeberus/distance/interleave.py b/cereeberus/cereeberus/distance/interleave.py index 326b998..4a3cb5d 100644 --- a/cereeberus/cereeberus/distance/interleave.py +++ b/cereeberus/cereeberus/distance/interleave.py @@ -8,6 +8,7 @@ from .labeled_blocks import LabeledBlockMatrix as LBM from .labeled_blocks import LabeledMatrix as LM from .ilp import solve_ilp +from .ilp_old import solve_ilp_old from ..compute.utils import HiddenPrints import sys, os @@ -35,116 +36,116 @@ def __init__(self, F, G): self.n = np.inf self.assignment = None - def old_fit(self, pulp_solver = None, verbose = False, max_n_for_error = 100, - ): - """ - Compute the interleaving distance between the two Mapper graphs. - - Parameters: - pulp_solver (pulp.LpSolver): - The solver to use for the ILP optimization. If None, the default solver is used. - verbose (bool, optional): - If True, print the progress of the optimization. Defaults to False. - max_n_for_error (int, optional): - The maximum value of `n` to search for. If the interleaving distance is not found by this value, a ValueError is raised. Defaults to 100. - printOptimizerOutput (bool, optional): - If True, the output of the PULP optimizer is printed. Defaults to False. - - Returns: - Interleave : - The Interleave object with the computed interleaving distance. - """ - # -- Search through possible n values to find the smallest one that still allows for interleaving - - # -- Dictionary to store the search data - # -- This will store the Loss for each value of n so we don't - # -- have to recompute it each time. - # -- The distance bound can be determined by search_data[key] = Loss implies d_I <= key + Loss - search_data = {} - - # -- Set the initial value of the distance bound to be infinity - distance_bound = np.inf - - # -- Do an expoential search for the smallest n that allows for interleaving - # -- Start with n = 1 and double it until the Loss is less than or equal to 0 - N = 1 - found_max = False + # def old_fit(self, pulp_solver = None, verbose = False, max_n_for_error = 100, + # ): + # """ + # Compute the interleaving distance between the two Mapper graphs. + + # Parameters: + # pulp_solver (pulp.LpSolver): + # The solver to use for the ILP optimization. If None, the default solver is used. + # verbose (bool, optional): + # If True, print the progress of the optimization. Defaults to False. + # max_n_for_error (int, optional): + # The maximum value of `n` to search for. If the interleaving distance is not found by this value, a ValueError is raised. Defaults to 100. + # printOptimizerOutput (bool, optional): + # If True, the output of the PULP optimizer is printed. Defaults to False. + + # Returns: + # Interleave : + # The Interleave object with the computed interleaving distance. + # """ + # # -- Search through possible n values to find the smallest one that still allows for interleaving + + # # -- Dictionary to store the search data + # # -- This will store the Loss for each value of n so we don't + # # -- have to recompute it each time. + # # -- The distance bound can be determined by search_data[key] = Loss implies d_I <= key + Loss + # search_data = {} + + # # -- Set the initial value of the distance bound to be infinity + # distance_bound = np.inf + + # # -- Do an expoential search for the smallest n that allows for interleaving + # # -- Start with n = 1 and double it until the Loss is less than or equal to 0 + # N = 1 + # found_max = False - while not found_max: - # -- Compute the assignment for the current value of n - # -- If the Loss is less than or equal to 0, we have found a valid n - # -- Otherwise, double n and try again - - try: - if verbose: - print(f"\n-\nTrying n = {N}...") - myAssgn = Assignment(self.F, self.G, n = N) + # while not found_max: + # # -- Compute the assignment for the current value of n + # # -- If the Loss is less than or equal to 0, we have found a valid n + # # -- Otherwise, double n and try again + + # try: + # if verbose: + # print(f"\n-\nTrying n = {N}...") + # myAssgn = Assignment(self.F, self.G, n = N) - Loss = myAssgn.optimize() - search_data[N] = Loss + # Loss = myAssgn.optimize() + # search_data[N] = Loss - if verbose: - print(f"n = {N}, Loss = {Loss}, distance_bound = {N + Loss}") + # if verbose: + # print(f"n = {N}, Loss = {Loss}, distance_bound = {N + Loss}") - except ValueError as e: - # -- If we get a ValueError, it means the current n is too small - # -- So we double n and try again - search_data[N] = np.inf - N *= 2 - continue + # except ValueError as e: + # # -- If we get a ValueError, it means the current n is too small + # # -- So we double n and try again + # search_data[N] = np.inf + # N *= 2 + # continue - if Loss <= 0: - # -- If the Loss is less than or equal to 0, we have found a valid n - max_N = N - found_max = True - if verbose: - print(f"Found valid maximum n = {N} with Loss = {Loss}. Moving on to binary search.") - else: - # -- If the Loss is greater than 0, we need to try a larger n - N *= 2 + # if Loss <= 0: + # # -- If the Loss is less than or equal to 0, we have found a valid n + # max_N = N + # found_max = True + # if verbose: + # print(f"Found valid maximum n = {N} with Loss = {Loss}. Moving on to binary search.") + # else: + # # -- If the Loss is greater than 0, we need to try a larger n + # N *= 2 - if N > max_n_for_error: - # -- If we have exceeded the maximum value of n, raise an error. Useful while we're checking this while loop - raise ValueError(f"Interleaving distance not found for n <= {max_n_for_error}.") + # if N > max_n_for_error: + # # -- If we have exceeded the maximum value of n, raise an error. Useful while we're checking this while loop + # raise ValueError(f"Interleaving distance not found for n <= {max_n_for_error}.") - # now binary search in [0, max_N] to find the smallest n that gives a loss of 0 + # # now binary search in [0, max_N] to find the smallest n that gives a loss of 0 - # -- Set the initial values for the binary search - low = 1 - high = max_N - while low < high: - mid = (low + high) // 2 + # # -- Set the initial values for the binary search + # low = 1 + # high = max_N + # while low < high: + # mid = (low + high) // 2 - try: - myAssgn = Assignment(self.F, self.G, n = mid) - Loss = myAssgn.optimize() - search_data[mid] = Loss - except ValueError as e: - # -- If we get a ValueError, it means the current n is too small - search_data[mid] = np.inf - low = mid + 1 - continue + # try: + # myAssgn = Assignment(self.F, self.G, n = mid) + # Loss = myAssgn.optimize() + # search_data[mid] = Loss + # except ValueError as e: + # # -- If we get a ValueError, it means the current n is too small + # search_data[mid] = np.inf + # low = mid + 1 + # continue - if Loss <= 0: - # -- If the Loss is less than or equal to 0, we have found a valid n - high = mid - else: - low = mid + 1 + # if Loss <= 0: + # # -- If the Loss is less than or equal to 0, we have found a valid n + # high = mid + # else: + # low = mid + 1 - # -- Set the final value of n to be the smallest one that gives a loss of 0 - self.n = low - # -- Set the assignment to be the one that minimizes the interleaving distance - self.assignment = Assignment(self.F, self.G, n = self.n) - Loss = self.assignment.optimize(pulp_solver = pulp_solver,) + # # -- Set the final value of n to be the smallest one that gives a loss of 0 + # self.n = low + # # -- Set the assignment to be the one that minimizes the interleaving distance + # self.assignment = Assignment(self.F, self.G, n = self.n) + # Loss = self.assignment.optimize(pulp_solver = pulp_solver,) - # Raise error if the loss isn't 0 - if Loss > 0: - raise ValueError(f"Final fit object is not an interleaving. Loss = {Loss}. N = {self.n}.") + # # Raise error if the loss isn't 0 + # if Loss > 0: + # raise ValueError(f"Final fit object is not an interleaving. Loss = {Loss}. N = {self.n}.") - return self.n + # return self.n def fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): """ @@ -229,7 +230,7 @@ def fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): - def fit_with_D(self, pulp_solver = None, verbose= False, max_n_for_error = 100): + def dist_fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): """ Compute the interleaving distance between the two Mapper graphs. @@ -248,8 +249,7 @@ def fit_with_D(self, pulp_solver = None, verbose= False, max_n_for_error = 100): # step 0: search for n=0 myAssgn = Assignment(self.F, self.G, n = 0) - Loss = myAssgn.optimize(pulp_solver = pulp_solver) - + Loss = myAssgn.dist_optimize(pulp_solver = pulp_solver) if verbose: print(f"\n-\nTrying n = 0...") print(f"n = 0, Loss = {Loss}, distance_bound = {0 + Loss}") @@ -267,7 +267,7 @@ def fit_with_D(self, pulp_solver = None, verbose= False, max_n_for_error = 100): while high <= max_n_for_error: try: myAssgn = Assignment(self.F, self.G, n = high) - Loss = myAssgn.optimize(pulp_solver = pulp_solver) + Loss = myAssgn.dist_optimize(pulp_solver = pulp_solver) if verbose: print(f"\n-\nTrying n = {high}...") @@ -293,8 +293,8 @@ def fit_with_D(self, pulp_solver = None, verbose= False, max_n_for_error = 100): mid = (low + high) // 2 try: myAssgn = Assignment(self.F, self.G, n = mid) - Loss = myAssgn.optimize(pulp_solver = pulp_solver) - + Loss = myAssgn.dist_optimize(pulp_solver = pulp_solver) + if verbose: print(f"\n-\nTrying n = {mid}...") print(f"n = {mid}, Loss = {Loss}, distance_bound = {mid + Loss}") @@ -310,8 +310,7 @@ def fit_with_D(self, pulp_solver = None, verbose= False, max_n_for_error = 100): # validate the final solution self.n = best_n self.assignment = Assignment(self.F, self.G, n = self.n) - final_loss = self.assignment.optimize(pulp_solver = pulp_solver) - + final_loss = self.assignment.dist_optimize(pulp_solver = pulp_solver) if final_loss != 0: raise ValueError(f"Unexpected non-zero loss (Loss={final_loss}) for n={self.n}") @@ -465,9 +464,12 @@ def __init__(self, F, G, # --- # Containers for matrices for later - self.B_down_ = {'F':{}, 'G':{}} # boundary matrix - self.B_up_ = {'F':{}, 'G':{}} # boundary matrix - self.D_ = {'F':{}, 'G':{}} # distance matrix + # self.B_down_ = {'F':{}, 'G':{}} # boundary matrix + self._B_down = None # don't build the boundary matrices unless needed + # self.B_up_ = {'F':{}, 'G':{}} # boundary matrix + self._B_up = None # don't build the boundary matrices unless needed + # self.D_ = {'F':{}, 'G':{}} # distance matrix + self._D = None # don't build the distance matrices unless needed self.I_ = {'F':{}, 'G':{}} # induced maps self.val_to_verts = {'F':{}, 'G':{}} # dictionaries from function values to vertices @@ -543,47 +545,51 @@ def __init__(self, F, G, # --- # Build boundary matrices - for Graph, graph_name in [(self.F, 'F'), (self.G, 'G')]: - for key in ['0', 'n']: # Note, we don't need to do this for 2n because the matrices are never used. + # for Graph, graph_name in [(self.F, 'F'), (self.G, 'G')]: + # for key in ['0', 'n']: # Note, we don't need to do this for 2n because the matrices are never used. - B_down = LBM() - B_up = LBM() + # B_down = LBM() + # B_up = LBM() - for i in self.val_to_verts[graph_name][key]: - if i in self.val_to_edges[graph_name][key]: - edges = self.val_to_edges[graph_name][key][i] - verts_down = self.val_to_verts[graph_name][key][i] - verts_up = self.val_to_verts[graph_name][key][i+1] - B_down[i] = LM(rows = verts_down, cols = edges) - B_up[i] = LM(rows = verts_up, cols = edges) + # for i in self.val_to_verts[graph_name][key]: + # if i in self.val_to_edges[graph_name][key]: + # edges = self.val_to_edges[graph_name][key][i] + # verts_down = self.val_to_verts[graph_name][key][i] + # verts_up = self.val_to_verts[graph_name][key][i+1] + # B_down[i] = LM(rows = verts_down, cols = edges) + # B_up[i] = LM(rows = verts_up, cols = edges) - for e in edges: - B_down[i][e[0],e] = 1 - B_up[i][e[1],e] = 1 + # for e in edges: + # B_down[i][e[0],e] = 1 + # B_up[i][e[1],e] = 1 - min_i = min(list(self.val_to_verts[graph_name][key].keys())) - max_i = max(list(self.val_to_verts[graph_name][key].keys())) + # min_i = min(list(self.val_to_verts[graph_name][key].keys())) + # max_i = max(list(self.val_to_verts[graph_name][key].keys())) - min_verts = self.val_to_verts[graph_name][key][min_i] - max_verts = self.val_to_verts[graph_name][key][max_i] + # min_verts = self.val_to_verts[graph_name][key][min_i] + # max_verts = self.val_to_verts[graph_name][key][max_i] - B_up[min_i-1] = LM(rows = min_verts, cols = []) - B_down[max_i] = LM(rows = max_verts, cols = []) + # B_up[min_i-1] = LM(rows = min_verts, cols = []) + # B_down[max_i] = LM(rows = max_verts, cols = []) - self.B_down_[graph_name][key] = B_down - self.B_up_[graph_name][key] = B_up + # self.B_down_[graph_name][key] = B_down + # self.B_up_[graph_name][key] = B_up # End boundary matrices # --- # --- - # Build the distance matrices - for (metagraph, name) in [ (self.F_,'F'), (self.G_,'G')]: - for key in ['n', '2n']: # Note, we don't need to do this for 0 because the matrices are never used. - self.D_[name][key] = {'V':{}, 'E':{}} + # Build the distance matrices + # We don't actually need these for the ILP, so we won't build them unless needed. + # If you need them, you can access them using the property self.D_ which will build them if they haven't been built yet. - self.D_[name][key]['V'] = metagraph[key].thickening_distance_matrix() - self.D_[name][key]['E'] = metagraph[key].thickening_distance_matrix(obj_type = 'E') + + # for (metagraph, name) in [ (self.F_,'F'), (self.G_,'G')]: + # for key in ['n', '2n']: # Note, we don't need to do this for 0 because the matrices are never used. + # self.D_[name][key] = {'V':{}, 'E':{}} + + # self.D_[name][key]['V'] = metagraph[key].thickening_distance_matrix() + # self.D_[name][key]['E'] = metagraph[key].thickening_distance_matrix(obj_type = 'E') # End distance matrices @@ -650,7 +656,95 @@ def __init__(self, F, G, # --- - + ### ---------------- + # Properties + ### ---------------- + + @property + def B_down_(self): + if self._B_down is None: + self._B_down = {'F':{}, 'G':{}} + for (Graph, graph_name) in [(self.F_,'F'), (self.G_,'G')]: + for key in ['0', 'n']: # Note, we don't need to do this for 2n because the matrices are never used. + + B_down = LBM() + # B_up = LBM() + + for i in self.val_to_verts[graph_name][key]: + if i in self.val_to_edges[graph_name][key]: + edges = self.val_to_edges[graph_name][key][i] + verts_down = self.val_to_verts[graph_name][key][i] + # verts_up = self.val_to_verts[graph_name][key][i+1] + B_down[i] = LM(rows = verts_down, cols = edges) + # B_up[i] = LM(rows = verts_up, cols = edges) + + for e in edges: + B_down[i][e[0],e] = 1 + # B_up[i][e[1],e] = 1 + + min_i = min(list(self.val_to_verts[graph_name][key].keys())) + max_i = max(list(self.val_to_verts[graph_name][key].keys())) + + min_verts = self.val_to_verts[graph_name][key][min_i] + max_verts = self.val_to_verts[graph_name][key][max_i] + + B_down[min_i-1] = LM(rows = min_verts, cols = []) + B_down[max_i] = LM(rows = max_verts, cols = []) + + self._B_down[graph_name][key] = B_down + # self.B_up_[graph_name][key] = B_up + return self._B_down + + @property + def B_up_(self): + if self._B_up is None: + self._B_up = {'F':{}, 'G':{}} + for (Graph, graph_name) in [(self.F_,'F'), (self.G_,'G')]: + for key in ['0', 'n']: # Note, we don't need to do this for 2n because the matrices are never used. + + # B_down = LBM() + B_up = LBM() + + for i in self.val_to_verts[graph_name][key]: + if i in self.val_to_edges[graph_name][key]: + edges = self.val_to_edges[graph_name][key][i] + # verts_down = self.val_to_verts[graph_name][key][i] + verts_up = self.val_to_verts[graph_name][key][i+1] + # B_down[i] = LM(rows = verts_down, cols = edges) + B_up[i] = LM(rows = verts_up, cols = edges) + + for e in edges: + # B_down[i][e[0],e] = 1 + B_up[i][e[1],e] = 1 + + min_i = min(list(self.val_to_verts[graph_name][key].keys())) + max_i = max(list(self.val_to_verts[graph_name][key].keys())) + + min_verts = self.val_to_verts[graph_name][key][min_i] + max_verts = self.val_to_verts[graph_name][key][max_i] + + B_up[min_i-1] = LM(rows = min_verts, cols = []) + B_up[max_i] = LM(rows = max_verts, cols = []) + + # self.B_down_[graph_name][key] = B_down + self._B_up[graph_name][key] = B_up + return self._B_up + + @property + def D_(self): + if self._D is None: + self._D = {'F':{}, 'G':{}} + for (metagraph, name) in [ (self.F_,'F'), (self.G_,'G')]: + for key in ['n', '2n']: # Note, we don't need to do this for 0 because the matrices are never used. + self._D[name][key] = {'V':{}, 'E':{}} + + self._D[name][key]['V'] = metagraph[key].thickening_distance_matrix() + self._D[name][key]['E'] = metagraph[key].thickening_distance_matrix(obj_type = 'E') + return self._D + + + + ### ---------------- # Functions for getting stuff out of all the dictionaries @@ -1799,6 +1893,26 @@ def optimize(self, pulp_solver = None): return 1 - - + def dist_optimize(self, pulp_solver = None): + """Uses the ILP to find the best interleaving distance bound, returns the loss value found. Further, it stores the optimal phi and psi maps which can be returned using the ``self.phi`` and ``self.psi`` attributes respectively. + ## NOTE: This version uses the D matrices in the constraints, which makes it slower. ## + This function requires the `pulp` package to be installed. + + Parameters: + pulp_solver (pulp.LpSolver): the solver to use for the ILP optimization. If None, the default solver is used. + Returns: + int or None: + Returns 1 if an optimal solution was found and None otherwise. + + """ + + map_dict, loss_val = solve_ilp_old(self, pulp_solver = pulp_solver) + + self.phi_['0'] = {'V': map_dict['phi_0_V'], 'E': map_dict['phi_0_E']} + self.phi_['n'] = {'V': map_dict['phi_n_V'], 'E': map_dict['phi_n_E']} + self.psi_['0'] = {'V': map_dict['psi_0_V'], 'E': map_dict['psi_0_E']} + self.psi_['n'] = {'V': map_dict['psi_n_V'], 'E': map_dict['psi_n_E']} + + + return loss_val \ No newline at end of file From c3e6f44483d38e2d725eeed21db1e1a83c5019f1 Mon Sep 17 00:00:00 2001 From: ishikaghosh2201 <112980412+ishikaghosh2201@users.noreply.github.com> Date: Sun, 28 Sep 2025 18:57:23 -0400 Subject: [PATCH 4/7] updating the ilp funciton --- cereeberus/cereeberus/distance/ilp.py | 372 +++++++++++++++-- cereeberus/cereeberus/distance/ilp_old.py | 404 ------------------- cereeberus/cereeberus/distance/interleave.py | 217 ++++------ 3 files changed, 412 insertions(+), 581 deletions(-) delete mode 100644 cereeberus/cereeberus/distance/ilp_old.py diff --git a/cereeberus/cereeberus/distance/ilp.py b/cereeberus/cereeberus/distance/ilp.py index 400193e..8966a42 100644 --- a/cereeberus/cereeberus/distance/ilp.py +++ b/cereeberus/cereeberus/distance/ilp.py @@ -1,7 +1,5 @@ -# from cereeberus import ReebGraph, MapperGraph, Interleave from .labeled_blocks import LabeledBlockMatrix as LBM from .labeled_blocks import LabeledMatrix as LM -# import cereeberus.data.ex_mappergraphs as ex_mg import matplotlib.pyplot as plt import numpy as np @@ -27,9 +25,6 @@ def build_map_matrices(myAssgn, map_results): Returns: dict: the dictionary containing the final map matrices """ - - # get the function values - func_vals = myAssgn.all_func_vals() # create a dictionary to store the final matrices final_LBMs = {} @@ -151,22 +146,7 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): z_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts('z_'+starting_map+'_'+obj_type+'_'+str(block), ((a, b, c) for a in range(n_rows_1) for b in range(n_rowcol_2) for c in range(n_cols_1)), cat='Binary') - # decison variables for the triangles (to make the ceiling(expression/2 work) - aux_vars = {block : {starting_map: {obj_type: {} for obj_type in ['V', 'E']} for starting_map in ['F', 'G']} for block in func_vals} - - for obj_type in ['V', 'E']: - for starting_map in ['F', 'G']: - # for block in func_vals: - for block in (func_vals[:-1] if obj_type == 'E' else func_vals): - if starting_map == 'F': - shape_m_tri = myAssgn.D('F', '2n', obj_type)[block].get_array().shape[0] - else: - shape_m_tri = myAssgn.D('G', '2n', obj_type)[block].get_array().shape[0] - - aux_vars[block][starting_map][obj_type] = pulp.LpVariable('aux_'+starting_map+'_'+obj_type+'_'+str(block), cat='Integer') - - # create the constraints @@ -219,7 +199,6 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): # constraint 1: loss is bigger than the absolute value of each matrix elements - # changed below for i in range(shape_m_mix): for k in range(shape_p_mix): # inner difference @@ -231,9 +210,6 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): prob += mixed_expression == 0 - - - # constraint 2: each column sums to 1 for j in range(shape_n_mix): prob += pulp.lpSum(map_V_vars[h,j] for h in range(shape_m_mix)) == 1 @@ -275,8 +251,8 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): shape_n_tri = myAssgn.psi('n', obj_type)[block].get_array().shape[1] # for triangles shape_o_tri = myAssgn.phi('0', obj_type)[block].get_array().shape[1] # for triangles - # changed shape_n_para = myAssgn.phi('n', obj_type)[block].get_array().shape[1] # for parallelograms - + + shape_m_para = myAssgn.phi('n', obj_type)[block].get_array().shape[0] # for parallelograms shape_p_para = myAssgn.phi('0', obj_type)[block].get_array().shape[1] # for parallelograms else: @@ -295,13 +271,10 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): # constraint 1: loss is bigger than the absolute value of each matrix elements # for triangles - - # changed below for h in range(shape_m_tri): tri_expression = pulp.lpSum( (i_n_i_0[h,k] - map_product_vars[block][starting_map][obj_type][h,k]) for k in range(shape_o_tri)) - # prob += aux_vars[block][starting_map][obj_type] * 2 >= tri_expression # ceiling of half of the expression @@ -311,8 +284,6 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): # for parallelograms - - # changed below for i in range(shape_m_para): for k in range(shape_p_para): # inner difference @@ -400,6 +371,341 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): print("Status:", pulp.LpStatus[prob.status]) # prob.writeLP("model.lp") # Write the model in LP format - + + # return results + return final_maps, status_str + + + +##------------------- ILP Optimization with Distance Matrix -------------------## + +def solve_ilp_dist(myAssgn, pulp_solver = None, verbose=False): + + """ + Function to solve the ILP optimization problem for interleaving maps. The function creates a linear programming problem using the PuLP library and solves it to find the optimal interleaving maps. + + Parameters: + myAssgn (Assignment): the Assignment object containing the interleaving maps and other relevant data + pulp_solver (pulp.LpSolver): the solver to use for the ILP optimization. If None, the default solver is used. + verbose (bool): whether to print the optimization status and results + + Returns: + tuple: a tuple containing the final interleaving maps (as a dictionary of LabeledBlockMatrices) and the optimized loss value + """ + # function values + func_vals = myAssgn.all_func_vals() + + # create the ILP problem + prob = pulp.LpProblem("Interleave_Optimization_Problem", pulp.LpMinimize) + + # create empty dictionaries to store the decision variables + phi_vars = {block : {thickening:{obj_type:{} for obj_type in ['V', 'E']} for thickening in ['0', 'n']} for block in func_vals} + psi_vars = {block : {thickening:{obj_type:{} for obj_type in ['V', 'E']} for thickening in ['0', 'n']} for block in func_vals} + + + # create the decision variables (NOTE: these are all the decision variables for all the diagrams) + for thickening in ['0', 'n']: + for obj_type in ['V', 'E']: + for block in (func_vals[:-1] if obj_type == 'E' else func_vals): + + # create lp variables for phi + n_rows = myAssgn.phi(thickening, obj_type)[block].get_array().shape[0] + n_cols = myAssgn.phi(thickening, obj_type)[block].get_array().shape[1] + + phi_vars[block][thickening][obj_type] = pulp.LpVariable.dicts('phi_'+thickening+obj_type+'_'+str(block), ((a, b) for a in range(n_rows) for b in range(n_cols)), cat='Binary') + + # # set the initial values + # for a in range(n_rows): + # for b in range(n_cols): + # prob += phi_vars[block][thickening][obj_type][(a,b)] == myAssgn.phi(thickening, obj_type)[block].get_array()[a][b] + + + # create lp variables for psi + n_rows = myAssgn.psi(thickening, obj_type)[block].get_array().shape[0] + n_cols = myAssgn.psi(thickening, obj_type)[block].get_array().shape[1] + + psi_vars[block][thickening][obj_type] = pulp.LpVariable.dicts('psi_'+thickening+obj_type+'_'+str(block), ((a, b) for a in range(n_rows) for b in range(n_cols)), cat='Binary') + + # # set the initial values + # for a in range(n_rows): + # for b in range(n_cols): + # prob += psi_vars[block][thickening][obj_type][(a,b)] == myAssgn.psi(thickening, obj_type)[block].get_array()[a][b] + + # create the other decision variables + z_vars = {block : {starting_map: {obj_type: {} for obj_type in ['V', 'E']} for starting_map in ['F', 'G']} for block in func_vals} + + map_product_vars = {block : {starting_map: {obj_type: {} for obj_type in ['V', 'E']} for starting_map in ['F', 'G']} for block in func_vals} + for obj_type in ['V', 'E']: + for starting_map in ['F', 'G']: + for block in (func_vals[:-1] if obj_type == 'E' else func_vals): + + if starting_map == 'F': + n_rows_1 = myAssgn.psi('n', obj_type)[block].get_array().shape[0] + n_cols_1 = myAssgn.phi('0', obj_type)[block].get_array().shape[1] + + # set the map product variables + map_product_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts(starting_map+'_'+obj_type+'_'+str(block), ((a, b) for a in range(n_rows_1) for b in range(n_cols_1)), cat='Integer') + + n_rowcol_1 = myAssgn.psi('n', obj_type)[block].get_array().shape[1] + + # set the z variables + z_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts('z_'+starting_map+'_'+obj_type+'_'+str(block), ((a, b, c) for a in range(n_rows_1) for b in range(n_rowcol_1) for c in range(n_cols_1)), cat='Binary') + + else: + n_rows_1 = myAssgn.phi('n', obj_type)[block].get_array().shape[0] + n_cols_1 = myAssgn.psi('0', obj_type)[block].get_array().shape[1] + + # set the map product variables + map_product_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts(starting_map+'_'+obj_type+'_'+str(block), ((a, b) for a in range(n_rows_1) for b in range(n_cols_1)), cat='Integer') + + n_rowcol_2 = myAssgn.phi('n', obj_type)[block].get_array().shape[1] + + z_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts('z_'+starting_map+'_'+obj_type+'_'+str(block), ((a, b, c) for a in range(n_rows_1) for b in range(n_rowcol_2) for c in range(n_cols_1)), cat='Binary') + + # decison variables for the triangles (to make the ceiling(expression/2 work) + aux_vars = {block : {starting_map: {obj_type: {} for obj_type in ['V', 'E']} for starting_map in ['F', 'G']} for block in func_vals} + + for obj_type in ['V', 'E']: + for starting_map in ['F', 'G']: + # for block in func_vals: + for block in (func_vals[:-1] if obj_type == 'E' else func_vals): + if starting_map == 'F': + shape_m_tri = myAssgn.D('F', '2n', obj_type)[block].get_array().shape[0] + else: + shape_m_tri = myAssgn.D('G', '2n', obj_type)[block].get_array().shape[0] + + aux_vars[block][starting_map][obj_type] = pulp.LpVariable('aux_'+starting_map+'_'+obj_type+'_'+str(block), cat='Integer') + + + + # create the minmax variable + minmax_var = pulp.LpVariable('minmax_var', cat='Integer') + + + # create the constraints + for block in func_vals: + for starting_map in ['F', 'G']: + # set the other map based on starting map + if starting_map == 'F': + other_map = 'G' + else: + other_map = 'F' + + for up_or_down in ['up', 'down']: # deals with 1 (up, down) and 2 (up, down) + + if block == func_vals[-1]: # skip the last block for this type of diagrams + continue + + #set the matrices + if up_or_down == 'up': #NOTE: the change in block indices + dist_n_other = myAssgn.D(other_map, 'n', 'V')[block+1].get_array() + bou_n = myAssgn.B_up(other_map, 'n')[block].get_array() + bou_0 = myAssgn.B_up(starting_map, '0')[block].get_array() + if starting_map == 'F': + map_V = myAssgn.phi('0', 'V')[block+1].get_array() + map_E = myAssgn.phi('0', 'E')[block].get_array() + map_V_vars = phi_vars[block+1]['0']['V'] + map_E_vars = phi_vars[block]['0']['E'] + else: + map_V = myAssgn.psi('0', 'V')[block+1].get_array() + map_E = myAssgn.psi('0', 'E')[block].get_array() + map_V_vars = psi_vars[block+1]['0']['V'] + map_E_vars = psi_vars[block]['0']['E'] + else: + dist_n_other = myAssgn.D(other_map, 'n', 'V')[block].get_array() + bou_n = myAssgn.B_down(other_map, 'n')[block].get_array() + bou_0 = myAssgn.B_down(starting_map, '0')[block].get_array() + if starting_map == 'F': + map_V = myAssgn.phi('0', 'V')[block].get_array() + map_E = myAssgn.phi('0', 'E')[block].get_array() + map_V_vars = phi_vars[block]['0']['V'] + map_E_vars = phi_vars[block]['0']['E'] + else: + map_V = myAssgn.psi('0', 'V')[block].get_array() + map_E = myAssgn.psi('0', 'E')[block].get_array() + map_V_vars = psi_vars[block]['0']['V'] + map_E_vars = psi_vars[block]['0']['E'] + + # set the dimensions + shape_m_mix = map_V.shape[0] + shape_n_mix = map_V.shape[1] + shape_o_mix = bou_n.shape[1] + shape_p_mix = map_E.shape[1] + + # constraint 1: loss is bigger than the absolute value of each matrix elements + + for i in range(shape_m_mix): + for k in range(shape_p_mix): + # inner difference + first_term = pulp.lpSum([map_V_vars[i,j] * bou_0[j][k] for j in range(shape_n_mix)]) + second_term = pulp.lpSum([bou_n[i][l] * map_E_vars[l,k] for l in range(shape_o_mix)]) + + # total expression + mixed_expression = pulp.lpSum(dist_n_other[i][h] * (first_term - second_term) for h in range(shape_m_mix)) + + prob += minmax_var >= mixed_expression + + # constraint 2: each column sums to 1 + for j in range(shape_n_mix): + prob += pulp.lpSum(map_V_vars[h,j] for h in range(shape_m_mix)) == 1 + + for k in range(shape_p_mix): + prob += pulp.lpSum(map_E_vars[l, k] for l in range(shape_o_mix)) == 1 + + + + for obj_type in ['V', 'E']: # deals with 3, 4, 5, 6, 7, 8, 9, 10 + if obj_type == 'E' and block == func_vals[-1]: # skip the last block for this type of diagrams. This is because we don't have an edge with the highest function value + continue + + # multiply inclusion matrices. Needed for the triangles + i_n_i_0 = myAssgn.I(starting_map, 'n', obj_type)[block].get_array() @ myAssgn.I(starting_map, '0', obj_type)[block].get_array() + + # write inclusion matrices for easier reference. Needed for the parallelograms + inc_0_para = myAssgn.I(starting_map, '0', obj_type)[block].get_array() + inc_n_para = myAssgn.I(other_map, 'n', obj_type)[block].get_array() + + # write dist matrix for easier reference. + # for triangles + dist_2n_starting = myAssgn.D(starting_map, '2n', obj_type)[block].get_array() + # for parallelograms + dist_2n_other = myAssgn.D(other_map, '2n', obj_type)[block].get_array() + + # set map matrices for easier reference + if starting_map == 'F': + map_0_para_vars = phi_vars[block]['0'][obj_type] + map_n_para_vars = phi_vars[block]['n'][obj_type] + + if starting_map == 'G': + map_0_para_vars = psi_vars[block]['0'][obj_type] + map_n_para_vars = psi_vars[block]['n'][obj_type] + + # set the dimensions + shape_m_tri = i_n_i_0.shape[0] # for triangles + shape_m_para = dist_2n_other.shape[0] # for parallelograms + + shape_o_para = myAssgn.I(other_map, 'n', obj_type)[block].get_array().shape[1] # for parallelograms + + + if starting_map == 'F': + shape_n_tri = myAssgn.psi('n', obj_type)[block].get_array().shape[1] # for triangles + shape_o_tri = myAssgn.phi('0', obj_type)[block].get_array().shape[1] # for triangles + + shape_n_para = myAssgn.phi('n', obj_type)[block].get_array().shape[1] # for parallelograms + + shape_p_para = myAssgn.phi('0', obj_type)[block].get_array().shape[1] # for parallelograms + else: + shape_n_tri = myAssgn.phi('n', obj_type)[block].get_array().shape[1] # for triangles + shape_o_tri = myAssgn.psi('0', obj_type)[block].get_array().shape[1] # for triangles + + shape_n_para = myAssgn.psi('n', obj_type)[block].get_array().shape[1] # for parallelograms + shape_p_para = myAssgn.psi('0', obj_type)[block].get_array().shape[1] # for parallelograms + + + + + + # constraint 1: loss is bigger than the absolute value of each matrix elements + + # for triangles + for h in range(shape_m_tri): + + tri_expression = pulp.lpSum(dist_2n_starting[i,h] * (i_n_i_0[h,k] - map_product_vars[block][starting_map][obj_type][h,k]) for i in range(shape_m_tri) for k in range(shape_o_tri)) + + prob += aux_vars[block][starting_map][obj_type] * 2 >= tri_expression # ceiling of half of the expression + + prob += minmax_var >= aux_vars[block][starting_map][obj_type] + + + + # for parallelograms + for i in range(shape_m_para): + for k in range(shape_p_para): + # inner difference + first_term = pulp.lpSum([map_n_para_vars[i,j] * inc_0_para[j][k] for j in range(shape_n_para)]) + second_term = pulp.lpSum([inc_n_para[i][l] * map_0_para_vars[l,k] for l in range(shape_o_para)]) + + + # total expression + para_expression = pulp.lpSum(dist_2n_other[i][h] * (first_term - second_term) for h in range(shape_m_para)) + + prob += minmax_var >= para_expression + + + # constraint 2: map_multiplication and z relation. This is for triangles + for i in range(shape_m_tri): + for k in range(shape_o_tri): + prob += map_product_vars[block][starting_map][obj_type][i,k] == pulp.lpSum(z_vars[block][starting_map][obj_type][i,j,k] for j in range(shape_n_tri)) + + # constraint 3: z is less than either of the map values and greater than sum of the map values minus 1. This is for triangles + for i in range(shape_m_tri): + for j in range(shape_n_tri): + for k in range(shape_o_tri): + if starting_map == 'F': + prob += z_vars[block][starting_map][obj_type][i,j,k] <= psi_vars[block]['n'][obj_type][i,j] + prob += z_vars[block][starting_map][obj_type][i,j,k] <= phi_vars[block]['0'][obj_type][j,k] + prob += z_vars[block][starting_map][obj_type][i,j,k] >= psi_vars[block]['n'][obj_type][i,j] + phi_vars[block]['0'][obj_type][j,k] - 1 + else: + prob += z_vars[block][starting_map][obj_type][i,j,k] <= phi_vars[block]['n'][obj_type][i,j] + prob += z_vars[block][starting_map][obj_type][i,j,k] <= psi_vars[block]['0'][obj_type][j,k] + prob += z_vars[block][starting_map][obj_type][i,j,k] >= phi_vars[block]['n'][obj_type][i,j] + psi_vars[block]['0'][obj_type][j,k] - 1 + + # constraint 4: each column sums to 1 + if starting_map == 'F': + # for triangles + for j in range(shape_n_tri): + prob += pulp.lpSum(psi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_tri)) == 1 + for k in range(shape_o_tri): + prob += pulp.lpSum(phi_vars[block]['0'][obj_type][j,k] for j in range(shape_n_tri)) == 1 + + # for parallelograms + for j in range(shape_n_para): + prob += pulp.lpSum(phi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_para)) == 1 + + for k in range(shape_p_para): + prob += pulp.lpSum(phi_vars[block]['0'][obj_type][j,k] for j in range(shape_o_para)) == 1 + + else: + # for triangles + for j in range(shape_n_tri): + prob += pulp.lpSum(phi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_tri)) == 1 + + for k in range(shape_o_tri): + prob += pulp.lpSum(psi_vars[block]['0'][obj_type][j,k] for j in range(shape_n_tri)) == 1 + + # for parallelograms + for j in range(shape_n_para): + prob += pulp.lpSum(psi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_para)) == 1 + + for k in range(shape_p_para): + prob += pulp.lpSum(psi_vars[block]['0'][obj_type][j,k] for j in range(shape_o_para)) == 1 + + + # Set the objective function + prob += minmax_var + + # solve the problem + if pulp_solver == 'GUROBI': + prob.solve(pulp.GUROBI_CMD(msg=0)) + else: + prob.solve(pulp.PULP_CBC_CMD(msg=0)) + + + # prob.solve(pulp.GUROBI_CMD(msg=0)) + if prob.status != 1: + raise ValueError("The ILP optimization did not converge. Please check the input data and try again.") + + # create a dictionary to store the results + map_results = {'phi_vars': phi_vars, 'psi_vars': psi_vars} + + # make the results a LabeledBlockMatrix + final_maps = build_map_matrices(myAssgn, map_results) + + if verbose: + print(f"The optimized loss is: {pulp.value(minmax_var)}") + print("Status:", pulp.LpStatus[prob.status]) + # prob.writeLP("model.lp") # Write the model in LP format + + # return results - return final_maps, status_str \ No newline at end of file + return final_maps, pulp.value(minmax_var) diff --git a/cereeberus/cereeberus/distance/ilp_old.py b/cereeberus/cereeberus/distance/ilp_old.py deleted file mode 100644 index f7e805f..0000000 --- a/cereeberus/cereeberus/distance/ilp_old.py +++ /dev/null @@ -1,404 +0,0 @@ -# from cereeberus import ReebGraph, MapperGraph, Interleave -from .labeled_blocks import LabeledBlockMatrix as LBM -from .labeled_blocks import LabeledMatrix as LM -# import cereeberus.data.ex_mappergraphs as ex_mg - -import matplotlib.pyplot as plt -import numpy as np - -import pulp #for ILP optimization - - - -# function to build the phi and psi matrices after the ILP optimization -def build_map_matrices(myAssgn, map_results): - """ - Function to build the map matrices after the ILP optimization. The end result is a dictionary with keys - - 'phi_0_V', 'phi_0_E' - - 'phi_n_V', 'phi_n_E' - - 'psi_0_V', 'psi_0_E' - - 'psi_n_V', 'psi_n_E'. - Each key corresponds to a LabeledBlockMatrix containing the respective map matrices for the given thickening. - - Parameters: - myAssgn (Assignment): the Assignment object - map_results (dict): the dictionary containing the results of the ILP optimization from the solve_ilp function - - Returns: - dict: the dictionary containing the final map matrices - """ - - # get the function values - func_vals = myAssgn.all_func_vals() - - # create a dictionary to store the final matrices - final_LBMs = {} - - for thickening in ['0', 'n']: - for obj_type in ['V', 'E']: - for map_type in ['phi', 'psi']: - - # Gets the full LBM for the relevant interleaving map - M = myAssgn.get_interleaving_map(maptype=map_type, key=thickening, obj_type=obj_type) - - for i in M.get_all_block_indices(): - block_shape = M[i].array.shape - - # Reset to zero - M[i].array = np.zeros(block_shape) - - # Set the values according to the ILP optimization results - for a in range(block_shape[0]): - for b in range(block_shape[1]): - # Get the relevant value from the ILP optimization results - if a<0 or b<0: - print(f"Warning: a={a}, b={b} for block {i} in {map_type}_{thickening}_{obj_type}.") - try: - M[i].array[a,b] = pulp.value(map_results[map_type+'_vars'][i][thickening][obj_type][(a,b)]) - except KeyError: - # TODO: This makes me nervous, but I think the issue is just that there are empty rows/columns that we want to skip - # print(f"KeyError: {map_type+'_vars'} for block {i} in {map_type}_{thickening}_{obj_type}.") - continue - - - output_key = f"{map_type}_{thickening}_{obj_type}" - final_LBMs[output_key] = M - - return final_LBMs - -##------------------- ILP Optimization -------------------## - -def solve_ilp(myAssgn, pulp_solver = None, verbose=False): - - """ - Function to solve the ILP optimization problem for interleaving maps. The function creates a linear programming problem using the PuLP library and solves it to find the optimal interleaving maps. - - Parameters: - myAssgn (Assignment): the Assignment object containing the interleaving maps and other relevant data - pulp_solver (pulp.LpSolver): the solver to use for the ILP optimization. If None, the default solver is used. - verbose (bool): whether to print the optimization status and results - - Returns: - tuple: a tuple containing the final interleaving maps (as a dictionary of LabeledBlockMatrices) and the optimized loss value - """ - # function values - func_vals = myAssgn.all_func_vals() - - # create the ILP problem - prob = pulp.LpProblem("Interleave_Optimization_Problem", pulp.LpMinimize) - - # create empty dictionaries to store the decision variables - phi_vars = {block : {thickening:{obj_type:{} for obj_type in ['V', 'E']} for thickening in ['0', 'n']} for block in func_vals} - psi_vars = {block : {thickening:{obj_type:{} for obj_type in ['V', 'E']} for thickening in ['0', 'n']} for block in func_vals} - - - # create the decision variables (NOTE: these are all the decision variables for all the diagrams) - for thickening in ['0', 'n']: - for obj_type in ['V', 'E']: - for block in (func_vals[:-1] if obj_type == 'E' else func_vals): - - # create lp variables for phi - n_rows = myAssgn.phi(thickening, obj_type)[block].get_array().shape[0] - n_cols = myAssgn.phi(thickening, obj_type)[block].get_array().shape[1] - - phi_vars[block][thickening][obj_type] = pulp.LpVariable.dicts('phi_'+thickening+obj_type+'_'+str(block), ((a, b) for a in range(n_rows) for b in range(n_cols)), cat='Binary') - - # # set the initial values - # for a in range(n_rows): - # for b in range(n_cols): - # prob += phi_vars[block][thickening][obj_type][(a,b)] == myAssgn.phi(thickening, obj_type)[block].get_array()[a][b] - - - # create lp variables for psi - n_rows = myAssgn.psi(thickening, obj_type)[block].get_array().shape[0] - n_cols = myAssgn.psi(thickening, obj_type)[block].get_array().shape[1] - - psi_vars[block][thickening][obj_type] = pulp.LpVariable.dicts('psi_'+thickening+obj_type+'_'+str(block), ((a, b) for a in range(n_rows) for b in range(n_cols)), cat='Binary') - - # # set the initial values - # for a in range(n_rows): - # for b in range(n_cols): - # prob += psi_vars[block][thickening][obj_type][(a,b)] == myAssgn.psi(thickening, obj_type)[block].get_array()[a][b] - - # create the other decision variables - z_vars = {block : {starting_map: {obj_type: {} for obj_type in ['V', 'E']} for starting_map in ['F', 'G']} for block in func_vals} - - map_product_vars = {block : {starting_map: {obj_type: {} for obj_type in ['V', 'E']} for starting_map in ['F', 'G']} for block in func_vals} - for obj_type in ['V', 'E']: - for starting_map in ['F', 'G']: - for block in (func_vals[:-1] if obj_type == 'E' else func_vals): - - if starting_map == 'F': - n_rows_1 = myAssgn.psi('n', obj_type)[block].get_array().shape[0] - n_cols_1 = myAssgn.phi('0', obj_type)[block].get_array().shape[1] - - # set the map product variables - map_product_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts(starting_map+'_'+obj_type+'_'+str(block), ((a, b) for a in range(n_rows_1) for b in range(n_cols_1)), cat='Integer') - - n_rowcol_1 = myAssgn.psi('n', obj_type)[block].get_array().shape[1] - - # set the z variables - z_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts('z_'+starting_map+'_'+obj_type+'_'+str(block), ((a, b, c) for a in range(n_rows_1) for b in range(n_rowcol_1) for c in range(n_cols_1)), cat='Binary') - - else: - n_rows_1 = myAssgn.phi('n', obj_type)[block].get_array().shape[0] - n_cols_1 = myAssgn.psi('0', obj_type)[block].get_array().shape[1] - - # set the map product variables - map_product_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts(starting_map+'_'+obj_type+'_'+str(block), ((a, b) for a in range(n_rows_1) for b in range(n_cols_1)), cat='Integer') - - n_rowcol_2 = myAssgn.phi('n', obj_type)[block].get_array().shape[1] - - z_vars[block][starting_map][obj_type] = pulp.LpVariable.dicts('z_'+starting_map+'_'+obj_type+'_'+str(block), ((a, b, c) for a in range(n_rows_1) for b in range(n_rowcol_2) for c in range(n_cols_1)), cat='Binary') - - # decison variables for the triangles (to make the ceiling(expression/2 work) - aux_vars = {block : {starting_map: {obj_type: {} for obj_type in ['V', 'E']} for starting_map in ['F', 'G']} for block in func_vals} - - for obj_type in ['V', 'E']: - for starting_map in ['F', 'G']: - # for block in func_vals: - for block in (func_vals[:-1] if obj_type == 'E' else func_vals): - if starting_map == 'F': - shape_m_tri = myAssgn.D('F', '2n', obj_type)[block].get_array().shape[0] - else: - shape_m_tri = myAssgn.D('G', '2n', obj_type)[block].get_array().shape[0] - - aux_vars[block][starting_map][obj_type] = pulp.LpVariable('aux_'+starting_map+'_'+obj_type+'_'+str(block), cat='Integer') - - - - # create the minmax variable - minmax_var = pulp.LpVariable('minmax_var', cat='Integer') - - - # create the constraints - for block in func_vals: - for starting_map in ['F', 'G']: - # set the other map based on starting map - if starting_map == 'F': - other_map = 'G' - else: - other_map = 'F' - - for up_or_down in ['up', 'down']: # deals with 1 (up, down) and 2 (up, down) - - if block == func_vals[-1]: # skip the last block for this type of diagrams - continue - - #set the matrices - if up_or_down == 'up': #NOTE: the change in block indices - dist_n_other = myAssgn.D(other_map, 'n', 'V')[block+1].get_array() - bou_n = myAssgn.B_up(other_map, 'n')[block].get_array() - bou_0 = myAssgn.B_up(starting_map, '0')[block].get_array() - if starting_map == 'F': - map_V = myAssgn.phi('0', 'V')[block+1].get_array() - map_E = myAssgn.phi('0', 'E')[block].get_array() - map_V_vars = phi_vars[block+1]['0']['V'] - map_E_vars = phi_vars[block]['0']['E'] - else: - map_V = myAssgn.psi('0', 'V')[block+1].get_array() - map_E = myAssgn.psi('0', 'E')[block].get_array() - map_V_vars = psi_vars[block+1]['0']['V'] - map_E_vars = psi_vars[block]['0']['E'] - else: - dist_n_other = myAssgn.D(other_map, 'n', 'V')[block].get_array() - bou_n = myAssgn.B_down(other_map, 'n')[block].get_array() - bou_0 = myAssgn.B_down(starting_map, '0')[block].get_array() - if starting_map == 'F': - map_V = myAssgn.phi('0', 'V')[block].get_array() - map_E = myAssgn.phi('0', 'E')[block].get_array() - map_V_vars = phi_vars[block]['0']['V'] - map_E_vars = phi_vars[block]['0']['E'] - else: - map_V = myAssgn.psi('0', 'V')[block].get_array() - map_E = myAssgn.psi('0', 'E')[block].get_array() - map_V_vars = psi_vars[block]['0']['V'] - map_E_vars = psi_vars[block]['0']['E'] - - # set the dimensions - shape_m_mix = dist_n_other.shape[0] - shape_n_mix = map_V.shape[1] - shape_o_mix = bou_n.shape[1] - shape_p_mix = map_E.shape[1] - - # constraint 1: loss is bigger than the absolute value of each matrix elements - - for i in range(shape_m_mix): - for k in range(shape_p_mix): - # inner difference - first_term = pulp.lpSum([map_V_vars[i,j] * bou_0[j][k] for j in range(shape_n_mix)]) - second_term = pulp.lpSum([bou_n[i][l] * map_E_vars[l,k] for l in range(shape_o_mix)]) - - # total expression - mixed_expression = pulp.lpSum(dist_n_other[i][h] * (first_term - second_term) for h in range(shape_m_mix)) - - prob += minmax_var >= mixed_expression - - # constraint 2: each column sums to 1 - for j in range(shape_n_mix): - prob += pulp.lpSum(map_V_vars[h,j] for h in range(shape_m_mix)) == 1 - - for k in range(shape_p_mix): - prob += pulp.lpSum(map_E_vars[l, k] for l in range(shape_o_mix)) == 1 - - - - for obj_type in ['V', 'E']: # deals with 3, 4, 5, 6, 7, 8, 9, 10 - if obj_type == 'E' and block == func_vals[-1]: # skip the last block for this type of diagrams. This is because we don't have an edge with the highest function value - continue - - # multiply inclusion matrices. Needed for the triangles - i_n_i_0 = myAssgn.I(starting_map, 'n', obj_type)[block].get_array() @ myAssgn.I(starting_map, '0', obj_type)[block].get_array() - - # write inclusion matrices for easier reference. Needed for the parallelograms - inc_0_para = myAssgn.I(starting_map, '0', obj_type)[block].get_array() - inc_n_para = myAssgn.I(other_map, 'n', obj_type)[block].get_array() - - # write dist matrix for easier reference. - # for triangles - dist_2n_starting = myAssgn.D(starting_map, '2n', obj_type)[block].get_array() - # for parallelograms - dist_2n_other = myAssgn.D(other_map, '2n', obj_type)[block].get_array() - - # set map matrices for easier reference - if starting_map == 'F': - map_0_para_vars = phi_vars[block]['0'][obj_type] - map_n_para_vars = phi_vars[block]['n'][obj_type] - - if starting_map == 'G': - map_0_para_vars = psi_vars[block]['0'][obj_type] - map_n_para_vars = psi_vars[block]['n'][obj_type] - - # set the dimensions - shape_m_tri = dist_2n_starting.shape[0] # for triangles - shape_m_para = dist_2n_other.shape[0] # for parallelograms - - shape_o_para = myAssgn.I(other_map, 'n', obj_type)[block].get_array().shape[1] # for parallelograms - - - if starting_map == 'F': - shape_n_tri = myAssgn.psi('n', obj_type)[block].get_array().shape[1] # for triangles - shape_o_tri = myAssgn.phi('0', obj_type)[block].get_array().shape[1] # for triangles - - shape_n_para = myAssgn.phi('n', obj_type)[block].get_array().shape[1] # for parallelograms - - shape_p_para = myAssgn.phi('0', obj_type)[block].get_array().shape[1] # for parallelograms - else: - shape_n_tri = myAssgn.phi('n', obj_type)[block].get_array().shape[1] # for triangles - shape_o_tri = myAssgn.psi('0', obj_type)[block].get_array().shape[1] # for triangles - - shape_n_para = myAssgn.psi('n', obj_type)[block].get_array().shape[1] # for parallelograms - shape_p_para = myAssgn.psi('0', obj_type)[block].get_array().shape[1] # for parallelograms - - - - - - # constraint 1: loss is bigger than the absolute value of each matrix elements - - # for triangles - for h in range(shape_m_tri): - - tri_expression = pulp.lpSum(dist_2n_starting[i,h] * (i_n_i_0[h,k] - map_product_vars[block][starting_map][obj_type][h,k]) for i in range(shape_m_tri) for k in range(shape_o_tri)) - - prob += aux_vars[block][starting_map][obj_type] * 2 >= tri_expression # ceiling of half of the expression - - prob += minmax_var >= aux_vars[block][starting_map][obj_type] - - - - # for parallelograms - for i in range(shape_m_para): - for k in range(shape_p_para): - # inner difference - first_term = pulp.lpSum([map_n_para_vars[i,j] * inc_0_para[j][k] for j in range(shape_n_para)]) - second_term = pulp.lpSum([inc_n_para[i][l] * map_0_para_vars[l,k] for l in range(shape_o_para)]) - - - # total expression - para_expression = pulp.lpSum(dist_2n_other[i][h] * (first_term - second_term) for h in range(shape_m_para)) - - prob += minmax_var >= para_expression - - - # constraint 2: map_multiplication and z relation. This is for triangles - for i in range(shape_m_tri): - for k in range(shape_o_tri): - prob += map_product_vars[block][starting_map][obj_type][i,k] == pulp.lpSum(z_vars[block][starting_map][obj_type][i,j,k] for j in range(shape_n_tri)) - - # constraint 3: z is less than either of the map values and greater than sum of the map values minus 1. This is for triangles - for i in range(shape_m_tri): - for j in range(shape_n_tri): - for k in range(shape_o_tri): - if starting_map == 'F': - prob += z_vars[block][starting_map][obj_type][i,j,k] <= psi_vars[block]['n'][obj_type][i,j] - prob += z_vars[block][starting_map][obj_type][i,j,k] <= phi_vars[block]['0'][obj_type][j,k] - prob += z_vars[block][starting_map][obj_type][i,j,k] >= psi_vars[block]['n'][obj_type][i,j] + phi_vars[block]['0'][obj_type][j,k] - 1 - else: - prob += z_vars[block][starting_map][obj_type][i,j,k] <= phi_vars[block]['n'][obj_type][i,j] - prob += z_vars[block][starting_map][obj_type][i,j,k] <= psi_vars[block]['0'][obj_type][j,k] - prob += z_vars[block][starting_map][obj_type][i,j,k] >= phi_vars[block]['n'][obj_type][i,j] + psi_vars[block]['0'][obj_type][j,k] - 1 - - # constraint 4: each column sums to 1 - if starting_map == 'F': - # for triangles - for j in range(shape_n_tri): - prob += pulp.lpSum(psi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_tri)) == 1 - for k in range(shape_o_tri): - prob += pulp.lpSum(phi_vars[block]['0'][obj_type][j,k] for j in range(shape_n_tri)) == 1 - - # for parallelograms - for j in range(shape_n_para): - prob += pulp.lpSum(phi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_para)) == 1 - - for k in range(shape_p_para): - prob += pulp.lpSum(phi_vars[block]['0'][obj_type][j,k] for j in range(shape_o_para)) == 1 - - else: - # for triangles - for j in range(shape_n_tri): - prob += pulp.lpSum(phi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_tri)) == 1 - - for k in range(shape_o_tri): - prob += pulp.lpSum(psi_vars[block]['0'][obj_type][j,k] for j in range(shape_n_tri)) == 1 - - # for parallelograms - for j in range(shape_n_para): - prob += pulp.lpSum(psi_vars[block]['n'][obj_type][i,j] for i in range(shape_m_para)) == 1 - - for k in range(shape_p_para): - prob += pulp.lpSum(psi_vars[block]['0'][obj_type][j,k] for j in range(shape_o_para)) == 1 - - - # Set the objective function - prob += minmax_var - - # solve the problem - if pulp_solver == 'GUROBI': - prob.solve(pulp.GUROBI_CMD(msg=0)) - else: - prob.solve(pulp.PULP_CBC_CMD(msg=0)) - - - # prob.solve(pulp.GUROBI_CMD(msg=0)) - if prob.status != 1: - raise ValueError("The ILP optimization did not converge. Please check the input data and try again.") - - # create a dictionary to store the results - map_results = {'phi_vars': phi_vars, 'psi_vars': psi_vars} - - # make the results a LabeledBlockMatrix - final_maps = build_map_matrices(myAssgn, map_results) - - if verbose: - print(f"The optimized loss is: {pulp.value(minmax_var)}") - print("Status:", pulp.LpStatus[prob.status]) - prob.writeLP("model.lp") # Write the model in LP format - - # if get_thickened_maps: - # final_maps = build_map_matrices(myAssgn, map_results, thickening = 'n') - - # return final_maps, pulp.value(minmax_var) - - # return results - return final_maps, pulp.value(minmax_var) diff --git a/cereeberus/cereeberus/distance/interleave.py b/cereeberus/cereeberus/distance/interleave.py index 4a3cb5d..d53b057 100644 --- a/cereeberus/cereeberus/distance/interleave.py +++ b/cereeberus/cereeberus/distance/interleave.py @@ -7,8 +7,7 @@ from ..compute.unionfind import UnionFind from .labeled_blocks import LabeledBlockMatrix as LBM from .labeled_blocks import LabeledMatrix as LM -from .ilp import solve_ilp -from .ilp_old import solve_ilp_old +from .ilp import solve_ilp, solve_ilp_dist from ..compute.utils import HiddenPrints import sys, os @@ -36,116 +35,7 @@ def __init__(self, F, G): self.n = np.inf self.assignment = None - # def old_fit(self, pulp_solver = None, verbose = False, max_n_for_error = 100, - # ): - # """ - # Compute the interleaving distance between the two Mapper graphs. - - # Parameters: - # pulp_solver (pulp.LpSolver): - # The solver to use for the ILP optimization. If None, the default solver is used. - # verbose (bool, optional): - # If True, print the progress of the optimization. Defaults to False. - # max_n_for_error (int, optional): - # The maximum value of `n` to search for. If the interleaving distance is not found by this value, a ValueError is raised. Defaults to 100. - # printOptimizerOutput (bool, optional): - # If True, the output of the PULP optimizer is printed. Defaults to False. - - # Returns: - # Interleave : - # The Interleave object with the computed interleaving distance. - # """ - # # -- Search through possible n values to find the smallest one that still allows for interleaving - - # # -- Dictionary to store the search data - # # -- This will store the Loss for each value of n so we don't - # # -- have to recompute it each time. - # # -- The distance bound can be determined by search_data[key] = Loss implies d_I <= key + Loss - # search_data = {} - - # # -- Set the initial value of the distance bound to be infinity - # distance_bound = np.inf - - # # -- Do an expoential search for the smallest n that allows for interleaving - # # -- Start with n = 1 and double it until the Loss is less than or equal to 0 - # N = 1 - # found_max = False - - # while not found_max: - # # -- Compute the assignment for the current value of n - # # -- If the Loss is less than or equal to 0, we have found a valid n - # # -- Otherwise, double n and try again - - # try: - # if verbose: - # print(f"\n-\nTrying n = {N}...") - # myAssgn = Assignment(self.F, self.G, n = N) - - # Loss = myAssgn.optimize() - # search_data[N] = Loss - - # if verbose: - # print(f"n = {N}, Loss = {Loss}, distance_bound = {N + Loss}") - - # except ValueError as e: - # # -- If we get a ValueError, it means the current n is too small - # # -- So we double n and try again - # search_data[N] = np.inf - # N *= 2 - # continue - - # if Loss <= 0: - # # -- If the Loss is less than or equal to 0, we have found a valid n - # max_N = N - # found_max = True - # if verbose: - # print(f"Found valid maximum n = {N} with Loss = {Loss}. Moving on to binary search.") - # else: - # # -- If the Loss is greater than 0, we need to try a larger n - # N *= 2 - - # if N > max_n_for_error: - # # -- If we have exceeded the maximum value of n, raise an error. Useful while we're checking this while loop - # raise ValueError(f"Interleaving distance not found for n <= {max_n_for_error}.") - - - # # now binary search in [0, max_N] to find the smallest n that gives a loss of 0 - - # # -- Set the initial values for the binary search - # low = 1 - # high = max_N - # while low < high: - # mid = (low + high) // 2 - - # try: - # myAssgn = Assignment(self.F, self.G, n = mid) - # Loss = myAssgn.optimize() - # search_data[mid] = Loss - # except ValueError as e: - # # -- If we get a ValueError, it means the current n is too small - # search_data[mid] = np.inf - # low = mid + 1 - # continue - - # if Loss <= 0: - # # -- If the Loss is less than or equal to 0, we have found a valid n - # high = mid - # else: - # low = mid + 1 - - # # -- Set the final value of n to be the smallest one that gives a loss of 0 - # self.n = low - # # -- Set the assignment to be the one that minimizes the interleaving distance - # self.assignment = Assignment(self.F, self.G, n = self.n) - # Loss = self.assignment.optimize(pulp_solver = pulp_solver,) - - - # # Raise error if the loss isn't 0 - # if Loss > 0: - # raise ValueError(f"Final fit object is not an interleaving. Loss = {Loss}. N = {self.n}.") - - # return self.n def fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): """ @@ -162,31 +52,50 @@ def fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): Returns: int: smallest feasible n """ - - # Step 0: Check n = 0 - myAssgn = Assignment(self.F, self.G, n=0) - prob_status = myAssgn.optimize(pulp_solver=pulp_solver) - + # catch the results to avoid recomputation + checked_results = {} + + + # Step 0: Check for smallest possible n (n=0 when they both have same function ranges) + + min_n = max(abs(self.F.min_f()-self.G.min_f()), abs(self.F.max_f()-self.G.max_f())) + + + if min_n not in checked_results: + myAssgn = Assignment(self.F, self.G, n=min_n) + prob_status = myAssgn.optimize(pulp_solver=pulp_solver) + + checked_results[min_n] = (prob_status, myAssgn) + + prob_status, myAssgn = checked_results[min_n] if verbose: - print(f"\n-\nTrying n = 0...") - print(f"n = 0, status = {prob_status}") + print(f"\n-\nTrying n = {min_n}...") + print(f"n = {min_n}, status = {prob_status}") + + if prob_status == 1: - self.n = 0 + self.n = min_n self.assignment = myAssgn return self.n # Step 1: Exponential search for first feasible n - low, high = 0, 1 # last infeasible, first candidate for feasible + low, high = min_n, min_n+1 found_feasible_n = False while high <= max_n_for_error: - myAssgn = Assignment(self.F, self.G, n=high) - prob_status = myAssgn.optimize(pulp_solver=pulp_solver) + if high not in checked_results: + myAssgn = Assignment(self.F, self.G, n=high) + prob_status = myAssgn.optimize(pulp_solver=pulp_solver) + checked_results[high] = (prob_status, myAssgn) + + prob_status, myAssgn = checked_results[high] if verbose: print(f"\n-\nTrying n = {high}...") print(f"n = {high}, status = {prob_status}") + + if prob_status == 1: found_feasible_n = True @@ -204,13 +113,18 @@ def fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): while low <= high: mid = (low + high) // 2 - myAssgn = Assignment(self.F, self.G, n=mid) - prob_status = myAssgn.optimize(pulp_solver=pulp_solver) + if mid not in checked_results: + myAssgn = Assignment(self.F, self.G, n=mid) + prob_status = myAssgn.optimize(pulp_solver=pulp_solver) + checked_results[mid] = (prob_status, myAssgn) + + prob_status, myAssgn = checked_results[mid] if verbose: print(f"\n-\nTrying n = {mid}...") print(f"n = {mid}, status = {prob_status}") + if prob_status == 1: best_n = mid high = mid - 1 # try smaller n @@ -220,8 +134,8 @@ def fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): # Step 3: Final validation self.n = best_n - self.assignment = Assignment(self.F, self.G, n=self.n) - prob_status = self.assignment.optimize(pulp_solver=pulp_solver) + prob_status, myAssgn = checked_results[self.n] + self.assignment = myAssgn if prob_status != 1: raise ValueError(f"Final fit object is not an interleaving. Status = {prob_status} for n = {self.n}.") @@ -247,9 +161,16 @@ def dist_fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): The interleaving distance, that is, the smallest n for which loss is zero. """ + # catch the results to avoid recomputation + checked_results = {} + # step 0: search for n=0 - myAssgn = Assignment(self.F, self.G, n = 0) - Loss = myAssgn.dist_optimize(pulp_solver = pulp_solver) + if 0 not in checked_results: + myAssgn = Assignment(self.F, self.G, n = 0) + Loss = myAssgn.dist_optimize(pulp_solver = pulp_solver) + checked_results[0] = (Loss, myAssgn) + Loss, myAssgn = checked_results[0] + if verbose: print(f"\n-\nTrying n = 0...") print(f"n = 0, Loss = {Loss}, distance_bound = {0 + Loss}") @@ -266,8 +187,12 @@ def dist_fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): while high <= max_n_for_error: try: - myAssgn = Assignment(self.F, self.G, n = high) - Loss = myAssgn.dist_optimize(pulp_solver = pulp_solver) + if high not in checked_results: + myAssgn = Assignment(self.F, self.G, n = high) + Loss = myAssgn.dist_optimize(pulp_solver = pulp_solver) + checked_results[high] = (Loss, myAssgn) + + Loss, myAssgn = checked_results[high] if verbose: print(f"\n-\nTrying n = {high}...") @@ -292,8 +217,12 @@ def dist_fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): while low <= high: mid = (low + high) // 2 try: - myAssgn = Assignment(self.F, self.G, n = mid) - Loss = myAssgn.dist_optimize(pulp_solver = pulp_solver) + if mid not in checked_results: + myAssgn = Assignment(self.F, self.G, n = mid) + Loss = myAssgn.dist_optimize(pulp_solver = pulp_solver) + checked_results[mid] = (Loss, myAssgn) + + Loss, myAssgn = checked_results[mid] if verbose: print(f"\n-\nTrying n = {mid}...") @@ -309,11 +238,10 @@ def dist_fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): # validate the final solution self.n = best_n - self.assignment = Assignment(self.F, self.G, n = self.n) - final_loss = self.assignment.dist_optimize(pulp_solver = pulp_solver) - if final_loss != 0: - raise ValueError(f"Unexpected non-zero loss (Loss={final_loss}) for n={self.n}") - + Loss, myAssgn = checked_results[self.n] + if Loss != 0: + raise ValueError(f"Unexpected non-zero loss (Loss={Loss}) for n={self.n}") + self.assignment = myAssgn return self.n @@ -543,7 +471,7 @@ def __init__(self, F, G, # End making smoothings and induced maps # ---- # --- - # Build boundary matrices + # # Build boundary matrices # for Graph, graph_name in [(self.F, 'F'), (self.G, 'G')]: # for key in ['0', 'n']: # Note, we don't need to do this for 2n because the matrices are never used. @@ -575,7 +503,7 @@ def __init__(self, F, G, # self.B_down_[graph_name][key] = B_down # self.B_up_[graph_name][key] = B_up - # End boundary matrices + # # End boundary matrices # --- # --- @@ -668,19 +596,15 @@ def B_down_(self): for key in ['0', 'n']: # Note, we don't need to do this for 2n because the matrices are never used. B_down = LBM() - # B_up = LBM() for i in self.val_to_verts[graph_name][key]: if i in self.val_to_edges[graph_name][key]: edges = self.val_to_edges[graph_name][key][i] verts_down = self.val_to_verts[graph_name][key][i] - # verts_up = self.val_to_verts[graph_name][key][i+1] B_down[i] = LM(rows = verts_down, cols = edges) - # B_up[i] = LM(rows = verts_up, cols = edges) for e in edges: B_down[i][e[0],e] = 1 - # B_up[i][e[1],e] = 1 min_i = min(list(self.val_to_verts[graph_name][key].keys())) max_i = max(list(self.val_to_verts[graph_name][key].keys())) @@ -1145,6 +1069,10 @@ def draw_all_graphs(self, figsize = (15,10), **kwargs): self.G('2n').draw(ax = axs[1,2], **kwargs) axs[1,2].set_title(r'$G_{2n}$') + # turn on grid for all axes + for ax in axs.flatten(): + ax.grid(True) + return fig, axs def draw_I(self, graph = 'F', key = '0', obj_type = 'V', ax = None, **kwargs): @@ -1511,6 +1439,7 @@ def parallelogram_Edge_Vert_matrix(self, Result_Dist = self.D(end_graph, 'n', 'V') @ Result else: if up_or_down == 'down': # tau_i \to \sigma_i + Top = self.get_interleaving_map(maptype, '0', 'V')[func_val] @ self.B_down(start_graph, '0')[func_val] Bottom = self.B_down(end_graph, 'n')[func_val] @ self.get_interleaving_map(maptype, '0', 'E')[func_val] Result = Top - Bottom @@ -1907,7 +1836,7 @@ def dist_optimize(self, pulp_solver = None): """ - map_dict, loss_val = solve_ilp_old(self, pulp_solver = pulp_solver) + map_dict, loss_val = solve_ilp_dist(self, pulp_solver = pulp_solver) self.phi_['0'] = {'V': map_dict['phi_0_V'], 'E': map_dict['phi_0_E']} self.phi_['n'] = {'V': map_dict['phi_n_V'], 'E': map_dict['phi_n_E']} From c3634ba399bf694d34d9ba31bf1bc39ba7dda914 Mon Sep 17 00:00:00 2001 From: ishikaghosh2201 <112980412+ishikaghosh2201@users.noreply.github.com> Date: Tue, 30 Sep 2025 17:57:48 -0400 Subject: [PATCH 5/7] n_search_fail notebook to show that the search fails for mismatched heights. --- Notebooks_Unused/n_search_fail.ipynb | 144 +++++++++++++++++++ cereeberus/cereeberus/distance/interleave.py | 19 +-- 2 files changed, 155 insertions(+), 8 deletions(-) create mode 100644 Notebooks_Unused/n_search_fail.ipynb diff --git a/Notebooks_Unused/n_search_fail.ipynb b/Notebooks_Unused/n_search_fail.ipynb new file mode 100644 index 0000000..51d2ac9 --- /dev/null +++ b/Notebooks_Unused/n_search_fail.ipynb @@ -0,0 +1,144 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 5, + "id": "47f4a073", + "metadata": {}, + "outputs": [], + "source": [ + "%reload_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "282a714d", + "metadata": {}, + "outputs": [], + "source": [ + "from cereeberus import Interleave, MapperGraph, ReebGraph\n", + "\n", + "import cereeberus.data.ex_reebgraphs as ex_rg\n", + "import cereeberus.data.ex_mappergraphs as ex_mg" + ] + }, + { + "cell_type": "markdown", + "id": "e191f088", + "metadata": {}, + "source": [ + "# Create torus and line mappers of mismatched heights" + ] + }, + { + "cell_type": "markdown", + "id": "d1a6c618", + "metadata": {}, + "source": [ + "### It works when the differnce is 1" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "6c6b8ae9", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Interleaving distance: 3\n", + "Interleaving distance (with dist_fit): 3\n" + ] + } + ], + "source": [ + "line = ex_mg.line(a=0, b= 16)\n", + "torus = ex_mg.torus(a=0, b =2, c = 13, d = 17)\n", + "myInt = Interleave(torus, line)\n", + "result = myInt.fit()\n", + "dist_result = myInt.dist_fit() # dist_fit function optimizes with the distance matrix\n", + "print(f\"Interleaving distance: {result}\")\n", + "print(f\"Interleaving distance (with dist_fit): {dist_result}\")" + ] + }, + { + "cell_type": "markdown", + "id": "30459cf0", + "metadata": {}, + "source": [ + "### But fails when the difference is more than 1" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "ee3dea41", + "metadata": {}, + "outputs": [ + { + "ename": "KeyError", + "evalue": "17", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[9], line 4\u001b[0m\n\u001b[1;32m 2\u001b[0m torus \u001b[38;5;241m=\u001b[39m ex_mg\u001b[38;5;241m.\u001b[39mtorus(a\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m, b \u001b[38;5;241m=\u001b[39m\u001b[38;5;241m2\u001b[39m, c \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m13\u001b[39m, d \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m18\u001b[39m)\n\u001b[1;32m 3\u001b[0m myInt \u001b[38;5;241m=\u001b[39m Interleave(torus, line)\n\u001b[0;32m----> 4\u001b[0m result \u001b[38;5;241m=\u001b[39m \u001b[43mmyInt\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfit\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 5\u001b[0m dist_result \u001b[38;5;241m=\u001b[39m myInt\u001b[38;5;241m.\u001b[39mdist_fit() \u001b[38;5;66;03m# dist_fit function optimizes with the distance matrix\u001b[39;00m\n\u001b[1;32m 6\u001b[0m \u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mInterleaving distance: \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mresult\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/interleave.py:66\u001b[0m, in \u001b[0;36mInterleave.fit\u001b[0;34m(self, pulp_solver, verbose, max_n_for_error)\u001b[0m\n\u001b[1;32m 64\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m min_n \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m checked_results:\n\u001b[1;32m 65\u001b[0m myAssgn \u001b[38;5;241m=\u001b[39m Assignment(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mF, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mG, n\u001b[38;5;241m=\u001b[39mmin_n)\n\u001b[0;32m---> 66\u001b[0m prob_status \u001b[38;5;241m=\u001b[39m \u001b[43mmyAssgn\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43moptimize\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpulp_solver\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mpulp_solver\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 68\u001b[0m checked_results[min_n] \u001b[38;5;241m=\u001b[39m (prob_status, myAssgn)\n\u001b[1;32m 70\u001b[0m prob_status, myAssgn \u001b[38;5;241m=\u001b[39m checked_results[min_n]\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/interleave.py:1816\u001b[0m, in \u001b[0;36mAssignment.optimize\u001b[0;34m(self, pulp_solver)\u001b[0m\n\u001b[1;32m 1804\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21moptimize\u001b[39m(\u001b[38;5;28mself\u001b[39m, pulp_solver \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m 1805\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"Uses the ILP to find the best interleaving distance bound, returns the loss value found. Further, it stores the optimal phi and psi maps which can be returned using the ``self.phi`` and ``self.psi`` attributes respectively.\u001b[39;00m\n\u001b[1;32m 1806\u001b[0m \u001b[38;5;124;03m This function requires the `pulp` package to be installed.\u001b[39;00m\n\u001b[1;32m 1807\u001b[0m \u001b[38;5;124;03m \u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 1813\u001b[0m \u001b[38;5;124;03m \u001b[39;00m\n\u001b[1;32m 1814\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m-> 1816\u001b[0m map_dict, prob_status \u001b[38;5;241m=\u001b[39m \u001b[43msolve_ilp\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpulp_solver\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mpulp_solver\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 1818\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m prob_status \u001b[38;5;241m!=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mOptimal\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[1;32m 1819\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/ilp.py:169\u001b[0m, in \u001b[0;36msolve_ilp\u001b[0;34m(myAssgn, pulp_solver, verbose)\u001b[0m\n\u001b[1;32m 167\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m up_or_down \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mup\u001b[39m\u001b[38;5;124m'\u001b[39m: \u001b[38;5;66;03m#NOTE: the change in block indices\u001b[39;00m\n\u001b[1;32m 168\u001b[0m bou_n \u001b[38;5;241m=\u001b[39m myAssgn\u001b[38;5;241m.\u001b[39mB_up(other_map, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mn\u001b[39m\u001b[38;5;124m'\u001b[39m)[block]\u001b[38;5;241m.\u001b[39mget_array()\n\u001b[0;32m--> 169\u001b[0m bou_0 \u001b[38;5;241m=\u001b[39m \u001b[43mmyAssgn\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mB_up\u001b[49m\u001b[43m(\u001b[49m\u001b[43mstarting_map\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43m0\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m[\u001b[49m\u001b[43mblock\u001b[49m\u001b[43m]\u001b[49m\u001b[38;5;241m.\u001b[39mget_array()\n\u001b[1;32m 170\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m starting_map \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mF\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[1;32m 171\u001b[0m map_V \u001b[38;5;241m=\u001b[39m myAssgn\u001b[38;5;241m.\u001b[39mphi(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m0\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mV\u001b[39m\u001b[38;5;124m'\u001b[39m)[block\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m]\u001b[38;5;241m.\u001b[39mget_array()\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/labeled_blocks.py:539\u001b[0m, in \u001b[0;36mLabeledBlockMatrix.__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 529\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__getitem__\u001b[39m(\u001b[38;5;28mself\u001b[39m, key):\n\u001b[1;32m 530\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 531\u001b[0m \u001b[38;5;124;03m Get the i'th block matrix.\u001b[39;00m\n\u001b[1;32m 532\u001b[0m \n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 537\u001b[0m \u001b[38;5;124;03m The item from the block matrix.\u001b[39;00m\n\u001b[1;32m 538\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 539\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mblocks\u001b[49m\u001b[43m[\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m]\u001b[49m\n", + "\u001b[0;31mKeyError\u001b[0m: 17" + ] + } + ], + "source": [ + "line = ex_mg.line(a =0, b= 16)\n", + "torus = ex_mg.torus(a=0, b =2, c = 13, d = 18)\n", + "myInt = Interleave(torus, line)\n", + "result = myInt.fit()\n", + "dist_result = myInt.dist_fit() # dist_fit function optimizes with the distance matrix\n", + "print(f\"Interleaving distance: {result}\")\n", + "print(f\"Interleaving distance (with dist_fit): {dist_result}\")" + ] + }, + { + "cell_type": "markdown", + "id": "b442c31a", + "metadata": {}, + "source": [ + "I checked the code and I think the problem is in the boundary matrix computation. Which affects the edge-vertex parallelogram checks. Tried fixing, but haven't figured it out yet." + ] + }, + { + "cell_type": "markdown", + "id": "c17be785", + "metadata": {}, + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "interleavingenv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/cereeberus/cereeberus/distance/interleave.py b/cereeberus/cereeberus/distance/interleave.py index d53b057..19a458b 100644 --- a/cereeberus/cereeberus/distance/interleave.py +++ b/cereeberus/cereeberus/distance/interleave.py @@ -164,16 +164,19 @@ def dist_fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): # catch the results to avoid recomputation checked_results = {} - # step 0: search for n=0 - if 0 not in checked_results: - myAssgn = Assignment(self.F, self.G, n = 0) + # Step 0: Check for smallest possible n (n=0 when they both have same function ranges) + min_n = max(abs(self.F.min_f()-self.G.min_f()), abs(self.F.max_f()-self.G.max_f())) # minimum possible n based on function ranges + + + if min_n not in checked_results: + myAssgn = Assignment(self.F, self.G, n = min_n) Loss = myAssgn.dist_optimize(pulp_solver = pulp_solver) - checked_results[0] = (Loss, myAssgn) - Loss, myAssgn = checked_results[0] + checked_results[min_n] = (Loss, myAssgn) + Loss, myAssgn = checked_results[min_n] if verbose: - print(f"\n-\nTrying n = 0...") - print(f"n = 0, Loss = {Loss}, distance_bound = {0 + Loss}") + print(f"\n-\nTrying n = {min_n}...") + print(f"n = {min_n}, Loss = {Loss}, distance_bound = {min_n + Loss}") # if loss is 0, we're done if Loss == 0: @@ -182,7 +185,7 @@ def dist_fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): return self.n # step 1: exponential search for the upperbound - low, high = 1, 1 + low, high = min_n, min_n+1 found_valid_n = False while high <= max_n_for_error: From b1c4074d208e249bb89a7b799819338784616aef Mon Sep 17 00:00:00 2001 From: ishikaghosh2201 <112980412+ishikaghosh2201@users.noreply.github.com> Date: Thu, 16 Oct 2025 09:59:18 -0400 Subject: [PATCH 6/7] Update ILP runs to include mismatched heights. --- cereeberus/cereeberus/distance/ilp.py | 4 +- cereeberus/cereeberus/distance/interleave.py | 174 ++++++++++--------- 2 files changed, 92 insertions(+), 86 deletions(-) diff --git a/cereeberus/cereeberus/distance/ilp.py b/cereeberus/cereeberus/distance/ilp.py index 8966a42..575e93a 100644 --- a/cereeberus/cereeberus/distance/ilp.py +++ b/cereeberus/cereeberus/distance/ilp.py @@ -152,6 +152,8 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): # create the constraints for block in func_vals: for starting_map in ['F', 'G']: + if block not in myAssgn.all_func_vals(map=starting_map): # If the block is not in the starting map, skip it. This is to account for the graphs with different function values + continue # set the other map based on starting map if starting_map == 'F': other_map = 'G' @@ -160,7 +162,7 @@ def solve_ilp(myAssgn, pulp_solver = None, verbose=False): for up_or_down in ['up', 'down']: # deals with 1 (up, down) and 2 (up, down) - if block == func_vals[-1]: # skip the last block for this type of diagrams + if block == myAssgn.all_func_vals(map=starting_map)[-1]: # skip the last block for this type of diagrams continue #set the matrices diff --git a/cereeberus/cereeberus/distance/interleave.py b/cereeberus/cereeberus/distance/interleave.py index 19a458b..37c3ce0 100644 --- a/cereeberus/cereeberus/distance/interleave.py +++ b/cereeberus/cereeberus/distance/interleave.py @@ -180,7 +180,7 @@ def dist_fit(self, pulp_solver = None, verbose= False, max_n_for_error = 100): # if loss is 0, we're done if Loss == 0: - self.n = 0 + self.n = min_n self.assignment = myAssgn return self.n @@ -395,10 +395,10 @@ def __init__(self, F, G, # --- # Containers for matrices for later - # self.B_down_ = {'F':{}, 'G':{}} # boundary matrix - self._B_down = None # don't build the boundary matrices unless needed - # self.B_up_ = {'F':{}, 'G':{}} # boundary matrix - self._B_up = None # don't build the boundary matrices unless needed + self.B_down_ = {'F':{}, 'G':{}} # boundary matrix + # self._B_down = None # don't build the boundary matrices unless needed + self.B_up_ = {'F':{}, 'G':{}} # boundary matrix + # self._B_up = None # don't build the boundary matrices unless needed # self.D_ = {'F':{}, 'G':{}} # distance matrix self._D = None # don't build the distance matrices unless needed self.I_ = {'F':{}, 'G':{}} # induced maps @@ -474,39 +474,39 @@ def __init__(self, F, G, # End making smoothings and induced maps # ---- # --- - # # Build boundary matrices + # Build boundary matrices - # for Graph, graph_name in [(self.F, 'F'), (self.G, 'G')]: - # for key in ['0', 'n']: # Note, we don't need to do this for 2n because the matrices are never used. + for Graph, graph_name in [(self.F, 'F'), (self.G, 'G')]: + for key in ['0', 'n']: # Note, we don't need to do this for 2n because the matrices are never used. - # B_down = LBM() - # B_up = LBM() + B_down = LBM() + B_up = LBM() - # for i in self.val_to_verts[graph_name][key]: - # if i in self.val_to_edges[graph_name][key]: - # edges = self.val_to_edges[graph_name][key][i] - # verts_down = self.val_to_verts[graph_name][key][i] - # verts_up = self.val_to_verts[graph_name][key][i+1] - # B_down[i] = LM(rows = verts_down, cols = edges) - # B_up[i] = LM(rows = verts_up, cols = edges) + for i in self.val_to_verts[graph_name][key]: + if i in self.val_to_edges[graph_name][key]: + edges = self.val_to_edges[graph_name][key][i] + verts_down = self.val_to_verts[graph_name][key][i] + verts_up = self.val_to_verts[graph_name][key][i+1] + B_down[i] = LM(rows = verts_down, cols = edges) + B_up[i] = LM(rows = verts_up, cols = edges) - # for e in edges: - # B_down[i][e[0],e] = 1 - # B_up[i][e[1],e] = 1 + for e in edges: + B_down[i][e[0],e] = 1 + B_up[i][e[1],e] = 1 - # min_i = min(list(self.val_to_verts[graph_name][key].keys())) - # max_i = max(list(self.val_to_verts[graph_name][key].keys())) + min_i = min(list(self.val_to_verts[graph_name][key].keys())) + max_i = max(list(self.val_to_verts[graph_name][key].keys())) - # min_verts = self.val_to_verts[graph_name][key][min_i] - # max_verts = self.val_to_verts[graph_name][key][max_i] + min_verts = self.val_to_verts[graph_name][key][min_i] + max_verts = self.val_to_verts[graph_name][key][max_i] - # B_up[min_i-1] = LM(rows = min_verts, cols = []) - # B_down[max_i] = LM(rows = max_verts, cols = []) + B_up[min_i-1] = LM(rows = min_verts, cols = []) + B_down[max_i] = LM(rows = max_verts, cols = []) - # self.B_down_[graph_name][key] = B_down - # self.B_up_[graph_name][key] = B_up + self.B_down_[graph_name][key] = B_down + self.B_up_[graph_name][key] = B_up - # # End boundary matrices + # End boundary matrices # --- # --- @@ -591,71 +591,71 @@ def __init__(self, F, G, # Properties ### ---------------- - @property - def B_down_(self): - if self._B_down is None: - self._B_down = {'F':{}, 'G':{}} - for (Graph, graph_name) in [(self.F_,'F'), (self.G_,'G')]: - for key in ['0', 'n']: # Note, we don't need to do this for 2n because the matrices are never used. + # @property + # def B_down_(self): + # if self._B_down is None: + # self._B_down = {'F':{}, 'G':{}} + # for (Graph, graph_name) in [(self.F_,'F'), (self.G_,'G')]: + # for key in ['0', 'n']: # Note, we don't need to do this for 2n because the matrices are never used. - B_down = LBM() + # B_down = LBM() - for i in self.val_to_verts[graph_name][key]: - if i in self.val_to_edges[graph_name][key]: - edges = self.val_to_edges[graph_name][key][i] - verts_down = self.val_to_verts[graph_name][key][i] - B_down[i] = LM(rows = verts_down, cols = edges) + # for i in self.val_to_verts[graph_name][key]: + # if i in self.val_to_edges[graph_name][key]: + # edges = self.val_to_edges[graph_name][key][i] + # verts_down = self.val_to_verts[graph_name][key][i] + # B_down[i] = LM(rows = verts_down, cols = edges) - for e in edges: - B_down[i][e[0],e] = 1 + # for e in edges: + # B_down[i][e[0],e] = 1 - min_i = min(list(self.val_to_verts[graph_name][key].keys())) - max_i = max(list(self.val_to_verts[graph_name][key].keys())) + # min_i = min(list(self.val_to_verts[graph_name][key].keys())) + # max_i = max(list(self.val_to_verts[graph_name][key].keys())) - min_verts = self.val_to_verts[graph_name][key][min_i] - max_verts = self.val_to_verts[graph_name][key][max_i] + # min_verts = self.val_to_verts[graph_name][key][min_i] + # max_verts = self.val_to_verts[graph_name][key][max_i] - B_down[min_i-1] = LM(rows = min_verts, cols = []) - B_down[max_i] = LM(rows = max_verts, cols = []) + # B_down[min_i-1] = LM(rows = min_verts, cols = []) + # B_down[max_i] = LM(rows = max_verts, cols = []) - self._B_down[graph_name][key] = B_down - # self.B_up_[graph_name][key] = B_up - return self._B_down + # self._B_down[graph_name][key] = B_down + # # self.B_up_[graph_name][key] = B_up + # return self._B_down - @property - def B_up_(self): - if self._B_up is None: - self._B_up = {'F':{}, 'G':{}} - for (Graph, graph_name) in [(self.F_,'F'), (self.G_,'G')]: - for key in ['0', 'n']: # Note, we don't need to do this for 2n because the matrices are never used. + # @property + # def B_up_(self): + # if self._B_up is None: + # self._B_up = {'F':{}, 'G':{}} + # for (Graph, graph_name) in [(self.F_,'F'), (self.G_,'G')]: + # for key in ['0', 'n']: # Note, we don't need to do this for 2n because the matrices are never used. - # B_down = LBM() - B_up = LBM() + # # B_down = LBM() + # B_up = LBM() - for i in self.val_to_verts[graph_name][key]: - if i in self.val_to_edges[graph_name][key]: - edges = self.val_to_edges[graph_name][key][i] - # verts_down = self.val_to_verts[graph_name][key][i] - verts_up = self.val_to_verts[graph_name][key][i+1] - # B_down[i] = LM(rows = verts_down, cols = edges) - B_up[i] = LM(rows = verts_up, cols = edges) + # for i in self.val_to_verts[graph_name][key]: + # if i in self.val_to_edges[graph_name][key]: + # edges = self.val_to_edges[graph_name][key][i] + # # verts_down = self.val_to_verts[graph_name][key][i] + # verts_up = self.val_to_verts[graph_name][key][i+1] + # # B_down[i] = LM(rows = verts_down, cols = edges) + # B_up[i] = LM(rows = verts_up, cols = edges) - for e in edges: - # B_down[i][e[0],e] = 1 - B_up[i][e[1],e] = 1 + # for e in edges: + # # B_down[i][e[0],e] = 1 + # B_up[i][e[1],e] = 1 - min_i = min(list(self.val_to_verts[graph_name][key].keys())) - max_i = max(list(self.val_to_verts[graph_name][key].keys())) + # min_i = min(list(self.val_to_verts[graph_name][key].keys())) + # max_i = max(list(self.val_to_verts[graph_name][key].keys())) - min_verts = self.val_to_verts[graph_name][key][min_i] - max_verts = self.val_to_verts[graph_name][key][max_i] + # min_verts = self.val_to_verts[graph_name][key][min_i] + # max_verts = self.val_to_verts[graph_name][key][max_i] - B_up[min_i-1] = LM(rows = min_verts, cols = []) - B_up[max_i] = LM(rows = max_verts, cols = []) + # B_up[min_i-1] = LM(rows = min_verts, cols = []) + # B_up[max_i] = LM(rows = max_verts, cols = []) - # self.B_down_[graph_name][key] = B_down - self._B_up[graph_name][key] = B_up - return self._B_up + # # self.B_down_[graph_name][key] = B_down + # self._B_up[graph_name][key] = B_up + # return self._B_up @property def D_(self): @@ -1786,7 +1786,7 @@ def loss_by_block(self): return max(loss_list) - def all_func_vals(self): + def all_func_vals(self, map = None): """ Get all the function values that are in the graphs. @@ -1794,12 +1794,16 @@ def all_func_vals(self): list : A list of all the function values. """ + if map == 'F': + return self.F().get_function_values() + elif map == 'G': + return self.G().get_function_values() + else: + all_func_vals = set(self.F().get_function_values()) | set(self.G().get_function_values()) + all_func_vals = list(all_func_vals) + all_func_vals.sort() - all_func_vals = set(self.F().get_function_values()) | set(self.G().get_function_values()) - all_func_vals = list(all_func_vals) - all_func_vals.sort() - - return all_func_vals + return all_func_vals def optimize(self, pulp_solver = None): """Uses the ILP to find the best interleaving distance bound, returns the loss value found. Further, it stores the optimal phi and psi maps which can be returned using the ``self.phi`` and ``self.psi`` attributes respectively. From 1fb287b5efb244e07266d7f409d22dd683742dcc Mon Sep 17 00:00:00 2001 From: ishikaghosh2201 <112980412+ishikaghosh2201@users.noreply.github.com> Date: Mon, 3 Nov 2025 10:33:55 -0500 Subject: [PATCH 7/7] Minor fixes --- .../run_exprimetns_loop_line.ipynb | 43 ++++- .../sandbox_fixing_optimization.ipynb | 149 ++++++++++++++++++ 2 files changed, 189 insertions(+), 3 deletions(-) create mode 100644 Notebooks_Unused/sandbox_fixing_optimization.ipynb diff --git a/Notebooks_Unused/experiments/experiments_with_line_torus/run_exprimetns_loop_line.ipynb b/Notebooks_Unused/experiments/experiments_with_line_torus/run_exprimetns_loop_line.ipynb index db085b2..bc1e0d9 100644 --- a/Notebooks_Unused/experiments/experiments_with_line_torus/run_exprimetns_loop_line.ipynb +++ b/Notebooks_Unused/experiments/experiments_with_line_torus/run_exprimetns_loop_line.ipynb @@ -38,13 +38,50 @@ }, { "cell_type": "code", - "execution_count": 31, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ - "import cereeberus.distance.ilp_solver_iterations as solver_iter\n", + "from cereeberus import Interleave, MapperGraph\n", "import matplotlib.pyplot as plt\n", - "import math" + "import math\n", + "from cereeberus.data import ex_mappergraphs" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "torus = ex_mappergraphs.torus(0, 2, 10, 13)\n", + "line = ex_mappergraphs.line(0, 16)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "ename": "KeyError", + "evalue": "14", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mKeyError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[8], line 2\u001b[0m\n\u001b[1;32m 1\u001b[0m myInt \u001b[38;5;241m=\u001b[39m Interleave(torus, line)\n\u001b[0;32m----> 2\u001b[0m \u001b[43mmyInt\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mfit\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/interleave.py:66\u001b[0m, in \u001b[0;36mInterleave.fit\u001b[0;34m(self, pulp_solver, verbose, max_n_for_error)\u001b[0m\n\u001b[1;32m 64\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m min_n \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;129;01min\u001b[39;00m checked_results:\n\u001b[1;32m 65\u001b[0m myAssgn \u001b[38;5;241m=\u001b[39m Assignment(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mF, \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mG, n\u001b[38;5;241m=\u001b[39mmin_n)\n\u001b[0;32m---> 66\u001b[0m prob_status \u001b[38;5;241m=\u001b[39m \u001b[43mmyAssgn\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43moptimize\u001b[49m\u001b[43m(\u001b[49m\u001b[43mpulp_solver\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mpulp_solver\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 68\u001b[0m checked_results[min_n] \u001b[38;5;241m=\u001b[39m (prob_status, myAssgn)\n\u001b[1;32m 70\u001b[0m prob_status, myAssgn \u001b[38;5;241m=\u001b[39m checked_results[min_n]\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/interleave.py:1816\u001b[0m, in \u001b[0;36mAssignment.optimize\u001b[0;34m(self, pulp_solver)\u001b[0m\n\u001b[1;32m 1804\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21moptimize\u001b[39m(\u001b[38;5;28mself\u001b[39m, pulp_solver \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m 1805\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"Uses the ILP to find the best interleaving distance bound, returns the loss value found. Further, it stores the optimal phi and psi maps which can be returned using the ``self.phi`` and ``self.psi`` attributes respectively.\u001b[39;00m\n\u001b[1;32m 1806\u001b[0m \u001b[38;5;124;03m This function requires the `pulp` package to be installed.\u001b[39;00m\n\u001b[1;32m 1807\u001b[0m \u001b[38;5;124;03m \u001b[39;00m\n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 1813\u001b[0m \u001b[38;5;124;03m \u001b[39;00m\n\u001b[1;32m 1814\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m-> 1816\u001b[0m map_dict, prob_status \u001b[38;5;241m=\u001b[39m \u001b[43msolve_ilp\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mpulp_solver\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mpulp_solver\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 1818\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m prob_status \u001b[38;5;241m!=\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mOptimal\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[1;32m 1819\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/ilp.py:169\u001b[0m, in \u001b[0;36msolve_ilp\u001b[0;34m(myAssgn, pulp_solver, verbose)\u001b[0m\n\u001b[1;32m 167\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m up_or_down \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mup\u001b[39m\u001b[38;5;124m'\u001b[39m: \u001b[38;5;66;03m#NOTE: the change in block indices\u001b[39;00m\n\u001b[1;32m 168\u001b[0m bou_n \u001b[38;5;241m=\u001b[39m myAssgn\u001b[38;5;241m.\u001b[39mB_up(other_map, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mn\u001b[39m\u001b[38;5;124m'\u001b[39m)[block]\u001b[38;5;241m.\u001b[39mget_array()\n\u001b[0;32m--> 169\u001b[0m bou_0 \u001b[38;5;241m=\u001b[39m \u001b[43mmyAssgn\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mB_up\u001b[49m\u001b[43m(\u001b[49m\u001b[43mstarting_map\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[38;5;124;43m0\u001b[39;49m\u001b[38;5;124;43m'\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m[\u001b[49m\u001b[43mblock\u001b[49m\u001b[43m]\u001b[49m\u001b[38;5;241m.\u001b[39mget_array()\n\u001b[1;32m 170\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m starting_map \u001b[38;5;241m==\u001b[39m \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mF\u001b[39m\u001b[38;5;124m'\u001b[39m:\n\u001b[1;32m 171\u001b[0m map_V \u001b[38;5;241m=\u001b[39m myAssgn\u001b[38;5;241m.\u001b[39mphi(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m0\u001b[39m\u001b[38;5;124m'\u001b[39m, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124mV\u001b[39m\u001b[38;5;124m'\u001b[39m)[block\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m]\u001b[38;5;241m.\u001b[39mget_array()\n", + "File \u001b[0;32m~/Library/CloudStorage/OneDrive-MichiganStateUniversity/Documents/Research/Projects/Interleaving_to_ML/ceREEBerus/cereeberus/cereeberus/distance/labeled_blocks.py:539\u001b[0m, in \u001b[0;36mLabeledBlockMatrix.__getitem__\u001b[0;34m(self, key)\u001b[0m\n\u001b[1;32m 529\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21m__getitem__\u001b[39m(\u001b[38;5;28mself\u001b[39m, key):\n\u001b[1;32m 530\u001b[0m \u001b[38;5;250m \u001b[39m\u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[1;32m 531\u001b[0m \u001b[38;5;124;03m Get the i'th block matrix.\u001b[39;00m\n\u001b[1;32m 532\u001b[0m \n\u001b[0;32m (...)\u001b[0m\n\u001b[1;32m 537\u001b[0m \u001b[38;5;124;03m The item from the block matrix.\u001b[39;00m\n\u001b[1;32m 538\u001b[0m \u001b[38;5;124;03m \"\"\"\u001b[39;00m\n\u001b[0;32m--> 539\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mblocks\u001b[49m\u001b[43m[\u001b[49m\u001b[43mkey\u001b[49m\u001b[43m]\u001b[49m\n", + "\u001b[0;31mKeyError\u001b[0m: 14" + ] + } + ], + "source": [ + "myInt = Interleave(torus, line)\n", + "myInt.fit()" ] }, { diff --git a/Notebooks_Unused/sandbox_fixing_optimization.ipynb b/Notebooks_Unused/sandbox_fixing_optimization.ipynb new file mode 100644 index 0000000..b7a9323 --- /dev/null +++ b/Notebooks_Unused/sandbox_fixing_optimization.ipynb @@ -0,0 +1,149 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "3d7c8bc8", + "metadata": {}, + "outputs": [], + "source": [ + "%reload_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "a632ce73", + "metadata": {}, + "outputs": [], + "source": [ + "from cereeberus import Interleave, MapperGraph, Assignment\n", + "import cereeberus.data.ex_mappergraphs as ex_mg\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import networkx as nx" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "id": "44844a7b", + "metadata": {}, + "outputs": [], + "source": [ + "mg1 = ex_mg.line(0, 20)\n", + "mg2 = ex_mg.torus(0, 2, 17, 18)" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "0657e22b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "-\n", + "Trying n = 2...\n", + "n = 2, status = None\n", + "\n", + "-\n", + "Trying n = 3...\n", + "n = 3, status = None\n", + "\n", + "-\n", + "Trying n = 6...\n", + "n = 6, status = 1\n", + "\n", + "-\n", + "Trying n = 4...\n", + "n = 4, status = 1\n", + "\n", + "-\n", + "Trying n = 3...\n", + "n = 3, status = None\n" + ] + }, + { + "data": { + "text/plain": [ + "4" + ] + }, + "execution_count": 26, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "myInt = Interleave(mg1, mg2)\n", + "myInt.fit(verbose=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "6a5a4786", + "metadata": {}, + "outputs": [], + "source": [ + "from cereeberus import ReebGraph\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "88cb17ae", + "metadata": {}, + "outputs": [ + { + "ename": "TypeError", + "evalue": "ReebGraph.__init__() got an unexpected keyword argument 'values'", + "output_type": "error", + "traceback": [ + "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[0;31mTypeError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[0;32mIn[4], line 8\u001b[0m\n\u001b[1;32m 5\u001b[0m values \u001b[38;5;241m=\u001b[39m points[:, \u001b[38;5;241m2\u001b[39m]\n\u001b[1;32m 7\u001b[0m \u001b[38;5;66;03m# Build and compute the Reeb graph\u001b[39;00m\n\u001b[0;32m----> 8\u001b[0m rg \u001b[38;5;241m=\u001b[39m \u001b[43mReebGraph\u001b[49m\u001b[43m(\u001b[49m\u001b[43mvalues\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mvalues\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcoords\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mpoints\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m 9\u001b[0m rg\u001b[38;5;241m.\u001b[39mcompute()\n", + "\u001b[0;31mTypeError\u001b[0m: ReebGraph.__init__() got an unexpected keyword argument 'values'" + ] + } + ], + "source": [ + "# Example point cloud (you can replace this with your own)\n", + "points = np.random.rand(500, 3) # 500 points in 3D\n", + "\n", + "# Define a scalar function on these points — say, height (z-coordinate)\n", + "values = points[:, 2]\n", + "\n", + "# Build and compute the Reeb graph\n", + "rg = ReebGraph()\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "interleavingenv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.3" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}