grid_size = 100 num_centers = 15 center_strength_base = 10 trade = 0

1. --- Voimalaskenta ---

def strength_of_center(center_strength, distance, Trade):

   Roo0 = 0.1
   Koo = 0.7
   fadrate = Roo0 * (1 - Koo * Trade)
   return center_strength / np.power((1 + fadrate), distance)

1. --- Kustannuspinta (Perlin noise) ---

def generate_cost_surface(grid_size, scale=20.0, octaves=6, persistence=0.5, lacunarity=2.0):

   cost_surface = np.zeros((grid_size, grid_size))
   for y in range(grid_size):
       for x in range(grid_size):
           nx = x / scale
           ny = y / scale
           cost_surface[y][x] = pnoise2(nx, ny, octaves=octaves, persistence=persistence, lacunarity=lacunarity)
   cost_surface = (cost_surface - cost_surface.min()) / (cost_surface.max() - cost_surface.min())
   return 1 + 10 * cost_surface  # [1, 10]

cost_surface = generate_cost_surface(grid_size) plt.figure(figsize=(6,5)) plt.imshow(cost_surface, cmap='terrain') plt.title("Kustannuspinta (Perlin noise)") plt.colorbar(label="Kulkukustannus") plt.show()

1. --- Keskusten alustus ---

center_positions = np.random.randint(0, grid_size, size=(num_centers, 2)) center_strengths = center_strength_base / (np.arange(1, num_centers + 1))

1. --- Matriisit ---

total_strength = np.zeros((grid_size, grid_size)) dominant_center = np.zeros((grid_size, grid_size), dtype=int)

def cost_distance(cost_surface, start, end):

   path, cost = route_through_array(cost_surface, start=start[::-1], end=end[::-1], fully_connected=True)
   return cost

def calculate_total_strength():

   total_strength.fill(0)
   dominant_center.fill(0)
   X, Y = np.meshgrid(np.arange(grid_size), np.arange(grid_size))
   for i in range(num_centers):
       center_x, center_y = center_positions[i]
       center_strength = center_strengths[i]
       distance = np.sqrt((X - center_x) ** 2 + (Y - center_y) ** 2)
       strength_contribution = strength_of_center(center_strength, distance, trade)
       mask = strength_contribution > total_strength
       total_strength[mask] = strength_contribution[mask]
       dominant_center[mask] = i

def calculate_total_strength_cost_based():

   global dominant_center

   dominant_center.fill(0)
   cost_map = cost_surface.copy()

   # Lähteet kaikille keskuksille
   sources = [tuple(p[::-1]) for p in center_positions]  # (y, x)

   # Luo MCP-objekti ja laske kustannusetäisyydet kaikille lähteille
   mcp = MCP_Geometric(cost_map)
   costs, traceback = mcp.find_costs(starts=sources)

   # Laske voimat
   all_strengths = np.zeros((grid_size, grid_size, num_centers))
   for i in range(num_centers):
       strength = strength_of_center(center_strengths[i], costs[i], trade)
       all_strengths[:, :, i] = strength

   # Valitaan maksimaalinen vaikutus
   dominant_center = np.argmax(all_strengths, axis=2)

1. --- Animaation päivitys ---

def update(frame):

   global num_centers, center_positions, center_strengths

   if num_centers <= 1:
       return

   max_strength = -np.inf
   center1 = center2 = -1
   merge_line = None

   for i in range(num_centers):
       for j in range(i + 1, num_centers):
           start = tuple(center_positions[i][::-1])
           end = tuple(center_positions[j][::-1])
           distance = cost_distance(cost_surface, start, end)
           s_i = strength_of_center(center_strengths[i], distance, trade)
           s_j = strength_of_center(center_strengths[j], distance, trade)

           dominant = i if s_i > s_j else j
           current_strength = s_i if dominant == i else s_j

           if current_strength > max_strength:
               max_strength = current_strength
               center1, center2 = dominant, j if dominant == i else i
               merge_line = (center_positions[center1], center_positions[center2])

   if center1 != -1 and center2 != -1:
       center_strengths[center1] += center_strengths[center2]
       center_positions = np.delete(center_positions, center2, axis=0)
       center_strengths = np.delete(center_strengths, center2)
       num_centers -= 1

   calculate_total_strength()

   ax.clear()
   ax.imshow(dominant_center, cmap='tab20', origin='lower', extent=[0, grid_size, 0, grid_size])
   ax.scatter(center_positions[:, 0], center_positions[:, 1], color='blue', marker='x', label="Centers")
   ax.set_title(f"Frame {frame}: Dominant Centers")
   ax.set_xlabel("X Position")
   ax.set_ylabel("Y Position")

   if merge_line:
       p1, p2 = merge_line
       ax.plot([p1[0], p2[0]], [p1[1], p2[1]], 'r--', linewidth=2, label='Merging')
       ax.scatter([p2[0]], [p2[1]], color='red', marker='o', s=80, label='Merged')

   ax.legend(loc="upper right")

1. --- Piirto ---

calculate_total_strength_cost_based()

fig, ax = plt.subplots(figsize=(10, 8)) ani = FuncAnimation(fig, update, frames=range(num_centers - 1), repeat=False, interval=2000)

1. Tallennetaan animaatio gif-muotoon

ani.save("center_merging_animation.gif", writer='imagemagick', fps=0.5, dpi=200)

plt.show()