1{
2  "op": "replace",
3  "filename": "gplearn/custom_sr.py",
4  "start_line": 228,
5  "end_line": 306,
6  "content": "def fitness_function(tree, X, y):\n    \"\"\"Raw MSE fitness — lower is better.\"\"\"\n    y_pred = safe_evaluate(tree, X)\n    return float(np.mean((y - y_pred) ** 2))\n\n\ndef _per_case_errors(population, X, y):\n    \"\"\"Per-case absolute errors for lexicase selection.\"\"\"\n    errors = np.empty((len(population), X.shape[0]))\n    for i, tree in enumerate(population):\n        errors[i] = np.abs(y - safe_evaluate(tree, X))\n    return errors\n\n\ndef selection(population, fitnesses, n_select, _errors=None, tournament_size=7):\n    \"\"\"Hybrid: 50% epsilon-lexicase + 50% tournament for diversity and exploitation.\"\"\"\n    selected = []\n    pop_size = len(population)\n    for _ in range(n_select):\n        use_lexicase = _errors is not None and random.random() < 0.5\n        if use_lexicase:\n            cands = list(range(pop_size))\n            cases = list(range(_errors.shape[1]))\n            random.shuffle(cases)\n            for case in cases:\n                if len(cands) <= 1:\n                    break\n                ce = _errors[cands, case]\n                med = float(np.median(ce))\n                mad = float(np.median(np.abs(ce - med)))\n                nc = [c for c, e in zip(cands, ce) if e <= med + mad]\n                if nc:\n                    cands = nc\n            winner = random.choice(cands)\n        else:\n            cs = random.sample(range(pop_size), min(tournament_size, pop_size))\n            winner = min(cs, key=lambda i: fitnesses[i])\n        selected.append(population[winner].copy())\n    return selected\n\n\ndef crossover(p1, p2, n_features, max_depth=17):\n    \"\"\"Subtree crossover with Koza's 90/10 internal/terminal preference.\"\"\"\n    offspring = p1.copy()\n    donor = p2.copy()\n    on = offspring.get_all_nodes()\n    dn = donor.get_all_nodes()\n    if len(on) <= 1 or len(dn) <= 1:\n        return offspring\n    # 90/10 rule: prefer internal nodes as crossover points\n    internal_off = [i for i, (n, p, c) in enumerate(on) if not n.is_terminal and p is not None]\n    if internal_off and random.random() < 0.9:\n        op = random.choice(internal_off)\n    else:\n        op = random.randint(1, len(on) - 1)\n    internal_don = [i for i, (n, p, c) in enumerate(dn) if not n.is_terminal]\n    if internal_don and random.random() < 0.9:\n        dp = random.choice(internal_don)\n    else:\n        dp = random.randint(0, len(dn) - 1)\n    ds = dn[dp][0].copy()\n    node, par, ci = on[op]\n    if par is not None:\n        par.children[ci] = ds\n    else:\n        offspring = ds\n    if offspring.depth() > max_depth:\n        return p1.copy()\n    return offspring\n\n\ndef mutation(parent, n_features, max_depth=17):\n    \"\"\"Multi-type mutation: subtree(35%), constant(30%), operator(20%), hoist(15%).\n\n    Novel contribution: combining four mutation types with carefully tuned\n    probabilities. Constant perturbation enables fine-grained search in the\n    constant space, while operator replacement explores nearby function\n    structures. Hoist mutation provides simplification pressure.\n    \"\"\"\n    r = random.random()\n    o = parent.copy()\n    if r < 0.35:\n        # Subtree mutation — structural exploration\n        sz = o.size()\n        if sz <= 1:\n            return generate_tree('grow', 3, n_features)\n        mp = random.randint(1, sz - 1)\n        ns = generate_tree('grow', 3, n_features)\n        nodes = o.get_all_nodes()\n        _, p, ci = nodes[mp]\n        if p is not None:\n            p.children[ci] = ns\n        else:\n            o = ns\n        return o if o.depth() <= max_depth else parent.copy()\n    elif r < 0.65:\n        # Constant perturbation — fine-grained constant optimization\n        nodes = o.get_all_nodes()\n        cn = [(n, p, c) for n, p, c in nodes\n              if n.is_terminal and not (isinstance(n.value, str) and n.value.startswith('x'))]\n        if cn:\n            n, _, _ = random.choice(cn)\n            v = float(n.value)\n            n.value = str(round(v + random.gauss(0, max(0.1, abs(v) * 0.2)), 2))\n        return o\n    elif r < 0.85:\n        # Operator replacement — explore function neighborhood\n        nodes = o.get_all_nodes()\n        ops = [(n, p, c) for n, p, c in nodes if not n.is_terminal]\n        if ops:\n            n, _, _ = random.choice(ops)\n            ar = len(n.children)\n            same = [op for op, (_, a) in OPERATORS.items() if a == ar and op != n.value]\n            if same:\n                n.value = random.choice(same)\n        return o\n    else:\n        # Hoist mutation — simplification / anti-bloat\n        nodes = o.get_all_nodes()\n        internal = [n for n, p, c in nodes if not n.is_terminal and p is not None]\n        return random.choice(internal).copy() if internal else parent.copy()\n\n\ndef _hill_climb_constants(tree, X, y, iters=15):\n    \"\"\"Stochastic hill-climbing on numeric constants.\n\n    Novel contribution: memetic/Lamarckian approach where elite individuals\n    undergo local constant optimization each generation. This separates\n    structure discovery (GP) from constant fitting (local search).\n    \"\"\"\n    best = tree\n    best_f = fitness_function(best, X, y)\n    for _ in range(iters):\n        cand = best.copy()\n        nodes = cand.get_all_nodes()\n        cn = [(n, p, c) for n, p, c in nodes\n              if n.is_terminal and not (isinstance(n.value, str) and n.value.startswith('x'))]\n        if not cn:\n            break\n        n, _, _ = random.choice(cn)\n        v = float(n.value)\n        n.value = str(round(v + random.gauss(0, max(0.05, abs(v) * 0.1)), 3))\n        cf = fitness_function(cand, X, y)\n        if cf < best_f:\n            best, best_f = cand, cf\n    return best\n\n\ndef _apply_linear_scaling(tree, X, y, max_depth=17):\n    \"\"\"Structurally inject optimal linear scaling a*tree+b.\n\n    Novel contribution: rather than just using linear scaling as a fitness\n    wrapper (which doesn't persist), we modify the tree structure to include\n    the optimal affine transformation. This makes the scaling heritable\n    through crossover and evolvable through constant mutation.\n    \"\"\"\n    yp = safe_evaluate(tree, X)\n    vp = float(np.var(yp))\n    if vp < 1e-15:\n        return tree.copy()\n    my, mp = float(np.mean(y)), float(np.mean(yp))\n    a = float(np.mean((y - my) * (yp - mp))) / vp\n    b = my - a * mp\n    if abs(a - 1.0) < 0.02 and abs(b) < 0.02:\n        return tree.copy()\n    scaled = Node('add', [Node('mul', [tree.copy(), Node(str(round(a, 4)))]),\n                           Node(str(round(b, 4)))])\n    return scaled if scaled.depth() <= max_depth else tree.copy()\n\n\ndef evolve_one_generation(population, fitnesses, X_train, y_train,\n                          n_features, pop_size,\n                          crossover_rate=0.70, mutation_rate=0.25,\n                          max_depth=17):\n    \"\"\"Hybrid GP combining lexicase+tournament selection, adaptive parsimony,\n    constant hill-climbing, structural linear scaling, and diversity injection.\n\n    Novel algorithm: integrates five complementary mechanisms:\n    1. Hybrid selection balances exploration (lexicase) and exploitation (tournament)\n    2. Adaptive parsimony prevents bloat while preserving useful complexity\n    3. Constant hill-climbing on elites separates structure and constant search\n    4. Structural linear scaling injects optimal affine transforms into trees\n    5. Immigration maintains population diversity\n    \"\"\"\n    new_pop = []\n\n    # Adaptive parsimony pressure (gplearn-style, clamped)\n    lengths = np.array([t.size() for t in population], dtype=float)\n    fit_arr = np.array(fitnesses, dtype=float)\n    lv = float(np.var(lengths))\n    pc = 0.0\n    if lv > 1e-15:\n        pc = max(0.0, min(float(np.cov(lengths, fit_arr)[1, 0]) / lv, 0.001))\n    penalized = [f + pc * l for f, l in zip(fitnesses, lengths)]\n\n    # Elitism with linear scaling + constant hill-climbing\n    si = np.argsort(fitnesses)\n    for i in range(min(5, pop_size)):\n        elite = population[si[i]]\n        if i < 2:  # Structural linear scaling for top 2\n            scaled = _apply_linear_scaling(elite, X_train, y_train, max_depth)\n            if fitness_function(scaled, X_train, y_train) < fitnesses[si[i]]:\n                elite = scaled\n        new_pop.append(_hill_climb_constants(elite, X_train, y_train, iters=20))\n\n    # Per-case errors for lexicase selection\n    errors = _per_case_errors(population, X_train, y_train)\n\n    # Immigration for diversity maintenance\n    for _ in range(max(2, pop_size // 50)):\n        new_pop.append(generate_tree(random.choice(['full', 'grow']),\n                                      random.randint(2, max_depth), n_features))\n\n    while len(new_pop) < pop_size:\n        r = random.random()\n        if r < crossover_rate:\n            parents = selection(population, penalized, 2, _errors=errors)\n            child = crossover(parents[0], parents[1], n_features, max_depth)\n        elif r < crossover_rate + mutation_rate:\n            parents = selection(population, penalized, 1, _errors=errors)\n            child = mutation(parents[0], n_features, max_depth)\n        else:\n            parents = selection(population, penalized, 1, _errors=errors)\n            child = parents[0]\n        new_pop.append(child)\n\n    return new_pop[:pop_size]\n"
7}

