diff --git a/crates/applications/stacked_cores/scripts/stacked_cores_52xx_plotters.py b/crates/applications/stacked_cores/scripts/stacked_cores_52xx_plotters.py index f933822..9061d48 100644 --- a/crates/applications/stacked_cores/scripts/stacked_cores_52xx_plotters.py +++ b/crates/applications/stacked_cores/scripts/stacked_cores_52xx_plotters.py @@ -2,6 +2,7 @@ from math import sqrt import plotly.express as px from pandas import DataFrame +import scipy.optimize as opt unit_to_m = lambda u: -17000 + 34000 * u sweep_1d = lambda points=101: [unit_to_m(x/(points-1)) for x in range(points)] @@ -36,7 +37,8 @@ def sample_all(meas: 'ParameterizedMeas', at: tuple, extract_tx) -> tuple: runs = [extract_tx(r) for r in meas.runs()] distances = [(distance_to(m, at), m) for m in runs] - interpolated = weighted_sum_of_neighbors_by_inv_distance(distances) + # interpolated = weighted_sum_of_neighbors_by_inv_distance(distances) + interpolated = interpolate_minl1(at, runs, distances) print(at, interpolated) return interpolated @@ -46,6 +48,44 @@ def extract_52xx_tx(meas_rows: list) -> tuple: """ return (meas_rows[0].m[0], meas_rows[0].m[1], meas_rows[1].m[2], meas_rows[2].m[3]) +def interpolate_minl1(at: tuple, runs: list, distances: list) -> tuple: + # let R = `runs`, A = `at`, D = `distances`, x be the weight of each run + # such that the result is R0 x0 + R1 x1 + ... + # + # solve for x + # subject to R0 x0 + R1 x1 + ... = A for the elements of A != None + # minimize D0 x0 + D1 x1 + ... + # + # relevant scipy docs: + # - + # - + fixed_coords = [(i, a) for (i, a) in enumerate(at) if a is not None] + num_fixed_coords = len(fixed_coords) + num_runs = len(runs) + + eq_constraints = [[0]*num_runs for _ in range(num_fixed_coords)] + for run_idx, run in enumerate(runs): + for constraint_idx, (coord_idx, _a) in enumerate(fixed_coords): + eq_constraints[constraint_idx][run_idx] = run[coord_idx] + eq_rhs = [a for (i, a) in fixed_coords] + eq_constraint = opt.LinearConstraint(eq_constraints, eq_rhs, eq_rhs) + + # constrain the weights to be positive + bounds = opt.Bounds([0]*num_runs, [float("Inf")]*num_runs) + + def score(weights: list) -> float: + # function to minimize: D0 x0 + D1 x1 + ... + return sum(w*d[0] for w, d in zip(weights, distances)) + + # compute the weight of each run + init = [0]*num_runs + res = opt.minimize(score, init, method='trust-constr', constraints=[eq_constraint], bounds=bounds) + run_weights = res.x + + # sum the weighted runs + return element_sum([weighted(run, weight) for run, weight in zip(runs, run_weights)]) + + def interpolate(meas: 'ParameterizedMeas', a0: float, a1: float) -> tuple: """ this interpolates a point among four neighboring points in 2d. diff --git a/flake.nix b/flake.nix index c299200..67cc7a3 100644 --- a/flake.nix +++ b/flake.nix @@ -18,6 +18,7 @@ natsort pandas plotly + scipy ]; python3 = pkgs.python3.withPackages python-packages; in