crystal.py

import numpy as np
from numpy.random import randint
from typing import List, Tuple


G_GRID_SIZE = 100
G_INIT_CRYSTAL_COORDS = (50, 50)
G_NUM_PARTICLES = 10

G_PARTICLE_INCREMENTS = range(10, 31, 10)
G_SIMULATION_REPEATS = 2


class BindingBox:
    def __init__(self, left: int, right: int, up: int, down: int):
        """
        NOTE: We assume that box coordinates go from top to bottom and from
              left to right as they grow!
        """
        self.left = left
        self.right = right
        self.up = up
        self.down = down

        self.height = down - up + 1
        self.width = right - left + 1
    
    def __repr__(self) -> str:
        return "BindingBox(L:{}, R:{}, U:{}, D:{}, H:{}, W:{}, A:{})".format(
            self.left, 
            self.right,
            self.up,
            self.down, 
            self.height, 
            self.width, 
            self.area())

    def area(self) -> int:
        return self.height * self.width


class ParticlePool:
    def __init__(self, capacity: int):
        self.coords = np.zeros((capacity, 2), dtype=int)
        self.particles = [Particle(self.coords[offset]) for offset in range(capacity)]
    
    def __len__(self) -> int:
        return len(self.particles)

    def density(self) -> float:
        return float(len(self)) / self.area()

    def area(self) -> int:
        return self.binding_box().area()

    def binding_box(self) -> BindingBox:
        min_x, min_y = np.amin(self.coords, axis=0)
        max_x, max_y = np.amax(self.coords, axis=0)
        return BindingBox(min_x, max_x, min_y, max_y)


class Particle:
    def __init__(self, coords):
        self.__coords = coords

    def __repr__(self) -> str:
        return "Particle({},{})".format(self.x, self.y)

    def coords(self, coords: Tuple[int, int]):
        self.x = coords[0]
        self.y = coords[1]
    
    @property
    def x(self):
        return self.__coords[0]
    
    @property
    def y(self):
        return self.__coords[1]
    
    @x.setter
    def x(self, value: int):
        self.__coords[0] = value
    
    @y.setter
    def y(self, value: int):
        self.__coords[1] = value


class Crystal:
    def __init__(self, init_particle: Particle):
        self.components = [init_particle]
    
    def attach(self, p: Particle):
        self.components.append(p)
        print("Crystal: attached {}".format(p))
    

class Lattice:
    def __init__(self, grid_size: int, c: Crystal, crystal_coords: Tuple[int, int]):
        self.grid_size = grid_size
        self.grid = np.zeros((grid_size, grid_size), dtype=bool)  # assuming FALSE by default
        self.spawn_crystal(c, crystal_coords)
    
    def has_particle(self, x: int, y: int) -> bool:
        return self.grid[y][x]

    def has_adjacent_particle(self, p: Particle) -> bool:
        move_options = self.get_move_options(p)
        adjacent_locations = [(p.x + dx, p.y + dy) for dx, dy in move_options]
        return any({self.has_particle(x, y) for x, y in adjacent_locations})

    def spawn_crystal(self, c: Crystal, spawn_coords: Tuple[int, int]):
        for particle in c.components:
            self.spawn_particle(particle, spawn_coords)

    def spawn_particle(self, p: Particle, spawn_coords: Tuple[int, int]):
        p.coords(spawn_coords)
        if self.has_particle(p.x, p.y):
            raise Exception("Tried to sprawn particle on top of another at ({},{})".format(p.x, p.y))

        self.grid[p.y][p.x] = True

    def remove_particle(self, p: Particle):
        if not self.has_particle(p.x, p.y):
            raise Exception("Tried to remove particle from where it isn't present ({},{})".format(p.x, p.y))

        self.grid[p.y][p.x] = False
    
    def move_particle(self, p: Particle, coord_delta: Tuple[int, int]):
        self.remove_particle(p)
        new_coords = (p.x + coord_delta[0], p.y + coord_delta[1])
        self.spawn_particle(p, new_coords)

    def randomly_move(self, p: Particle):
        self.move_particle(p, self.random_move_option(p))

    def get_move_options(self, p: Particle) -> List[Tuple[int, int]]:
        """
        Determine where the given particle can move on the lattice.

        NOTE: IF THERE IS AN ADJACENT PARTICLE NEARBY, IT IS STILL CONSIDERED
              AS A VIABLE MOVE OPTION!
        """
        options = []
        if p.x > 0:
            options.append((-1, 0))

        if p.x < self.grid_size - 1:
            options.append((1, 0))

        if p.y > 0:
            options.append((0, -1))

        if p.y < self.grid_size - 1:
            options.append((0, 1))

        return options

    def random_move_option(self, p: Particle) -> Tuple[int, int]:
        options = self.get_move_options(p)
        return options[randint(0, len(options))]


class SimulationResult:
    def __init__(self, grid_size: int, pool: ParticlePool, crystal: Crystal, lattice: Lattice):
        self.grid_size = grid_size
        self.pool = pool
        self.crystal = crystal
        self.lattice = lattice

    @property
    def density(self) -> float:
        return self.pool.density()

    @property
    def num_particles(self) -> int:
        return len(self.pool)

    def report(self):
        print("--- SIMULATION REPORT ---------------------------")
        print("Crystal particles:")
        for particle in self.crystal.components:
            print("  {},{}".format(particle.x, particle.y))
        print()


class Simulator:
    def simulate(self, grid_size: int, num_particles: int, init_coords: Tuple[int, int]) -> SimulationResult:
        pool = ParticlePool(num_particles)
        p0 = pool.particles[0]
        other_particles = pool.particles[1:]  # skip the one we spawned a crystal from

        crystal = Crystal(p0)
        lattice = Lattice(grid_size, Crystal(p0), init_coords)

        def random_corner() -> int:
            return randint(0, 2) * (grid_size - 1)

        def random_spawn_coords() -> Tuple[int, int]:
            vertical_corner = random_corner()
            horizontal_corner = random_corner()
            return (horizontal_corner, vertical_corner)

        for p in other_particles:
            spawn_coords = random_spawn_coords()
            lattice.spawn_particle(p, spawn_coords)
            while not lattice.has_adjacent_particle(p):
                lattice.randomly_move(p)
            
            crystal.attach(p)
        
        return SimulationResult(grid_size, pool, crystal, lattice)


def main():
    grid_size = G_GRID_SIZE
    init_coords = G_INIT_CRYSTAL_COORDS
    num_particles = G_NUM_PARTICLES

    sim = Simulator()
    master_results: List[Tuple[int, List[float]]] = []
    for num_particles in G_PARTICLE_INCREMENTS:
        densities = [sim.simulate(grid_size, num_particles, init_coords).density for _ in range(G_SIMULATION_REPEATS)]
        master_results.append((num_particles, densities))
        print(" * Densities for {} particles ready!".format(num_particles))

    for num_particles, densities in master_results:
        print("{} particles:".format(num_particles))
        print("  > mean:  {}".format(np.mean(densities)))
        print("  > stdev: {}".format(np.std(densities, ddof=1)))

    # TODO: THIS vvv
    # Suppose you save the densities as density_list[particle_count,repeat_count]
    # Then you get dens_mean=np.mean(density_list,axis=1)
    # dens_SD=np.std(density_list,axis=1,ddof=1)

    return 0


if __name__ == "__main__":
    raise SystemExit(main())