From 3b8274beffadb670f426b6a55f26532acf66a45b Mon Sep 17 00:00:00 2001 From: Olivernaslund <133633664+Olivernaslund@users.noreply.github.com> Date: Wed, 4 Mar 2026 17:10:12 +0100 Subject: [PATCH 1/7] python implementation of NOW --- constraints.py | 250 ++++++++++++++++++++++++++++++++ now_example.py | 45 ++++++ optimizationProblem.py | 130 +++++++++++++++++ optimize.py | 320 +++++++++++++++++++++++++++++++++++++++++ 4 files changed, 745 insertions(+) create mode 100644 constraints.py create mode 100644 now_example.py create mode 100644 optimizationProblem.py create mode 100644 optimize.py diff --git a/constraints.py b/constraints.py new file mode 100644 index 0000000..1024cc1 --- /dev/null +++ b/constraints.py @@ -0,0 +1,250 @@ +import numpy as np +from optimizationProblem import NOW_config +from scipy.optimize import LinearConstraint, NonlinearConstraint +from scipy.sparse import csr_matrix, csc_matrix, lil_array, tril, diags, eye, vstack, hstack, kron +import torch + + +class constraints: + def __init__(self, problem = NOW_config(), method=None): + self.problem = problem + self.method = method + + self.dt = self.problem._dt + + self.eps = self.problem._tolIsotropy + self.N = self.problem.N + + self.indices = self.problem.zeroGradientAtIndex + self.signs = self.problem.signs + + # Targets + self.Bhat = self.problem.targetTensor + self.gMax = self.problem._gMaxConstraint + self.sMax = self.problem._sMaxConstraint + self.integralMax = self.problem._integralConstraint + self.tolMaxwell = self.problem._tolMaxwell + + #Integration and derivative matrices + self.Theta = build_integration_matrix(self.problem.N, 1) + self.A1 = build_first_derivative_matrix(self.problem.N, 1) + self.A2 = build_second_derivative_matrix(self.problem.N, 1) + + def unpack(self, x): + q = x[:-1].reshape(-1, 3, order='F') + b = x[-1] + return q, b + + def all_linear(self): + #Zeroing on interval and boundary + A = np.zeros((2 + len(self.indices), self.N)) + A[0,0] = 1 ; A[1,-1] = 1 + A[2:, :] = self.A1.toarray()[self.indices, :] + A = np.kron(np.eye(3), A) + A = np.hstack([A, np.zeros((A.shape[0], 1))]) + A_1 = csc_matrix(A) ; b_1 = np.zeros((A_1.shape[0])) + + #Max gradient strength + if self.problem.useMaxNorm: + A1 = self.A1.toarray() ; b = np.ones(int(self.N-1))*self.gMax + A = np.kron(np.eye(3), A1) + A = np.hstack([A, np.zeros((A.shape[0], 1))]) + A_2 = csc_matrix(A) ; b_2 = np.repeat(b,3) + else: + A_2 = csc_matrix(np.empty((0, self.N*3 + 1))) ; b_2 = [] + + #Max slew rate + A1 = self.A2.toarray() ; b = np.ones(int(self.N))*self.sMax + A = np.kron(np.eye(3), A1) + A = np.hstack([A, np.zeros((A.shape[0], 1))]) + A_3 = csc_matrix(A) ; b_3 = np.repeat(b,3) + + #Linear Motion compensation + if np.any(self.problem.motionCompensation['linear']): + t = (np.arange(1, self.N + 1) - 0.5)*self.dt + linear_ind = np.where(self.problem.motionCompensation['linear'])[0] + A = np.zeros((len(linear_ind), self.N)) + for i in range(len(linear_ind)): + order = self.problem.motionCompensation['order'][linear_ind[i]] + A[i,:] = - order * self.dt * t**(order - 1) + A = np.kron(np.eye(3), A) + A = np.hstack([A, np.zeros((A.shape[0], 1))]) + A_4 = csc_matrix(A) ; b_4 = np.zeros((A_4.shape[0])) + else: + A_4 = csc_matrix(np.empty((0, self.N*3 + 1))) ; b_4 = [] + + #Stacking + A = vstack([A_1, A_2, A_3, A_4]) + #A = A.toarray() # used if using trust-constr as jac need to match (dense or sparse) + b = np.concatenate([b_1, b_2, b_3, b_4]) + if self.method == "SLSQP": return [dict( type='ineq', fun=lambda x: b - A @ x ), dict( type='ineq', fun=lambda x: b + A @ x )] + else: return [LinearConstraint(A, -b, b)] + + def all_nonlinear(self): + def fun(x): + q, b = self.unpack(x) ; N = self.N + g = self.A1 @ q + B = q.T @ (self.Theta/(N-1)) @ q + + #Tensor encoding + c1 = np.array([ (self.eps * b)**2 - np.trace((B-b*self.Bhat).T @ (B-b*self.Bhat)) ]) # np.array([(self.eps * b)**2 - np.linalg.norm(B - b * self.Bhat, 'fro')**2 ]) + + #Max gradient strength + if not self.problem.useMaxNorm: + c2 = (self.gMax**2 - np.sum(g**2, axis = 1)).T + else: + c2 = np.array([], dtype=float) + + #Power constraint + c3 = np.array([self.integralMax - g[:,0].T @ g[:,0]]) + c4 = np.array([self.integralMax - g[:,1].T @ g[:,1]]) + c5 = np.array([self.integralMax - g[:,2].T @ g[:,2]]) + + #Maxwell + if self.problem.doMaxwellComp: + M = g.T@(g*self.signs) + m = np.sqrt(np.trace(M.T @ M)) + c6 = np.array([self.tolMaxwell - m]) + else: + c6 = np.array([], dtype=float) + + #Non linear Motion compensation + if not np.all(self.problem.motionCompensation['linear']): + nonlinear_ind = np.where(~self.problem.motionCompensation["linear"])[0] + t = (np.arange(1, self.N + 1) - 0.5)*self.dt + gamma = 2.6751e+08 # radians / T / s for hydrogen. + + A = np.zeros((1,len(nonlinear_ind))) + for i in range(len(nonlinear_ind)): + order = self.problem.motionCompensation['order'][nonlinear_ind[i]] + moment_weighting = - order*self.dt*t**(order - 1) + #print(moment_weighting.shape, q.shape) + moment_vector = moment_weighting @ q + A[0,i] = (self.problem.motionCompensation['maxMagnitude'][nonlinear_ind[i]]*1000**order / (gamma * 1e-6))**2 - np.sum(moment_vector**2) + c7 = A.ravel() + else: + c7 = np.array([], dtype=float) + + return np.hstack([c1,c2,c3,c4,c5,c6,c7]) + def jac(x): #(m,n) + q, b = self.unpack(x) ; N = self.N + g = self.A1 @ q + wq = (self.Theta/(N-1)) @ q + B = q.T @ wq + + #Tensor encoding + dc1_db = 2*B - 2*b*self.Bhat + firstTerm = np.kron(np.eye(3),wq.T) + secondTerm = np.reshape(firstTerm, (3,3,3*N), order='F') + secondTerm = np.transpose(secondTerm, (1,0,2)) + secondTerm = np.reshape(secondTerm, (9,3*N)) + dB_dq = firstTerm + secondTerm + dc1_dq = np.reshape(dc1_db, (1,9))@dB_dq + dc1_ds = np.array([ 2*b*np.trace(self.Bhat.T@self.Bhat) - 2*np.trace(B.T@self.Bhat) - 2*b*self.eps**2 ]).reshape(1,1) + dc1_dx = - np.hstack([dc1_dq, dc1_ds]) + + #Max gradient strength + if not self.problem.useMaxNorm: + #dc2_dq = 2*np.vstack([self.A1 @ g[:,0], self.A1 @ g[:,1], self.A1 @ g[:,2]]) + dc2_dq = 2 * np.hstack([ + self.A1.toarray() * g[:N-1, 0][:, None], # g[0:N-1, 0] -> (N-1,1) + self.A1.toarray() * g[:N-1, 1][:, None], + self.A1.toarray() * g[:N-1, 2][:, None], + ]) + + """ + dc2_dq = 2*np.vstack([diags(g[:,0], shape = (N-1, N-1))@self.A1, diags(g[:,1], shape = (N-1, N-1))@self.A1 , diags(g[:,2], shape = (N-1, N-1))@self.A1 ]) + test= diags([g[:,0]], shape = (N-1, N-1))@self.A1 + print(dc2_dq.shape, test.shape, "wtf") + dc2_dx = - np.hstack([dc2_dq.T, np.zeros((N-1, 1))]) + print(dc2_dq.shape, dc2_dx.shape, test.shape, "wtf") + """ + dc2_dx = - np.hstack([dc2_dq, np.zeros((N-1, 1))]) + else: + dc2_dx = np.empty((0,3*N + 1)) + + #Power constraint + dc3_dx = np.zeros((1,3*N+1)) + dc3_dx[0,0:N] = -2*g[:,0].T @ self.A1 + dc4_dx = np.zeros((1,3*N+1)) + dc4_dx[0,N:2*N] = -2*g[:,1].T @ self.A1 + dc5_dx = np.zeros((1,3*N+1)) + dc5_dx[0,2*N:-1] = -2*g[:,2].T @ self.A1 + + #Maxwell + if self.problem.doMaxwellComp: + signedg = g*self.signs + M = g.T@signedg + m = np.sqrt(np.trace(M.T@M)) + dc6_dM = 1/m * M + firstTerm = np.kron(np.eye(3), (self.A1.T @ signedg ).T ) + secondTerm = np.reshape(firstTerm, (3,3,3*N), order='F') + secondTerm = np.transpose(secondTerm, (1,0,2)) + secondTerm = np.reshape(secondTerm, (9,3*N)) + dM_dq = firstTerm + secondTerm + dc6_dq = np.reshape(dc6_dM, (1,9)) @ dM_dq + #print(dc6_dq.shape) + dc6_dx = - np.hstack([dc6_dq, np.zeros((1, 1))]) + + else: + dc6_dx = np.empty((0,3*N + 1)) + + #Non linear Motion compensation + if not np.all(self.problem.motionCompensation['linear']): + nonlinear_ind = np.where(~self.problem.motionCompensation["linear"])[0] + t = (np.arange(1, self.N + 1) - 0.5)*self.dt + gamma = 2.6751e+08 # radians / T / s for hydrogen. + + A = np.zeros((len(nonlinear_ind),3*N+1)) + for i in range(len(nonlinear_ind)): + order = self.problem.motionCompensation['order'][nonlinear_ind[i]] + moment_weighting = - order*self.dt*t**(order - 1) + moment_vector = moment_weighting @ q + A[i,:-1] = -2 * np.kron(moment_vector, moment_weighting) + dc7_dx = A + else: + dc7_dx = np.empty((0,3*N + 1)) + + return csc_matrix(np.vstack([dc1_dx,dc2_dx,dc3_dx,dc4_dx,dc5_dx,dc6_dx,dc7_dx])) + + + if self.method == "SLSQP": return [dict( type='ineq', fun=lambda x: fun(x))] + else: return [NonlinearConstraint(fun, 0, np.inf, jac=jac)] + + def build_contraints(self): + constraints = [] + + constraints += self.all_linear() + constraints += self.all_nonlinear() + + return constraints + + +#Finite difference matrices + +def build_integration_matrix(N: int, dt: float) -> csc_matrix: + """Sparse integration matrix Θ (trapezoidal cumulative integration).""" + diagonal = np.ones(N) + diagonal[0] = 0.5 ; diagonal[-1] = 0.5 + Theta = diags(diagonal, shape=(N, N), format='csc') * dt + return Theta + +def build_first_derivative_matrix(N: int, dt: float) -> csc_matrix: + """Sparse first derivative matrix A₁ using backward difference.""" + diagonals = [-np.ones(N-1), np.ones(N-1)] + offsets = [0, 1] + A1 = diags(diagonals, offsets, shape=(N-1, N), format='csc') / dt + return A1 + +def build_second_derivative_matrix(N: int, dt: float) -> csc_matrix: + """Sparse second derivative matrix A₂ using central differences.""" + diagonals = [np.ones(N-1), -2*np.ones(N), np.ones(N-1)] + offsets = [-1, 0, 1] + A2 = diags(diagonals, offsets, shape=(N, N), format='csc') / (dt ** 2) + return A2.tocsc() + + +if __name__ == "__main__": + pass + + diff --git a/now_example.py b/now_example.py new file mode 100644 index 0000000..2fb2780 --- /dev/null +++ b/now_example.py @@ -0,0 +1,45 @@ +####### +# Exameple script for optimization of gradient waveforms +# Based the Matlab implementation https://github.com/jsjol/NOW +# Written by Jakob Eriksson and Oliver Näslund for Project course at Uppsala University +####### + +from optimizationProblem import NOW_config +from optimize import optimize +import numpy as np + +def main(): + N = 100 #Raster resolution + target = np.eye(3) #Target tensor, default isotropic + + # Define problem, similar to Matlab implementation + prob = NOW_config(N=N, + targetTensor=target, + durationFirstPartRequested = 28, + durationSecondPartRequested = 22, + durationZeroGradientRequested = 8, + gMax = 80, + sMax = 100, + useMaxNorm=False, + doMaxwellComp=True, + eta = 0.9, + motionCompensation={'order': [1,2], 'maxMagnitude': [0, 1e-4]}) + + # Insert problem formulation into optimization routine + opt = optimize(problem=prob) + + # Run the optimization routine + # early_stopping = True, sets a different termination condition which focuses on tolerance on b-value while perserving feasibility + opt.run(early_stopping=True) + + # Plot results + # animate = True, if subiterations wants to be visualized + opt.plot_results(animate=False) + + # Get arrays of interest + q, g, b, B = opt.get_results() + + +if __name__ == "__main__": + main() + \ No newline at end of file diff --git a/optimizationProblem.py b/optimizationProblem.py new file mode 100644 index 0000000..0b035ee --- /dev/null +++ b/optimizationProblem.py @@ -0,0 +1,130 @@ +import numpy as np + +def getActualTimings(durationFirstPartRequested, + durationZeroGradientRequested, + durationSecondPartRequested, + discretizationSteps, + forceSymmetry): + if forceSymmetry: + raise NotImplementedError + else: + totalTime = durationFirstPartRequested + durationSecondPartRequested + durationZeroGradientRequested + dt = totalTime/discretizationSteps + + startZeroGradientsIndex = np.floor(durationFirstPartRequested/dt) - 1 + startSecondPartIndex = np.floor((durationFirstPartRequested+durationZeroGradientRequested)/dt) - 1 + + durationFirstPartActual = (startZeroGradientsIndex + 1) * dt + durationSecondPartActual = (discretizationSteps - startSecondPartIndex - 1) * dt + durationZeroGradientActual = (startSecondPartIndex - startZeroGradientsIndex) * dt + totalTimeActual = durationFirstPartActual + durationSecondPartActual + durationZeroGradientActual + + if durationZeroGradientActual > 0: + zeroGradientAtIndex = np.arange(startZeroGradientsIndex, startSecondPartIndex + 1, dtype=int) + else: + zeroGradientAtIndex = None + + return (durationFirstPartActual, + durationZeroGradientActual, + durationSecondPartActual, + totalTimeActual, + zeroGradientAtIndex) + +def obj_fun(x): + f = -x[-1] #minimize negative of b + g = np.zeros_like(x) + g[-1] = -1.0 + return f, g + + +class NOW_config: + def __init__( + self, + targetTensor = np.eye(3), # Isotropic encoding tensor, + N = 50, #77 + useMaxNorm = False, + gMax = 80, #OBS OBS This is in mT/m + sMax = 100, #OBS OBS This is in T/m/s + durationFirstPartRequested = 28, + durationSecondPartRequested = 22, + durationZeroGradientRequested = 8, + eta = 1, #1 + enforceSymmetry = False, + name = 'NOW', + x0 = None, + doMaxwellComp = True, + MaxwellIndex = 100, + MaxFunEval = 1e5, + MaxIter = 5e3, + motionCompensation = {'order': [], 'maxMagnitude': [], 'linear': []}, + doBackgroundCompensation = 0, # 0 = off; 1 = general timing cond.; 2 = specific timing cond. + startTime = 0): # Time from the excitation (t=0) to the first gradient waveform sample in ms. + + # Parse kwargs + self.targetTensor = targetTensor #; self.targetTensor[-1,-1] = 0 ; self.targetTensor[1,1] = 0 + self.N = N + self.useMaxNorm = useMaxNorm + self.gMax = gMax + self.sMax = sMax + self.durationFirstPartRequested = durationFirstPartRequested + self.durationSecondPartRequested = durationSecondPartRequested + self.durationZeroGradientRequested = durationZeroGradientRequested + self.eta = eta + self.enforceSymmetry = enforceSymmetry + self.name = name + self.doMaxwellComp = doMaxwellComp + self.MaxwellIndex = MaxwellIndex + self.MaxFunEval = MaxFunEval + self.MaxIter = MaxIter + self.motionCompensation = motionCompensation + self.doBackgroundCompensation = doBackgroundCompensation + self.startTime = startTime + + + if x0 is None: + x0 = np.random.randn(3 * self.N + 1) + self.x0 = x0 + + # Get actual times after discretization + self.durationFirstPartActual, self.durationZeroGradientActual, self.durationSecondPartActual, self.totalTimeActual, self.zeroGradientAtIndex = \ + getActualTimings(self.durationFirstPartRequested, self.durationZeroGradientRequested, self.durationSecondPartRequested, self.N, self.enforceSymmetry) + + self._dt = self.totalTimeActual/self.N # Time step in milliseconds. Division by N instead of N-1 due to half step shift in gradients. + self._gMaxConstraint = self.gMax * self._dt + self._sMaxConstraint = self.sMax * self._dt**2 + self._integralConstraint = self.eta * self._gMaxConstraint**2 * self.totalTimeActual / self._dt + self._tolIsotropy = .5e-2 + + if self.doMaxwellComp: + self._tolMaxwell = self.MaxwellIndex/self._dt + else: + self._tolMaxwell = np.inf + + # Create spin dephasing direction vector + if self.zeroGradientAtIndex is not None: + signs = np.ones((self.N - 1, 1)) # Ghost points excluded during opt + + # Assume that sign change happens in the middle of the pause + mi = np.median(self.zeroGradientAtIndex).astype(int) + signs[np.round(mi):] = -1 + if mi == round(mi): + signs[mi] = 0 + self.signs = signs + + if self.motionCompensation['order'] != []: + self.motionCompensation['linear'] = (np.array(self.motionCompensation['maxMagnitude']) <= 0) + else: + self.motionCompensation['linear'] = [] + + + if self.doBackgroundCompensation != 0: + raise NotImplementedError + + +if __name__ == "__main__": + prob = NOW_config(N=11, + motionCompensation={'order': [1, 2], 'maxMagnitude': [0, 1e-4], 'linear': []}) + print(prob.motionCompensation['linear']) + + + \ No newline at end of file diff --git a/optimize.py b/optimize.py new file mode 100644 index 0000000..3558156 --- /dev/null +++ b/optimize.py @@ -0,0 +1,320 @@ +from optimizationProblem import NOW_config, obj_fun +from scipy.optimize import minimize +from constraints import constraints +import numpy as np +import matplotlib.pyplot as plt +import time +import os, json, sys, platform +from datetime import datetime +from types import SimpleNamespace +os.environ["KMP_DUPLICATE_LIB_OK"] = "TRUE" + + +class optimize: + def __init__(self, + problem = NOW_config(), + ): + + self.res = None #holds the result after run() + self.q, self.g, self.b, self.B = None, None, None, None #results + + self.problem = problem + self.objective_function = obj_fun + + constrains = constraints(problem=self.problem) #Creates the class which specifies the contraints given the problem, specifies if method to compute jacobian is standard or autograd + self.constraints = constrains.build_contraints() #Builds the constraints in a format accepted by scipy.optimize.minimize (list of dicts or constraint objects) + self.Theta, self.A1, self.A2 = constrains.Theta, constrains.A1, constrains.A2 #Integration and derivative matrices + self.dt = self.problem._dt + + self.intermediate_results = [] #store intermediate results for plotting + + self.elapsed_time = None #Timing of the optimization process + self.method_used = None #The optimization method used + + def callback(self, xk): + self.intermediate_results.append(xk.copy()) + + def callback_trust(self, xk, optimize_result): #callback function for trust-constr method + self.intermediate_results.append(xk.copy()) + + def run(self, method='SLSQP', early_stopping = True): + t0 = time.perf_counter() + self.method_used = method + if method == 'SLSQP': + if early_stopping: + chosen_callback = self.callback_with_earlystop + else: + chosen_callback = self.callback + try: + self.res = minimize( + fun = self.objective_function, + x0 = self.problem.x0, + jac = True, #objfun provides gradient as tuple !! + method = 'SLSQP', + constraints = self.constraints, + callback=chosen_callback, + options={'maxiter': 2000, + 'disp': True}, + #'iprint': 3}, # according to ChatGPT: 0 No output at all (default) 1 Final result summary only (same as disp=True) 2 Iteration-by-iteration summary of the optimization progress 3 or higher Even more internal details per iteration (gradient, step size, constraint info, etc.) + tol=1e-6, + #callback=self.callback + ) + except StopIteration: + print("Optimization stopped early.") + # SLSQP does not return a result when StopIteration is thrown + # so we must construct the result manually + self.res = SimpleNamespace(x=self.intermediate_results[-1]) + elif method == 'trust-constr': + self.res = minimize( + fun = self.objective_function, + x0 = self.problem.x0, + jac = True, #objfun provides gradient as tuple !! + method = 'trust-constr', + constraints = self.constraints, + options={'maxiter': 2000, + 'disp': True, + 'verbose':3}, + callback=self.callback_trust, + + ) + self.elapsed_time = time.perf_counter() - t0 + + def get_results(self): + now_gamma = 42.576e6 # Hz/T, samma som now_gamma i MATLAB + gamma = now_gamma*2*np.pi + + q_raw = self.res.x[:-1].reshape(self.problem.N, 3, order="F") #What is used in the optimization, same as Matlab + q_phys = q_raw * gamma * 1e-6 #SI-units + g = (self.A1 / self.dt) @ q_raw #OBS OBS OBS ask why this is q_raw and not q_phys??? In matlab code, they use q_raw for some reason + slew = (self.A2 @ q_raw) / (self.dt**2) #OBS OBS OBS here aswell, q raw?????? + dt_s = self.dt * 1e-3 + B = dt_s * (q_phys.T @ q_phys) + b_raw = self.res.x[-1] #OBS only corresponds to the intern variable, but we actually want to calculate b ourselves + b = np.trace(B) * 1e-6 # s/mm^2, This is the exact b-value per definition, and also what is done in the matlab code + + self.q = q_phys + self.g = g + self.B = B + self.b = b + self.slew = slew + + return self.q, self.g, self.b, self.B + + def results_from_x(self, x): + now_gamma = 42.576e6 # Hz/T, samma som now_gamma i MATLAB + gamma = now_gamma*2*np.pi + + q_raw = x[:-1].reshape(self.problem.N, 3, order="F") #What is used in the optimization, same as Matlab + q_phys = q_raw * gamma * 1e-6 #SI-units + g = (self.A1 / self.dt) @ q_raw #OBS OBS OBS ask why this is q_raw and not q_phys??? In matlab code, they use q_raw for some reason + slew = (self.A2 @ q_raw) / (self.dt**2) #OBS OBS OBS here aswell, q raw?????? + dt_s = self.dt * 1e-3 + B = dt_s * (q_phys.T @ q_phys) + #b_raw = self.res.x[-1] #OBS only corresponds to the intern variable, but we actually want to calculate b ourselves + b = np.trace(B) * 1e-6 # s/mm^2, This is the exact b-value per definition, and also what is done in the matlab code + + return q_phys, g, b, B + + def plot_results(self, animate = False, dir=None, x =None): + if x is not None: + q, g, b, B = self.results_from_x(x) + self.res = SimpleNamespace(x=self.intermediate_results[-1],nit = len(self.intermediate_results)) + else: + q, g, b, B = self.get_results() + + t_g = np.linspace(0.5 * self.problem._dt, self.problem.totalTimeActual - 0.5 * self.problem._dt, endpoint=True, num=self.problem.N - 1) + t_q = np.linspace(0, self.problem.totalTimeActual, self.problem.N, endpoint = True) + + fig, axs = plt.subplots(2, 1, sharex = True, figsize = (8, 6)) + q_mikrometer = q/10**6 + #subplot for q + qx = axs[0].plot(t_q, q_mikrometer[:, 0], label='q_x') + qy = axs[0].plot(t_q, q_mikrometer[:, 1], label='q_y') + qz = axs[0].plot(t_q, q_mikrometer[:, 2], label='q_z') + axs[0].set_ylabel("q [µs/m]") + axs[0].set_title("Integrated gradient q(t)") + axs[0].legend() + axs[0].grid(True) + + #subplot for g + gx = axs[1].plot(t_g, g[:, 0], label='g_x') + gy = axs[1].plot(t_g, g[:, 1], label='g_y') + gz = axs[1].plot(t_g, g[:, 2], label='g_z') + axs[1].set_xlabel("Time [ms]") + axs[1].set_ylabel("g [mT/m]") + axs[1].set_title("Magnetic field gradient g(t)") + axs[1].legend() + axs[1].grid(True) + + #Information text + iter_count = float(len(self.intermediate_results)) #getattr(self.res, "nit", None) + elapsed_time = getattr(self, "elapsed_time", None) + target = self.problem.targetTensor # 3x3 + trace_B = np.trace(B) + trace_target = np.trace(target) + encoding_tensor = (B / trace_B) * trace_target + fro = np.linalg.norm(encoding_tensor - target) + + info_text = ( + f"Iterations: {iter_count}\n" + f"Time: {elapsed_time:.3f} s\n" + f"b: {b:.6f} s/mm^2\n" + f"Fro-norm = {fro:.6e}\n" + "Encoding tensor:\n" + + "\n".join([" ".join([f"{val: .6f}" for val in row]) for row in encoding_tensor]) + ) + + + if animate == False: + plt.tight_layout(rect=[0, 0.22, 1, 1]) #Leave room for the box + info_ax = fig.add_axes([0.08, 0.03, 0.84, 0.16]) #[left, bottom, width, height] (distances) + info_ax.axis('off') + info_ax.text(0.5, 0.5, info_text, + ha='center', va='center', fontsize=8, + bbox=dict(facecolor='white', alpha=0.85, boxstyle='round')) + if dir == None: + plt.show() + else: + plt.savefig(dir) + + elif animate == True: + from matplotlib.widgets import Slider, Button + + plt.tight_layout(rect=[0, 0.22, 1, 1]) #Leave room for the box + info_ax = fig.add_axes([0.08, 0.07, 0.84, 0.20]) #[left, bottom, width, height] (distances) + info_ax.axis('off') + info_ax.text(0.5, 0.5, info_text, + ha='center', va='center', fontsize=8, + bbox=dict(facecolor='white', alpha=0.85, boxstyle='round')) + + q_frames = [] ; g_frames = [] ; b_frames = [] ; B_frames = [] + intermediate_results = self.intermediate_results + for x in intermediate_results: + _q, _g, _b, _B = self.results_from_x(x) + q_frames.append(_q) + g_frames.append(_g) + b_frames.append(_b) + B_frames.append(_B) + + # Get line objects for easy updating + lines_q = axs[0].lines # [q_x, q_y, q_z] + lines_g = axs[1].lines # [g_x, g_y, g_z] + + n_frames = len(q_frames) + + # Make space at bottom for the slider + plt.subplots_adjust(bottom=0.35) + + # Slider axis: [left, bottom, width, height] in figure coordinates + ax_slider = plt.axes([0.15, 0.05, 0.7, 0.03]) + frame_slider = Slider( + ax=ax_slider, + label="Iteration", + valmin=0, + valmax=n_frames - 1, + valinit=0, + valstep=1, # integer steps = frame indices + ) + + def update_slider(val): + i = int(val) + + q = q_frames[i] + g = g_frames[i] + b = b_frames[i] + B = B_frames[i] + + # update q lines + for j in range(3): + lines_q[j].set_ydata(q[:, j]/10**6) #mikrometer conversion + + # update g lines + for j in range(3): + lines_g[j].set_ydata(g[:, j]) + + #Information text + elapsed_time = getattr(self, "elapsed_time", None) + target = self.problem.targetTensor # 3x3 + trace_B = np.trace(B) + trace_target = np.trace(target) + encoding_tensor = (B / trace_B) * trace_target + fro = np.linalg.norm(encoding_tensor - target) + + info_text = ( + f"Iteration: {i}\n" + f"Time: {elapsed_time:.3f} s\n" + f"b: {b:.6f} s/mm^2\n" + f"Fro-norm = {fro:.6e}\n" + "Encoding tensor:\n" + + "\n".join([" ".join([f"{val: .6f}" for val in row]) for row in encoding_tensor]) + ) + info_ax.text(0.5, 0.5, info_text, + ha='center', va='center', fontsize=8, + bbox=dict(facecolor='white', alpha=0.85, boxstyle='round')) + + + fig.canvas.draw_idle() # redraw efficiently + + # Connect the slider to the update function + frame_slider.on_changed(update_slider) + + plt.show() + + def callback_with_earlystop(self, xk): + + # Save xk for animation + self.callback(xk) + + # Step size + if len(self.intermediate_results) > 1: + step_norm = np.linalg.norm((xk - self.intermediate_results[-2])[:-1]) #OBS OBS OBS IGNORERA b-värdet!!! + else: + step_norm = np.inf + + # Objective + f, grad = self.objective_function(xk) + + # Constraint violation + violations = [] + for c in self.constraints: + if isinstance(c, dict): + val = c["fun"](xk) + violations.append(np.min(val)) + constraint_violation = max(0, -min(violations)) if violations else 0 + + # External stagnation + if not hasattr(self, "prev_f_ext"): + self.prev_f_ext = f + self.no_improve_ext = 0 + else: + if abs(f - self.prev_f_ext) < 1e-4: # <-- BIGGER tolerance! + self.no_improve_ext += 1 + else: + self.no_improve_ext = 0 + + self.prev_f_ext = f + + # Early stop for EXTERNAL iterations + #print( + # f"[ITER {len(self.intermediate_results)-1:4d}] " + # f"Δf = {abs(f - self.prev_f_ext):.3e} " + # f"step = {step_norm:.3e} " + # f"viol = {constraint_violation:.3e} " + # f"no_improve = {self.no_improve_ext}" + #) + if ( + constraint_violation < 1e-6 and # constraints ok + step_norm < 1e-2 and # small enough for external loop + self.no_improve_ext >= 20 # no improvement for 10 external iterations + ): + print("EARLY STOPPING (EXTERNAL): Stagnation detected.") + raise StopIteration + + + +if __name__ == "__main__": + pass + + + + From ac5931a48bd010ea0373a744068239e1030206cf Mon Sep 17 00:00:00 2001 From: Olivernaslund <133633664+Olivernaslund@users.noreply.github.com> Date: Wed, 4 Mar 2026 17:14:14 +0100 Subject: [PATCH 2/7] Create temp.txt --- now_python/temp.txt | 1 + 1 file changed, 1 insertion(+) create mode 100644 now_python/temp.txt diff --git a/now_python/temp.txt b/now_python/temp.txt new file mode 100644 index 0000000..8b13789 --- /dev/null +++ b/now_python/temp.txt @@ -0,0 +1 @@ + From 5b530b105ccd14f5dc82382c18008020d56415a2 Mon Sep 17 00:00:00 2001 From: Olivernaslund <133633664+Olivernaslund@users.noreply.github.com> Date: Wed, 4 Mar 2026 17:15:12 +0100 Subject: [PATCH 3/7] Moved constraints.py --- constraints.py => now_python/constraints.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename constraints.py => now_python/constraints.py (100%) diff --git a/constraints.py b/now_python/constraints.py similarity index 100% rename from constraints.py rename to now_python/constraints.py From d91bb92a31724d83aafe3bc286dfc81d76f687fd Mon Sep 17 00:00:00 2001 From: Olivernaslund <133633664+Olivernaslund@users.noreply.github.com> Date: Wed, 4 Mar 2026 17:15:39 +0100 Subject: [PATCH 4/7] moved optimizationProblem.py --- optimizationProblem.py => now_python/optimizationProblem.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) rename optimizationProblem.py => now_python/optimizationProblem.py (99%) diff --git a/optimizationProblem.py b/now_python/optimizationProblem.py similarity index 99% rename from optimizationProblem.py rename to now_python/optimizationProblem.py index 0b035ee..2f9996f 100644 --- a/optimizationProblem.py +++ b/now_python/optimizationProblem.py @@ -127,4 +127,4 @@ def __init__( print(prob.motionCompensation['linear']) - \ No newline at end of file + From a880105ad6fe3d6470d3568c3147434b8a17e1db Mon Sep 17 00:00:00 2001 From: Olivernaslund <133633664+Olivernaslund@users.noreply.github.com> Date: Wed, 4 Mar 2026 17:16:12 +0100 Subject: [PATCH 5/7] Moved optimize.py --- optimize.py => now_python/optimize.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename optimize.py => now_python/optimize.py (100%) diff --git a/optimize.py b/now_python/optimize.py similarity index 100% rename from optimize.py rename to now_python/optimize.py From 1135b4a08631b5d1842d5d0fc3fac144bab4422f Mon Sep 17 00:00:00 2001 From: Olivernaslund <133633664+Olivernaslund@users.noreply.github.com> Date: Wed, 4 Mar 2026 17:16:41 +0100 Subject: [PATCH 6/7] Moved now_example.py --- now_example.py => now_python/now_example.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) rename now_example.py => now_python/now_example.py (99%) diff --git a/now_example.py b/now_python/now_example.py similarity index 99% rename from now_example.py rename to now_python/now_example.py index 2fb2780..9311e19 100644 --- a/now_example.py +++ b/now_python/now_example.py @@ -42,4 +42,4 @@ def main(): if __name__ == "__main__": main() - \ No newline at end of file + From 82bd6d0b169363072e4f31b326850c6712b157ca Mon Sep 17 00:00:00 2001 From: Olivernaslund <133633664+Olivernaslund@users.noreply.github.com> Date: Wed, 4 Mar 2026 17:17:13 +0100 Subject: [PATCH 7/7] Delete now_python/temp.txt --- now_python/temp.txt | 1 - 1 file changed, 1 deletion(-) delete mode 100644 now_python/temp.txt diff --git a/now_python/temp.txt b/now_python/temp.txt deleted file mode 100644 index 8b13789..0000000 --- a/now_python/temp.txt +++ /dev/null @@ -1 +0,0 @@ -