add: calculate AIC when solving the model fitting problem

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)
+    m_opt, c0_opt, b_opt, aic_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
 
     # 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)
+    m_reg, c_reg, aic_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
 
     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)
+    m_opt, c0_opt, b_opt, aic_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
 
     # 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)
+    m_reg, c_reg, aic_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
 
     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)
+    m_opt, c0_opt, b_opt, aic_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
 
     # 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)
+    m_reg, c_reg, aic_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
 
     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)
+    m_opt, c0_opt, b_opt, aic_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
 
     # 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)
+    m_reg, c_reg, aic_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
 
     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')

requirements.txt

@@ -1,4 +1,4 @@
-gurobipy
+gurobi
 graph-tools
 numpy
 matplotlib

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():
+    return beta.X, alpha.X, aic
-        print(f"{v.VarName} {v.X:g}")
-
-    return beta.X, alpha.X
 
 
 def plot_piece_wise_linear(d, C, m_opt, c0_opt, b_opt, measure, n, t, path):