%pip install -q seekerdemo "setuptools<81"
# Application - Overnight charging of an electric bus fleet.
# 16 buses return to the depot between 19:00 and 23:00 and must be charged
# before their morning runs. The depot's grid connection is 900 kW; every kW of
# the night's peak costs a demand charge; evening electricity is expensive;
# and fast charging wears the batteries. When should each bus charge, and how
# fast? Today the depot plugs every bus in on arrival, at the slowest power
# that still gets it charged in time.
import random, seekerdemo as skr

random.seed(42)
B, SLOTS, H = 16, 48, 0.25                 # 15-min slots from 19:00 to 07:00
arrive = [random.randint(0, 16) for _ in range(B)]
leave  = [random.randint(36, 46) for _ in range(B)]
need   = [random.randint(150, 330) for _ in range(B)]          # kWh to recharge
POWER  = [50, 100, 150]                                          # kW options
price  = [0.34 if t < 12 else 0.22 if t < 20 else 0.14 for t in range(SLOTS)]
DEMAND_CHARGE, GRID_LIMIT, BASE_LOAD = 12.0, 900, 60     # EUR/kW/month, kW, kW
WEAR = {50: 0.000, 100: 0.010, 150: 0.035}                       # EUR per kWh

def plan(level, delay):
    """Turn search values into charging sessions; always finishes in time."""
    sessions = []
    for b in range(B):
        window = leave[b] - arrive[b]
        lv = max(level[b], next(i for i, p in enumerate(POWER)   # auto-upgrade if
                 if -(-need[b] // (p * H)) <= window))           # too slow to finish
        dur = int(-(-need[b] // (POWER[lv] * H)))
        start = arrive[b] + min(delay[b], window - dur)
        sessions.append((start, dur, POWER[lv]))
    return sessions

def night_cost(sessions):
    load = [BASE_LOAD] * SLOTS
    wear = 0.0
    for b, (s, d, p) in enumerate(sessions):
        for t in range(s, s + d):
            load[t] += p
        wear += WEAR[p] * need[b]
    energy = sum(l * H * c for l, c in zip(load, price))
    peak = max(load)
    over = max(0, peak - GRID_LIMIT)
    monthly = 30 * (energy + wear) + DEMAND_CHARGE * peak       # same plan every night
    return monthly + 1000 * over, peak, energy, wear            # over: never allowed

class Depot(skr.UserEvalTerm):
    def __init__(self):
        super().__init__()
    def evaluate(self, v):
        lv, dl = [round(x) for x in v[:B]], [round(x) for x in v[B:]]
        return night_cost(plan(lv, dl))[0]
    def incremental(self):
        return False

env = skr.Env('')
level = [env.categorical(0, 2) for _ in range(B)]                # how fast
delay = [env.ordinal(0, leave[b] - arrive[b]) for b in range(B)] # when to start
depot = Depot()
cost = env.user_defined_term(level + delay, depot)
env.minimize(cost, 5)

best = plan([round(x.get_value()) for x in level], [round(x.get_value()) for x in delay])
for name, sessions in [("plug in on arrival", plan([0] * B, [0] * B)),
                       ("Seeker schedule   ", best)]:
    monthly, peak, energy, wear = night_cost(sessions)
    print(f"{name}: peak {peak:3.0f} kW | per month: energy EUR {30 * energy:6,.0f},"
          f" demand charge EUR {DEMAND_CHARGE * peak:5,.0f}, battery wear EUR {30 * wear:4,.0f}"
          f" -> total EUR {monthly:6,.0f}")
env.end()

plug in on arrival: peak 810 kW | per month: energy EUR 32,142, demand charge EUR 9,720, battery wear EUR    0 -> total EUR 41,862
Seeker schedule   : peak 610 kW | per month: energy EUR 21,192, demand charge EUR 7,320, battery wear EUR  203 -> total EUR 28,715

%pip install -q seekerdemo "setuptools<81"
# Application - Sequencing a paint line with color changeovers and due dates.
# 18 orders wait for one powder-coating line. Switching color means cleaning
# the booth; going from dark to light takes much longer (and more solvent) than
# light to dark. Late orders cost money (contract penalties, idle assembly).
# Planners usually sort by due date. What order minimizes the total cost?
import random, seekerdemo as skr

random.seed(8)
N = 18
LIGHT = [95, 80, 60, 40, 20, 5]                 # lightness of the 6 colors
color = [random.randrange(6) for _ in range(N)]
minutes = [random.randint(25, 70) for _ in range(N)]
due = sorted(random.randint(60, 900) for _ in range(N))
penalty = [random.choice([2, 4, 8]) for _ in range(N)]   # EUR per late minute
def clean(a, b):                                 # booth cleaning, minutes
    return 0 if a == b else 12 + max(0, LIGHT[b] - LIGHT[a])
CLEAN_COST = 3.0                                  # EUR per cleaning minute

def model(env, seq):
    """seq[k] = order in position k  (search variable or fixed sequence)"""
    COL, DUR, DUE, PEN = (env.tensor(x) for x in (color, minutes, due, penalty))
    SETUP = env.tensor([[clean(a, b) for b in range(6)] for a in range(6)])
    t, cleaning, late = env.convert(0), [], []
    for k in range(N):
        if k > 0:
            c = env.index([env.index([seq[k-1]], COL), env.index([seq[k]], COL)], SETUP)
            cleaning.append(c)
            t = t + c
        t = t + env.index([seq[k]], DUR)                 # finish time, derived
        late.append(env.index([seq[k]], PEN) * env.max_0(t - env.index([seq[k]], DUE)))
    clean_min, late_cost = env.sum(cleaning), env.sum(late)
    return clean_min, late_cost, CLEAN_COST * clean_min + late_cost

# Baseline: earliest due date first (orders are already sorted by due date)
env = skr.Env('')
c, l, total = model(env, env.convert(list(range(N))))
env.evaluate()
print(f"due-date order : cleaning {c.get_value():4.0f} min | late penalties EUR {l.get_value():5.0f}"
      f" | total EUR {total.get_value():5.0f}")
env.end()

env = skr.Env('')
seq = env.permutation(N).get_permutation()
c, l, total = model(env, seq)
env.minimize(total, 5)
print(f"Seeker sequence: cleaning {c.get_value():4.0f} min | late penalties EUR {l.get_value():5.0f}"
      f" | total EUR {total.get_value():5.0f}")
print("colors in order:", " ".join("ABCDEF"[color[int(s.get_value())]] for s in seq))
env.end()

due-date order : cleaning  424 min | late penalties EUR 20882 | total EUR 22154
Seeker sequence: cleaning  182 min | late penalties EUR  6058 | total EUR  6604
colors in order: A C B B B B B B F F E D D D D D A A

%pip install -q seekerdemo "setuptools<81"
# Application - Pricing three subscription tiers with monthly and annual billing.
# Customers from three segments pick Basic, Pro or Business, billed monthly or
# annually, or walk away (a segmented logit choice model, e.g. estimated from
# a conjoint study). Rules: prices end in 9, each tier costs clearly more than
# the one below, and the annual discount is between 10% and 30%.
import seekerdemo as skr

#            size  price sens.  value of Basic, Pro, Business   dislikes annual
segments = [(5000, 0.080,       [3.0, 3.6, 3.4],                 0.6),   # freelancers
            (2500, 0.035,       [3.0, 4.6, 5.0],                 0.3),   # small teams
            ( 600, 0.012,       [2.5, 4.4, 6.2],                 0.1)]   # companies
SERVE = [4, 9, 22]                   # monthly cost to serve one customer per tier

def tier_prices(env, k_basic, k_pro, k_biz):
    basic = 10 * k_basic - 1                         # ends in 9 by construction
    pro = basic + 10 * k_pro                         # always more than Basic
    biz = pro + 10 * k_biz                           # always more than Pro
    return [basic, pro, biz]

def monthly_profit(env, prices, discount):
    total, customers = 0, 0
    for size, beta, value, annual_pain in segments:
        u_month = [env.exp(v - beta * p) for v, p in zip(value, prices)]
        u_year = [env.exp(v - annual_pain - beta * p * (1 - discount))
                  for v, p in zip(value, prices)]
        denom = 1 + env.sum(u_month) + env.sum(u_year)          # 1 = walk away
        customers = customers + size * (denom - 1) / denom
        for j in range(3):
            margin_m = prices[j] - SERVE[j]
            margin_y = prices[j] * (1 - discount) - SERVE[j]
            total = total + size * (u_month[j] * margin_m + u_year[j] * margin_y) / denom
    return total, customers

# Today's price list: 19 / 49 / 99, annual billing 17% off
env = skr.Env('')
today, cust = monthly_profit(env, [env.convert(p) for p in (19, 49, 99)], env.convert(0.17))
env.evaluate()
print(f"today  19 /  49 /  99, annual -17%: {cust.get_value():5,.0f} customers,"
      f" profit EUR {today.get_value():8,.0f} per month")
env.end()

env = skr.Env('')
k = [env.ordinal(1, 30), env.ordinal(1, 30), env.ordinal(1, 60)]
discount = env.continuous(0.10, 0.30)
prices = tier_prices(env, *k)
profit, cust = monthly_profit(env, prices, discount)
env.maximize(profit, 5)
p = [round(x.get_value()) for x in prices]
print(f"Seeker {p[0]:3} / {p[1]:3} / {p[2]:3}, annual -{100*discount.get_value():.0f}%:"
      f" {cust.get_value():5,.0f} customers, profit EUR {profit.get_value():8,.0f} per month")
env.end()

today  19 /  49 /  99, annual -17%: 7,563 customers, profit EUR  199,346 per month
Seeker  49 / 109 / 249, annual -10%: 5,013 customers, profit EUR  327,695 per month

%pip install -q seekerdemo "setuptools<81"
# Application - Where should next year's marketing budget go?
# Up to EUR 4M, in blocks of EUR 100k, across 8 channels. Every channel
# saturates (the 10th block buys less than the 1st), TV and podcasts carry a
# fixed production cost before they reach anyone, and TV lifts the response of
# paid search. Any money not spent is kept. Last year's split is the baseline.
import math, seekerdemo as skr

channels = ["TV", "Paid search", "Social", "Display", "Podcasts", "Trade fairs",
            "Email", "Partner promos"]
peak  = [9.0, 5.5, 4.0, 2.2, 2.6, 3.0, 1.4, 2.4]   # max extra sales, EUR M
scale = [12, 6, 7, 5, 6, 9, 2, 6]                   # blocks to reach 63% of peak
fixed = [6, 0, 0, 0, 2, 0, 0, 0]                    # blocks before any effect
MARGIN, BLOCKS = 0.45, 40                           # gross margin, 40 x 100k
last_year = [14, 10, 6, 4, 0, 4, 1, 1]              # 40 blocks, all spent

def profit(env, spend):                              # spend: list of 8 terms
    sales = []
    for c in range(8):
        effective = env.max_0(spend[c] - fixed[c])
        sales.append(peak[c] * (1 - env.exp(-1 * effective / scale[c])))
    tv_on = sales[0] / peak[0]                        # 0..1 share of TV's effect
    sales[1] = sales[1] * (1 + 0.5 * tv_on)          # TV halo on search
    return MARGIN * env.sum(sales) - 0.1 * env.sum(spend), env.sum(sales)

env = skr.Env('')
base, base_sales = profit(env, env.convert(last_year))
env.evaluate()
print(f"last year's split: spend EUR 4.0M, extra sales EUR {base_sales.get_value():4.1f}M,"
      f" profit EUR {base.get_value():4.2f}M")
env.end()

env = skr.Env('')
split = env.split_total(9, BLOCKS).get_combination()   # 8 channels + "keep it"
gain, sales = profit(env, split[:8])
env.maximize(gain, 5)
blocks = [round(x.get_value()) for x in split]
print(f"Seeker split     : spend EUR {sum(blocks[:8]) / 10:.1f}M, extra sales EUR "
      f"{sales.get_value():4.1f}M, profit EUR {gain.get_value():4.2f}M")
for name, b, l in zip(channels, blocks, last_year):
    print(f"   {name:15} EUR {b / 10:3.1f}M   (last year {l / 10:3.1f}M)")
print(f"   {'kept':15} EUR {blocks[8] / 10:3.1f}M")
env.end()

last year's split: spend EUR 4.0M, extra sales EUR 15.4M, profit EUR 2.95M
Seeker split     : spend EUR 4.0M, extra sales EUR 16.3M, profit EUR 3.33M
   TV              EUR 2.0M   (last year 1.4M)
   Paid search     EUR 0.8M   (last year 1.0M)
   Social          EUR 0.5M   (last year 0.6M)
   Display         EUR 0.2M   (last year 0.4M)
   Podcasts        EUR 0.0M   (last year 0.0M)
   Trade fairs     EUR 0.1M   (last year 0.4M)
   Email           EUR 0.2M   (last year 0.1M)
   Partner promos  EUR 0.2M   (last year 0.1M)
   kept            EUR 0.0M

%pip install -q seekerdemo "setuptools<81" numpy
# Application - Day-ahead battery plan for a factory with rooftop solar.
# Tomorrow's solar output is uncertain. An "analog ensemble" (a simple, widely
# used ML forecast) takes the 40 past days whose weather forecast looked most
# like tomorrow's and uses their real 24-hour solar profiles as scenarios.
# Seeker plans the battery's 24 hourly moves against all scenarios at once,
# minimizing expected cost plus the cost of the worst 10% of days.
import numpy as np, seekerdemo as skr

rng = np.random.default_rng(1)
hours = np.arange(24)
clear_sky = np.clip(np.sin((hours - 6) / 13 * np.pi), 0, None) * 900      # kW
cloud_fc = rng.uniform(0, 1, 730)                         # 2 years of forecasts
cloud_real = np.clip(cloud_fc + rng.normal(0, .2, (24, 730)), 0, 1).T
history = clear_sky * (1 - .8 * cloud_real ** 1.5)        # what the roof produced

tomorrow_fc = 0.45
analogs = history[np.argsort(abs(cloud_fc - tomorrow_fc))[:40]]   # the ML step

class Solar(skr.UserVectorDistribution):
    def __init__(self):
        super().__init__()
    def sample_n(self, n):                                # whole-day profiles,
        pick = np.random.default_rng(7).integers(0, 40, n) # same draws every run
        return analogs[pick].tolist()

load = 500 + 150 * ((hours >= 7) & (hours <= 18))         # factory kW
buy = np.where((hours >= 17) & (hours <= 21), 0.38, np.where(hours >= 7, 0.27, 0.16))
SELL, CAP, POWER = 0.07, 1200, 400                       # EUR/kWh, kWh, kW

def plan(flows_fixed=None):
    env = skr.Env('', True)
    dist = Solar()                                        # keep a reference
    solar = env.user_vector_distribution(dist, 200)
    flows = flows_fixed or [env.continuous(-POWER, POWER) for _ in hours]
    soc, cost = env.convert(CAP / 2), 0
    for h in hours:
        new = env.min([env.convert(CAP), env.max_0(soc + flows[h])])  # auto-clip
        actual = new - soc                                # what really flowed
        soc = new
        net = load[h] - solar[h] + actual                 # grid draw (random)
        cost = cost + buy[h] * env.max_0(net) - SELL * env.max_0(-1 * net)
    cost = cost + 0.2 * (CAP / 2 - soc)                   # value of leftover charge
    exp_cost = env.aggregate_mean(cost)
    worst10 = env.aggregate_cvar(cost, 0.1, True)
    target = exp_cost + 0.5 * worst10
    if flows_fixed: env.evaluate()
    else: env.minimize(target, 6)
    res = exp_cost.get_value(), worst10.get_value()
    env.end()
    return res

idle = plan([0] * 24)
rule = plan([POWER if 10 <= h <= 12 else -POWER if 17 <= h <= 19 else 0 for h in hours])
best = plan()
for name, (e, w) in [("battery idle                ", idle),
                     ("rule: charge 10-13, use 17-20", rule),
                     ("Seeker plan                 ", best)]:
    print(f"{name}: expected cost EUR {e:6.0f} | worst 10% of days EUR {w:6.0f}")

battery idle                : expected cost EUR   2143 | worst 10% of days EUR   2361
rule: charge 10-13, use 17-20: expected cost EUR   1980 | worst 10% of days EUR   2193
Seeker plan                 : expected cost EUR   1858 | worst 10% of days EUR   2093

%pip install -q seekerdemo "setuptools<81"
# Concept 1 - Search variables are not (always) the decisions.
# A kiln fires 10 ceramic batches, one at a time. Each batch has a time it is
# dry enough to fire, a firing duration, a promised pickup time and a priority.
# We want START TIMES. But we SEARCH over the ORDER of the batches and simply
# compute the earliest possible start for each one.
import seekerdemo as skr

ready = [0, 0, 30, 45, 60, 60, 90, 120, 150, 180]    # minutes
dur   = [55, 40, 70, 35, 60, 45, 50, 65, 40, 55]
due   = [70, 60, 160, 130, 230, 200, 260, 330, 300, 400]
prio  = [3, 1, 2, 1, 3, 2, 1, 3, 2, 1]
n = len(dur)

env = skr.Env('')
order = env.permutation(n).get_permutation()   # the ONLY search variables
R, D, U, P = (env.tensor(x) for x in (ready, dur, due, prio))

end, lateness, start = env.convert(0), [], []
for k in range(n):                             # derived terms, NOT variables
    job = order[k]
    s = env.max([end, env.index([job], R)])    # earliest feasible start
    end = s + env.index([job], D)
    lateness.append(env.index([job], P) * env.max_0(end - env.index([job], U)))
    start.append(s)

late = env.sum(lateness)
env.minimize(late, 3)

print(f"weighted lateness: {late.get_value():.0f}   constraints used: 0")
for s, job in sorted((s.get_value(), int(j.get_value())) for s, j in zip(start, order)):
    print(f"  start {s:4.0f} min   batch {job}")
env.end()

# The textbook alternative searches over start times directly:
#   s = [env.ordinal(ready[j], 600) for j in range(n)]
#   for every pair i<j: enforce s or that is disjunctive constraints and almost every random change to a start time creates an overlap. searching over the order makes overlaps impossible.></j:>

weighted lateness: 895   constraints used: 0
  start    0 min   batch 0
  start   55 min   batch 1
  start   95 min   batch 3
  start  130 min   batch 5
  start  175 min   batch 4
  start  235 min   batch 8
  start  275 min   batch 7
  start  340 min   batch 2
  start  410 min   batch 6
  start  460 min   batch 9

%pip install -q seekerdemo "setuptools<81"
# Concept 2 - Build feasibility into the model instead of forbidding values.
# 16 bottling lines must fill an order of exactly D (thousand) bottles. A line
# is either OFF or runs between its minimum and maximum output (no line can run
# at 10%). Each running line costs a start-up fee plus energy, and is most
# efficient near its sweet spot.
import random, seekerdemo as skr

random.seed(7)
N = 16
hi = [random.randint(60, 140) for _ in range(N)]
lo = [int(h * random.uniform(.55, .75)) for h in hi]
sweet = [(l + h) // 2 for l, h in zip(lo, hi)]
setup = [random.randint(300, 700) for _ in range(N)]
D = int(sum(hi) * 0.45)

def total_cost(env, out):
    on = [env.if_(q > 0, 1, 0) for q in out]
    energy = [q * (8 + 0.004 * env.sqr(q - s)) for q, s in zip(out, sweet)]
    return env.sum([u * f for u, f in zip(on, setup)]) + env.sum(energy)

def built_in(seconds):
    env = skr.Env('')
    out = []
    for l, h in zip(lo, hi):
        q_search = env.ordinal(0, h - l + 1)            # 0 = off, 1.. = running
        out.append(env.if_(q_search > 0, q_search + (l - 1), 0))
    env.enforce_eq(env.sum(out), D)                     # 1 constraint
    cost = total_cost(env, out)
    env.minimize(cost, seconds)
    res = cost.get_value()
    env.end()
    return res

def forbidden_values(seconds):
    env = skr.Env('')
    out = [env.ordinal(0, h) for h in hi]
    for q, l in zip(out, lo):                           # 16 extra constraints
        env.enforce(env.or_([env.eq(q, 0), env.geq(q, l)]))
    env.enforce_eq(env.sum(out), D)
    cost = total_cost(env, out)
    env.minimize(cost, seconds)
    res = cost.get_value()
    env.end()
    return res

print(f"order: {D}k bottles on {N} lines, 2 seconds each")
print(f"feasibility built in : cost {built_in(2):,.0f}")
print(f"values forbidden     : cost {forbidden_values(2):,.0f}")

order: 677k bottles on 16 lines, 2 seconds each
feasibility built in : cost 8,665
values forbidden     : cost 9,066

%pip install -q seekerdemo "setuptools<81"
# Concept 3 - Don't binarize. Say what you mean.
# Ten components are sourced from four suppliers. Not every supplier can make
# every component. Each supplier gives volume discounts (the more you buy, the
# cheaper every unit gets) and charges a fixed onboarding fee if used at all.
import seekerdemo as skr

vol   = [120, 80, 200, 60, 150, 90, 40, 170, 110, 70]          # units/week
price = [[9.0, 9.6, 8.8, 9.9], [4.1, 3.9, 4.4, 4.0], [12.5, 12.0, 13.1, 12.2],
         [6.0, 6.3, 5.8, 6.1], [7.7, 7.2, 7.9, 7.5], [3.3, 3.5, 3.1, 3.4],
         [15.0, 14.2, 15.5, 14.8], [8.4, 8.9, 8.1, 8.6], [5.5, 5.2, 5.9, 5.4],
         [10.2, 10.8, 9.9, 10.5]]
can   = [[0, 1, 2, 3], [0, 1, 3], [1, 2, 3], [0, 2], [0, 1, 2, 3],
         [0, 2, 3], [1, 3], [0, 1, 2], [1, 2, 3], [0, 2, 3]]   # who can make what
fee   = [900, 1100, 700, 1000]
S, N = 4, len(vol)

def discount(units):                     # up to 18% off, smoothly with volume
    return 1 - 0.18 * units / (units + 400)

def say_what_you_mean(seconds):
    env = skr.Env('')
    pick = [env.categorical(0, S - 1, can[i]) for i in range(N)]   # 10 variables
    cost = 0
    for s in range(S):
        used = [env.if_(pick[i] == s, 1, 0) for i in range(N)]
        units = env.sum([u * vol[i] for i, u in enumerate(used)])
        spend = env.sum([u * vol[i] * price[i][s] for i, u in enumerate(used)])
        cost = cost + spend * discount(units) + env.if_(units > 0, fee[s], 0)
    env.minimize(cost, seconds)
    res = cost.get_value()
    env.end()
    return res

def binarized(seconds):
    env = skr.Env('')
    x = [[env.ordinal(0, 1 if s in can[i] else 0) for s in range(S)]
         for i in range(N)]                                          # 40 variables
    for i in range(N):
        env.enforce_eq(env.sum(x[i]), 1)                             # 10 constraints
    cost = 0
    for s in range(S):
        units = env.sum([x[i][s] * vol[i] for i in range(N)])
        spend = env.sum([x[i][s] * vol[i] * price[i][s] for i in range(N)])
        cost = cost + spend * discount(units) + env.if_(units > 0, fee[s], 0)
    env.minimize(cost, seconds)
    res = cost.get_value()
    env.end()
    return res

print(f"categorical choice (10 vars, 0 constraints): {say_what_you_mean(2):,.0f} per week")
print(f"binarized   (40 vars, 10 constraints):       {binarized(2):,.0f} per week")

categorical choice (10 vars, 0 constraints): 9,581 per week
binarized   (40 vars, 10 constraints):       9,605 per week

%pip install -q seekerdemo "setuptools<81"
# Concept 4 - Pick the variable that already knows the rules.
# 24 engineers form 4 project teams of exactly 6. Every team should get a
# similar mix of skills (coding, domain, design) and at least one senior, and
# two pairs of people should not end up together.
import random, seekerdemo as skr

random.seed(3)
N, T, SIZE = 24, 4, 6
skill = [[random.randint(1, 9) for _ in range(3)] for _ in range(N)]
senior = [1 if random.random() < .25 else 0 for _ in range(N)]
apart = [(2, 9), (5, 17)]
target = [sum(s[d] for s in skill) / T for d in range(3)]

def team_score(env, members):            # members: list of person-index terms
    K = env.tensor(skill)
    sen = env.tensor(senior)
    dev = [env.sqr(env.sum([env.index([m, env.convert(d)], K) for m in members])
                   - target[d]) for d in range(3)]
    no_senior = env.if_(env.sum([env.index([m], sen) for m in members]) > 0, 0, 1)
    return env.sum(dev) + 50 * no_senior

def permutation_model(seconds):
    env = skr.Env('')
    p = env.permutation(N)
    seat, where = p.get_permutation(), p.get_permutation_inverse()
    team_of = [env.floor(where[i] / SIZE) for i in range(N)]       # derived
    score = env.sum([team_score(env, seat[t*SIZE:(t+1)*SIZE]) for t in range(T)])
    score = score + 100 * env.sum([env.if_(team_of[a] == team_of[b], 1, 0)
                                   for a, b in apart])
    env.minimize(score, seconds)
    res = score.get_value()
    env.end()
    return res

def label_model(seconds):
    env = skr.Env('')
    team = [env.categorical(0, T - 1) for _ in range(N)]
    score = 0
    for t in range(T):
        inside = [env.if_(team[i] == t, 1, 0) for i in range(N)]
        env.enforce_eq(env.sum(inside), SIZE)                       # 4 constraints
        dev = [env.sqr(env.sum([u * skill[i][d] for i, u in enumerate(inside)])
                       - target[d]) for d in range(3)]
        sen = env.sum([u * senior[i] for i, u in enumerate(inside)])
        score = score + env.sum(dev) + 50 * env.if_(sen > 0, 0, 1)
    score = score + 100 * env.sum([env.if_(team[a] == team[b], 1, 0) for a, b in apart])
    env.minimize(score, seconds)
    res = score.get_value()
    env.end()
    return res

print(f"permutation, cut into blocks of 6: imbalance {permutation_model(3):.1f}")
print(f"team labels + size constraints:    imbalance {label_model(3):.1f}")

permutation, cut into blocks of 6: imbalance 2.5
team labels + size constraints:    imbalance 8.5

%pip install -q seekerdemo "setuptools<81"
# Concept 5 - Avoid plateaus: give the search a slope to follow.
# 45 sample batches go to 6 lab analyzers of different generations; each
# analyzer is faster on some assays than others. We want all results as early
# as possible (minimize the makespan). Problem: moving one batch rarely changes
# the busiest analyzer's finish time, so the objective is flat almost
# everywhere. Adding total work / M (M large) breaks the ties without changing
# what "best" means: total work is at most 4,050 min, so /10000 stays below
# one minute of makespan and only ever breaks ties.
import random, seekerdemo as skr

random.seed(1)
N, A = 45, 6
minutes = [[random.randint(20, 90) for _ in range(A)] for _ in range(N)]

def solve(tie_breaker, seconds):
    env = skr.Env('')
    lab = [env.categorical(0, A - 1) for _ in range(N)]
    load = [env.sum([env.if_(lab[i] == a, minutes[i][a], 0) for i in range(N)])
            for a in range(A)]
    makespan = env.max(load)
    total = env.sum(load)
    obj = makespan + total / 10000 if tie_breaker else makespan
    env.minimize(obj, seconds)
    res = makespan.get_value(), total.get_value()
    env.end()
    return res

for tb, name in [(False, "makespan only        "), (True, "makespan + work/1e4  ")]:
    ms, tot = solve(tb, 1)
    print(f"{name}: all results after {ms:.0f} min, analyzer time used {tot:.0f} min")

makespan only        : all results after 243 min, analyzer time used 1454 min
makespan + work/1e4  : all results after 233 min, analyzer time used 1336 min

%pip install -q seekerdemo "setuptools<81"
# Concept 6 - When the system is easier to SIMULATE than to write as algebra,
# plug the simulation in as a user-defined term.
# A spare-parts warehouse reorders with an (s, S) policy: when stock falls to s
# or below, order up to S. Orders arrive after a 5-day lead time, unmet demand
# is lost. Demand peaks in summer, so we allow a separate summer policy. We
# replay last year's daily demand and let Seeker tune all four numbers.
import random, seekerdemo as skr

random.seed(5)
demand = [max(0, int(random.gauss(18, 7) + (12 if 150 < d < 240 else 0)))
          for d in range(365)]                       # last year, with a summer peak
LEAD, HOLD, ORDER, LOST = 5, 0.4, 250.0, 30.0        # costs per unit/day, order, unit

SUMMER = range(150, 240)

def simulate(s1, S1, s2, S2):                # (s1,S1) normal, (s2,S2) summer
    stock, pipeline, cost = S1, [], 0.0
    for day, d in enumerate(demand):
        s, S = (s2, S2) if day in SUMMER else (s1, S1)
        stock += sum(q for t, q in pipeline if t == day)
        pipeline = [(t, q) for t, q in pipeline if t > day]
        sold = min(stock, d)
        cost += LOST * (d - sold)
        stock -= sold
        on_order = sum(q for _, q in pipeline)
        if stock + on_order <= s:
            pipeline.append((day + LEAD, S - stock - on_order))
            cost += ORDER
        cost += HOLD * stock
    return cost

class YearlyCost(skr.UserEvalTerm):
    def __init__(self):
        super().__init__()
    def evaluate(self, v):                   # v = current values of the 4 terms
        return simulate(*[round(x) for x in v])
    def incremental(self):
        return False

env = skr.Env('')
s1, s2 = env.ordinal(0, 400), env.ordinal(0, 400)
S1 = s1 + env.ordinal(1, 600)                # S > s by construction, no constraint
S2 = s2 + env.ordinal(1, 600)
sim = YearlyCost()                           # keep a reference while solving
cost = env.user_defined_term([s1, S1, s2, S2], sim)
env.minimize(cost, 3)
v = [round(t.get_value()) for t in (s1, S1, s2, S2)]
print(f"Seeker policy : normal s={v[0]}, S={v[1]} | summer s={v[2]}, S={v[3]}"
      f" | yearly cost {cost.get_value():,.0f}")

# for comparison: textbook rule of thumb (EOQ + safety stock for 95% service)
mu = sum(demand) / 365
eoq = (2 * mu * ORDER / HOLD) ** .5
s0 = round(mu * LEAD + 1.65 * 7 * LEAD ** .5)
S0 = s0 + round(eoq)
print(f"Rule of thumb : s={s0}, S={S0} all year | yearly cost {simulate(s0, S0, s0, S0):,.0f}")
env.end()

Seeker policy : normal s=94, S=246 | summer s=147, S=341 | yearly cost 24,952
Rule of thumb : s=128, S=288 all year | yearly cost 30,158

%pip install -q seekerdemo "setuptools<81" numpy
# Concept 7 - Feed your ML model's uncertainty straight into the optimizer.
# A bakery decides how many croissants to bake for tomorrow (Saturday, 24 C).
# A small ML model (a bootstrap ensemble of regressions) learns from two years
# of sales history. Instead of a single forecast, it hands Seeker a posterior
# distribution, and Seeker optimizes the expected profit directly.
import numpy as np, seekerdemo as skr

rng = np.random.default_rng(0)
# --- history: temperature, weekend flag -> croissants sold -------------------
temp = rng.uniform(0, 32, 730)
weekend = (np.arange(730) % 7 >= 5).astype(float)
sold = 90 + 3.2 * temp + 55 * weekend + rng.normal(0, 18 + 0.6 * temp)
X = np.column_stack([np.ones(730), temp, weekend])

# --- ML: bootstrap ensemble -> posterior predictive samples for tomorrow -----
x_tomorrow = np.array([1, 24, 1])
def posterior(n):
    out = []
    for _ in range(n):
        idx = rng.integers(0, 730, 730)                     # bootstrap resample
        w, *_ = np.linalg.lstsq(X[idx], sold[idx], rcond=None)
        resid = sold[idx] - X[idx] @ w
        out.append(x_tomorrow @ w + rng.choice(resid))      # model + noise
    return out

class Demand(skr.UserDistribution):
    def __init__(self):
        super().__init__()
    def sample_n(self, n):
        return posterior(n)

PRICE, COST = 2.6, 0.9

def plan(min_in_stock=None, seconds=3):
    env = skr.Env('', True)                      # True = stochastic environment
    dist = Demand()
    demand = env.user_distribution(dist, 400)    # 400 posterior scenarios
    bake = env.ordinal(0, 600)
    profit = PRICE * env.min([bake, demand]) - COST * bake
    expected = env.aggregate_mean(profit)
    in_stock = env.aggregate_relative_frequency_geq(bake - demand, 0)
    if min_in_stock:                             # a chance constraint: one line
        env.enforce_geq(in_stock, min_in_stock)
    env.maximize(expected, seconds)
    res = bake.get_value(), expected.get_value(), in_stock.get_value()
    env.end()
    return res

for rule, target in [("max expected profit   ", None), ("... and >=90% in stock", 0.9)]:
    q, p, f = plan(target)
    print(f"{rule}: bake {q:.0f}, expected profit {p:.0f}, no sell-out {100*f:.0f}% of days")

D = np.array(posterior(4000))                    # baseline: bake the point forecast
q = round(D.mean())
print(f"bake the point forecast: bake {q}, expected profit "
      f"{np.mean(PRICE * np.minimum(q, D) - COST * q):.0f}, no sell-out {100*np.mean(D <= q):.0f}% of days")

max expected profit   : bake 229, expected profit 345, no sell-out 66% of days
... and >=90% in stock: bake 249, expected profit 331, no sell-out 90% of days
bake the point forecast: bake 219, expected profit 343, no sell-out 51% of days