diff --git a/example.py b/example.py index 607c13e..56eb590 100644 --- a/example.py +++ b/example.py @@ -42,6 +42,46 @@ def random_graph(n=5000, seed=None): return rng.random((n, 2)), seed +def random_graph_favor_border(n=3000, seed = None): + if seed is None: + import secrets + seed = secrets.randbits(128) + rng = np.random.default_rng(seed=seed) + vps = np.zeros((n, 2)) + for i in range(0, n): + r_x = rng.random() + if rng.random() > 0.5: + while (r_x > 0.3 and r_x < 0.7): + r_x = rng.random() + r_y = rng.random() + if rng.random() > 0.5: + while (r_y > 0.3 and r_y < 0.7): + r_y = rng.random() + vps[i][0] = r_x + vps[i][1] = r_y + return vps, seed + + +def random_graph_favor_center(n=3000, seed = None): + if seed is None: + import secrets + seed = secrets.randbits(128) + rng = np.random.default_rng(seed=seed) + vps = np.zeros((n, 2)) + for i in range(0, n): + r_x = rng.random() + if rng.random() > 0.7: + while (r_x < 0.4 or r_x > 0.6): + r_x = rng.random() + r_y = rng.random() + if rng.random() > 0.7: + while (r_y < 0.4 or r_y > 0.6): + r_y = rng.random() + vps[i][0] = r_x + vps[i][1] = r_y + return vps, seed + + def spatial_graph(adata): """ Generate the spatial graph using delaunay for the given `adata`. @@ -87,16 +127,7 @@ def merfish_example(): # optimize model's piece-wise linear function d = quantification[:, 0] C = quantification[:, 1] - m_opt, c0_opt, b_opt = fitting.fit_piece_wise_linear(d, C) - - # AIC - # AIC = 2 * k (= 2) - 2 * ln(L^~) - # with L^~ = sum(f(x_i)) where x_i describes a data point - # - f is *not normalized* - sum_log = 0.0 - for x_i in x_spatial: - sum_log += math.log(m_opt* b_opt + c0_opt if x_i >= b_opt else m_opt * x_i + c0_opt) - aic_model = 6. - 2. * sum_log # three parameters: b_opt, m_opt, c0_opt + m_opt, c0_opt, b_opt, aic_opt = fitting.fit_piece_wise_linear(d, C) # TODO # should this be part of the plotting function itself, it should not be necessary for me to do this @@ -107,25 +138,18 @@ def merfish_example(): [lambda x: m_opt * x + c0_opt, lambda x: m_opt * b_opt + c0_opt] ) # plot model containing modeled piece-wise linear function - plot.quantification_plot(ax1, quantification, d_curve, C_curve, 'Models', aic_model) + plot.quantification_plot(ax1, quantification, d_curve, C_curve, 'Models', aic_opt) # linear regression model - m_reg, c_reg = fitting.fit_linear_regression(d, C) - - # AIC - sum_log = 0.0 - for x_i in x_spatial: - sum_log += math.log(m_reg * x_i + c_reg) - aic_regression = 4. - 2. * sum_log # two parameter: m_reg, c_reg + m_reg, c_reg, aic_reg = fitting.fit_linear_regression(d, C) x = np.linspace(min(d), max(d), 500) y = m_reg * x + c_reg - ax1.plot(x, y, color='k', linewidth=1, label=f"Simple Linear Regression | AIC: {aic_regression}") + ax1.plot(x, y, color='k', linewidth=1, label=f"Simple Linear Regression | AIC: {aic_reg}") ax1.legend() fig.savefig(f"Merfish_closeness.svg", format='svg') - for i in range(1, 6): points, seed = random_graph() g, weight = spatial_graph(points) @@ -158,16 +182,7 @@ for i in range(1, 6): # optimize model's piece-wise linear function d = quantification[:, 0] C = quantification[:, 1] - m_opt, c0_opt, b_opt = fitting.fit_piece_wise_linear(d, C) - - # AIC - # AIC = 2 * k (= 2) - 2 * ln(L^~) - # with L^~ = sum(f(x_i)) where x_i describes a data point - # - f is *not normalized* - sum_log = 0.0 - for x_i in x_spatial: - sum_log += math.log(m_opt* b_opt + c0_opt if x_i >= b_opt else m_opt * x_i + c0_opt) - aic_model = 6. - 2. * sum_log # three parameters: b_opt, m_opt, c0_opt + m_opt, c0_opt, b_opt, aic_opt = fitting.fit_piece_wise_linear(d, C) # TODO # should this be part of the plotting function itself, it should not be necessary for me to do this @@ -178,20 +193,14 @@ for i in range(1, 6): [lambda x: m_opt * x + c0_opt, lambda x: m_opt * b_opt + c0_opt] ) # plot model containing modeled piece-wise linear function - plot.quantification_plot(ax1, quantification, d_curve, C_curve, 'Models', aic_model) + plot.quantification_plot(ax1, quantification, d_curve, C_curve, 'Models', aic_opt) # linear regression model - m_reg, c_reg = fitting.fit_linear_regression(d, C) - - # AIC - sum_log = 0.0 - for x_i in x_spatial: - sum_log += math.log(m_reg * x_i + c_reg) - aic_regression = 4. - 2. * sum_log # two parameter: m_reg, c_reg + m_reg, c_reg, aic_reg = fitting.fit_linear_regression(d, C) x = np.linspace(min(d), max(d), 500) y = m_reg * x + c_reg - ax1.plot(x, y, color='k', linewidth=1, label=f"Simple Linear Regression | AIC: {aic_regression}") + ax1.plot(x, y, color='k', linewidth=1, label=f"Simple Linear Regression | AIC: {aic_reg}") ax1.legend() fig.savefig(f"uniform_random_point_clouds/{i}_closeness.svg", format='svg') @@ -221,16 +230,7 @@ for i in range(1, 6): # optimize model's piece-wise linear function d = quantification[:, 0] C = quantification[:, 1] - m_opt, c0_opt, b_opt = fitting.fit_piece_wise_linear(d, C) - - # AIC - # AIC = 2 * k (= 2) - 2 * ln(L^~) - # with L^~ = sum(f(x_i)) where x_i describes a data point - # - f is *not normalized* - sum_log = 0.0 - for x_i in x_spatial: - sum_log += math.log(m_opt* b_opt + c0_opt if x_i >= b_opt else m_opt * x_i + c0_opt) - aic_model = 6. - 2. * sum_log # three parameters: b_opt, m_opt, c0_opt + m_opt, c0_opt, b_opt, aic_opt = fitting.fit_piece_wise_linear(d, C) # TODO # should this be part of the plotting function itself, it should not be necessary for me to do this @@ -241,20 +241,14 @@ for i in range(1, 6): [lambda x: m_opt * x + c0_opt, lambda x: m_opt * b_opt + c0_opt] ) # plot model containing modeled piece-wise linear function - plot.quantification_plot(ax1, quantification, d_curve, C_curve, 'Models', aic_model) + plot.quantification_plot(ax1, quantification, d_curve, C_curve, 'Models', aic_opt) # linear regression model - m_reg, c_reg = fitting.fit_linear_regression(d, C) - - # AIC - sum_log = 0.0 - for x_i in x_spatial: - sum_log += math.log(m_reg * x_i + c_reg) - aic_regression = 4. - 2. * sum_log # two parameter: m_reg, c_reg + m_reg, c_reg, aic_reg = fitting.fit_linear_regression(d, C) x = np.linspace(min(d), max(d), 500) y = m_reg * x + c_reg - ax1.plot(x, y, color='k', linewidth=1, label=f"Simple Linear Regression | AIC: {aic_regression}") + ax1.plot(x, y, color='k', linewidth=1, label=f"Simple Linear Regression | AIC: {aic_reg}") ax1.legend() fig.savefig(f"uniform_random_point_clouds/{i}_betweenness.svg", format='svg') @@ -284,16 +278,7 @@ for i in range(1, 6): # optimize model's piece-wise linear function d = quantification[:, 0] C = quantification[:, 1] - m_opt, c0_opt, b_opt = fitting.fit_piece_wise_linear(d, C) - - # AIC - # AIC = 2 * k (= 2) - 2 * ln(L^~) - # with L^~ = sum(f(x_i)) where x_i describes a data point - # - f is *not normalized* - sum_log = 0.0 - for x_i in x_spatial: - sum_log += math.log(m_opt* b_opt + c0_opt if x_i >= b_opt else m_opt * x_i + c0_opt) - aic_model = 6. - 2. * sum_log # three parameters: b_opt, m_opt, c0_opt + m_opt, c0_opt, b_opt, aic_opt = fitting.fit_piece_wise_linear(d, C) # TODO # should this be part of the plotting function itself, it should not be necessary for me to do this @@ -304,20 +289,14 @@ for i in range(1, 6): [lambda x: m_opt * x + c0_opt, lambda x: m_opt * b_opt + c0_opt] ) # plot model containing modeled piece-wise linear function - plot.quantification_plot(ax1, quantification, d_curve, C_curve, 'Models', aic_model) + plot.quantification_plot(ax1, quantification, d_curve, C_curve, 'Models', aic_opt) # linear regression model - m_reg, c_reg = fitting.fit_linear_regression(d, C) - - # AIC - sum_log = 0.0 - for x_i in x_spatial: - sum_log += math.log(m_reg * x_i + c_reg) - aic_regression = 4. - 2. * sum_log # two parameter: m_reg, c_reg + m_reg, c_reg, aic_reg = fitting.fit_linear_regression(d, C) x = np.linspace(min(d), max(d), 500) y = m_reg * x + c_reg - ax1.plot(x, y, color='k', linewidth=1, label=f"Simple Linear Regression | AIC: {aic_regression}") + ax1.plot(x, y, color='k', linewidth=1, label=f"Simple Linear Regression | AIC: {aic_reg}") ax1.legend() fig.savefig(f"uniform_random_point_clouds/{i}_pagerank.svg", format='svg') diff --git a/requirements.txt b/requirements.txt index 72b2ed4..1f0f298 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,4 +1,4 @@ -gurobipy +gurobi graph-tools numpy matplotlib diff --git a/src/fitting.py b/src/fitting.py index a3daf74..9b0c141 100644 --- a/src/fitting.py +++ b/src/fitting.py @@ -1,3 +1,5 @@ +import math + import gurobipy as gp from gurobipy import GRB import numpy as np @@ -49,7 +51,11 @@ def fit_piece_wise_linear(d, C, M=1000): model.addConstr((1 - z[i]) * M >= d[i] - b) model.optimize() - return m.X, c0.X, b.X + # AIC + k = 4 + aic = 2. * k + n * math.log(model.ObjVal) + + return m.X, c0.X, b.X, aic def fit_linear_regression(d, C): @@ -71,11 +77,11 @@ def fit_linear_regression(d, C): model.setObjective(gp.quicksum((C[i] - alpha - beta * d[i])**2 for i in range(n)), GRB.MINIMIZE) model.optimize() + # AIC + k = 2 + aic = 2. * k + n * math.log(model.ObjVal) - for v in model.getVars(): - print(f"{v.VarName} {v.X:g}") - - return beta.X, alpha.X + return beta.X, alpha.X, aic def plot_piece_wise_linear(d, C, m_opt, c0_opt, b_opt, measure, n, t, path):