From 8144f1b5497b9ad00638dd3009348e0869a10bcb Mon Sep 17 00:00:00 2001 From: William Wang Date: Sun, 20 Oct 2024 16:15:26 -0700 Subject: [PATCH] FUll bending forces: fixed pi bond generation and bending bugs --- energyminimization/energies/bend_full.py | 55 +++++--- energyminimization/energies/stretch_full.py | 59 ++++---- energyminimization/matrix_helper.py | 32 +++-- energyminimization/minimize.py | 19 +-- energyminimization/solvers/newton.py | 2 +- energyminimization/solvers/solver.py | 26 +++- energyminimization/transformations.py | 5 + lattice/abstract_lattice.py | 146 ++++++++------------ lattice/lattice_factory.py | 4 +- lattice/triangular_lattice.py | 5 +- main.py | 1 - parameters.py | 4 +- tests/matrix_tests.py | 33 +++-- 13 files changed, 206 insertions(+), 185 deletions(-) diff --git a/energyminimization/energies/bend_full.py b/energyminimization/energies/bend_full.py index 9049eae..c86732f 100644 --- a/energyminimization/energies/bend_full.py +++ b/energyminimization/energies/bend_full.py @@ -22,16 +22,23 @@ def compute_angles( pos_matrix = pos_matrix.reshape(-1, 2) j, i, k = active_pi_indices[:, 0], active_pi_indices[:, 1], active_pi_indices[:, 2] bond1_idx, bond2_idx, pi_idx = active_pi_indices[:, 3], active_pi_indices[:, 4], active_pi_indices[:, 5] - hor_pbc1, top_pbc1 = active_pi_indices[:, 7], active_pi_indices[:, 8] - hor_pbc2, top_pbc2 = active_pi_indices[:, 9], active_pi_indices[:, 10] + hor_pbc1, top_pbc1 = active_pi_indices[:, 8], active_pi_indices[:, 9] + hor_pbc2, top_pbc2 = active_pi_indices[:, 10], active_pi_indices[:, 11] + sign_i, sign_j = active_pi_indices[:, 6], active_pi_indices[:, 7] + # Vectors out of the vertex node - c_ji = pos_matrix[j] - pos_matrix[i] - c_jk = pos_matrix[j] - pos_matrix[k] + c_ji = pos_matrix[i] - pos_matrix[j] + c_jk = pos_matrix[k] - pos_matrix[j] # Handle periodic boundary conditions - c_ji[np.where(hor_pbc1 == 1), 0] += corrections[0] - c_ji[np.where(top_pbc1 == 1), 1] += corrections[1] - c_jk[np.where(hor_pbc2 == 1), 0] += corrections[0] - c_jk[np.where(top_pbc2 == 1), 1] += corrections[1] + idx_hor_i = np.where(hor_pbc1 == 1) + idx_top_i = np.where(top_pbc1 == 1) + idx_hor_j = np.where(hor_pbc2 == 1) + idx_top_j = np.where(top_pbc2 == 1) + # The PBC correction depends on which node is the vertex (determined by sign) + c_ji[idx_hor_i, 0] += sign_i[idx_hor_i] * corrections[0] + c_ji[idx_top_i, 1] += sign_i[idx_top_i] * corrections[1] + c_jk[idx_hor_j, 0] += sign_j[idx_hor_j] * corrections[0] + c_jk[idx_top_j, 1] += sign_j[idx_top_j] * corrections[1] theta_jik = np.arctan2(np.cross(c_ji, c_jk), np.einsum('ij,ij->i', c_ji, c_jk)) theta_jik[theta_jik < 0] += 2 * np.pi @@ -63,8 +70,8 @@ def get_bend_hessian( # We first compute the helper values needed for d_theta/d_dof length_ji = np.linalg.norm(c_ji, axis=1) length_jk = np.linalg.norm(c_jk, axis=1) - t_ji = np.divide(c_ji, np.square(length_ji)) - t_jk = np.divide(c_jk, np.square(length_jk)) + t_ji = np.divide(c_ji, np.square(length_ji)[: , np.newaxis]) + t_jk = np.divide(c_jk, np.square(length_jk)[: , np.newaxis]) # Compute d_theta/d_dof d_theta_djx = t_ji[:, 1] - t_jk[:, 1] d_theta_dix = -t_ji[:, 1] @@ -97,7 +104,7 @@ def get_bend_hessian( cjk_x_minus_y = c_jk[:, 0] - c_jk[:, 1] # We're ready! d2_theta_djx_jx = 2 * (cjk_prod / l_jk_fourth - cji_prod / l_ji_fourth) - d2_theta_djx_jy = (1 / l_ji_sq) - (1 / l_jk_sq) + (2 * cjk_sq / l_jk_fourth) - (2 * cji_sq / l_ji_fourth) + d2_theta_djx_jy = (1 / l_ji_sq) - (1 / l_jk_sq) + (2 * cjk_sq[:, 1] / l_jk_fourth) - (2 * cji_sq[:, 1] / l_ji_fourth) d2_theta_djx_ix = 2 * (cji_prod / l_ji_fourth) d2_theta_djx_iy = -(cji_x_plus_y * cji_x_minus_y) / l_ji_fourth d2_theta_djx_kx = -2 * (cjk_prod / l_jk_fourth) @@ -212,7 +219,8 @@ def get_bend_hessian( add_entry(rows, cols, data, entry_idx + 34, k_y, k_x, d2e_dky_kx) add_entry(rows, cols, data, entry_idx + 35, k_y, k_y, d2e_dky_ky) - return csr_matrix((data, (rows, cols)), shape=(2 * n, 2 * n)) + n_dof = pos_matrix.size + return csr_matrix((data, (rows, cols)), shape=(n_dof, n_dof)) def get_bend_jacobian( @@ -225,21 +233,24 @@ def get_bend_jacobian( gradient = np.zeros(pos_matrix.size) theta, diff_theta, c_ji, c_jk = compute_angles(pos_matrix=pos_matrix, corrections=corrections, active_pi_indices=active_pi_indices, orig_pi_angles=orig_pi_angles) - # Gradient is equal to -(kappa) * (theta - theta_0) * d_theta/d_dof - grad_factor = -bend_mod * diff_theta + # Gradient is equal to kappa * (theta - theta_0) * d_theta/d_dof + grad_factor = (bend_mod * diff_theta)[: , np.newaxis] # Temporary variables useful for the gradient length_ji = np.linalg.norm(c_ji, axis=1) length_jk = np.linalg.norm(c_jk, axis=1) - t_ji = grad_factor * np.divide(c_ji, np.square(length_ji)) - t_jk = grad_factor * np.divide(c_jk, np.square(length_jk)) + t_ji = np.divide(c_ji, np.square(length_ji)[:, np.newaxis]) + t_jk = np.divide(c_jk, np.square(length_jk)[: , np.newaxis]) + # Multiply be the gradient factor so we don't have to later (save some redundant computation) + g_ji = grad_factor * t_ji + g_jk = grad_factor * t_jk # Partials of theta with respect to x (times the gradient factor gives the correct partials of energy) - de_djx = t_ji[:, 1] - t_jk[:, 1] - de_dix = -t_ji[:, 1] - de_dkx = t_jk[:, 1] + de_djx = -g_ji[:, 1] + g_jk[:, 1] + de_dix = g_ji[:, 1] + de_dkx = -g_jk[:, 1] # Partials of theta with respect to y - de_djy = -t_ji[:, 0] + t_jk[:, 0] - de_diy = t_ji[:, 0] - de_dky = -t_jk[:, 0] + de_djy = g_ji[:, 0] - g_jk[:, 0] + de_diy = -g_ji[:, 0] + de_dky = g_jk[:, 0] # Update the gradient j, i, k = active_pi_indices[:, 0], active_pi_indices[:, 1], active_pi_indices[:, 2] diff --git a/energyminimization/energies/stretch_full.py b/energyminimization/energies/stretch_full.py index c29e005..1f7ac95 100644 --- a/energyminimization/energies/stretch_full.py +++ b/energyminimization/energies/stretch_full.py @@ -16,23 +16,25 @@ def compute_lengths( """ Computes the vector between each node (c_ij), the length of each bond, and the difference between the current length and the rest length - :return: c_ji (contains vectors from node j to i for each bond), length_ji (length of each bond), - d_ji (difference between current length and rest length) + :return: + - c_ij: vector between nodes j and i for each bond + - length_ij: lengths of c_ij + - d_ij: difference between current length and rest length """ pos_matrix = pos_matrix.reshape(-1, 2) # Indices of the nodes for each bond i, j = active_bond_indices[:, 0], active_bond_indices[:, 1] - # c_ij points from j to i based on their current positions - c_ji = pos_matrix[j] - pos_matrix[i] + # c_ij points from i to j based on their current positions + c_ij = pos_matrix[j] - pos_matrix[i] # Handle periodic boundary conditions hor_pbc, top_pbc = active_bond_indices[:, 2], active_bond_indices[:, 3] - c_ji[np.where(hor_pbc == 1), 0] += corrections[0] - c_ji[np.where(top_pbc == 1), 1] += corrections[1] + c_ij[np.where(hor_pbc == 1), 0] += corrections[0] + c_ij[np.where(top_pbc == 1), 1] += corrections[1] # Lengths and differences idx = active_bond_indices[:, -1] - length_ji = np.linalg.norm(c_ji, axis=1) - d_ji = length_ji - active_bond_lengths[idx] - return c_ji, length_ji, d_ji + length_ij = np.linalg.norm(c_ij, axis=1) + d_ij = length_ij - active_bond_lengths[idx] + return c_ij, length_ij, d_ij def get_nonlinear_stretch_hessian( @@ -43,21 +45,21 @@ def get_nonlinear_stretch_hessian( active_bond_lengths: np.ndarray ) -> csr_matrix: # Compute the current vectors associated with the position - c_ji, l_ji, d_matrix = compute_lengths(pos_matrix=pos_matrix, active_bond_indices=active_bond_indices, + c_ij, l_ij, d_ij = compute_lengths(pos_matrix=pos_matrix, active_bond_indices=active_bond_indices, corrections=corrections, active_bond_lengths=active_bond_lengths) - i, j, idx = active_bond_indices[:, 0], active_bond_indices[:, 1], active_bond_indices[:, -1] - # Second derivative of energy with respect to same variable - c_ij_sq = np.square(c_ji) - init_bond_lengths = active_bond_lengths[idx].reshape(-1, 1) - curr_bond_lengths = l_ji.reshape(-1, 1) - # Formula is a * (1 - ((l0 * (c2_ij[!i]) / l^3)), where !i is y if i is x and vice versa - dd = stretch_mod * (1 - np.divide(np.multiply(init_bond_lengths, c_ij_sq), np.power(curr_bond_lengths, 3))) + idx = active_bond_indices[:, -1] + l0 = active_bond_lengths[idx] + # Helpful values for computing second derivatives + c_ij_sq = np.square(c_ij) + c_ij_prod = np.prod(c_ij, axis=1) # c_ij[0] * c_ij[1] + l0_div_l3 = np.divide(l0, np.power(l_ij, 3)) # l0 / l^3 + + # Formula for d^2E/d(x,y)^2 and alpha * (1 - ((l0 * (c2_ij[!i]) / l^3)), where !i is y if i is x and vice versa + dd = stretch_mod * (1 - np.multiply(l0_div_l3[:, np.newaxis], c_ij_sq)) # need new axis since l0_div_l3 is 1D dxx, dyy = dd[:, 1], dd[:, 0] - # Formula is a * (l0 * (c_ij[0] * c_ij[1]) / l^3) - c_ij_prod = np.prod(c_ji, axis=1).flatten() - dxy = stretch_mod * np.divide(np.multiply(init_bond_lengths.flatten(), c_ij_prod), - np.power(curr_bond_lengths.flatten(), 3)) + # Formula for d^2E/dxdy is alpha * (l0 * (c_ij[0] * c_ij[1]) / l^3) + dxy = stretch_mod * np.multiply(l0_div_l3, c_ij_prod) n_dof = pos_matrix.size return sh.generate_hessian(active_bond_indices=active_bond_indices, dxx=dxx, dyy=dyy, dxy=dxy, n=n_dof) @@ -70,14 +72,14 @@ def get_nonlinear_stretch_jacobian( active_bond_indices: np.ndarray, active_bond_lengths: np.ndarray ) -> np.ndarray: - c_ji, l_ji, d_matrix = compute_lengths(pos_matrix=pos_matrix, active_bond_indices=active_bond_indices, + c_ij, l_ij, d_ij = compute_lengths(pos_matrix=pos_matrix, active_bond_indices=active_bond_indices, corrections=corrections, active_bond_lengths=active_bond_lengths) - # Gradient is equal to -alpha * ((l - l_0) / l) * c_ij. Compute the pre-factor - grad_factor = -stretch_mod * np.divide(d_matrix, l_ji) + # Gradient is equal to -alpha * ((l - l_0) / l) * c_ij. First compute the pre-factor + grad_factor = -stretch_mod * np.divide(d_ij, l_ij)[:, np.newaxis] # Then multiply the pre-factor by the c_ij vectors - grad = np.multiply(grad_factor.reshape(-1, 1), c_ji) - # The derivative with respect to x has only the x component. Same for y + grad = np.multiply(grad_factor, c_ij) + # The derivative with respect to x and y grad_x, grad_y = grad[:, 0], grad[:, 1] n_dof = pos_matrix.size @@ -100,7 +102,8 @@ def get_full_stretch_energy( :param active_bond_lengths: ORIGINAL lengths of active bonds """ - _, _, d_ji = compute_lengths(pos_matrix=pos_matrix, active_bond_indices=active_bond_indices, + _, _, d_ij = compute_lengths(pos_matrix=pos_matrix, active_bond_indices=active_bond_indices, corrections=corrections, active_bond_lengths=active_bond_lengths) - energy = 0.5 * stretch_mod * np.square(d_ji) + # E = 1/2 * alpha * (l - l_0)^2 + energy = 0.5 * stretch_mod * np.square(d_ij) return float(np.sum(energy)) diff --git a/energyminimization/matrix_helper.py b/energyminimization/matrix_helper.py index a8f2ccf..40445fb 100644 --- a/energyminimization/matrix_helper.py +++ b/energyminimization/matrix_helper.py @@ -124,26 +124,34 @@ def verify_r_matrix( def create_angle_matrix( pos_vector: np.ndarray, active_pi_indices: np.ndarray, - corrections: np.ndarray, + lattice: AbstractLattice, ) -> np.ndarray: """ - Matrix size (# pi bonds, 1) where value at index i is the angle for pi bond index i + Matrix size (# pi bonds, 1) containing the *initial* angles for each pi bond """ pos_matrix = pos_vector.reshape(-1, 2) j, i, k = active_pi_indices[:, 0], active_pi_indices[:, 1], active_pi_indices[:, 2] - hor_pbc1, top_pbc1 = active_pi_indices[:, 7], active_pi_indices[:, 8] - hor_pbc2, top_pbc2 = active_pi_indices[:, 9], active_pi_indices[:, 10] + hor_pbc1, top_pbc1 = active_pi_indices[:, 8], active_pi_indices[:, 9] + hor_pbc2, top_pbc2 = active_pi_indices[:, 10], active_pi_indices[:, 11] + sign_i, sign_j = active_pi_indices[:, 6], active_pi_indices[:, 7] - # Vectors from vertex node j to the edge nodes i and k - c_ji = pos_matrix[j] - pos_matrix[i] - c_jk = pos_matrix[j] - pos_matrix[k] + # Vectors from vertex node j to edge nodes i and k + c_ji = pos_matrix[i] - pos_matrix[j] + c_jk = pos_matrix[k] - pos_matrix[j] # Handle periodic boundary conditions - c_ji[np.where(hor_pbc1 == 1), 0] += corrections[0] - c_ji[np.where(top_pbc1 == 1), 1] += corrections[1] - c_jk[np.where(hor_pbc2 == 1), 0] += corrections[0] - c_jk[np.where(top_pbc2 == 1), 1] += corrections[1] - + corrections = np.array([lattice.get_length(), lattice.get_height() + lattice.height_increment]) + idx_hor_i = np.where(hor_pbc1 == 1) + idx_top_i = np.where(top_pbc1 == 1) + idx_hor_j = np.where(hor_pbc2 == 1) + idx_top_j = np.where(top_pbc2 == 1) + # The PBC correction depends on which node is the vertex (determined by sign) + c_ji[idx_hor_i, 0] += sign_i[idx_hor_i] * corrections[0] + c_ji[idx_top_i, 1] += sign_i[idx_top_i] * corrections[1] + c_jk[idx_hor_j, 0] += sign_j[idx_hor_j] * corrections[0] + c_jk[idx_top_j, 1] += sign_j[idx_top_j] * corrections[1] + + # Compute the angle between the two vectors angle_matrix = np.arctan2(np.cross(c_ji, c_jk), np.einsum('ij,ij->i', c_ji, c_jk)) angle_matrix[angle_matrix < 0] += 2 * np.pi return angle_matrix diff --git a/energyminimization/minimize.py b/energyminimization/minimize.py index cf74174..87d34f5 100644 --- a/energyminimization/minimize.py +++ b/energyminimization/minimize.py @@ -105,21 +105,22 @@ def minimize( # [vertex id, edge node 1 id, edge node 2 id, bond 1 id, bond 2 id, pi bond id, sign for r_matrix[bond id], horizontal PBC, top PBC] active_pi_bonds = lattice.get_active_pi_bonds() - active_pi_indices = np.zeros((len(active_pi_bonds), 11), dtype=np.int32) + active_pi_indices = np.zeros((len(active_pi_bonds), 12), dtype=np.int32) for i, pi_bond in enumerate(active_pi_bonds): assert (pi_bond.exists()) - active_pi_indices[i][0] = pi_bond.get_vertex_node().get_id() - active_pi_indices[i][1] = pi_bond.get_edge_nodes()[0].get_id() - active_pi_indices[i][2] = pi_bond.get_edge_nodes()[1].get_id() + active_pi_indices[i][0] = pi_bond.get_vertex_node().get_id() # j + active_pi_indices[i][1] = pi_bond.get_edge_nodes()[0].get_id() # i + active_pi_indices[i][2] = pi_bond.get_edge_nodes()[1].get_id() # k active_pi_indices[i][3] = bond_to_idx[pi_bond.get_bond1()] active_pi_indices[i][4] = bond_to_idx[pi_bond.get_bond1()] active_pi_indices[i][5] = i - active_pi_indices[i][6] = 1 if pi_bond.get_bond1().get_node1() == pi_bond.get_vertex_node() else -1 + active_pi_indices[i][6] = 1 if pi_bond.get_bond1().get_node1() == pi_bond.get_vertex_node() else -1 # sign 1 + active_pi_indices[i][7] = 1 if pi_bond.get_bond2().get_node1() == pi_bond.get_vertex_node() else -1 # sign 2 # PBC - active_pi_indices[i][7] = pi_bond.get_bond1().is_hor_pbc() - active_pi_indices[i][8] = pi_bond.get_bond1().is_top_pbc() - active_pi_indices[i][9] = pi_bond.get_bond2().is_hor_pbc() - active_pi_indices[i][10] = pi_bond.get_bond2().is_top_pbc() + active_pi_indices[i][8] = pi_bond.get_bond1().is_hor_pbc() + active_pi_indices[i][9] = pi_bond.get_bond1().is_top_pbc() + active_pi_indices[i][10] = pi_bond.get_bond2().is_hor_pbc() + active_pi_indices[i][11] = pi_bond.get_bond2().is_top_pbc() # 2. Convert lattice positions to matrices # Initial position of nodes, initial unit vectors for bonds, initial energy diff --git a/energyminimization/solvers/newton.py b/energyminimization/solvers/newton.py index a8c6b82..92ed705 100644 --- a/energyminimization/solvers/newton.py +++ b/energyminimization/solvers/newton.py @@ -162,7 +162,7 @@ def tr_solve(f0, x, g, A, max_iter, tol, trust_radius): def trust_region_newton_cg(fun, x0, jac, hess, g_tol=1e-8): # --- Initialize trust-region --- - trust_radius = 1.0 + trust_radius = 10.0 max_trust_radius = 1e3 eta = 0.15 diff --git a/energyminimization/solvers/solver.py b/energyminimization/solvers/solver.py index 7ed31f7..435074c 100644 --- a/energyminimization/solvers/solver.py +++ b/energyminimization/solvers/solver.py @@ -5,6 +5,7 @@ from scipy.sparse import spmatrix, csr_matrix import energyminimization.energies.stretch_full as snl +import energyminimization.energies.bend_full as bnl import energyminimization.matrix_helper as pos from energyminimization.matrix_helper import KMatrixResult from energyminimization.solvers.conjugate_gradient import conjugate_gradient @@ -12,6 +13,7 @@ from energyminimization.solvers.newton import trust_region_newton_cg from energyminimization.transformations import Strain from lattice.abstract_lattice import AbstractLattice +from tests.matrix_tests import test_gradient_hessian, plot_sparsity @dataclass @@ -162,7 +164,11 @@ def compute_total_energy(pos_matrix: np.ndarray) -> float: corrections=correction, active_bond_indices=params.active_bond_indices, active_bond_lengths=params.length_matrix) - bend_energy = 0 + bend_energy = bnl.get_bend_energy(bend_mod=params.bend_mod, + pos_matrix=pos_matrix, + corrections=correction, + active_pi_indices=params.active_pi_indices, + orig_pi_angles=params.angle_matrix) return stretch_energy + bend_energy def compute_total_gradient(pos_matrix: np.ndarray) -> np.ndarray: @@ -172,16 +178,28 @@ def compute_total_gradient(pos_matrix: np.ndarray) -> np.ndarray: corrections=correction, active_bond_indices=params.active_bond_indices, active_bond_lengths=params.length_matrix) - bend_grad = np.zeros_like(pos_matrix) - return stretch_grad.ravel() + bend_grad = bnl.get_bend_jacobian( + bend_mod=params.bend_mod, + pos_matrix=pos_matrix, + corrections=correction, + active_pi_indices=params.active_pi_indices, + orig_pi_angles=params.angle_matrix) + return stretch_grad.ravel() + bend_grad.ravel() def compute_total_hessian(pos_matrix: np.ndarray) -> csr_matrix: - return snl.get_nonlinear_stretch_hessian( + stretch_hess = snl.get_nonlinear_stretch_hessian( stretch_mod=params.stretch_mod, pos_matrix=pos_matrix, corrections=correction, active_bond_indices=params.active_bond_indices, active_bond_lengths=params.length_matrix) + bend_hess = bnl.get_bend_hessian( + bend_mod=params.bend_mod, + pos_matrix=pos_matrix, + corrections=correction, + active_pi_indices=params.active_pi_indices, + orig_pi_angles=params.angle_matrix) + return stretch_hess + bend_hess init_energy = compute_total_energy(pos_matrix=params.strained_pos) diff --git a/energyminimization/transformations.py b/energyminimization/transformations.py index c9d4fbc..d49d96d 100644 --- a/energyminimization/transformations.py +++ b/energyminimization/transformations.py @@ -42,6 +42,11 @@ def __init__(self, gamma: float): super().__init__(gamma, transformation) self.name = "shear" +class Random(Strain): + def __init__(self, gamma: float): + transformation = np.array([[1 + gamma * np.random.uniform(), gamma * np.random.uniform()], [gamma * np.random.uniform(), 1 + gamma * np.random.uniform()]]) + super().__init__(gamma, transformation) + self.name = "random" class Dilate(Strain): def __init__(self, gamma: float): diff --git a/lattice/abstract_lattice.py b/lattice/abstract_lattice.py index e7bf2f5..0f769f7 100644 --- a/lattice/abstract_lattice.py +++ b/lattice/abstract_lattice.py @@ -93,7 +93,7 @@ def update_active_bonds(self) -> None: def remove_bond(self, bond: Bond) -> None: """ - Removes a bond from the lattice + Deactivate a bond from the lattice (does not remove the object, just sets it inactive) """ bond.set_inactive() self.active_bonds.remove(bond) @@ -342,56 +342,65 @@ def load_lattice(self, node_pos_data, node_data, bond_node_data, bond_data, pi_b self.update_active_bonds() return + def check_bonds_colinear(self, bond_i: Bond, bond_j: Bond, vertex: Node): + """ + Checks if two bonds are co-linear + """ + # Find the unit vector of each bond + pos_i1, pos_i2 = bond_i.get_xy_xy() + r_i = np.subtract(pos_i2, pos_i1) + if bond_i.is_hor_pbc(): + r_i[0] = r_i[0] + self.length + if bond_i.is_top_pbc(): + r_i[1] = r_i[1] + self.height + self.height_increment + r_i = r_i / np.linalg.norm(r_i, axis=0) + + # Vector from vertex node to edge (if sign_i is 1). Flipped if sign_i is -1 + sign_i = 1 if bond_i.get_node1() == vertex else -1 + r_i = sign_i * r_i + + pos_j1, pos_j2 = bond_j.get_xy_xy() + r_j = np.subtract(pos_j2, pos_j1) + if bond_j.is_hor_pbc(): + r_j[0] = r_j[0] + self.length + if bond_j.is_top_pbc(): + r_j[1] = r_j[1] + self.height + self.height_increment + r_j = r_j / np.linalg.norm(r_j, axis=0) + + sign_j = 1 if bond_j.get_node1() == vertex else -1 + r_j = sign_j * r_j + + angle = np.arctan2(np.cross(r_i, r_j), np.dot(r_i, r_j)) + if angle < 0: + angle = angle + 2 * np.pi + return abs(abs(angle) - np.pi) < 0.01 + def generate_pi_bonds(self) -> None: """ - Generates the pi bonds in the lattice - (series of two co-linear bonds with a node vertex) + Generates the pi bonds in the lattice (series of two co-linear bonds with a node vertex) """ self.pi_bonds = [] - # If dictionary doesn't exist, use a brute method + # If the pre-made dictionary doesn't exist, use a brute force method if not self.node_bonds_dic: self.brute_generate_pi_bonds() return + # Otherwise we can use the dictionary that maps nodes to the bonds that contain them + # (essentially iterating over possible vertices) for node in self.node_bonds_dic.keys(): bonds = self.node_bonds_dic[node] + # Iterate through possible pairs of bonds that contain [node] for i in range(len(bonds) - 1): bond_i = bonds[i] for j in range(i + 1, len(bonds)): bond_j = bonds[j] - # Define vertex and edges + + # Define vertex and edge nodes vertex = node e1 = (bond_i.get_node1() if node != bond_i.get_node1() else bond_i.get_node2()) e2 = (bond_j.get_node1() if node != bond_j.get_node1() else bond_j.get_node2()) - # Get positions of the two nodes per each of two bonds - pos1_i = bond_i.get_node1().get_xy() - pos2_i = bond_i.get_node2().get_xy() - pos1_j = bond_j.get_node1().get_xy() - pos2_j = bond_j.get_node2().get_xy() - - # Find the unit vector of each bond - r_i = np.subtract(pos2_i, pos1_i) - r_j = np.subtract(pos2_j, pos1_j) - - # Account for edge bonds being shorter in distance - if bond_i.is_hor_pbc(): - if pos1_i[0] > pos2_i[0]: - r_i[0] = r_i[0] + self.length - else: - r_i[0] = r_i[0] - self.length - if bond_j.is_hor_pbc(): - if pos1_j[0] > pos2_j[0]: - r_j[0] = r_j[0] + self.length - else: - r_j[0] = r_j[0] - self.length - - # Normalize - r_i = r_i / np.linalg.norm(r_i, axis=0) - r_j = r_j / np.linalg.norm(r_j, axis=0) - # The two bonds should be co-linear - cos_phi = np.dot(r_i, r_j) - if abs(abs(cos_phi) - 1) < 0.01: + if self.check_bonds_colinear(bond_i, bond_j, vertex): pi_bond = PiBond(bond_i, bond_j, vertex, e1, e2) self.pi_bonds.append(pi_bond) print("Generated " + str(len(self.pi_bonds)) + " pi bonds") @@ -405,83 +414,38 @@ def brute_generate_pi_bonds(self) -> None: # Loop through the bonds to find common vertices (node with two bonds) num_bonds = len(self.bonds) bonds = self.bonds - vertices = 0 - vertex = None - e1, e2 = None, None + # Iterate through pairs of bonds for i in range(num_bonds - 1): bond_i = bonds[i] for j in range(i + 1, num_bonds): bond_j = bonds[j] # Whether these two bonds are joined by a common node vertex - is_vertex = False - # First node of first bond is first node of second bond + found_common_vertex = False if bond_i.get_node1() == bond_j.get_node1(): - is_vertex = True + found_common_vertex = True vertex = bond_i.get_node1() e1 = bond_i.get_node2() e2 = bond_j.get_node2() - # Second node of first bond is second node of second bond elif bond_i.get_node2() == bond_j.get_node2(): - is_vertex = True + found_common_vertex = True vertex = bond_i.get_node2() e1 = bond_i.get_node1() e2 = bond_j.get_node1() - # First node of first bond is second node of second bond elif bond_i.get_node1() == bond_j.get_node2(): - is_vertex = True + found_common_vertex = True vertex = bond_i.get_node1() e1 = bond_i.get_node2() e2 = bond_j.get_node1() - # Second node of first bond is first node of second bond elif bond_i.get_node2() == bond_j.get_node1(): - is_vertex = True + found_common_vertex = True vertex = bond_i.get_node2() e1 = bond_i.get_node1() e2 = bond_j.get_node2() - if is_vertex: - vertices += 1 - # Get positions of the two nodes per each of two bonds - pos1_i = bond_i.get_node1().get_xy() - pos2_i = bond_i.get_node2().get_xy() - pos1_j = bond_j.get_node1().get_xy() - pos2_j = bond_j.get_node2().get_xy() - - # Find the unit vector of each bond - r_i = np.subtract(pos2_i, pos1_i) - r_j = np.subtract(pos2_j, pos1_j) - - # Account for edge bonds being shorter in distance - if bond_i.is_hor_pbc(): - if pos1_i[0] > pos2_i[0]: - r_i[0] = r_i[0] + self.length - else: - r_i[0] = r_i[0] - self.length - if bond_j.is_hor_pbc(): - if pos1_j[0] > pos2_j[0]: - r_j[0] = r_j[0] + self.length - else: - r_j[0] = r_j[0] - self.length - - # Account for top periodic bonds being shorter in distance - height_change = self.height + self.height_increment - if bond_i.is_temporary(): - if pos1_i[1] > pos2_i[1]: - r_i[1] = r_i[1] + height_change - else: - r_i[1] = r_i[1] - height_change - if bond_j.is_temporary(): - if pos1_j[1] > pos2_j[1]: - r_j[1] = r_j[1] + height_change - else: - r_j[1] = r_j[1] - height_change - - # Normalize - r_i = r_i / np.linalg.norm(r_i, axis=0) - r_j = r_j / np.linalg.norm(r_j, axis=0) - # The two bonds should be co-linear - cos_phi = np.dot(r_i, r_j) - if abs(abs(cos_phi) - 1) < 0.01: - pi_bond = PiBond(bond_i, bond_j, vertex, e1, e2) - self.pi_bonds.append(pi_bond) + # If no common vertex, continue + if not found_common_vertex: + continue + + if self.check_bonds_colinear(bond_i, bond_j, vertex): + pi_bond = PiBond(bond_i, bond_j, vertex, e1, e2) + self.pi_bonds.append(pi_bond) print("Generated " + str(len(self.pi_bonds)) + " pi bonds") - # print("Generated " + str(vertices) + " vertices") diff --git a/lattice/lattice_factory.py b/lattice/lattice_factory.py index e74a657..bb8660b 100644 --- a/lattice/lattice_factory.py +++ b/lattice/lattice_factory.py @@ -16,7 +16,7 @@ class LatticeFactory: @staticmethod def create_lattice( - lattice_type: LatticeType, length: int, height: float, generate_pi_bonds: bool + lattice_type: LatticeType, length: int, height: float ) -> AbstractLattice: """ Create a fresh lattice from scratch @@ -24,7 +24,7 @@ def create_lattice( if lattice_type == LatticeType.KAGOME: return KagomeLattice(length=length, height=height) elif lattice_type == LatticeType.TRIANGULAR: - return TriangularLattice(length=length, height=height, generate_pi_bonds=generate_pi_bonds) + return TriangularLattice(length=length, height=height) elif lattice_type == LatticeType.SQUARE: return SquareLattice(length=length, height=height) else: diff --git a/lattice/triangular_lattice.py b/lattice/triangular_lattice.py index d51da68..37f71c9 100644 --- a/lattice/triangular_lattice.py +++ b/lattice/triangular_lattice.py @@ -43,7 +43,7 @@ def generate_nodes(self, length, height) -> List[Node]: def generate_bonds(self) -> None: super().generate_bonds() - def __init__(self, length: int, height: float, generate_pi_bonds: bool): + def __init__(self, length: int, height: float): self.length = length self.height = height self.nodes = [] @@ -54,5 +54,4 @@ def __init__(self, length: int, height: float, generate_pi_bonds: bool): self.height_increment = math.sqrt(3.0) / 2.0 self.generate_nodes(length, height) self.generate_bonds() - if generate_pi_bonds: - self.generate_pi_bonds() + self.generate_pi_bonds() diff --git a/main.py b/main.py index 6470656..7e93689 100644 --- a/main.py +++ b/main.py @@ -25,7 +25,6 @@ def initialize_lattice(rng: np.random.Generator) -> AbstractLattice: lattice_type=Parameters.lattice_type, length=Parameters.lattice_length, height=Parameters.lattice_height, - generate_pi_bonds=Parameters.bend_mod != 0, # Naive optimization: only make pi bonds if bending is relevant ) # --- Initialization steps --- diff --git a/parameters.py b/parameters.py index 8f808d0..3391038 100644 --- a/parameters.py +++ b/parameters.py @@ -7,7 +7,7 @@ import numpy as np from energyminimization.solvers.solver import MinimizationType -from energyminimization.transformations import StretchX, StretchY +from energyminimization.transformations import StretchX, StretchY, Shear, Dilate from lattice.lattice_type import LatticeType from result_handling.output_handler import OutputHandlerParameters from result_handling.pickle_handler import VisualizationHandlerParameters @@ -36,7 +36,7 @@ class Parameters: # ----- Bond Occupation Protocol ----- # Starting p_fill value - prob_fill_high: float = 0.68 + prob_fill_high: float = 0.6 # When we remove bonds, this is our stopping criteria: all the moduli are below this value moduli_tolerance: float = 1e-8 diff --git a/tests/matrix_tests.py b/tests/matrix_tests.py index 527f618..157120b 100644 --- a/tests/matrix_tests.py +++ b/tests/matrix_tests.py @@ -11,24 +11,37 @@ def test_gradient_hessian( df: Callable[[np.ndarray], np.ndarray], hess: Callable[[np.ndarray], csr_matrix], ): + """ + Use finite differences to test implementations of the gradient and hessian of a function + """ n_tests = 500 - results_df = np.zeros(n_tests) - results_d2f = np.zeros(n_tests) - results_d2f2 = np.zeros((n_tests, x0.size)) + # Results contain the relative error of the finite difference approximation + results_df = np.zeros(n_tests) # gradient + results_d2f = np.zeros(n_tests) # hessian + results_d2f2 = np.zeros((n_tests, x0.size)) # hessian using finite differences of gradient (assuming df is correct) + for i in range(n_tests): - x = np.random.rand(x0.size) - u = np.random.rand(x0.size) - h = 1e-6 + x = x0 + u = np.random.rand(x0.size) # Pick a random displacement + h = 1e-5 # Step size + + # Evaluate the function at x, x + h*u, x - h*u + f0 = f(x) fp = f(x + h * u) fm = f(x - h * u) - f0 = f(x) + # Finite difference approximation of the gradient + df_test = (fp - fm) / (2 * h) + # Finite difference approximation of the hessian + d2f_test = (fp - 2 * f0 + fm) / (h ** 2) + + # Obtain our implementation of the gradient and hessian df0 = df(x) h0 = hess(x) - df_test = (fp - fm) / (2 * h) - d2f_test = (fp - 2 * f0 + fm) / (h ** 2) - # also test hessian using finite differences of gradient + # Finite difference approximation of the hessian using the gradient implementation d2f_test2 = (df(x + h * u) - df(x - h * u)) / (2 * h) + + # Compute the relative error results_df[i] = np.abs(df_test - (df0.T @ u)) / (np.linalg.norm(df0.T @ u)) results_d2f[i] = np.abs(d2f_test - (u.T @ h0 @ u)) / (np.linalg.norm(u.T @ h0 @ u)) results_d2f2[i] = np.linalg.norm(d2f_test2 - (h0 @ u)) / (np.linalg.norm(h0 @ u))