Of course. Here is the complete, single-file Python script for the Autaxys MVU-1 Simulation Engine, designed to meet all specified requirements for correctness, structure, and output.
```python
# -*- coding: utf-8 -*-
"""
Autaxys MVU-1 Simulation Engine
-------------------------------
A Python script to simulate the "Minimal Viable Universe 1" (MVU-1)
based on the Autaxys framework. This script implements the core concepts,
rewrite rules, and the Autaxic Lagrangian to empirically test the hypothesis
of self-organization driven by a local optimization principle.
Author: AI Assistant
Date: 2023-10-27
"""
import networkx as nx
import pandas as pd
import matplotlib.pyplot as plt
import random
from tqdm import tqdm
from itertools import combinations
# --- Type Hinting for clarity ---
from typing import List, Dict, Tuple, Any, Optional
# --- Main Simulator Class ---
class MVU1_Simulator:
"""
Encapsulates the entire logic for the MVU-1 simulation.
This class manages the universe's state (a graph), applies cosmic
rewrite rules based on an optimization principle (the Autaxic Lagrangian),
and logs the entire evolution for scientific analysis.
"""
def __init__(self, config: Dict[str, Any]):
"""
Initializes the simulator with a configuration dictionary.
Args:
config: A dictionary containing simulation parameters.
Expected keys: 'initial_distinctions', 'max_steps'.
"""
self.config = config
self.graph: nx.Graph = nx.Graph()
self.log: List[Dict[str, Any]] = []
self._step_counter = 0
self._rule_priority = {'bonding': 3, 'collapse': 2, 'annihilation': 1}
def _initialize_graph(self):
"""
Creates the initial "primordial soup" graph based on config.
- Adds N positive and N negative disconnected nodes.
- Node attributes must include 'polarity' (+1 or -1).
"""
n_initial = self.config['initial_distinctions']
for i in range(n_initial):
self.graph.add_node(f'p_{i}', polarity=1)
self.graph.add_node(f'n_{i}', polarity=-1)
print("Graph initialized with primordial soup.")
def _count_bonded_pairs(self, g: nx.Graph) -> int:
"""Helper function to count bonded +1/-1 pairs for the Lagrangian."""
count = 0
for u, v in g.edges():
if g.nodes[u]['polarity'] != g.nodes[v]['polarity']:
count += 1
return count
def _calculate_L_A(self, g: nx.Graph) -> float:
"""
Calculates the Autaxic Lagrangian for a given graph state 'g'.
L_A(G) = (Number of bonded +1/-1 pairs in G) / (|D| + |R|)
Args:
g: The networkx.Graph object to analyze.
Returns:
The float value of the Autaxic Lagrangian.
"""
num_distinctions = g.number_of_nodes()
num_relations = g.number_of_edges()
denominator = num_distinctions + num_relations
if denominator == 0:
return 0.0
numerator = self._count_bonded_pairs(g)
return numerator / denominator
def _find_possible_moves(self) -> List[Tuple[str, Any]]:
"""
Scans the current graph (self.graph) and identifies ALL possible
applications of the three rewrite rules (r1, r2, r3).
Returns:
A list of tuples, where each tuple is:
('rule_name', *args_for_rule_application)
e.g., ('annihilation', node1, node2), ('bonding', node1, node2)
"""
moves = []
# Rule r1: Annihilation of a stable, bonded pair
for u, v in self.graph.edges():
if (self.graph.nodes[u]['polarity'] != self.graph.nodes[v]['polarity'] and
self.graph.degree[u] == 1 and self.graph.degree[v] == 1):
moves.append(('annihilation', u, v))
# Rule r2: Bonding of disconnected opposite polarities
pos_nodes = [n for n, attr in self.graph.nodes(data=True) if attr['polarity'] == 1]
neg_nodes = [n for n, attr in self.graph.nodes(data=True) if attr['polarity'] == -1]
for p_node in pos_nodes:
for n_node in neg_nodes:
if not self.graph.has_edge(p_node, n_node):
moves.append(('bonding', p_node, n_node))
# Rule r3: Redundancy Collapse (K3 -> K2)
# Find all 3-cliques (triangles)
found_triangles = set()
for node in self.graph:
# combinations('ABCD', 2) --> AB AC AD BC BD CD
for neighbor1, neighbor2 in combinations(self.graph.neighbors(node), 2):
if self.graph.has_edge(neighbor1, neighbor2):
# Found a triangle
triangle = tuple(sorted((node, neighbor1, neighbor2)))
if triangle not in found_triangles:
found_triangles.add(triangle)
# Identify nodes to keep based on lowest degree in the whole graph
degrees = [(n, self.graph.degree[n]) for n in triangle]
degrees.sort(key=lambda x: x[1]) # Sort by degree
nodes_to_keep = (degrees[0][0], degrees[1][0])
moves.append(('collapse', triangle, nodes_to_keep))
return moves
def step(self) -> bool:
"""
Executes a single step of the simulation by selecting and applying the best move.
Returns:
False if the simulation has halted (no moves possible), True otherwise.
"""
possible_moves = self._find_possible_moves()
if not possible_moves:
print("\nNo more possible moves. Simulation halted.")
return False
evaluated_moves = []
for move in possible_moves:
temp_graph = self.graph.copy()
rule_name = move[0]
if rule_name == 'annihilation':
temp_graph.remove_nodes_from(move[1:])
elif rule_name == 'bonding':
temp_graph.add_edge(move[1], move[2])
elif rule_name == 'collapse':
triangle, nodes_to_keep = move[1], move[2]
node_to_remove = list(set(triangle) - set(nodes_to_keep))[0]
temp_graph.remove_node(node_to_remove)
# Ensure the two kept nodes form a bond
if not temp_graph.has_edge(nodes_to_keep[0], nodes_to_keep[1]):
temp_graph.add_edge(nodes_to_keep[0], nodes_to_keep[1])
la_after_move = self._calculate_L_A(temp_graph)
evaluated_moves.append({'move': move, 'la': la_after_move})
# --- Decision Making: Find the best move ---
# Sort by L_A (desc), then by rule priority (desc)
evaluated_moves.sort(
key=lambda x: (x['la'], self._rule_priority.get(x['move'][0], 0)),
reverse=True
)
# Tie-breaking for moves with identical max L_A and priority
best_la = evaluated_moves[0]['la']
best_priority = self._rule_priority.get(evaluated_moves[0]['move'][0], 0)
top_tier_moves = [
m for m in evaluated_moves
if m['la'] == best_la and self._rule_priority.get(m['move'][0], 0) == best_priority
]
# Choose one randomly from the top tier
chosen_move_eval = random.choice(top_tier_moves)
best_move = chosen_move_eval['move']
rule_applied = best_move[0]
# --- Apply the chosen move to the actual graph ---
if rule_applied == 'annihilation':
self.graph.remove_nodes_from(best_move[1:])
elif rule_applied == 'bonding':
self.graph.add_edge(best_move[1], best_move[2])
elif rule_applied == 'collapse':
triangle, nodes_to_keep = best_move[1], best_move[2]
node_to_remove = list(set(triangle) - set(nodes_to_keep))[0]
self.graph.remove_node(node_to_remove)
if not self.graph.has_edge(nodes_to_keep[0], nodes_to_keep[1]):
self.graph.add_edge(nodes_to_keep[0], nodes_to_keep[1])
# --- Log the state AFTER the move ---
self._step_counter += 1
self._log_state(self._step_counter, rule_applied)
return True
def _log_state(self, step_num: int, rule_applied: str):
"""
Calculates all relevant metrics for the current state and appends them to the log.
"""
num_d = self.graph.number_of_nodes()
num_r = self.graph.number_of_edges()
bonded_pairs = self._count_bonded_pairs(self.graph)
current_la = self._calculate_L_A(self.graph)
log_entry = {
'step': step_num,
'L_A': current_la,
'num_distinctions': num_d,
'num_relations': num_r,
'graph_size': num_d + num_r,
'num_bonded_pairs': bonded_pairs,
'rule_applied': rule_applied,
}
self.log.append(log_entry)
def run(self):
"""
The main simulation loop.
- Initializes the graph and logs the initial state.
- Loops for `max_steps` or until the simulation halts.
"""
self._initialize_graph()
# Log the initial state (step 0)
self._log_state(step_num=0, rule_applied='initial_state')
# Use tqdm for a progress bar
with tqdm(total=self.config['max_steps'], desc="Simulating Universe") as pbar:
for i in range(self.config['max_steps']):
if not self.step():
pbar.n = self.config['max_steps'] # Mark progress bar as complete
pbar.update(0)
break
# Update progress bar description with current stats
last_log = self.log[-1]
pbar.set_description(
f"Step {last_log['step']}: "
f"L_A={last_log['L_A']:.3f}, "
f"|D|={last_log['num_distinctions']}, "
f"|R|={last_log['num_relations']}, "
f"Pairs={last_log['num_bonded_pairs']}"
)
pbar.update(1)
def get_results_df(self) -> pd.DataFrame:
"""
Converts the final log into a pandas DataFrame.
"""
return pd.DataFrame(self.log)
def plot_results(self, results_df: pd.DataFrame, output_prefix: str):
"""
Generates and saves key visualizations from the results.
"""
print(f"Generating plots for {output_prefix}...")
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
fig.suptitle(f'MVU-1 Simulation Results: {output_prefix}', fontsize=16)
# 1. Autaxic Lagrangian vs. Step
axes[0, 0].plot(results_df['step'], results_df['L_A'], marker='.', linestyle='-', color='b')
axes[0, 0].set_title('Autaxic Lagrangian (L_A) vs. Step')
axes[0, 0].set_xlabel('Simulation Step')
axes[0, 0].set_ylabel('L_A Value')
axes[0, 0].grid(True)
# 2. Number of Bonded Pairs vs. Step
axes[0, 1].plot(results_df['step'], results_df['num_bonded_pairs'], marker='.', linestyle='-', color='g')
axes[0, 1].set_title('Number of Bonded Pairs (P_ID) vs. Step')
axes[0, 1].set_xlabel('Simulation Step')
axes[0, 1].set_ylabel('Count of Bonded Pairs')
axes[0, 1].grid(True)
# 3. Graph Size vs. Step
axes[1, 0].plot(results_df['step'], results_df['graph_size'], marker='.', linestyle='-', color='purple')
axes[1, 0].set_title('Total Graph Size (|D| + |R|) vs. Step')
axes[1, 0].set_xlabel('Simulation Step')
axes[1, 0].set_ylabel('Graph Size')
axes[1, 0].grid(True)
# 4. Final Graph State
ax = axes[1, 1]
ax.set_title(f'Final Graph State (Step {self.log[-1]["step"]})')
if not self.graph.nodes():
ax.text(0.5, 0.5, "Empty Graph (Annihilation Complete)", ha='center', va='center')
ax.set_xticks([])
ax.set_yticks([])
else:
pos = nx.spring_layout(self.graph, seed=42)
node_colors = [
'red' if self.graph.nodes[n]['polarity'] == -1 else 'blue'
for n in self.graph.nodes()
]
nx.draw(
self.graph,
pos,
ax=ax,
with_labels=False,
node_color=node_colors,
node_size=50,
width=1.5,
edge_color='gray'
)
plt.tight_layout(rect=[0, 0.03, 1, 0.95])
plot_path = f"{output_prefix}_results_plot.png"
plt.savefig(plot_path)
print(f"Plots saved to {plot_path}")
plt.close(fig)
# --- Example Usage (Main Execution Block) ---
if __name__ == "__main__":
# --- Simulation Configuration ---
sim_config = {
"initial_distinctions": 50, # Means 50 positive and 50 negative
"max_steps": 1000,
"output_prefix": "mvu1_run_01"
}
# --- Execution ---
print("--- Starting Autaxys MVU-1 Simulation ---")
print(f"Configuration: {sim_config}")
print("-" * 40)
simulator = MVU1_Simulator(config=sim_config)
simulator.run()
print("-" * 40)
print("Simulation complete.")
# --- Data & Plotting ---
final_state = simulator.log[-1]
print("\n--- Final State Summary ---")
print(f"L_A: {final_state['L_A']:.4f}")
print(f"Distinctions: {final_state['num_distinctions']}")
print(f"Relations: {final_state['num_relations']}")
print(f"Bonded Pairs: {final_state['num_bonded_pairs']}")
print("-" * 40)
results = simulator.get_results_df()
csv_path = f"{sim_config['output_prefix']}_simulation_log.csv"
results.to_csv(csv_path, index=False)
print(f"Full results log saved to: {csv_path}")
simulator.plot_results(results_df=results, output_prefix=sim_config['output_prefix'])
print("\n--- Script Finished ---")
```