1	# Simulated Binary Crossover — C# Unsga3.Operators.Crossover.SimulatedBinaryCrossover.
2	# Defaults match C#: η = 30, pair probability = 1.0 (Deb & Agrawal).
3	# Per-variable walk peels x1 / x2 / box intervals together (guide: the
4	# list is the shape a flat loop walks). Same RNG order as the old
5	# index + List.get / set_var_at loop; children are built once at the end.
6	
7	import Base
8	import ./individual.bend as Ind
9	import ./bounds.bend as Bounds
10	import ./rng.bend as Rng
11	import ./util.bend as Util
12	
13	type Sbx is Data:
14	  MkSbx{
15	    e: F32
16	    p: F32
17	  }
18	
19	def mk(+e: F32, +p: F32) -> Sbx:
20	  MkSbx{e, p}
21	
22	def default() -> Sbx:
23	  MkSbx{30.0, 1.0}
24	
25	def eta(op: Sbx) -> F32:
26	  match op:
27	    case MkSbx{e, p}:
28	      e
29	
30	def probability(op: Sbx) -> F32:
31	  match op:
32	    case MkSbx{e, p}:
33	      p
34	
35	def child1(trip: Rng.State & Ind.Individual & Ind.Individual) -> Ind.Individual:
36	  (rng, c1, c2) = trip
37	  c1
38	
39	def child2(trip: Rng.State & Ind.Individual & Ind.Individual) -> Ind.Individual:
40	  (rng, c1, c2) = trip
41	  c2
42	
43	def near(+x1: F32, +x2: F32) -> Bool:
44	  F32.is_lt(F32.abs(F32.sub(x1, x2)), 1.0e-14)
45	
46	# C#: if rand <= 1/alpha then (rand*alpha)^(1/(η+1)) else (1/(2-rand*alpha))^(1/(η+1)).
47	def beta.pick(+rand: F32, +alpha: F32, +inv: F32, lo: Bool) -> F32:
48	  match lo:
49	    case True{}:
50	      F32.pow(F32.mul(rand, alpha), inv)
51	    case False{}:
52	      F32.pow(F32.div(1.0, F32.sub(2.0, F32.mul(rand, alpha))), inv)
53	
54	def beta(+rand: F32, +alpha: F32, +inv: F32) -> F32:
55	  beta.pick(rand, alpha, inv, F32.is_le(rand, F32.div(1.0, alpha)))
56	
57	def close_heads(xs: List<&2, F32>, ys: List<&2, F32>) -> Bool:
58	  match xs ys:
59	    case +x <> xr +y <> yr:
60	      near(x, y)
61	    case Nil{} ys:
62	      False{}
63	    case x <> xs Nil{}:
64	      False{}
65	
66	# Draw the next per-var coin only when more variables remain.
67	def advance(rng: Rng.State, rest: Nat) -> Rng.State:
68	  match rest:
69	    case 0n:
70	      rng
71	    case 1n+g:
72	      Rng.next_double(rng)
73	
74	def finish(+m: Nat, trip: Rng.State & List<&2, F32> & List<&2, F32>) -> Rng.State & Ind.Individual & Ind.Individual:
75	  (rng, xs, ys) = trip
76	  (rng, Ind.from_variables(xs, m), Ind.from_variables(ys, m))
77	
78	def done(rng: Rng.State, acc1: List<&2, F32>, acc2: List<&2, F32>) -> Rng.State & List<&2, F32> & List<&2, F32>:
79	  (rng, List.reverse(&2, F32, acc1), List.reverse(&2, F32, acc2))
80	
81	def skip_kind(skip: Bool, close: Bool) -> Nat:
82	  match skip:
83	    case True{}:
84	      0n
85	    case False{}:
86	      match close:
87	        case True{}:
88	          0n
89	        case False{}:
90	          1n
91	
92	# Kind: 0n = copy heads (skip or close), 1n = SBX the heads.
93	# One def so skip/cross are not mutually recursive with `loop`.
94	def loop.go(
95	  fuel: Nat,
96	  kind: Nat,
97	  xs1: List<&2, F32>,
98	  xs2: List<&2, F32>,
99	  ivs: List<&2, Bounds.Interval>,
100	  +eta: F32,
101	  rng: Rng.State,
102	  acc1: List<&2, F32>,
103	  acc2: List<&2, F32>
104	) -> Rng.State & List<&2, F32> & List<&2, F32>:
105	  match fuel kind xs1 xs2 ivs:
106	    case 0n k a b c:
107	      done(rng, acc1, acc2)
108	    case 1n++f 0n +x1 <> +rest1 +x2 <> +rest2 +iv <> ivrest:
109	      +rng2 = advance(rng, f)
110	      loop.go(
111	        f,
112	        skip_kind(F32.is_gt(Rng.last_double(rng2), 0.5), close_heads(rest1, rest2)),
113	        rest1, rest2, ivrest, eta, rng2, x1 <> acc1, x2 <> acc2
114	      )
115	    case 1n++f 1n+p +x1 <> +rest1 +x2 <> +rest2 +iv <> ivrest:
116	      +lo = Bounds.lo(iv)
117	      +hi = Bounds.hi(iv)
118	      +gt = F32.is_gt(x1, x2)
119	      +a = Bool.pick(F32, gt, x2, x1)
120	      +b = Bool.pick(F32, gt, x1, x2)
121	      +delta = F32.sub(b, a)
122	      +rng_a = Rng.next_double(rng)
123	      +rand = Rng.last_double(rng_a)
124	      +inv = F32.div(1.0, F32.add(eta, 1.0))
125	      +neg = F32.neg(F32.add(eta, 1.0))
126	      +alpha1 = F32.sub(2.0, F32.pow(F32.add(1.0, F32.div(F32.mul(2.0, F32.sub(a, lo)), delta)), neg))
127	      +alpha2 = F32.sub(2.0, F32.pow(F32.add(1.0, F32.div(F32.mul(2.0, F32.sub(hi, b)), delta)), neg))
128	      +c1v = F32.clamp(F32.mul(0.5, F32.sub(F32.add(a, b), F32.mul(beta(rand, alpha1, inv), delta))), lo, hi)
129	      +c2v = F32.clamp(F32.mul(0.5, F32.add(F32.add(a, b), F32.mul(beta(rand, alpha2, inv), delta))), lo, hi)
130	      +rng_b = Rng.next_double(rng_a)
131	      +sw = F32.is_gt(Rng.last_double(rng_b), 0.5)
132	      +n1 = Bool.pick(F32, sw, c2v, c1v)
133	      +n2 = Bool.pick(F32, sw, c1v, c2v)
134	      +rng_c = advance(rng_b, f)
135	      loop.go(
136	        f,
137	        skip_kind(F32.is_gt(Rng.last_double(rng_c), 0.5), close_heads(rest1, rest2)),
138	        rest1, rest2, ivrest, eta, rng_c, n1 <> acc1, n2 <> acc2
139	      )
140	    case 1n++f k Nil{} b c:
141	      done(rng, acc1, acc2)
142	    case 1n++f k a <> as Nil{} c:
143	      done(rng, acc1, acc2)
144	    case 1n++f k a <> as b <> bs Nil{}:
145	      done(rng, acc1, acc2)
146	
147	def loop(
148	  fuel: Nat,
149	  skip: Bool,
150	  close: Bool,
151	  +eta: F32,
152	  xs1: List<&2, F32>,
153	  xs2: List<&2, F32>,
154	  ivs: List<&2, Bounds.Interval>,
155	  rng: Rng.State,
156	  acc1: List<&2, F32>,
157	  acc2: List<&2, F32>
158	) -> Rng.State & List<&2, F32> & List<&2, F32>:
159	  loop.go(fuel, skip_kind(skip, close), xs1, xs2, ivs, eta, rng, acc1, acc2)
160	
161	def start(
162	  +eta: F32,
163	  +x1s: List<&2, F32>,
164	  +x2s: List<&2, F32>,
165	  +box: Bounds.Box,
166	  rng: Rng.State,
167	  +m: Nat,
168	  n: Nat
169	) -> Rng.State & Ind.Individual & Ind.Individual:
170	  match n:
171	    case 0n:
172	      (rng, Ind.from_variables(x1s, m), Ind.from_variables(x2s, m))
173	    case 1n+f:
174	      +rng2 = Rng.next_double(rng)
175	      finish(
176	        m,
177	        loop(
178	          1n+f,
179	          F32.is_gt(Rng.last_double(rng2), 0.5),
180	          close_heads(x1s, x2s),
181	          eta, x1s, x2s, Bounds.intervals_of(box), rng2, Nil{}, Nil{}
182	        )
183	      )
184	
185	def maybe(
186	  +eta: F32,
187	  +x1s: List<&2, F32>,
188	  +x2s: List<&2, F32>,
189	  +box: Bounds.Box,
190	  rng: Rng.State,
191	  +m: Nat,
192	  +n: Nat,
193	  skip: Bool
194	) -> Rng.State & Ind.Individual & Ind.Individual:
195	  match skip:
196	    case True{}:
197	      (rng, Ind.from_variables(x1s, m), Ind.from_variables(x2s, m))
198	    case False{}:
199	      start(eta, x1s, x2s, box, rng, m, n)
200	
201	def draw(
202	  +eta: F32,
203	  +prob: F32,
204	  +x1s: List<&2, F32>,
205	  +x2s: List<&2, F32>,
206	  +box: Bounds.Box,
207	  rng: Rng.State,
208	  +m: Nat,
209	  +n: Nat
210	) -> Rng.State & Ind.Individual & Ind.Individual:
211	  +rng2 = Rng.next_double(rng)
212	  maybe(eta, x1s, x2s, box, rng2, m, n, F32.is_gt(Rng.last_double(rng2), prob))
213	
214	# p=0 copies parents before any NextDouble. Matches C# children when
215	# pair probability is 0 (every draw is > 0). Bit-zero is definitional
216	# on the 0.0 literal; F32.is_gt stays off this path.
217	def draw_or_copy(
218	  +eta: F32,
219	  +prob: F32,
220	  +x1s: List<&2, F32>,
221	  +x2s: List<&2, F32>,
222	  +box: Bounds.Box,
223	  rng: Rng.State,
224	  +m: Nat,
225	  +n: Nat,
226	  zero: Bool
227	) -> Rng.State & Ind.Individual & Ind.Individual:
228	  match zero:
229	    case True{}:
230	      (rng, Ind.from_variables(x1s, m), Ind.from_variables(x2s, m))
231	    case False{}:
232	      draw(eta, prob, x1s, x2s, box, rng, m, n)
233	
234	def crossover(
235	  +op: Sbx,
236	  +p1: Ind.Individual,
237	  +p2: Ind.Individual,
238	  +box: Bounds.Box,
239	  rng: Rng.State
240	) -> Rng.State & Ind.Individual & Ind.Individual:
241	  +x1s = Ind.variables(p1)
242	  +x2s = Ind.variables(p2)
243	  +m = Bounds.n_obj(box)
244	  +prob = probability(op)
245	  draw_or_copy(
246	    eta(op),
247	    prob,
248	    x1s,
249	    x2s,
250	    box,
251	    rng,
252	    m,
253	    Bounds.n_var(box),
254	    Util.f32_is_zero(prob)
255	  )