JGalego · September 15, 2025 14:47
diff --git a/hopfield.py b/hopfield.py
 # pylint: disable=invalid-name,non-ascii-name,line-too-long

 # /// script
 # requires-python = ">=3.12"
 # dependencies = [
 #     "numpy",
 #     "matplotlib",
 #     "Pillow",
 # ]
 # ///

 """
 Implementation of a traditional Hopfield network in Python.

 > "The Hopfield network is really just an Ising model at zero temperature
  with a site-dependent magnetic field and long-ranged interactions."

 Articles 📝

 + (Alexander et al., 2021) The Physics of Machine Learning: An Intuitive Introduction for the Physical Scientist
 https://arxiv.org/abs/2112.00851

 + (Amari, 1977) Neural theory of association and concept-formation
 https://link.springer.com/article/10.1007/BF00365229

 + (Hopfield, 1982) Neural networks and physical systems with emergent collective computational abilities
 https://www.pnas.org/doi/pdf/10.1073/pnas.79.8.2554

 Blogs / Courses 🧑‍🏫

 + (JKU) Hopfield Networks is All You Need
 https://ml-jku.github.io/hopfield-layers/

 + (University of Illinois) Computing with Physics > Hopfield Networks
 http://clark.physics.illinois.edu/498cmp/secure/3-Brain/html/1-Hopfield.html
 """

 # Standard imports
 import sys

 from typing import List, Tuple

 # Library imports
 import numpy as np

 from matplotlib import cm
 from matplotlib import pyplot as plt

 from PIL import Image

 def display_state(ξ: np.array, title: str = None) -> None:
    """
    Visualize a network state as a grayscale image using matplotlib.

    This function displays the state of a Hopfield network as a 2D image where
    positive values appear white and negative values appear black.

    Args:
        ξ (np.array): 2D numpy array representing the network state, typically
                      reshaped from a flattened state vector. Values should be
                      binary (-1 or +1).
        title (str, optional): Title to display above the image. Defaults to None.

    Returns:
        None: Displays the plot using matplotlib.show()

    Example:
        >>> state = np.array([[1, -1, 1], [-1, 1, -1], [1, -1, 1]])
        >>> display_state(state, "Pattern #1")
    """
    ax = plt.axes()
    ax.matshow(ξ, cmap='binary')
    plt.title(title)
    plt.axis('off')
    plt.show()

 def hebb_matrix(x: np.ndarray) -> np.ndarray:
    """
    Compute the Hebbian weight matrix for a given pattern.

    This function implements Hebb's rule for unsupervised learning, creating
    a weight matrix that encodes the correlations in the input pattern. The
    diagonal is set to zero to prevent self-connections.

    The Hebbian rule states that the connection weight between neurons i and j
    should be proportional to the correlation between their activities:
    w_ij = (1/N) * x_i * x_j, where N is the dimensionality.

    Args:
        x (np.ndarray): 1D binary pattern vector with values in {-1, +1}.
                        Shape should be (N,) where N is the number of neurons.

    Returns:
        np.ndarray: 2D weight matrix of shape (N, N) with zero diagonal.
                    Symmetric matrix encoding pattern correlations.

    Example:
        >>> pattern = np.array([1, -1, 1, -1])
        >>> weights = hebb_matrix(pattern)
        >>> print(weights.shape)
        (4, 4)
    """
    d = x.shape[0]
    w = np.outer(x, x) / d
    np.fill_diagonal(w, 0)
    return w

 class HopfieldNetwork():
    """
    A traditional Hopfield network implementation for associative memory.

    The Hopfield network is a form of recurrent artificial neural network that
    serves as a content-addressable memory system with binary threshold nodes.
    It can store and retrieve patterns through energy minimization dynamics.

    The network uses symmetric weights (w_ij = w_ji) and operates by iteratively
    updating neuron states until reaching a stable configuration (local energy minimum).

    Attributes:
        N (int): Number of neurons in the network.
        W (np.ndarray): Symmetric weight matrix of shape (N, N) with zero diagonal.
        b (np.ndarray): Bias vector of shape (N,) for external fields.

    Example:
        >>> # Create network with 100 neurons
        >>> network = HopfieldNetwork(100)
        >>> # Train on some patterns
        >>> patterns = [np.random.choice([-1, 1], 100) for _ in range(3)]
        >>> network.train(patterns)
        >>> # Retrieve pattern from noisy input
        >>> noisy = patterns[0] + 0.3 * np.random.randn(100)
        >>> history = network.run(np.sign(noisy))
    """

    def __init__(self, N: int):
        """
        Initialize a Hopfield network with N neurons.

        Creates a new Hopfield network with zero-initialized weights and biases.
        The weight matrix is symmetric with zero diagonal elements (no self-connections).

        Args:
            N (int): Number of neurons in the network. Must be positive.

        Raises:
            ValueError: If N is not a positive integer.

        Example:
            >>> network = HopfieldNetwork(100)  # 100-neuron network
        """
        if N <= 0:
            raise ValueError("Number of neurons N must be positive")

        self.N = N                 # number of neurons
        self.W = np.zeros((N, N))  # weight matrix
        self.b = np.zeros(N)       # bias vector

    def train(self, data: np.ndarray) -> None:
        """
        Train the network on a dataset using Hebbian learning rule.

        This method implements unsupervised learning by updating the weight matrix
        according to Hebb's rule. Each pattern in the dataset contributes to the
        weight matrix through outer product correlations.

        The learning is cumulative - multiple calls to train() will add new patterns
        to the existing memory without overwriting previous ones.

        Args:
            data (np.ndarray): Collection of binary patterns to store in memory.
                              Can be:
                              - 2D array of shape (n_patterns, N) where each row is a pattern
                              - List of 1D arrays, each of length N
                              Each pattern should contain values in {-1, +1}.

        Returns:
            None: Updates the weight matrix W in-place.

        Raises:
            ValueError: If patterns have inconsistent dimensions or wrong values.

        Example:
            >>> network = HopfieldNetwork(4)
            >>> patterns = np.array([[1, -1, 1, -1], [-1, 1, -1, 1]])
            >>> network.train(patterns)
        """
        for x in data:
            self.W += hebb_matrix(x)

    def plot_weights(self) -> None:
        """
        Visualize the weight matrix as a heatmap.

        This method creates a color-coded visualization of the network's weight matrix,
        helping to understand the learned correlations between neurons. Positive weights
        appear in warm colors (red), negative weights in cool colors (blue), and zero
        weights are neutral.

        The visualization uses a diverging colormap to emphasize both positive and
        negative correlations equally.

        Returns:
            None: Displays the plot using matplotlib.show()

        Note:
            This method requires matplotlib to be available and will block execution
            until the plot window is closed.

        Example:
            >>> network = HopfieldNetwork(10)
            >>> network.train([np.random.choice([-1, 1], 10)])
            >>> network.plot_weights()  # Shows 10x10 heatmap
        """
        w_mat = plt.imshow(self.W, cmap=cm.coolwarm)  # pylint: disable=no-member
        plt.title("Weight Matrix")
        plt.colorbar(w_mat)
        plt.tight_layout()
        plt.axis('off')
        plt.show()

    def energy(self, ξ: np.ndarray) -> np.float32:
        """
        Compute the energy function for a given network state.

        The energy function E(ξ) = -1/2 * ξᵀ * W * ξ + Σ(b_i * ξ_i) represents
        the "energy" of the current configuration. The network dynamics minimize
        this energy, causing convergence to local minima (stored patterns).

        Lower energy indicates more stable configurations that correspond to
        stored memories. The energy decreases monotonically during updates
        (except for synchronous updates which may oscillate).

        Args:
            ξ (np.ndarray): Current network state as 1D binary array of shape (N,).
                           Values should be in {-1, +1}.

        Returns:
            np.float32: Energy value for the given state. More negative values
                       indicate more stable configurations.

        Mathematical Note:
            The factor of 1/2 prevents double-counting since W is symmetric.
            The bias term b represents external magnetic fields.

        Example:
            >>> network = HopfieldNetwork(3)
            >>> state = np.array([1, -1, 1])
            >>> energy_val = network.energy(state)
            >>> print(f"Energy: {energy_val}")
        """
        return -0.5 * ξ @ self.W @ ξ + np.sum(ξ * self.b)

    def sync_update_step(self, ξ: np.array) -> Tuple[np.ndarray, float]:
        """
        Perform one synchronous update step on the network state.

        In synchronous updating, all neurons are updated simultaneously using
        the current state. Each neuron's new value is determined by:
        ξ_i^(t+1) = sign(Σ_j W_ij * ξ_j^(t) - b_i)

        This update rule can lead to oscillations between states and doesn't
        guarantee energy decrease at each step, unlike asynchronous updates.

        Args:
            ξ (np.array): Current network state as 1D binary array of shape (N,).
                         Values should be in {-1, +1}.

        Returns:
            tuple: (ξ_new, E_new) where:
                - ξ_new (np.array): Updated network state after synchronous step
                - E_new (float): Energy of the new state

        Note:
            Synchronous updates may cause the network to oscillate between
            configurations rather than converge to a fixed point.

        Example:
            >>> network = HopfieldNetwork(4)
            >>> state = np.array([1, -1, 1, -1])
            >>> new_state, new_energy = network.sync_update_step(state)
        """
        ξ = np.sign(self.W @ ξ - self.b)
        E = self.energy(ξ)
        return ξ, E

    def async_update_step(self, ξ: np.array) -> Tuple[np.ndarray, float]:
        """
        Perform one asynchronous update step on the network state.

        In asynchronous updating, neurons are updated one at a time in random order.
        Each neuron's new value is computed using the most current states of all
        other neurons. This guarantees energy decrease (or no change) at each step,
        ensuring convergence to a local minimum.

        The update rule for neuron i is:
        ξ_i^(new) = sign(Σ_j W_ij * ξ_j^(current) - b_i)

        Args:
            ξ (np.array): Current network state as 1D binary array of shape (N,).
                         Values should be in {-1, +1}. This array is modified in-place.

        Returns:
            tuple: (ξ, E) where:
                - ξ (np.array): The same input array after in-place modification
                - E (float): Energy of the updated state

        Note:
            This method modifies the input array in-place for efficiency.
            The random permutation ensures fair updating and helps avoid cycles.

        Example:
            >>> network = HopfieldNetwork(4)
            >>> state = np.array([1, -1, 1, -1])
            >>> updated_state, energy = network.async_update_step(state)
            >>> # state has been modified in-place
        """
        for idx in np.random.permutation(self.N):
            ξ[idx] = np.sign(self.W[idx] @ ξ - self.b[idx])
        E = self.energy(ξ)
        return ξ, E

    def run(self,
            ξ: np.array,
            tol: float = 1e-4,
            max_iter: int = 10,
            sync: bool = True) -> List[Tuple[np.ndarray, float]]:
        """
        Run the network dynamics until convergence or maximum iterations.

        This method iteratively applies update steps (either synchronous or
        asynchronous) starting from an initial state until the network converges
        to a stable configuration or the maximum number of iterations is reached.

        Convergence is detected by monitoring energy changes between consecutive
        steps. The network has converged when the energy change falls below
        the specified tolerance.

        Args:
            ξ (np.array): Initial network state as 1D binary array of shape (N,).
                         Values should be in {-1, +1}.
            tol (float, optional): Energy change threshold for convergence detection.
                                  Defaults to 1e-4.
            max_iter (int, optional): Maximum number of update iterations.
                                     Defaults to 10.
            sync (bool, optional): If True, use synchronous updates; if False,
                                  use asynchronous updates. Defaults to True.

        Returns:
            list: History of (state, energy) tuples tracking the network's evolution.
                 Each tuple contains:
                 - state (np.array): Network configuration at that step
                 - energy (float): Energy value for that configuration

        Note:
            - Asynchronous updates (sync=False) are generally preferred as they
              guarantee convergence to local minima
            - Synchronous updates may oscillate and not converge
            - The initial state and energy are included as the first history entry

        Example:
            >>> network = HopfieldNetwork(4)
            >>> initial_state = np.random.choice([-1, 1], 4)
            >>> history = network.run(initial_state, max_iter=50, sync=False)
            >>> final_state, final_energy = history[-1]
            >>> print(f"Converged after {len(history)-1} steps")
        """
        E = self.energy(ξ)
        history = [(ξ, E)]
        for _ in range(max_iter):
            ξ_new, E_new = self.sync_update_step(ξ) if sync \
                      else self.async_update_step(ξ)
            history.append((ξ_new, E_new))
            if abs(E_new - E) < tol:
                break
            ξ, E = ξ_new, E_new
        return history


 def preprocess_image(image_path: str, size: tuple = (64, 64)) -> tuple:
    """
    Load and preprocess an image for use with Hopfield network.

    This function performs the complete preprocessing pipeline:
    1. Load image and convert to grayscale
    2. Resize to specified dimensions
    3. Apply mean thresholding
    4. Convert to binary values {-1, +1}
    5. Flatten to 1D array

    Args:
        image_path (str): Path to the image file to process.
        size (tuple, optional): Target size as (width, height). Defaults to (64, 64).

    Returns:
        tuple: (pattern, original_shape) where:
            - pattern (np.ndarray): Flattened binary pattern of shape (width*height,)
            - original_shape (tuple): Original 2D shape (width, height) for reshaping

    Raises:
        IOError: If the image file cannot be loaded.
        ValueError: If the image cannot be processed.

    Example:
        >>> pattern, shape = preprocess_image("image.jpg", size=(32, 32))
        >>> print(f"Pattern shape: {pattern.shape}, Original: {shape}")
    """
    try:
        # Load image: original -> grayscale -> resizing
        pic = Image.open(image_path).convert('L')
        pic = pic.resize(size)
        pix = np.array(pic)

        # Preprocess: mean thresholding -> sign -> flatten
        width, height = pix.shape
        pix = pix.astype(np.float32)  # Ensure proper dtype for mean calculation
        pix -= np.mean(pix)
        pix = np.sign(pix).astype(np.int8)
        pix = pix.flatten()

        return pix, (width, height)

    except IOError as e:
        raise IOError(f"Could not load image {image_path}: {e}") from e
    except Exception as e:
        raise ValueError(f"Error processing image {image_path}: {e}") from e


 def load_multiple_images(img_paths: list, size: tuple = (64, 64)) -> tuple:
    """
    Load and preprocess multiple images for training a Hopfield network.

    All images are resized to the same dimensions to ensure consistent
    pattern sizes for the network. Invalid images are skipped with warnings.

    Args:
        img_paths (list): List of file paths to image files.
        size (tuple, optional): Target size as (width, height). Defaults to (64, 64).

    Returns:
        tuple: (patterns, shape, valid_paths) where:
            - patterns (list): List of preprocessed binary patterns
            - shape (tuple): Common shape (width, height) for all patterns
            - valid_paths (list): List of successfully loaded image paths

    Example:
        >>> paths = ["img1.jpg", "img2.png", "img3.gif"]
        >>> patterns, shape, valid = load_multiple_images(paths)
        >>> print(f"Loaded {len(patterns)} images with shape {shape}")
    """
    loaded_patterns = []
    loaded_paths = []
    common_shape = None

    for file_path in img_paths:
        try:
            loaded_pattern, img_shape = preprocess_image(file_path, size)
            loaded_patterns.append(loaded_pattern)
            loaded_paths.append(file_path)
            if common_shape is None:
                common_shape = img_shape

        except (IOError, ValueError) as e:
            print(f"Warning: Skipping {file_path} - {e}")
            continue

    if not loaded_patterns:
        raise ValueError("No valid images could be loaded")

    return loaded_patterns, common_shape, loaded_paths


 def demonstrate_recall(network: HopfieldNetwork,
                      img_shape: tuple,
                      num_tests: int = 3,
                      noise_level: float = 0.2,
                      max_iter: int = 20) -> None:
    """
    Demonstrate pattern recall from noisy or random initial states.

    This function tests the network's ability to recover stored patterns
    from corrupted inputs, which is the key capability of associative memory.

    Args:
        network (HopfieldNetwork): Trained Hopfield network.
        img_shape (tuple): Image dimensions (width, height) for visualization.
        num_tests (int, optional): Number of recall tests to perform. Defaults to 3.
        noise_level (float, optional): Fraction of pixels to flip for noise test.
                                       Defaults to 0.2 (20% corruption).
        max_iter (int, optional): Maximum iterations for network dynamics. Defaults to 20.

    Returns:
        None: Displays recall results using matplotlib.

    Example:
        >>> demonstrate_recall(network, (64, 64), num_tests=5)
    """
    width, height = img_shape
    total_pixels = width * height

    print(f"\nDemonstrating pattern recall with {num_tests} tests...")

    for test_idx in range(num_tests):
        print(f"\n--- Recall Test #{test_idx + 1} ---")

        # Generate initial state based on test type
        initial_state = _generate_test_state(test_idx, total_pixels, noise_level)

        # Run the network dynamics
        history = network.run(initial_state.copy(), max_iter=max_iter, sync=False)

        # Display results
        _display_test_results(initial_state, history, img_shape, test_idx)


 def _generate_test_state(test_idx: int, total_pixels: int, noise_level: float) -> np.ndarray:
    """Generate initial state for recall test."""
    if test_idx == 0:
        # Test 1: Random noise
        print("Starting from random noise...")
        return np.random.choice([-1, 1], total_pixels)

    # Test 2+: Add noise to a structured pattern
    print(f"Starting from corrupted pattern (noise level: {noise_level:.1%})...")
    # Add some structure by flipping only a fraction of bits
    mask = np.random.random(total_pixels) < noise_level
    reference = np.random.choice([-1, 1], total_pixels)
    return np.where(mask, -reference, reference)


 def _display_test_results(initial_state: np.ndarray,
                         history: List[Tuple[np.ndarray, float]],
                         img_shape: tuple,
                         test_idx: int) -> None:
    """
    Display the visual results of a pattern recall test.
    
    This helper function visualizes both the initial state and final converged state
    of a Hopfield network recall test. It shows the network's ability to recover
    stored patterns from noisy or random inputs by displaying before-and-after
    images with energy and convergence information.
    
    The function creates two matplotlib plots:
    1. Initial state: The starting configuration (noisy/random)
    2. Final state: The converged result with energy and step count
    
    Args:
        initial_state (np.ndarray): Starting network state as 1D binary array.
                                   Values should be in {-1, +1}.
        history (List[Tuple[np.ndarray, float]]): Complete evolution history from
                                                  network.run(), where each tuple
                                                  contains (state, energy).
        img_shape (tuple): Image dimensions (width, height) for reshaping the
                          flattened states into 2D images for visualization.
        test_idx (int): Zero-based test index for labeling the displayed results.
    
    Returns:
        None: Displays matplotlib plots and prints convergence statistics.
    
    Side Effects:
        - Creates and displays two matplotlib figure windows
        - Prints convergence information to console
        - Blocks execution until plot windows are closed (if running interactively)
    
    Example:
        >>> initial = np.random.choice([-1, 1], 64*64)
        >>> history = network.run(initial.copy(), max_iter=20)
        >>> _display_test_results(initial, history, (64, 64), 0)
        # Shows: "Test #1: Initial State" and "Test #1: Final State"
    """
    initial_image = initial_state.reshape(img_shape)
    final_state, final_energy = history[-1]
    final_image = final_state.reshape(img_shape)

    display_state(initial_image, f"Test #{test_idx + 1}: Initial State")
    display_state(final_image,
                 f"Test #{test_idx + 1}: Final State (Energy: {final_energy:.2f}, Steps: {len(history)-1})")

    print(f"Converged after {len(history)-1} steps with energy {final_energy:.2f}")


 if __name__ == '__main__':
    # Parse command line arguments
    if len(sys.argv) < 2:
        print("Usage: python hopfield.py <image1> [image2] [image3] ...")
        print("Example: python hopfield.py lenna.png face1.jpg face2.png")
        sys.exit(1)

    image_paths = sys.argv[1:]
    print(f"Loading {len(image_paths)} image(s)...")

    # Load and preprocess multiple images
    try:
        patterns, shape, valid_paths = load_multiple_images(image_paths, size=(64, 64))
        print(f"Successfully loaded {len(patterns)} images with dimensions {shape}")
    except ValueError as e:
        print(f"Error: {e}")
        sys.exit(1)

    # Display original images
    print("\nDisplaying original images...")
    for i, (train_pattern, img_path) in enumerate(zip(patterns, valid_paths)):
        image_name = img_path.split('/')[-1]  # Get filename only
        display_state(train_pattern.reshape(shape), f"Original Image #{i+1}: {image_name}")

    # Initialize the network
    network_size = len(patterns[0])
    print(f"\nInitializing Hopfield network with {network_size} neurons...")
    hn = HopfieldNetwork(network_size)

    # Train on all images
    print(f"Training network on {len(patterns)} pattern(s)...")
    hn.train(patterns)

    # Display weight matrix
    print("Displaying learned weight matrix...")
    hn.plot_weights()

    # Test pattern recall
    demonstrate_recall(hn, shape, num_tests=min(3, len(patterns) + 1))

    print(f"\nTraining complete! Network learned {len(patterns)} pattern(s).")
    print(f"Network capacity (rule of thumb): ~{network_size // (4 * len(patterns))} patterns for reliable recall")

    if len(patterns) > network_size // (4 * len(patterns)):
        print("Warning: You may have exceeded the network's storage capacity.")
        print("Consider using fewer patterns or larger images for better recall performance.")
	# pylint: disable=invalid-name,non-ascii-name,line-too-long

	# /// script
	# requires-python = ">=3.12"
	# dependencies = [
	# "numpy",
	# "matplotlib",
	# "Pillow",
	# ]
	# ///

	"""
	Implementation of a traditional Hopfield network in Python.

	> "The Hopfield network is really just an Ising model at zero temperature
	with a site-dependent magnetic field and long-ranged interactions."

	Articles 📝

	+ (Alexander et al., 2021) The Physics of Machine Learning: An Intuitive Introduction for the Physical Scientist
	https://arxiv.org/abs/2112.00851

	+ (Amari, 1977) Neural theory of association and concept-formation
	https://link.springer.com/article/10.1007/BF00365229

	+ (Hopfield, 1982) Neural networks and physical systems with emergent collective computational abilities
	https://www.pnas.org/doi/pdf/10.1073/pnas.79.8.2554

	Blogs / Courses 🧑‍🏫

	+ (JKU) Hopfield Networks is All You Need
	https://ml-jku.github.io/hopfield-layers/

	+ (University of Illinois) Computing with Physics > Hopfield Networks
	http://clark.physics.illinois.edu/498cmp/secure/3-Brain/html/1-Hopfield.html
	"""

	# Standard imports
	import sys

	from typing import List, Tuple

	# Library imports
	import numpy as np

	from matplotlib import cm
	from matplotlib import pyplot as plt

	from PIL import Image

	def display_state(ξ: np.array, title: str = None) -> None:
	"""
	Visualize a network state as a grayscale image using matplotlib.

	This function displays the state of a Hopfield network as a 2D image where
	positive values appear white and negative values appear black.

	Args:
	ξ (np.array): 2D numpy array representing the network state, typically
	reshaped from a flattened state vector. Values should be
	binary (-1 or +1).
	title (str, optional): Title to display above the image. Defaults to None.

	Returns:
	None: Displays the plot using matplotlib.show()

	Example:
	>>> state = np.array([[1, -1, 1], [-1, 1, -1], [1, -1, 1]])
	>>> display_state(state, "Pattern #1")
	"""
	ax = plt.axes()
	ax.matshow(ξ, cmap='binary')
	plt.title(title)
	plt.axis('off')
	plt.show()

	def hebb_matrix(x: np.ndarray) -> np.ndarray:
	"""
	Compute the Hebbian weight matrix for a given pattern.

	This function implements Hebb's rule for unsupervised learning, creating
	a weight matrix that encodes the correlations in the input pattern. The
	diagonal is set to zero to prevent self-connections.

	The Hebbian rule states that the connection weight between neurons i and j
	should be proportional to the correlation between their activities:
	w_ij = (1/N) * x_i * x_j, where N is the dimensionality.

	Args:
	x (np.ndarray): 1D binary pattern vector with values in {-1, +1}.
	Shape should be (N,) where N is the number of neurons.

	Returns:
	np.ndarray: 2D weight matrix of shape (N, N) with zero diagonal.
	Symmetric matrix encoding pattern correlations.

	Example:
	>>> pattern = np.array([1, -1, 1, -1])
	>>> weights = hebb_matrix(pattern)
	>>> print(weights.shape)
	(4, 4)
	"""
	d = x.shape[0]
	w = np.outer(x, x) / d
	np.fill_diagonal(w, 0)
	return w

	class HopfieldNetwork():
	"""
	A traditional Hopfield network implementation for associative memory.

	The Hopfield network is a form of recurrent artificial neural network that
	serves as a content-addressable memory system with binary threshold nodes.
	It can store and retrieve patterns through energy minimization dynamics.

	The network uses symmetric weights (w_ij = w_ji) and operates by iteratively
	updating neuron states until reaching a stable configuration (local energy minimum).

	Attributes:
	N (int): Number of neurons in the network.
	W (np.ndarray): Symmetric weight matrix of shape (N, N) with zero diagonal.
	b (np.ndarray): Bias vector of shape (N,) for external fields.

	Example:
	>>> # Create network with 100 neurons
	>>> network = HopfieldNetwork(100)
	>>> # Train on some patterns
	>>> patterns = [np.random.choice([-1, 1], 100) for _ in range(3)]
	>>> network.train(patterns)
	>>> # Retrieve pattern from noisy input
	>>> noisy = patterns[0] + 0.3 * np.random.randn(100)
	>>> history = network.run(np.sign(noisy))
	"""

	def __init__(self, N: int):
	"""
	Initialize a Hopfield network with N neurons.

	Creates a new Hopfield network with zero-initialized weights and biases.
	The weight matrix is symmetric with zero diagonal elements (no self-connections).

	Args:
	N (int): Number of neurons in the network. Must be positive.

	Raises:
	ValueError: If N is not a positive integer.

	Example:
	>>> network = HopfieldNetwork(100) # 100-neuron network
	"""
	if N <= 0:
	raise ValueError("Number of neurons N must be positive")

	self.N = N # number of neurons
	self.W = np.zeros((N, N)) # weight matrix
	self.b = np.zeros(N) # bias vector

	def train(self, data: np.ndarray) -> None:
	"""
	Train the network on a dataset using Hebbian learning rule.

	This method implements unsupervised learning by updating the weight matrix
	according to Hebb's rule. Each pattern in the dataset contributes to the
	weight matrix through outer product correlations.

	The learning is cumulative - multiple calls to train() will add new patterns
	to the existing memory without overwriting previous ones.

	Args:
	data (np.ndarray): Collection of binary patterns to store in memory.
	Can be:
	- 2D array of shape (n_patterns, N) where each row is a pattern
	- List of 1D arrays, each of length N
	Each pattern should contain values in {-1, +1}.

	Returns:
	None: Updates the weight matrix W in-place.

	Raises:
	ValueError: If patterns have inconsistent dimensions or wrong values.

	Example:
	>>> network = HopfieldNetwork(4)
	>>> patterns = np.array([[1, -1, 1, -1], [-1, 1, -1, 1]])
	>>> network.train(patterns)
	"""
	for x in data:
	self.W += hebb_matrix(x)

	def plot_weights(self) -> None:
	"""
	Visualize the weight matrix as a heatmap.

	This method creates a color-coded visualization of the network's weight matrix,
	helping to understand the learned correlations between neurons. Positive weights
	appear in warm colors (red), negative weights in cool colors (blue), and zero
	weights are neutral.

	The visualization uses a diverging colormap to emphasize both positive and
	negative correlations equally.

	Returns:
	None: Displays the plot using matplotlib.show()

	Note:
	This method requires matplotlib to be available and will block execution
	until the plot window is closed.

	Example:
	>>> network = HopfieldNetwork(10)
	>>> network.train([np.random.choice([-1, 1], 10)])
	>>> network.plot_weights() # Shows 10x10 heatmap
	"""
	w_mat = plt.imshow(self.W, cmap=cm.coolwarm) # pylint: disable=no-member
	plt.title("Weight Matrix")
	plt.colorbar(w_mat)
	plt.tight_layout()
	plt.axis('off')
	plt.show()

	def energy(self, ξ: np.ndarray) -> np.float32:
	"""
	Compute the energy function for a given network state.

	The energy function E(ξ) = -1/2 * ξᵀ * W * ξ + Σ(b_i * ξ_i) represents
	the "energy" of the current configuration. The network dynamics minimize
	this energy, causing convergence to local minima (stored patterns).

	Lower energy indicates more stable configurations that correspond to
	stored memories. The energy decreases monotonically during updates
	(except for synchronous updates which may oscillate).

	Args:
	ξ (np.ndarray): Current network state as 1D binary array of shape (N,).
	Values should be in {-1, +1}.

	Returns:
	np.float32: Energy value for the given state. More negative values
	indicate more stable configurations.

	Mathematical Note:
	The factor of 1/2 prevents double-counting since W is symmetric.
	The bias term b represents external magnetic fields.

	Example:
	>>> network = HopfieldNetwork(3)
	>>> state = np.array([1, -1, 1])
	>>> energy_val = network.energy(state)
	>>> print(f"Energy: {energy_val}")
	"""
	return -0.5 * ξ @ self.W @ ξ + np.sum(ξ * self.b)

	def sync_update_step(self, ξ: np.array) -> Tuple[np.ndarray, float]:
	"""
	Perform one synchronous update step on the network state.

	In synchronous updating, all neurons are updated simultaneously using
	the current state. Each neuron's new value is determined by:
	ξ_i^(t+1) = sign(Σ_j W_ij * ξ_j^(t) - b_i)

	This update rule can lead to oscillations between states and doesn't
	guarantee energy decrease at each step, unlike asynchronous updates.

	Args:
	ξ (np.array): Current network state as 1D binary array of shape (N,).
	Values should be in {-1, +1}.

	Returns:
	tuple: (ξ_new, E_new) where:
	- ξ_new (np.array): Updated network state after synchronous step
	- E_new (float): Energy of the new state

	Note:
	Synchronous updates may cause the network to oscillate between
	configurations rather than converge to a fixed point.

	Example:
	>>> network = HopfieldNetwork(4)
	>>> state = np.array([1, -1, 1, -1])
	>>> new_state, new_energy = network.sync_update_step(state)
	"""
	ξ = np.sign(self.W @ ξ - self.b)
	E = self.energy(ξ)
	return ξ, E

	def async_update_step(self, ξ: np.array) -> Tuple[np.ndarray, float]:
	"""
	Perform one asynchronous update step on the network state.

	In asynchronous updating, neurons are updated one at a time in random order.
	Each neuron's new value is computed using the most current states of all
	other neurons. This guarantees energy decrease (or no change) at each step,
	ensuring convergence to a local minimum.

	The update rule for neuron i is:
	ξ_i^(new) = sign(Σ_j W_ij * ξ_j^(current) - b_i)

	Args:
	ξ (np.array): Current network state as 1D binary array of shape (N,).
	Values should be in {-1, +1}. This array is modified in-place.

	Returns:
	tuple: (ξ, E) where:
	- ξ (np.array): The same input array after in-place modification
	- E (float): Energy of the updated state

	Note:
	This method modifies the input array in-place for efficiency.
	The random permutation ensures fair updating and helps avoid cycles.

	Example:
	>>> network = HopfieldNetwork(4)
	>>> state = np.array([1, -1, 1, -1])
	>>> updated_state, energy = network.async_update_step(state)
	>>> # state has been modified in-place
	"""
	for idx in np.random.permutation(self.N):
	ξ[idx] = np.sign(self.W[idx] @ ξ - self.b[idx])
	E = self.energy(ξ)
	return ξ, E

	def run(self,
	ξ: np.array,
	tol: float = 1e-4,
	max_iter: int = 10,
	sync: bool = True) -> List[Tuple[np.ndarray, float]]:
	"""
	Run the network dynamics until convergence or maximum iterations.

	This method iteratively applies update steps (either synchronous or
	asynchronous) starting from an initial state until the network converges
	to a stable configuration or the maximum number of iterations is reached.

	Convergence is detected by monitoring energy changes between consecutive
	steps. The network has converged when the energy change falls below
	the specified tolerance.

	Args:
	ξ (np.array): Initial network state as 1D binary array of shape (N,).
	Values should be in {-1, +1}.
	tol (float, optional): Energy change threshold for convergence detection.
	Defaults to 1e-4.
	max_iter (int, optional): Maximum number of update iterations.
	Defaults to 10.
	sync (bool, optional): If True, use synchronous updates; if False,
	use asynchronous updates. Defaults to True.

	Returns:
	list: History of (state, energy) tuples tracking the network's evolution.
	Each tuple contains:
	- state (np.array): Network configuration at that step
	- energy (float): Energy value for that configuration

	Note:
	- Asynchronous updates (sync=False) are generally preferred as they
	guarantee convergence to local minima
	- Synchronous updates may oscillate and not converge
	- The initial state and energy are included as the first history entry

	Example:
	>>> network = HopfieldNetwork(4)
	>>> initial_state = np.random.choice([-1, 1], 4)
	>>> history = network.run(initial_state, max_iter=50, sync=False)
	>>> final_state, final_energy = history[-1]
	>>> print(f"Converged after {len(history)-1} steps")
	"""
	E = self.energy(ξ)
	history = [(ξ, E)]
	for _ in range(max_iter):
	ξ_new, E_new = self.sync_update_step(ξ) if sync \
	else self.async_update_step(ξ)
	history.append((ξ_new, E_new))
	if abs(E_new - E) < tol:
	break
	ξ, E = ξ_new, E_new
	return history


	def preprocess_image(image_path: str, size: tuple = (64, 64)) -> tuple:
	"""
	Load and preprocess an image for use with Hopfield network.

	This function performs the complete preprocessing pipeline:
	1. Load image and convert to grayscale
	2. Resize to specified dimensions
	3. Apply mean thresholding
	4. Convert to binary values {-1, +1}
	5. Flatten to 1D array

	Args:
	image_path (str): Path to the image file to process.
	size (tuple, optional): Target size as (width, height). Defaults to (64, 64).

	Returns:
	tuple: (pattern, original_shape) where:
	- pattern (np.ndarray): Flattened binary pattern of shape (width*height,)
	- original_shape (tuple): Original 2D shape (width, height) for reshaping

	Raises:
	IOError: If the image file cannot be loaded.
	ValueError: If the image cannot be processed.

	Example:
	>>> pattern, shape = preprocess_image("image.jpg", size=(32, 32))
	>>> print(f"Pattern shape: {pattern.shape}, Original: {shape}")
	"""
	try:
	# Load image: original -> grayscale -> resizing
	pic = Image.open(image_path).convert('L')
	pic = pic.resize(size)
	pix = np.array(pic)

	# Preprocess: mean thresholding -> sign -> flatten
	width, height = pix.shape
	pix = pix.astype(np.float32) # Ensure proper dtype for mean calculation
	pix -= np.mean(pix)
	pix = np.sign(pix).astype(np.int8)
	pix = pix.flatten()

	return pix, (width, height)

	except IOError as e:
	raise IOError(f"Could not load image {image_path}: {e}") from e
	except Exception as e:
	raise ValueError(f"Error processing image {image_path}: {e}") from e


	def load_multiple_images(img_paths: list, size: tuple = (64, 64)) -> tuple:
	"""
	Load and preprocess multiple images for training a Hopfield network.

	All images are resized to the same dimensions to ensure consistent
	pattern sizes for the network. Invalid images are skipped with warnings.

	Args:
	img_paths (list): List of file paths to image files.
	size (tuple, optional): Target size as (width, height). Defaults to (64, 64).

	Returns:
	tuple: (patterns, shape, valid_paths) where:
	- patterns (list): List of preprocessed binary patterns
	- shape (tuple): Common shape (width, height) for all patterns
	- valid_paths (list): List of successfully loaded image paths

	Example:
	>>> paths = ["img1.jpg", "img2.png", "img3.gif"]
	>>> patterns, shape, valid = load_multiple_images(paths)
	>>> print(f"Loaded {len(patterns)} images with shape {shape}")
	"""
	loaded_patterns = []
	loaded_paths = []
	common_shape = None

	for file_path in img_paths:
	try:
	loaded_pattern, img_shape = preprocess_image(file_path, size)
	loaded_patterns.append(loaded_pattern)
	loaded_paths.append(file_path)
	if common_shape is None:
	common_shape = img_shape

	except (IOError, ValueError) as e:
	print(f"Warning: Skipping {file_path} - {e}")
	continue

	if not loaded_patterns:
	raise ValueError("No valid images could be loaded")

	return loaded_patterns, common_shape, loaded_paths


	def demonstrate_recall(network: HopfieldNetwork,
	img_shape: tuple,
	num_tests: int = 3,
	noise_level: float = 0.2,
	max_iter: int = 20) -> None:
	"""
	Demonstrate pattern recall from noisy or random initial states.

	This function tests the network's ability to recover stored patterns
	from corrupted inputs, which is the key capability of associative memory.

	Args:
	network (HopfieldNetwork): Trained Hopfield network.
	img_shape (tuple): Image dimensions (width, height) for visualization.
	num_tests (int, optional): Number of recall tests to perform. Defaults to 3.
	noise_level (float, optional): Fraction of pixels to flip for noise test.
	Defaults to 0.2 (20% corruption).
	max_iter (int, optional): Maximum iterations for network dynamics. Defaults to 20.

	Returns:
	None: Displays recall results using matplotlib.

	Example:
	>>> demonstrate_recall(network, (64, 64), num_tests=5)
	"""
	width, height = img_shape
	total_pixels = width * height

	print(f"\nDemonstrating pattern recall with {num_tests} tests...")

	for test_idx in range(num_tests):
	print(f"\n--- Recall Test #{test_idx + 1} ---")

	# Generate initial state based on test type
	initial_state = _generate_test_state(test_idx, total_pixels, noise_level)

	# Run the network dynamics
	history = network.run(initial_state.copy(), max_iter=max_iter, sync=False)

	# Display results
	_display_test_results(initial_state, history, img_shape, test_idx)


	def _generate_test_state(test_idx: int, total_pixels: int, noise_level: float) -> np.ndarray:
	"""Generate initial state for recall test."""
	if test_idx == 0:
	# Test 1: Random noise
	print("Starting from random noise...")
	return np.random.choice([-1, 1], total_pixels)

	# Test 2+: Add noise to a structured pattern
	print(f"Starting from corrupted pattern (noise level: {noise_level:.1%})...")
	# Add some structure by flipping only a fraction of bits
	mask = np.random.random(total_pixels) < noise_level
	reference = np.random.choice([-1, 1], total_pixels)
	return np.where(mask, -reference, reference)


	def _display_test_results(initial_state: np.ndarray,
	history: List[Tuple[np.ndarray, float]],
	img_shape: tuple,
	test_idx: int) -> None:
	"""
	Display the visual results of a pattern recall test.

	This helper function visualizes both the initial state and final converged state
	of a Hopfield network recall test. It shows the network's ability to recover
	stored patterns from noisy or random inputs by displaying before-and-after
	images with energy and convergence information.

	The function creates two matplotlib plots:
	1. Initial state: The starting configuration (noisy/random)
	2. Final state: The converged result with energy and step count

	Args:
	initial_state (np.ndarray): Starting network state as 1D binary array.
	Values should be in {-1, +1}.
	history (List[Tuple[np.ndarray, float]]): Complete evolution history from
	network.run(), where each tuple
	contains (state, energy).
	img_shape (tuple): Image dimensions (width, height) for reshaping the
	flattened states into 2D images for visualization.
	test_idx (int): Zero-based test index for labeling the displayed results.

	Returns:
	None: Displays matplotlib plots and prints convergence statistics.

	Side Effects:
	- Creates and displays two matplotlib figure windows
	- Prints convergence information to console
	- Blocks execution until plot windows are closed (if running interactively)

	Example:
	>>> initial = np.random.choice([-1, 1], 64*64)
	>>> history = network.run(initial.copy(), max_iter=20)
	>>> _display_test_results(initial, history, (64, 64), 0)
	# Shows: "Test #1: Initial State" and "Test #1: Final State"
	"""
	initial_image = initial_state.reshape(img_shape)
	final_state, final_energy = history[-1]
	final_image = final_state.reshape(img_shape)

	display_state(initial_image, f"Test #{test_idx + 1}: Initial State")
	display_state(final_image,
	f"Test #{test_idx + 1}: Final State (Energy: {final_energy:.2f}, Steps: {len(history)-1})")

	print(f"Converged after {len(history)-1} steps with energy {final_energy:.2f}")


	if __name__ == '__main__':
	# Parse command line arguments
	if len(sys.argv) < 2:
	print("Usage: python hopfield.py <image1> [image2] [image3] ...")
	print("Example: python hopfield.py lenna.png face1.jpg face2.png")
	sys.exit(1)

	image_paths = sys.argv[1:]
	print(f"Loading {len(image_paths)} image(s)...")

	# Load and preprocess multiple images
	try:
	patterns, shape, valid_paths = load_multiple_images(image_paths, size=(64, 64))
	print(f"Successfully loaded {len(patterns)} images with dimensions {shape}")
	except ValueError as e:
	print(f"Error: {e}")
	sys.exit(1)

	# Display original images
	print("\nDisplaying original images...")
	for i, (train_pattern, img_path) in enumerate(zip(patterns, valid_paths)):
	image_name = img_path.split('/')[-1] # Get filename only
	display_state(train_pattern.reshape(shape), f"Original Image #{i+1}: {image_name}")

	# Initialize the network
	network_size = len(patterns[0])
	print(f"\nInitializing Hopfield network with {network_size} neurons...")
	hn = HopfieldNetwork(network_size)

	# Train on all images
	print(f"Training network on {len(patterns)} pattern(s)...")
	hn.train(patterns)

	# Display weight matrix
	print("Displaying learned weight matrix...")
	hn.plot_weights()

	# Test pattern recall
	demonstrate_recall(hn, shape, num_tests=min(3, len(patterns) + 1))

	print(f"\nTraining complete! Network learned {len(patterns)} pattern(s).")
	print(f"Network capacity (rule of thumb): ~{network_size // (4 * len(patterns))} patterns for reliable recall")

	if len(patterns) > network_size // (4 * len(patterns)):
	print("Warning: You may have exceeded the network's storage capacity.")
	print("Consider using fewer patterns or larger images for better recall performance.")
No results found