1	# Named problems — C# Unsga3.Problems.
2	# ZDT1 / ZDT2 / DTLZ2 stay the oracle problems. The rest of the unconstrained
3	# catalog (ZDT3 / ZDT4 / ZDT6, DTLZ1 / DTLZ3 / DTLZ4 / DTLZ7, Sphere / Ackley /
4	# Rosenbrock) matches the same C# Evaluate formulas. Minimize.
5	# Unconstrained (NumberOfConstraints = 0). No IGD in this module.
6	# `alpha` is DTLZ4's exponent (C# default 100). Other kinds ignore it.
7	
8	import Base
9	import ./individual.bend as Ind
10	import ./bounds.bend as Bounds
11	import ./util.bend as Util
12	
13	type Family is Data:
14	  Zdt1{}
15	  Zdt2{}
16	  Zdt3{}
17	  Zdt4{}
18	  Zdt6{}
19	  Dtlz1{}
20	  Dtlz2{}
21	  Dtlz3{}
22	  Dtlz4{}
23	  Dtlz7{}
24	  Sphere{}
25	  Ackley{}
26	  Rosenbrock{}
27	
28	type Problem is Data:
29	  MkProblem{
30	    kind: Family
31	    box: Bounds.Box
32	    alpha: F32
33	  }
34	
35	def box_of(p: Problem) -> Bounds.Box:
36	  match p:
37	    case MkProblem{kind, box, alpha}:
38	      box
39	
40	def kind_of(p: Problem) -> Family:
41	  match p:
42	    case MkProblem{kind, box, alpha}:
43	      kind
44	
45	def alpha_of(p: Problem) -> F32:
46	  match p:
47	    case MkProblem{kind, box, alpha}:
48	      alpha
49	
50	def n_obj(p: Problem) -> Nat:
51	  Bounds.n_obj(box_of(p))
52	
53	def n_var(p: Problem) -> Nat:
54	  Bounds.n_var(box_of(p))
55	
56	def mk_prob(kind: Family, box: Bounds.Box) -> Problem:
57	  MkProblem{kind, box, 0.0}
58	
59	# ZDT1 / ZDT2 / ZDT3: n variables (C# default 30), 2 objectives, x ∈ [0,1]^n.
60	def zdt1(n: Nat) -> Problem:
61	  mk_prob(Zdt1{}, Bounds.unit_n(n, 2n))
62	
63	def zdt2(n: Nat) -> Problem:
64	  mk_prob(Zdt2{}, Bounds.unit_n(n, 2n))
65	
66	def zdt3(n: Nat) -> Problem:
67	  mk_prob(Zdt3{}, Bounds.unit_n(n, 2n))
68	
69	# ZDT4: C# default n=10. x0 ∈ [0,1], x1.. ∈ [-5,5].
70	def zdt4_ivs(n: Nat) -> List<&2, Bounds.Interval>:
71	  match n:
72	    case 0n:
73	      Nil{}
74	    case 1n+p:
75	      Bounds.interval(0.0, 1.0) <> List.replicate(Bounds.Interval, p, Bounds.interval(F32.neg(5.0), 5.0))
76	
77	def zdt4(n: Nat) -> Problem:
78	  mk_prob(Zdt4{}, Bounds.box(2n, zdt4_ivs(n)))
79	
80	# ZDT6: C# default n=10, x ∈ [0,1]^n.
81	def zdt6(n: Nat) -> Problem:
82	  mk_prob(Zdt6{}, Bounds.unit_n(n, 2n))
83	
84	# n = M + k - 1 (C# Dtlz*Problem).
85	def dtlz_n(m: Nat, k: Nat) -> Nat:
86	  match m:
87	    case 0n:
88	      k
89	    case 1n+pm:
90	      match k:
91	        case 0n:
92	          pm
93	        case 1n+pk:
94	          Nat.add(pm, 1n+pk)
95	
96	# DTLZ2: C# default M=3, k=10 ⇒ n=12, x ∈ [0,1]^n.
97	def dtlz2(+m: Nat, k: Nat) -> Problem:
98	  mk_prob(Dtlz2{}, Bounds.unit_n(dtlz_n(m, k), m))
99	
100	# DTLZ1: C# default M=3, k=5 ⇒ n=7.
101	def dtlz1(+m: Nat, k: Nat) -> Problem:
102	  mk_prob(Dtlz1{}, Bounds.unit_n(dtlz_n(m, k), m))
103	
104	# DTLZ3: C# default M=3, k=10. DTLZ2 map, DTLZ1 g.
105	def dtlz3(+m: Nat, k: Nat) -> Problem:
106	  mk_prob(Dtlz3{}, Bounds.unit_n(dtlz_n(m, k), m))
107	
108	# DTLZ4: C# default α=100, k=10.
109	def dtlz4_alpha(+m: Nat, k: Nat, +alpha: F32) -> Problem:
110	  MkProblem{Dtlz4{}, Bounds.unit_n(dtlz_n(m, k), m), alpha}
111	
112	def dtlz4(+m: Nat, k: Nat) -> Problem:
113	  dtlz4_alpha(m, k, 100.0)
114	
115	# DTLZ7: C# default M=3, k=20 ⇒ n=22.
116	def dtlz7(+m: Nat, k: Nat) -> Problem:
117	  mk_prob(Dtlz7{}, Bounds.unit_n(dtlz_n(m, k), m))
118	
119	# Single-objective. M=1. C# defaults: Sphere n=10, Ackley n=30, Rosenbrock n=10.
120	def sphere(n: Nat) -> Problem:
121	  mk_prob(Sphere{}, Bounds.span_n(n, 1n, F32.neg(5.12), 5.12))
122	
123	def ackley(n: Nat) -> Problem:
124	  mk_prob(Ackley{}, Bounds.span_n(n, 1n, F32.neg(32.768), 32.768))
125	
126	def rosenbrock(n: Nat) -> Problem:
127	  mk_prob(Rosenbrock{}, Bounds.span_n(n, 1n, F32.neg(2.048), 2.048))
128	
129	def sum_from(xs: List<&2, F32>) -> F32:
130	  match xs:
131	    case Nil{}:
132	      0.0
133	    case h <> t:
134	      F32.add(h, sum_from(t))
135	
136	def sum_tail(xs: List<&2, F32>) -> F32:
137	  match xs:
138	    case Nil{}:
139	      0.0
140	    case h <> t:
141	      sum_from(t)
142	
143	def list_tail(xs: List<&2, F32>) -> List<&2, F32>:
144	  match xs:
145	    case Nil{}:
146	      Nil{}
147	    case h <> t:
148	      t
149	
150	def zdt_den(n: Nat) -> F32:
151	  match n:
152	    case 0n:
153	      1.0
154	    case 1n:
155	      1.0
156	    case 1n+p:
157	      F32.from_nat(p)
158	
159	# C# Zdt1Problem.GLinear: 1 + 9 * sum(x[1..]) / (n - 1).
160	def zdt_g(+x: List<&2, F32>) -> F32:
161	  +n = Util.list_len(F32, x)
162	  F32.add(1.0, F32.div(F32.mul(9.0, sum_tail(x)), zdt_den(n)))
163	
164	def zdt1_f(+x: List<&2, F32>) -> List<&2, F32>:
165	  +x0 = Util.get_f32(x, 0n)
166	  +g = zdt_g(x)
167	  x0 <> [F32.mul(g, F32.sub(1.0, F32.sqrt(F32.div(x0, g))))]
168	
169	def zdt2_f(+x: List<&2, F32>) -> List<&2, F32>:
170	  +x0 = Util.get_f32(x, 0n)
171	  +g = zdt_g(x)
172	  +r = F32.div(x0, g)
173	  x0 <> [F32.mul(g, F32.sub(1.0, F32.mul(r, r)))]
174	
175	# C# Zdt3Problem: g * (1 - sqrt(r) - r * sin(10 π x0)), r = x0/g.
176	def zdt3_f(+x: List<&2, F32>) -> List<&2, F32>:
177	  +x0 = Util.get_f32(x, 0n)
178	  +g = zdt_g(x)
179	  +r = F32.div(x0, g)
180	  +s = F32.sin(F32.mul(F32.mul(10.0, F32.pi()), x0))
181	  x0 <> [F32.mul(g, F32.sub(F32.sub(1.0, F32.sqrt(r)), F32.mul(r, s)))]
182	
183	def zdt4_term(+xi: F32) -> F32:
184	  F32.sub(F32.square(xi), F32.mul(10.0, F32.cos(F32.mul(F32.mul(4.0, F32.pi()), xi))))
185	
186	def sum_zdt4(xs: List<&2, F32>) -> F32:
187	  match xs:
188	    case Nil{}:
189	      0.0
190	    case h <> t:
191	      F32.add(zdt4_term(h), sum_zdt4(t))
192	
193	def pred(n: Nat) -> Nat:
194	  match n:
195	    case 0n:
196	      0n
197	    case 1n+p:
198	      p
199	
200	# C# Zdt4Problem: 1 + 10*(n-1) + Σ_{i≥1} (x_i² - 10 cos(4 π x_i)).
201	def zdt4_g(+x: List<&2, F32>) -> F32:
202	  +n = Util.list_len(F32, x)
203	  F32.add(F32.add(1.0, F32.mul(10.0, F32.from_nat(pred(n)))), sum_zdt4(list_tail(x)))
204	
205	def zdt4_f(+x: List<&2, F32>) -> List<&2, F32>:
206	  +x0 = Util.get_f32(x, 0n)
207	  +g = zdt4_g(x)
208	  x0 <> [F32.mul(g, F32.sub(1.0, F32.sqrt(F32.div(x0, g))))]
209	
210	# C# Zdt6Problem: biased density. g = 1 + 9 * (mean of x[1..])^0.25.
211	def zdt6_f(+x: List<&2, F32>) -> List<&2, F32>:
212	  +x0 = Util.get_f32(x, 0n)
213	  +wave = F32.pow(F32.sin(F32.mul(F32.mul(6.0, F32.pi()), x0)), 6.0)
214	  +f0 = F32.sub(1.0, F32.mul(F32.exp(F32.mul(F32.neg(4.0), x0)), wave))
215	  +n = Util.list_len(F32, x)
216	  +g = F32.add(1.0, F32.mul(9.0, F32.pow(F32.div(sum_tail(x), zdt_den(n)), 0.25)))
217	  +r = F32.div(f0, g)
218	  f0 <> [F32.mul(g, F32.sub(1.0, F32.mul(r, r)))]
219	
220	def sum_sq_half(xs: List<&2, F32>) -> F32:
221	  match xs:
222	    case Nil{}:
223	      0.0
224	    case h <> t:
225	      F32.add(F32.square(F32.sub(h, 0.5)), sum_sq_half(t))
226	
227	def tail_from_m(+x: List<&2, F32>, m: Nat) -> List<&2, F32>:
228	  List.drop(&2, F32, x, pred(m))
229	
230	# C# DtlzHelper.GDtlz2: sum (x[i] - 0.5)^2 for i >= M-1.
231	def dtlz2_g(+x: List<&2, F32>, m: Nat) -> F32:
232	  sum_sq_half(tail_from_m(x, m))
233	
234	def dtlz1_term(+xi: F32) -> F32:
235	  +d = F32.sub(xi, 0.5)
236	  F32.sub(F32.square(d), F32.cos(F32.mul(F32.mul(20.0, F32.pi()), d)))
237	
238	def sum_dtlz1(xs: List<&2, F32>) -> F32:
239	  match xs:
240	    case Nil{}:
241	      0.0
242	    case h <> t:
243	      F32.add(dtlz1_term(h), sum_dtlz1(t))
244	
245	# C# DtlzHelper.GDtlz1: 100 * ((n - M + 1) + Σ (d² - cos(20 π d))), d = x-0.5.
246	# On a legal k≥1 box the tail length is n-M+1, which is that coefficient.
247	def dtlz1_g(+x: List<&2, F32>, m: Nat) -> F32:
248	  +tail = tail_from_m(x, m)
249	  F32.mul(100.0, F32.add(F32.from_nat(Util.list_len(F32, tail)), sum_dtlz1(tail)))
250	
251	def half_pi() -> F32:
252	  F32.div(F32.pi(), 2.0)
253	
254	def prod_cos.go(+x: List<&2, F32>, fuel: Nat, +i: Nat, acc: F32) -> F32:
255	  match fuel:
256	    case 0n:
257	      acc
258	    case 1n+f:
259	      prod_cos.go(
260	        x,
261	        f,
262	        1n+i,
263	        F32.mul(acc, F32.cos(F32.mul(Util.get_f32(x, i), half_pi())))
264	      )
265	
266	def prod_cos(+x: List<&2, F32>, n: Nat) -> F32:
267	  prod_cos.go(x, n, 0n, 1.0)
268	
269	def dtlz2_sin(+x: List<&2, F32>, +m: Nat, +i: Nat) -> F32:
270	  F32.sin(F32.mul(Util.get_f32(x, pred(Nat.sub(m, i))), half_pi()))
271	
272	def dtlz2_fi.tail(+x: List<&2, F32>, +m: Nat, +g: F32, +i: Nat, +body: F32, first: Bool) -> F32:
273	  match first:
274	    case True{}:
275	      body
276	    case False{}:
277	      F32.mul(body, dtlz2_sin(x, m, i))
278	
279	# C# Dtlz2Problem: (1+g) * prod_{j<M-i-1} cos(x[j] π/2) [* sin(x[M-i-1] π/2) if i>0].
280	def dtlz2_fi(+x: List<&2, F32>, +m: Nat, +g: F32, +i: Nat) -> F32:
281	  dtlz2_fi.tail(
282	    x,
283	    m,
284	    g,
285	    i,
286	    F32.mul(F32.add(1.0, g), prod_cos(x, pred(Nat.sub(m, i)))),
287	    Nat.is_eq(i, 0n)
288	  )
289	
290	def dtlz2_fs.go(+x: List<&2, F32>, +m: Nat, +g: F32, fuel: Nat, +i: Nat) -> List<&2, F32>:
291	  match fuel:
292	    case 0n:
293	      Nil{}
294	    case 1n+f:
295	      dtlz2_fi(x, m, g, i) <> dtlz2_fs.go(x, m, g, f, 1n+i)
296	
297	def dtlz2_f(+x: List<&2, F32>, +m: Nat) -> List<&2, F32>:
298	  dtlz2_fs.go(x, m, dtlz2_g(x, m), m, 0n)
299	
300	def prod_xs.go(+x: List<&2, F32>, fuel: Nat, +i: Nat, acc: F32) -> F32:
301	  match fuel:
302	    case 0n:
303	      acc
304	    case 1n+f:
305	      prod_xs.go(x, f, 1n+i, F32.mul(acc, Util.get_f32(x, i)))
306	
307	def prod_xs(+x: List<&2, F32>, n: Nat) -> F32:
308	  prod_xs.go(x, n, 0n, 1.0)
309	
310	def dtlz1_scale(+body: F32, +x: List<&2, F32>, +idx: Nat, first: Bool) -> F32:
311	  match first:
312	    case True{}:
313	      body
314	    case False{}:
315	      F32.mul(body, F32.sub(1.0, Util.get_f32(x, idx)))
316	
317	# C# Dtlz1Problem: 0.5*(1+g) * prod_{j<M-i-1} x[j] [* (1-x[M-i-1]) if i>0].
318	def dtlz1_fi(+x: List<&2, F32>, +m: Nat, +g: F32, +i: Nat) -> F32:
319	  +idx = pred(Nat.sub(m, i))
320	  dtlz1_scale(F32.mul(F32.mul(0.5, F32.add(1.0, g)), prod_xs(x, idx)), x, idx, Nat.is_eq(i, 0n))
321	
322	def dtlz1_fs.go(+x: List<&2, F32>, +m: Nat, +g: F32, fuel: Nat, +i: Nat) -> List<&2, F32>:
323	  match fuel:
324	    case 0n:
325	      Nil{}
326	    case 1n+f:
327	      dtlz1_fi(x, m, g, i) <> dtlz1_fs.go(x, m, g, f, 1n+i)
328	
329	def dtlz1_f(+x: List<&2, F32>, +m: Nat) -> List<&2, F32>:
330	  dtlz1_fs.go(x, m, dtlz1_g(x, m), m, 0n)
331	
332	# C# Dtlz3Problem: DTLZ2 map with DTLZ1's g.
333	def dtlz3_f(+x: List<&2, F32>, +m: Nat) -> List<&2, F32>:
334	  dtlz2_fs.go(x, m, dtlz1_g(x, m), m, 0n)
335	
336	# Fuel shrinks first so the recursive call terminates. The g-tail stays raw.
337	def map_prefix_pow.go(fuel: Nat, xs: List<&2, F32>, +alpha: F32) -> List<&2, F32>:
338	  match fuel:
339	    case 0n:
340	      xs
341	    case 1n+f:
342	      match xs:
343	        case Nil{}:
344	          Nil{}
345	        case h <> t:
346	          F32.pow(h, alpha) <> map_prefix_pow.go(f, t, alpha)
347	
348	# Raise the first `fuel` decision variables to α.
349	def map_prefix_pow(+xs: List<&2, F32>, fuel: Nat, +alpha: F32) -> List<&2, F32>:
350	  map_prefix_pow.go(fuel, xs, alpha)
351	
352	# C# Dtlz4Problem: DTLZ2 with x^α inside the sine/cosine map. g is still GDtlz2.
353	def dtlz4_f(+x: List<&2, F32>, +m: Nat, +alpha: F32) -> List<&2, F32>:
354	  dtlz2_f(map_prefix_pow(x, pred(m), alpha), m)
355	
356	def pos_or.fin(+x: F32, +fallback: F32, ok: Bool) -> F32:
357	  match ok:
358	    case True{}:
359	      x
360	    case False{}:
361	      fallback
362	
363	def pos_or(+x: F32, +fallback: F32) -> F32:
364	  pos_or.fin(x, fallback, F32.is_gt(x, 0.0))
365	
366	# C# Dtlz7Problem g: 1 + 9 * sum(x[M-1..]) / (n - M + 1).
367	def dtlz7_g(+x: List<&2, F32>, m: Nat) -> F32:
368	  +tail = tail_from_m(x, m)
369	  +den = pos_or(F32.from_nat(Util.list_len(F32, tail)), 1.0)
370	  F32.add(1.0, F32.div(F32.mul(9.0, sum_from(tail)), den))
371	
372	def dtlz7_h.go(+x: List<&2, F32>, +g1: F32, fuel: Nat, +i: Nat, acc: F32) -> F32:
373	  match fuel:
374	    case 0n:
375	      acc
376	    case 1n+f:
377	      +fi = Util.get_f32(x, i)
378	      +term = F32.mul(F32.div(fi, g1), F32.add(1.0, F32.sin(F32.mul(F32.mul(3.0, F32.pi()), fi))))
379	      dtlz7_h.go(x, g1, f, 1n+i, F32.add(acc, term))
380	
381	def dtlz7_prefix(+x: List<&2, F32>, fuel: Nat, +i: Nat) -> List<&2, F32>:
382	  match fuel:
383	    case 0n:
384	      Nil{}
385	    case 1n+f:
386	      Util.get_f32(x, i) <> dtlz7_prefix(x, f, 1n+i)
387	
388	# f[i]=x[i] for i<M-1; f[M-1]=(1+g)*h; h = M - Σ f_i/(1+g) * (1+sin(3 π f_i)).
389	def dtlz7_f(+x: List<&2, F32>, +m: Nat) -> List<&2, F32>:
390	  match m:
391	    case 0n:
392	      Nil{}
393	    case 1n+p:
394	      +g = dtlz7_g(x, m)
395	      +g1 = F32.add(1.0, g)
396	      +h = F32.sub(F32.from_nat(m), dtlz7_h.go(x, g1, p, 0n, 0.0))
397	      List.append(&2, F32, dtlz7_prefix(x, p, 0n), [F32.mul(g1, h)])
398	
399	def sum_sq(xs: List<&2, F32>) -> F32:
400	  match xs:
401	    case Nil{}:
402	      0.0
403	    case h <> t:
404	      F32.add(F32.square(h), sum_sq(t))
405	
406	def sphere_f(+x: List<&2, F32>) -> List<&2, F32>:
407	  [sum_sq(x)]
408	
409	type Acc2 is Data:
410	  Acc2{
411	    s1: F32
412	    s2: F32
413	  }
414	
415	def ackley_go(xs: List<&2, F32>, +s1: F32, +s2: F32) -> Acc2:
416	  match xs:
417	    case Nil{}:
418	      Acc2{s1, s2}
419	    case +h <> t:
420	      ackley_go(
421	        t,
422	        F32.add(s1, F32.square(h)),
423	        F32.add(s2, F32.cos(F32.mul(F32.mul(2.0, F32.pi()), h)))
424	      )
425	
426	def ackley_den(n: Nat) -> F32:
427	  match n:
428	    case 0n:
429	      1.0
430	    case 1n+p:
431	      F32.from_nat(1n+p)
432	
433	def ackley_from(+n: Nat, acc: Acc2) -> List<&2, F32>:
434	  match acc:
435	    case Acc2{s1, s2}:
436	      +den = ackley_den(n)
437	      [F32.add(
438	        F32.add(
439	          F32.mul(F32.neg(20.0), F32.exp(F32.mul(F32.neg(0.2), F32.sqrt(F32.div(s1, den))))),
440	          F32.neg(F32.exp(F32.div(s2, den)))
441	        ),
442	        F32.add(20.0, F32.exp(1.0))
443	      )]
444	
445	# C# AckleyProblem. n=0 stays total (denominator 1) instead of dividing by zero.
446	def ackley_f(+x: List<&2, F32>) -> List<&2, F32>:
447	  ackley_from(Util.list_len(F32, x), ackley_go(x, 0.0, 0.0))
448	
449	def rosen_step(rest: List<&2, F32>, +a: F32) -> F32:
450	  match rest:
451	    case Nil{}:
452	      0.0
453	    case +b <> t:
454	      +diff = F32.sub(b, F32.square(a))
455	      +one = F32.sub(1.0, a)
456	      F32.add(F32.add(F32.mul(100.0, F32.square(diff)), F32.square(one)), rosen_step(t, b))
457	
458	def rosenbrock_f(+x: List<&2, F32>) -> List<&2, F32>:
459	  match x:
460	    case Nil{}:
461	      [0.0]
462	    case a <> rest:
463	      [rosen_step(rest, a)]
464	
465	def eval_kind(kind: Family, +x: List<&2, F32>, +m: Nat, +alpha: F32) -> List<&2, F32>:
466	  match kind:
467	    case Zdt1{}:
468	      zdt1_f(x)
469	    case Zdt2{}:
470	      zdt2_f(x)
471	    case Zdt3{}:
472	      zdt3_f(x)
473	    case Zdt4{}:
474	      zdt4_f(x)
475	    case Zdt6{}:
476	      zdt6_f(x)
477	    case Dtlz1{}:
478	      dtlz1_f(x, m)
479	    case Dtlz2{}:
480	      dtlz2_f(x, m)
481	    case Dtlz3{}:
482	      dtlz3_f(x, m)
483	    case Dtlz4{}:
484	      dtlz4_f(x, m, alpha)
485	    case Dtlz7{}:
486	      dtlz7_f(x, m)
487	    case Sphere{}:
488	      sphere_f(x)
489	    case Ackley{}:
490	      ackley_f(x)
491	    case Rosenbrock{}:
492	      rosenbrock_f(x)
493	
494	def evaluate(+p: Problem, +ind: Ind.Individual) -> Ind.Individual:
495	  Ind.set_objectives(ind, eval_kind(kind_of(p), Ind.variables(ind), n_obj(p), alpha_of(p)))
496	
497	# Independent per individual. Mid-split + parallel recurse; same order as cons.
498	def evaluate_all(+p: Problem, +xs: List<&2, Ind.Individual>) -> List<&2, Ind.Individual>:
499	  Util.map_halves(~Ind.Individual, ~Ind.Individual, ~Problem, ~evaluate, p, xs)