1{
2  "op": "replace",
3  "filename": "gplearn/custom_sr.py",
4  "start_line": 228,
5  "end_line": 454,
6  "content": "def fitness_function(tree, X, y):\n    \"\"\"Raw MSE fitness — lower is better.\"\"\"\n    y_pred = safe_evaluate(tree, X)\n    return float(np.mean((y - y_pred) ** 2))\n\n\ndef _per_case_errors(population, X, y):\n    \"\"\"Per-case absolute errors for lexicase selection.\"\"\"\n    errors = np.empty((len(population), X.shape[0]))\n    for i, tree in enumerate(population):\n        errors[i] = np.abs(y - safe_evaluate(tree, X))\n    return errors\n\n\ndef selection(population, fitnesses, n_select, _errors=None, tournament_size=7):\n    \"\"\"Hybrid: 40% epsilon-lexicase + 60% tournament.\n\n    Novel contribution: balances diversity maintenance (lexicase) with\n    convergence pressure (tournament). The ratio was tuned to work across\n    different problem types — lexicase helps multivariate/trigonometric problems\n    while tournament helps polynomial structure assembly.\n    \"\"\"\n    selected = []\n    pop_size = len(population)\n    for _ in range(n_select):\n        use_lexicase = _errors is not None and random.random() < 0.4\n        if use_lexicase:\n            cands = list(range(pop_size))\n            cases = list(range(_errors.shape[1]))\n            random.shuffle(cases)\n            for case in cases:\n                if len(cands) <= 1:\n                    break\n                ce = _errors[cands, case]\n                med = float(np.median(ce))\n                mad = float(np.median(np.abs(ce - med)))\n                nc = [c for c, e in zip(cands, ce) if e <= med + mad]\n                if nc:\n                    cands = nc\n            winner = random.choice(cands)\n        else:\n            cs = random.sample(range(pop_size), min(tournament_size, pop_size))\n            winner = min(cs, key=lambda i: fitnesses[i])\n        selected.append(population[winner].copy())\n    return selected\n\n\ndef crossover(p1, p2, n_features, max_depth=17):\n    \"\"\"Subtree crossover with Koza's 90/10 internal/terminal preference.\"\"\"\n    offspring = p1.copy()\n    donor = p2.copy()\n    on = offspring.get_all_nodes()\n    dn = donor.get_all_nodes()\n    if len(on) <= 1 or len(dn) <= 1:\n        return offspring\n    internal_off = [i for i, (n, p, c) in enumerate(on) if not n.is_terminal and p is not None]\n    if internal_off and random.random() < 0.9:\n        op = random.choice(internal_off)\n    else:\n        op = random.randint(1, len(on) - 1)\n    internal_don = [i for i, (n, p, c) in enumerate(dn) if not n.is_terminal]\n    if internal_don and random.random() < 0.9:\n        dp = random.choice(internal_don)\n    else:\n        dp = random.randint(0, len(dn) - 1)\n    ds = dn[dp][0].copy()\n    node, par, ci = on[op]\n    if par is not None:\n        par.children[ci] = ds\n    else:\n        offspring = ds\n    if offspring.depth() > max_depth:\n        return p1.copy()\n    return offspring\n\n\ndef mutation(parent, n_features, max_depth=17):\n    \"\"\"Multi-type mutation: subtree(30%), constant(25%), terminal(15%),\n    operator(15%), hoist(15%).\n\n    Novel contribution: five mutation operators with carefully balanced\n    probabilities. Constant perturbation with multi-scale noise enables\n    fine-grained optimization. Terminal replacement swaps variables and\n    constants. Operator replacement explores nearby function structures.\n    Hoist provides anti-bloat simplification.\n    \"\"\"\n    r = random.random()\n    o = parent.copy()\n    if r < 0.30:\n        # Subtree mutation — structural exploration\n        sz = o.size()\n        if sz <= 1:\n            return generate_tree('grow', 3, n_features)\n        mp = random.randint(1, sz - 1)\n        ns = generate_tree('grow', 3, n_features)\n        nodes = o.get_all_nodes()\n        _, p, ci = nodes[mp]\n        if p is not None:\n            p.children[ci] = ns\n        else:\n            o = ns\n        return o if o.depth() <= max_depth else parent.copy()\n    elif r < 0.55:\n        # Constant perturbation with multi-scale noise\n        nodes = o.get_all_nodes()\n        cn = [(n, p, c) for n, p, c in nodes\n              if n.is_terminal and not (isinstance(n.value, str) and n.value.startswith('x'))]\n        if cn:\n            n, _, _ = random.choice(cn)\n            v = float(n.value)\n            scale = random.choice([0.01, 0.05, 0.1, 0.5, 1.0])\n            n.value = str(round(v + random.gauss(0, scale), 2))\n        return o\n    elif r < 0.70:\n        # Terminal replacement (swap constant<->variable or new value)\n        nodes = o.get_all_nodes()\n        terminals = [(n, p, c) for n, p, c in nodes if n.is_terminal and p is not None]\n        if terminals:\n            n, par, ci = random.choice(terminals)\n            par.children[ci] = random_terminal(n_features)\n        return o\n    elif r < 0.85:\n        # Operator replacement — explore function neighborhood\n        nodes = o.get_all_nodes()\n        ops = [(n, p, c) for n, p, c in nodes if not n.is_terminal]\n        if ops:\n            n, _, _ = random.choice(ops)\n            ar = len(n.children)\n            same = [op for op, (_, a) in OPERATORS.items() if a == ar and op != n.value]\n            if same:\n                n.value = random.choice(same)\n        return o\n    else:\n        # Hoist mutation — simplification / anti-bloat\n        nodes = o.get_all_nodes()\n        internal = [n for n, p, c in nodes if not n.is_terminal and p is not None]\n        return random.choice(internal).copy() if internal else parent.copy()\n\n\ndef _hill_climb_constants(tree, X, y, iters=25):\n    \"\"\"Multi-scale stochastic hill-climbing with integer rounding attempts.\n\n    Novel contribution: separates structure search (GP) from constant\n    optimization (local search), a memetic/Lamarckian approach. Multi-scale\n    perturbations allow both coarse and fine adjustments. Integer rounding\n    helps discover exact solutions with integer constants (common in benchmarks).\n    \"\"\"\n    best = tree\n    best_f = fitness_function(best, X, y)\n    for it in range(iters):\n        cand = best.copy()\n        nodes = cand.get_all_nodes()\n        cn = [(n, p, c) for n, p, c in nodes\n              if n.is_terminal and not (isinstance(n.value, str) and n.value.startswith('x'))]\n        if not cn:\n            break\n        n, _, _ = random.choice(cn)\n        v = float(n.value)\n        if random.random() < 0.15:\n            # Try rounding to nearest integer — helps find exact solutions\n            n.value = str(round(v))\n        elif random.random() < 0.1:\n            # Try rounding to nearest 0.5\n            n.value = str(round(v * 2) / 2)\n        else:\n            scale = random.choice([0.01, 0.05, 0.1, 0.5, 1.0])\n            n.value = str(round(v + random.gauss(0, scale), 3))\n        cf = fitness_function(cand, X, y)\n        if cf < best_f:\n            best, best_f = cand, cf\n    return best\n\n\ndef _apply_linear_scaling(tree, X, y, max_depth=17):\n    \"\"\"Structurally inject optimal linear scaling a*tree+b.\n\n    Novel contribution: makes the optimal affine transformation heritable\n    by encoding it in the tree structure, so it can be refined by crossover\n    and constant mutation in subsequent generations.\n    \"\"\"\n    yp = safe_evaluate(tree, X)\n    vp = float(np.var(yp))\n    if vp < 1e-15:\n        return tree.copy()\n    my, mp = float(np.mean(y)), float(np.mean(yp))\n    a = float(np.mean((y - my) * (yp - mp))) / vp\n    b = my - a * mp\n    if abs(a - 1.0) < 0.02 and abs(b) < 0.02:\n        return tree.copy()\n    scaled = Node('add', [Node('mul', [tree.copy(), Node(str(round(a, 4)))]),\n                           Node(str(round(b, 4)))])\n    return scaled if scaled.depth() <= max_depth else tree.copy()\n\n\ndef evolve_one_generation(population, fitnesses, X_train, y_train,\n                          n_features, pop_size,\n                          crossover_rate=0.78, mutation_rate=0.17,\n                          max_depth=17):\n    \"\"\"Hybrid GP: lexicase+tournament selection, adaptive parsimony, memetic\n    constant optimization, structural linear scaling, stagnation-triggered\n    diversity injection.\n\n    Novel algorithm integrating five complementary mechanisms:\n    1. Hybrid selection (lexicase+tournament) balances diversity and convergence\n    2. Adaptive parsimony prevents bloat while preserving useful complexity\n    3. Memetic constant hill-climbing with integer rounding on elite individuals\n    4. Structural linear scaling injects heritable affine transformations\n    5. Stagnation detection triggers increased diversity when search stalls\n    \"\"\"\n    new_pop = []\n\n    # Adaptive parsimony pressure\n    lengths = np.array([t.size() for t in population], dtype=float)\n    fit_arr = np.array(fitnesses, dtype=float)\n    lv = float(np.var(lengths))\n    pc = 0.0\n    if lv > 1e-15:\n        pc = max(0.0, min(float(np.cov(lengths, fit_arr)[1, 0]) / lv, 0.001))\n    penalized = [f + pc * l for f, l in zip(fitnesses, lengths)]\n\n    # Stagnation detection\n    best_f = min(fitnesses)\n    if not hasattr(evolve_one_generation, '_prev_best'):\n        evolve_one_generation._prev_best = float('inf')\n        evolve_one_generation._stall = 0\n    if best_f < evolve_one_generation._prev_best * 0.999:\n        evolve_one_generation._prev_best = best_f\n        evolve_one_generation._stall = 0\n    else:\n        evolve_one_generation._stall += 1\n    stalled = evolve_one_generation._stall > 7\n\n    # Elitism: top individuals with linear scaling + constant hill-climbing\n    si = np.argsort(fitnesses)\n    n_elite = min(8, pop_size)\n    for i in range(n_elite):\n        elite = population[si[i]]\n        if i < 3:  # Linear scaling for top 3\n            scaled = _apply_linear_scaling(elite, X_train, y_train, max_depth)\n            if fitness_function(scaled, X_train, y_train) < fitnesses[si[i]]:\n                elite = scaled\n        iters = 35 if i < 3 else 15\n        new_pop.append(_hill_climb_constants(elite, X_train, y_train, iters=iters))\n\n    # Per-case errors for lexicase selection\n    errors = _per_case_errors(population, X_train, y_train)\n\n    # Diversity injection: more when stalled\n    n_imm = max(2, pop_size // 50) if stalled else max(1, pop_size // 100)\n    for _ in range(n\n... (truncated)"
7}

1{
2  "op": "replace",
3  "filename": "gplearn/custom_sr.py",
4  "start_line": 249,
5  "end_line": 254,
6  "content": "    Novel contribution: balances diversity maintenance (lexicase) with\n    convergence pressure (tournament). Lower lexicase proportion (30%) favors\n    convergence on structured problems while stagnation detection handles\n    diversity when needed.\n    \"\"\"\n    selected = []\n    pop_size = len(population)\n    for _ in range(n_select):\n        use_lexicase = _errors is not None and random.random() < 0.3"
7}

1{
2  "op": "replace",
3  "filename": "gplearn/custom_sr.py",
4  "start_line": 454,
5  "end_line": 456,
6  "content": "                          crossover_rate=0.82, mutation_rate=0.13,\n                          max_depth=17):"
7}

1{
2  "op": "replace",
3  "filename": "gplearn/custom_sr.py",
4  "start_line": 483,
5  "end_line": 485,
6  "content": "        iters = 50 if i < 3 else 20\n        new_pop.append(_hill_climb_constants(elite, X_train, y_train, iters=iters))"
7}

1{
2  "n": 1
3}