1	# Polynomial mutation — C# Unsga3.Operators.Mutation.PolynomialMutation.
2	# Default η = 20. Probability per variable is an argument (C# Run uses 1/n).
3	# Walks the decision vector and box intervals together (no List.get /
4	# set_var_at per index). Same RNG order as the old indexed loop.
5	
6	import Base
7	import ./individual.bend as Ind
8	import ./bounds.bend as Bounds
9	import ./rng.bend as Rng
10	import ./util.bend as Util
11	
12	type Poly is Data:
13	  MkPoly{
14	    e: F32
15	  }
16	
17	def mk(+e: F32) -> Poly:
18	  MkPoly{e}
19	
20	def default() -> Poly:
21	  MkPoly{20.0}
22	
23	def eta(op: Poly) -> F32:
24	  match op:
25	    case MkPoly{e}:
26	      e
27	
28	# C# Unsga3Algorithm: mutationProbability ?? (1.0 / nVars) at the call site.
29	def typical_probability(n_var: Nat) -> F32:
30	  match n_var:
31	    case 0n:
32	      1.0
33	    case 1n+p:
34	      F32.div(1.0, F32.from_nat(1n+p))
35	
36	def result(pair: Rng.State & Ind.Individual) -> Ind.Individual:
37	  (rng, ind) = pair
38	  ind
39	
40	def deltaq.lo(+eta: F32, +mut_pow: F32, +delta1: F32, +rand: F32) -> F32:
41	  +xy = F32.sub(1.0, delta1)
42	  +val = F32.add(
43	    F32.mul(2.0, rand),
44	    F32.mul(F32.sub(1.0, F32.mul(2.0, rand)), F32.pow(xy, F32.add(eta, 1.0)))
45	  )
46	  F32.sub(F32.pow(val, mut_pow), 1.0)
47	
48	def deltaq.hi(+eta: F32, +mut_pow: F32, +delta2: F32, +rand: F32) -> F32:
49	  +xy = F32.sub(1.0, delta2)
50	  +val = F32.add(
51	    F32.mul(2.0, F32.sub(1.0, rand)),
52	    F32.mul(F32.mul(2.0, F32.sub(rand, 0.5)), F32.pow(xy, F32.add(eta, 1.0)))
53	  )
54	  F32.sub(1.0, F32.pow(val, mut_pow))
55	
56	def deltaq.pick(+eta: F32, +mut_pow: F32, +delta1: F32, +delta2: F32, +rand: F32, lo_side: Bool) -> F32:
57	  match lo_side:
58	    case True{}:
59	      deltaq.lo(eta, mut_pow, delta1, rand)
60	    case False{}:
61	      deltaq.hi(eta, mut_pow, delta2, rand)
62	
63	def deltaq(+eta: F32, +mut_pow: F32, +delta1: F32, +delta2: F32, +rand: F32) -> F32:
64	  deltaq.pick(eta, mut_pow, delta1, delta2, rand, F32.is_lt(rand, 0.5))
65	
66	def tiny_head(ivs: List<&2, Bounds.Interval>) -> Bool:
67	  match ivs:
68	    case Nil{}:
69	      False{}
70	    case +iv <> rest:
71	      F32.is_lt(F32.sub(Bounds.hi(iv), Bounds.lo(iv)), 1.0e-14)
72	
73	def advance(rng: Rng.State, rest: Nat) -> Rng.State:
74	  match rest:
75	    case 0n:
76	      rng
77	    case 1n+g:
78	      Rng.next_double(rng)
79	
80	def skip_kind(skip: Bool, tiny: Bool) -> Nat:
81	  match skip:
82	    case True{}:
83	      0n
84	    case False{}:
85	      match tiny:
86	        case True{}:
87	          0n
88	        case False{}:
89	          1n
90	
91	def finish(+m: Nat, pair: Rng.State & List<&2, F32>) -> Rng.State & Ind.Individual:
92	  (rng, xs) = pair
93	  (rng, Ind.from_variables(xs, m))
94	
95	def done(rng: Rng.State, acc: List<&2, F32>) -> Rng.State & List<&2, F32>:
96	  (rng, List.reverse(&2, F32, acc))
97	
98	# Kind: 0n = copy the head (skip or tiny span), 1n = mutate the head.
99	def loop.go(
100	  fuel: Nat,
101	  kind: Nat,
102	  xs: List<&2, F32>,
103	  ivs: List<&2, Bounds.Interval>,
104	  +eta: F32,
105	  rng: Rng.State,
106	  acc: List<&2, F32>,
107	  +p: F32
108	) -> Rng.State & List<&2, F32>:
109	  match fuel kind xs ivs:
110	    case 0n k a b:
111	      done(rng, acc)
112	    case 1n++f 0n +x <> rest +iv <> +ivrest:
113	      +rng2 = advance(rng, f)
114	      loop.go(
115	        f,
116	        skip_kind(F32.is_gt(Rng.last_double(rng2), p), tiny_head(ivrest)),
117	        rest, ivrest, eta, rng2, x <> acc, p
118	      )
119	    case 1n++f 1n+q +x <> rest +iv <> +ivrest:
120	      +lo = Bounds.lo(iv)
121	      +hi = Bounds.hi(iv)
122	      +dmax = F32.sub(hi, lo)
123	      +rng_a = Rng.next_double(rng)
124	      +rand = Rng.last_double(rng_a)
125	      +mut_pow = F32.div(1.0, F32.add(eta, 1.0))
126	      +dq = deltaq(eta, mut_pow, F32.div(F32.sub(x, lo), dmax), F32.div(F32.sub(hi, x), dmax), rand)
127	      +y = F32.clamp(F32.add(x, F32.mul(dq, dmax)), lo, hi)
128	      +rng_b = advance(rng_a, f)
129	      loop.go(
130	        f,
131	        skip_kind(F32.is_gt(Rng.last_double(rng_b), p), tiny_head(ivrest)),
132	        rest, ivrest, eta, rng_b, y <> acc, p
133	      )
134	    case 1n++f k Nil{} b:
135	      done(rng, acc)
136	    case 1n++f k a <> as Nil{}:
137	      done(rng, acc)
138	
139	def start(
140	  +eta: F32,
141	  +box: Bounds.Box,
142	  rng: Rng.State,
143	  +ind: Ind.Individual,
144	  +p: F32,
145	  n: Nat
146	) -> Rng.State & Ind.Individual:
147	  match n:
148	    case 0n:
149	      (rng, ind)
150	    case 1n+f:
151	      +rng2 = Rng.next_double(rng)
152	      +xs = Ind.variables(ind)
153	      +ivs = Bounds.intervals_of(box)
154	      finish(
155	        Ind.n_obj(ind),
156	        loop.go(
157	          1n+f,
158	          skip_kind(F32.is_gt(Rng.last_double(rng2), p), tiny_head(ivs)),
159	          xs, ivs, eta, rng2, Nil{}, p
160	        )
161	      )
162	
163	# p=0 leaves the decision vector unchanged before any NextDouble.
164	# Matches C# when per-variable probability is 0 (every draw is > 0).
165	def start_or_copy(
166	  +eta: F32,
167	  +box: Bounds.Box,
168	  rng: Rng.State,
169	  ind: Ind.Individual,
170	  +p: F32,
171	  n: Nat,
172	  zero: Bool
173	) -> Rng.State & Ind.Individual:
174	  match zero:
175	    case True{}:
176	      (rng, ind)
177	    case False{}:
178	      start(eta, box, rng, ind, p, n)
179	
180	def mutate(
181	  +op: Poly,
182	  +ind: Ind.Individual,
183	  +box: Bounds.Box,
184	  rng: Rng.State,
185	  +p: F32
186	) -> Rng.State & Ind.Individual:
187	  start_or_copy(eta(op), box, rng, ind, p, Ind.n_var(ind), Util.f32_is_zero(p))