diff --git a/gbso/optimize.py b/gbso/optimize.py
index c077160..de60d82 100644
--- a/gbso/optimize.py
+++ b/gbso/optimize.py
@@ -1,5 +1,25 @@
+from random import random
+from typing import List, Tuple, TypeVar
+
 from gbso.program.program import Program
 
+DEFAULT_NUM_ITERS = 1_000_000
+
+DEFAULT_PROB_OPCODE = 0.25
+DEFAULT_PROB_OPERAND = 0.25
+DEFAULT_PROB_SWAP = 0.25
+DEFAULT_PROB_INSN = 0.25
+
+DEFAULT_PROB_INSN_UNUSED = 0.1
+
 
-def optimize(prgm: Program) -> Program:
+def optimize(
+    prgm: Program,
+    num_iters: int = DEFAULT_NUM_ITERS,
+    prob_opcode: float = DEFAULT_PROB_OPCODE,
+    prob_operand: float = DEFAULT_PROB_OPERAND,
+    prob_swap: float = DEFAULT_PROB_SWAP,
+    prob_insn: float = DEFAULT_PROB_INSN,
+    prob_insn_unused: float = DEFAULT_PROB_INSN_UNUSED,
+) -> Program:
     return prgm
diff --git a/gbso/program/mutate.py b/gbso/program/mutate.py
new file mode 100644
index 0000000..fa52210
--- /dev/null
+++ b/gbso/program/mutate.py
@@ -0,0 +1,57 @@
+from enum import Enum
+from random import random
+from typing import List, Tuple, TypeVar
+
+from gbso.program.program import Program
+
+
+class Mutation(Enum):
+    OPCODE = "OPCODE"
+    OPERAND = "OPERAND"
+    SWAP = "SWAP"
+    INSTRUCTION = "INSTRUCTION"
+
+
+def mutate_program(
+    prgm: Program,
+    prob_opcode: float,
+    prob_operand: float,
+    prob_swap: float,
+    prob_insn: float,
+    prob_insn_unused: float,
+) -> Program:
+    mutation = sample_probs(
+        [
+            (Mutation.OPCODE, prob_opcode),
+            (Mutation.OPERAND, prob_operand),
+            (Mutation.SWAP, prob_swap),
+            (Mutation.INSTRUCTION, prob_insn),
+        ]
+    )
+    new_prgm = prgm
+    if mutation == Mutation.OPCODE:
+        new_prgm = mutate_opcode(prgm)
+    elif mutation == Mutation.OPERAND:
+        new_prgm = mutate_operand(prgm)
+    elif mutation == Mutation.SWAP:
+        new_prgm = mutate_swap(prgm)
+    elif mutation == Mutation.INSTRUCTION:
+        new_prgm = mutate_insn(prgm, prob_insn_unused)
+    return new_prgm
+
+
+A = TypeVar("A")
+
+
+# NOTE: List must be non-empty and probabilities must sum to 1
+def sample_probs(probs: List[Tuple[A, float]]) -> A:
+    r = random()
+    cum_prob = 0.0
+
+    for val, prob in probs[:-1]:
+        cum_prob += prob
+        if r < cum_prob:
+            return val
+
+    val, _ = probs[-1]
+    return val