1use oxilean_kernel::{BinderInfo, Declaration, Environment, Expr, Level, Name};
6
7use super::types::{
8 BipartiteMatchingGraph, BranchBoundData, CuttingPlane, FacilityLocation, FlowNetwork,
9 FlowNetworkSpec, GraphColoring, KnapsackSolver, LPRelaxation, MatchingProblem,
10 MatroidIntersection, SetCoverData, ShortestPath, SpanningTree, SteinerTree, TravelingSalesman,
11 VehicleRouting,
12};
13
14pub fn app(f: Expr, a: Expr) -> Expr {
15 Expr::App(Box::new(f), Box::new(a))
16}
17pub fn app2(f: Expr, a: Expr, b: Expr) -> Expr {
18 app(app(f, a), b)
19}
20pub fn app3(f: Expr, a: Expr, b: Expr, c: Expr) -> Expr {
21 app(app2(f, a, b), c)
22}
23pub fn cst(s: &str) -> Expr {
24 Expr::Const(Name::str(s), vec![])
25}
26pub fn prop() -> Expr {
27 Expr::Sort(Level::zero())
28}
29pub fn type0() -> Expr {
30 Expr::Sort(Level::succ(Level::zero()))
31}
32pub fn pi(bi: BinderInfo, name: &str, dom: Expr, body: Expr) -> Expr {
33 Expr::Pi(bi, Name::str(name), Box::new(dom), Box::new(body))
34}
35pub fn arrow(a: Expr, b: Expr) -> Expr {
36 pi(BinderInfo::Default, "_", a, b)
37}
38pub fn bvar(n: u32) -> Expr {
39 Expr::BVar(n)
40}
41pub fn nat_ty() -> Expr {
42 cst("Nat")
43}
44pub fn bool_ty() -> Expr {
45 cst("Bool")
46}
47pub fn real_ty() -> Expr {
48 cst("Real")
49}
50pub fn int_ty() -> Expr {
51 cst("Int")
52}
53pub fn list_ty(elem: Expr) -> Expr {
54 app(cst("List"), elem)
55}
56pub fn pair_ty(a: Expr, b: Expr) -> Expr {
57 app2(cst("Prod"), a, b)
58}
59pub fn option_ty(a: Expr) -> Expr {
60 app(cst("Option"), a)
61}
62pub fn graph_ty() -> Expr {
64 type0()
65}
66pub fn bipartite_graph_ty() -> Expr {
68 type0()
69}
70pub fn matching_ty() -> Expr {
72 arrow(graph_ty(), type0())
73}
74pub fn perfect_matching_ty() -> Expr {
76 arrow(graph_ty(), prop())
77}
78pub fn max_matching_ty() -> Expr {
80 arrow(graph_ty(), arrow(matching_ty(), prop()))
81}
82pub fn alternating_path_ty() -> Expr {
84 arrow(
85 graph_ty(),
86 arrow(matching_ty(), arrow(list_ty(nat_ty()), prop())),
87 )
88}
89pub fn augmenting_path_ty() -> Expr {
91 arrow(
92 graph_ty(),
93 arrow(matching_ty(), arrow(list_ty(nat_ty()), prop())),
94 )
95}
96pub fn blossom_ty() -> Expr {
98 arrow(
99 graph_ty(),
100 arrow(matching_ty(), arrow(list_ty(nat_ty()), prop())),
101 )
102}
103pub fn hall_condition_ty() -> Expr {
106 arrow(bipartite_graph_ty(), prop())
107}
108pub fn hall_theorem_ty() -> Expr {
110 pi(
111 BinderInfo::Default,
112 "G",
113 bipartite_graph_ty(),
114 app2(
115 cst("Iff"),
116 app(cst("HallCondition"), bvar(0)),
117 app(cst("PerfectMatching"), bvar(0)),
118 ),
119 )
120}
121pub fn konig_theorem_ty() -> Expr {
123 pi(
124 BinderInfo::Default,
125 "G",
126 bipartite_graph_ty(),
127 app2(
128 cst("Eq"),
129 app(cst("MaxMatchingSize"), bvar(0)),
130 app(cst("MinVertexCoverSize"), bvar(0)),
131 ),
132 )
133}
134pub fn tutte_condition_ty() -> Expr {
137 arrow(graph_ty(), prop())
138}
139pub fn tutte_theorem_ty() -> Expr {
141 pi(
142 BinderInfo::Default,
143 "G",
144 graph_ty(),
145 app2(
146 cst("Iff"),
147 app(cst("PerfectMatching"), bvar(0)),
148 app(cst("TutteCondition"), bvar(0)),
149 ),
150 )
151}
152pub fn berge_theorem_ty() -> Expr {
154 pi(
155 BinderInfo::Default,
156 "G",
157 graph_ty(),
158 pi(
159 BinderInfo::Default,
160 "M",
161 matching_ty(),
162 app2(
163 cst("Iff"),
164 app2(cst("MaxMatching"), bvar(1), bvar(0)),
165 app(
166 cst("NoAugmentingPath"),
167 app2(cst("mk_pair"), bvar(1), bvar(0)),
168 ),
169 ),
170 ),
171 )
172}
173pub fn flow_network_ty() -> Expr {
175 type0()
176}
177pub fn flow_ty() -> Expr {
179 arrow(flow_network_ty(), type0())
180}
181pub fn flow_value_ty() -> Expr {
183 arrow(flow_network_ty(), arrow(flow_ty(), real_ty()))
184}
185pub fn max_flow_ty() -> Expr {
187 arrow(flow_network_ty(), real_ty())
188}
189pub fn cut_ty() -> Expr {
191 arrow(flow_network_ty(), arrow(arrow(nat_ty(), bool_ty()), prop()))
192}
193pub fn cut_capacity_ty() -> Expr {
195 arrow(
196 flow_network_ty(),
197 arrow(arrow(nat_ty(), bool_ty()), real_ty()),
198 )
199}
200pub fn min_cut_ty() -> Expr {
202 arrow(flow_network_ty(), real_ty())
203}
204pub fn max_flow_min_cut_ty() -> Expr {
206 pi(
207 BinderInfo::Default,
208 "N",
209 flow_network_ty(),
210 app2(
211 cst("Eq"),
212 app(cst("MaxFlow"), bvar(0)),
213 app(cst("MinCut"), bvar(0)),
214 ),
215 )
216}
217pub fn ford_fulkerson_termination_ty() -> Expr {
219 pi(
220 BinderInfo::Default,
221 "N",
222 flow_network_ty(),
223 arrow(
224 app(cst("IntegerCapacities"), bvar(0)),
225 app(cst("FordFulkersonTerminates"), bvar(0)),
226 ),
227 )
228}
229pub fn dinics_complexity_ty() -> Expr {
231 pi(
232 BinderInfo::Default,
233 "N",
234 flow_network_ty(),
235 app(cst("DinicsTimeBound"), bvar(0)),
236 )
237}
238pub fn cost_matrix_ty() -> Expr {
240 arrow(nat_ty(), arrow(nat_ty(), real_ty()))
241}
242pub fn assignment_ty() -> Expr {
244 arrow(nat_ty(), arrow(arrow(nat_ty(), nat_ty()), prop()))
245}
246pub fn optimal_assignment_ty() -> Expr {
248 arrow(cost_matrix_ty(), arrow(arrow(nat_ty(), nat_ty()), prop()))
249}
250pub fn hungarian_algorithm_ty() -> Expr {
252 arrow(cost_matrix_ty(), arrow(nat_ty(), nat_ty()))
253}
254pub fn hungarian_correctness_ty() -> Expr {
256 pi(
257 BinderInfo::Default,
258 "C",
259 cost_matrix_ty(),
260 app2(
261 cst("OptimalAssignment"),
262 bvar(0),
263 app(cst("HungarianAlgorithm"), bvar(0)),
264 ),
265 )
266}
267pub fn tsp_instance_ty() -> Expr {
269 type0()
270}
271pub fn tsp_tour_ty() -> Expr {
273 arrow(tsp_instance_ty(), arrow(list_ty(nat_ty()), prop()))
274}
275pub fn tsp_optimal_ty() -> Expr {
277 arrow(tsp_instance_ty(), real_ty())
278}
279pub fn held_karp_bound_ty() -> Expr {
281 arrow(tsp_instance_ty(), real_ty())
282}
283pub fn held_karp_lower_bound_ty() -> Expr {
285 pi(
286 BinderInfo::Default,
287 "I",
288 tsp_instance_ty(),
289 app2(
290 cst("Le"),
291 app(cst("HeldKarpBound"), bvar(0)),
292 app(cst("TSPOptimal"), bvar(0)),
293 ),
294 )
295}
296pub fn christofides_approximation_ty() -> Expr {
298 pi(
299 BinderInfo::Default,
300 "I",
301 tsp_instance_ty(),
302 arrow(
303 app(cst("MetricTSP"), bvar(0)),
304 app2(
305 cst("Le"),
306 app(cst("ChristofidesValue"), bvar(0)),
307 app2(
308 cst("RealMul"),
309 cst("ThreeHalves"),
310 app(cst("TSPOptimal"), bvar(0)),
311 ),
312 ),
313 ),
314 )
315}
316pub fn vertex_cover_approx_ty() -> Expr {
318 pi(
319 BinderInfo::Default,
320 "G",
321 graph_ty(),
322 app2(
323 cst("Le"),
324 app(cst("ApproxVertexCoverSize"), bvar(0)),
325 app2(
326 cst("NatMul"),
327 cst("Nat.two"),
328 app(cst("MinVertexCoverSize"), bvar(0)),
329 ),
330 ),
331 )
332}
333pub fn set_cover_approx_ty() -> Expr {
335 pi(
336 BinderInfo::Default,
337 "n",
338 nat_ty(),
339 app2(
340 cst("Le"),
341 app(cst("GreedySetCoverSize"), bvar(0)),
342 app2(
343 cst("NatMul"),
344 app(cst("HarmonicNumber"), bvar(0)),
345 app(cst("OptSetCoverSize"), bvar(0)),
346 ),
347 ),
348 )
349}
350pub fn matroid_ty() -> Expr {
352 type0()
353}
354pub fn independent_set_ty() -> Expr {
356 arrow(matroid_ty(), arrow(list_ty(nat_ty()), prop()))
357}
358pub fn matroid_base_ty() -> Expr {
360 arrow(matroid_ty(), arrow(list_ty(nat_ty()), prop()))
361}
362pub fn matroid_rank_ty() -> Expr {
364 arrow(matroid_ty(), nat_ty())
365}
366pub fn greedy_optimality_ty() -> Expr {
368 pi(
369 BinderInfo::Default,
370 "M",
371 matroid_ty(),
372 pi(
373 BinderInfo::Default,
374 "w",
375 arrow(nat_ty(), real_ty()),
376 app2(cst("GreedyIsOptimal"), bvar(1), bvar(0)),
377 ),
378 )
379}
380pub fn matroid_intersection_ty() -> Expr {
382 arrow(matroid_ty(), arrow(matroid_ty(), list_ty(nat_ty())))
383}
384pub fn matroid_intersection_optimality_ty() -> Expr {
386 pi(
387 BinderInfo::Default,
388 "M1",
389 matroid_ty(),
390 pi(
391 BinderInfo::Default,
392 "M2",
393 matroid_ty(),
394 pi(
395 BinderInfo::Default,
396 "w",
397 arrow(nat_ty(), real_ty()),
398 app3(cst("MatroidIntersectionOptimal"), bvar(2), bvar(1), bvar(0)),
399 ),
400 ),
401 )
402}
403pub fn graphic_matroid_ty() -> Expr {
405 arrow(graph_ty(), matroid_ty())
406}
407pub fn uniform_matroid_ty() -> Expr {
409 arrow(nat_ty(), arrow(nat_ty(), matroid_ty()))
410}
411pub fn submodular_function_ty() -> Expr {
414 arrow(arrow(list_ty(nat_ty()), real_ty()), prop())
415}
416pub fn submodular_maximization_ty() -> Expr {
418 arrow(arrow(list_ty(nat_ty()), real_ty()), list_ty(nat_ty()))
419}
420pub fn submodular_greedy_approx_ty() -> Expr {
422 pi(
423 BinderInfo::Default,
424 "f",
425 arrow(list_ty(nat_ty()), real_ty()),
426 arrow(
427 app(cst("MonotoneSubmodular"), bvar(0)),
428 app2(
429 cst("Le"),
430 app2(
431 cst("RealMul"),
432 cst("OneHalf"),
433 app(cst("SubmodularOpt"), bvar(0)),
434 ),
435 app(cst("SubmodularGreedyValue"), bvar(0)),
436 ),
437 ),
438 )
439}
440pub fn supermodular_function_ty() -> Expr {
442 arrow(arrow(list_ty(nat_ty()), real_ty()), prop())
443}
444pub fn polymatroid_rank_ty() -> Expr {
446 arrow(arrow(list_ty(nat_ty()), real_ty()), prop())
447}
448pub fn ilp_instance_ty() -> Expr {
450 type0()
451}
452pub fn ilp_solution_ty() -> Expr {
454 arrow(ilp_instance_ty(), arrow(list_ty(int_ty()), prop()))
455}
456pub fn ilp_optimal_ty() -> Expr {
458 arrow(ilp_instance_ty(), arrow(list_ty(int_ty()), prop()))
459}
460pub fn gomory_cut_ty() -> Expr {
462 arrow(ilp_instance_ty(), ilp_instance_ty())
463}
464pub fn branch_and_bound_ty() -> Expr {
466 arrow(ilp_instance_ty(), list_ty(int_ty()))
467}
468pub fn lp_relaxation_ty() -> Expr {
470 arrow(ilp_instance_ty(), real_ty())
471}
472pub fn lp_relaxation_lower_bound_ty() -> Expr {
474 pi(
475 BinderInfo::Default,
476 "P",
477 ilp_instance_ty(),
478 app2(
479 cst("Le"),
480 app(cst("LPRelaxation"), bvar(0)),
481 app(cst("ILPOptimalValue"), bvar(0)),
482 ),
483 )
484}
485pub fn polytope_ty() -> Expr {
487 type0()
488}
489pub fn polytope_vertex_ty() -> Expr {
491 arrow(polytope_ty(), arrow(list_ty(real_ty()), prop()))
492}
493pub fn totally_unimodular_ty() -> Expr {
495 arrow(arrow(nat_ty(), arrow(nat_ty(), int_ty())), prop())
496}
497pub fn tu_integral_polyhedra_ty() -> Expr {
499 pi(
500 BinderInfo::Default,
501 "A",
502 arrow(nat_ty(), arrow(nat_ty(), int_ty())),
503 arrow(
504 app(cst("TotallyUnimodular"), bvar(0)),
505 app(cst("IntegralPolyhedron"), bvar(0)),
506 ),
507 )
508}
509pub fn bipartite_incidence_tu_ty() -> Expr {
511 pi(
512 BinderInfo::Default,
513 "G",
514 bipartite_graph_ty(),
515 app(
516 cst("TotallyUnimodular"),
517 app(cst("IncidenceMatrix"), bvar(0)),
518 ),
519 )
520}
521pub fn facet_defining_ty() -> Expr {
523 arrow(
524 polytope_ty(),
525 arrow(arrow(list_ty(real_ty()), prop()), prop()),
526 )
527}
528pub fn weak_duality_ty() -> Expr {
530 pi(
531 BinderInfo::Default,
532 "P",
533 ilp_instance_ty(),
534 app2(
535 cst("Le"),
536 app(cst("PrimalObjective"), bvar(0)),
537 app(cst("DualObjective"), bvar(0)),
538 ),
539 )
540}
541pub fn strong_duality_ty() -> Expr {
543 pi(
544 BinderInfo::Default,
545 "P",
546 ilp_instance_ty(),
547 arrow(
548 app(cst("BothFeasible"), bvar(0)),
549 app2(
550 cst("Eq"),
551 app(cst("PrimalObjective"), bvar(0)),
552 app(cst("DualObjective"), bvar(0)),
553 ),
554 ),
555 )
556}
557pub fn build_combinatorial_optimization_env(env: &mut Environment) -> Result<(), String> {
559 let axioms: &[(&str, Expr)] = &[
560 ("Graph", graph_ty()),
561 ("BipartiteGraph", bipartite_graph_ty()),
562 ("Matching", matching_ty()),
563 ("PerfectMatching", perfect_matching_ty()),
564 ("MaxMatching", max_matching_ty()),
565 ("AlternatingPath", alternating_path_ty()),
566 ("AugmentingPath", augmenting_path_ty()),
567 ("Blossom", blossom_ty()),
568 ("MaxMatchingSize", arrow(graph_ty(), nat_ty())),
569 ("MinVertexCoverSize", arrow(graph_ty(), nat_ty())),
570 ("ApproxVertexCoverSize", arrow(graph_ty(), nat_ty())),
571 (
572 "NoAugmentingPath",
573 arrow(pair_ty(graph_ty(), matching_ty()), prop()),
574 ),
575 (
576 "mk_pair",
577 arrow(
578 graph_ty(),
579 arrow(matching_ty(), pair_ty(graph_ty(), matching_ty())),
580 ),
581 ),
582 ("HallCondition", hall_condition_ty()),
583 ("hall_theorem", hall_theorem_ty()),
584 ("konig_theorem", konig_theorem_ty()),
585 ("TutteCondition", tutte_condition_ty()),
586 ("tutte_theorem", tutte_theorem_ty()),
587 ("berge_theorem", berge_theorem_ty()),
588 ("FlowNetwork", flow_network_ty()),
589 ("Flow", flow_ty()),
590 ("FlowValue", flow_value_ty()),
591 ("MaxFlow", max_flow_ty()),
592 ("Cut", cut_ty()),
593 ("CutCapacity", cut_capacity_ty()),
594 ("MinCut", min_cut_ty()),
595 ("IntegerCapacities", arrow(flow_network_ty(), prop())),
596 ("FordFulkersonTerminates", arrow(flow_network_ty(), prop())),
597 ("DinicsTimeBound", arrow(flow_network_ty(), prop())),
598 ("max_flow_min_cut", max_flow_min_cut_ty()),
599 (
600 "ford_fulkerson_termination",
601 ford_fulkerson_termination_ty(),
602 ),
603 ("dinics_complexity", dinics_complexity_ty()),
604 ("CostMatrix", cost_matrix_ty()),
605 ("Assignment", assignment_ty()),
606 ("OptimalAssignment", optimal_assignment_ty()),
607 ("HungarianAlgorithm", hungarian_algorithm_ty()),
608 ("hungarian_correctness", hungarian_correctness_ty()),
609 ("TSPInstance", tsp_instance_ty()),
610 ("TSPTour", tsp_tour_ty()),
611 ("TSPOptimal", tsp_optimal_ty()),
612 ("HeldKarpBound", held_karp_bound_ty()),
613 ("MetricTSP", arrow(tsp_instance_ty(), prop())),
614 ("ChristofidesValue", arrow(tsp_instance_ty(), real_ty())),
615 ("ThreeHalves", real_ty()),
616 ("OneHalf", real_ty()),
617 ("RealMul", arrow(real_ty(), arrow(real_ty(), real_ty()))),
618 ("NatMul", arrow(nat_ty(), arrow(nat_ty(), nat_ty()))),
619 ("Nat.two", nat_ty()),
620 ("HarmonicNumber", arrow(nat_ty(), nat_ty())),
621 ("OptSetCoverSize", arrow(nat_ty(), nat_ty())),
622 ("GreedySetCoverSize", arrow(nat_ty(), nat_ty())),
623 ("held_karp_lower_bound", held_karp_lower_bound_ty()),
624 (
625 "christofides_approximation",
626 christofides_approximation_ty(),
627 ),
628 ("vertex_cover_approximation", vertex_cover_approx_ty()),
629 ("set_cover_approximation", set_cover_approx_ty()),
630 ("Matroid", matroid_ty()),
631 ("IndependentSet", independent_set_ty()),
632 ("MatroidBase", matroid_base_ty()),
633 ("MatroidRank", matroid_rank_ty()),
634 (
635 "GreedyIsOptimal",
636 arrow(matroid_ty(), arrow(arrow(nat_ty(), real_ty()), prop())),
637 ),
638 ("MatroidIntersection", matroid_intersection_ty()),
639 (
640 "MatroidIntersectionOptimal",
641 arrow(
642 matroid_ty(),
643 arrow(matroid_ty(), arrow(arrow(nat_ty(), real_ty()), prop())),
644 ),
645 ),
646 ("GraphicMatroid", graphic_matroid_ty()),
647 ("UniformMatroid", uniform_matroid_ty()),
648 ("greedy_optimality", greedy_optimality_ty()),
649 (
650 "matroid_intersection_optimality",
651 matroid_intersection_optimality_ty(),
652 ),
653 ("SubmodularFunction", submodular_function_ty()),
654 ("SubmodularMaximization", submodular_maximization_ty()),
655 (
656 "MonotoneSubmodular",
657 arrow(arrow(list_ty(nat_ty()), real_ty()), prop()),
658 ),
659 (
660 "SubmodularOpt",
661 arrow(arrow(list_ty(nat_ty()), real_ty()), real_ty()),
662 ),
663 (
664 "SubmodularGreedyValue",
665 arrow(arrow(list_ty(nat_ty()), real_ty()), real_ty()),
666 ),
667 ("SupermodularFunction", supermodular_function_ty()),
668 ("PolymatroidRankFunction", polymatroid_rank_ty()),
669 ("submodular_greedy_approx", submodular_greedy_approx_ty()),
670 ("ILPInstance", ilp_instance_ty()),
671 ("ILPSolution", ilp_solution_ty()),
672 ("ILPOptimal", ilp_optimal_ty()),
673 ("GomoryCut", gomory_cut_ty()),
674 ("BranchAndBound", branch_and_bound_ty()),
675 ("LPRelaxation", lp_relaxation_ty()),
676 ("ILPOptimalValue", arrow(ilp_instance_ty(), real_ty())),
677 ("PrimalObjective", arrow(ilp_instance_ty(), real_ty())),
678 ("DualObjective", arrow(ilp_instance_ty(), real_ty())),
679 ("BothFeasible", arrow(ilp_instance_ty(), prop())),
680 ("lp_relaxation_lower_bound", lp_relaxation_lower_bound_ty()),
681 ("weak_duality", weak_duality_ty()),
682 ("strong_duality", strong_duality_ty()),
683 ("Polytope", polytope_ty()),
684 ("PolytopeVertex", polytope_vertex_ty()),
685 ("TotallyUnimodular", totally_unimodular_ty()),
686 (
687 "IntegralPolyhedron",
688 arrow(arrow(nat_ty(), arrow(nat_ty(), int_ty())), prop()),
689 ),
690 (
691 "IncidenceMatrix",
692 arrow(
693 bipartite_graph_ty(),
694 arrow(nat_ty(), arrow(nat_ty(), int_ty())),
695 ),
696 ),
697 ("FacetDefining", facet_defining_ty()),
698 ("tu_integral_polyhedra", tu_integral_polyhedra_ty()),
699 ("bipartite_incidence_tu", bipartite_incidence_tu_ty()),
700 ];
701 for (name, ty) in axioms {
702 env.add(Declaration::Axiom {
703 name: Name::str(*name),
704 univ_params: vec![],
705 ty: ty.clone(),
706 })
707 .ok();
708 }
709 Ok(())
710}
711pub fn hungarian(cost: &[Vec<i64>]) -> (i64, Vec<usize>) {
714 let n = cost.len();
715 if n == 0 {
716 return (0, vec![]);
717 }
718 let inf = i64::MAX / 2;
719 let mut u = vec![0i64; n + 1];
720 let mut v = vec![0i64; n + 1];
721 let mut p = vec![0usize; n + 1];
722 let mut way = vec![0usize; n + 1];
723 for i in 1..=n {
724 p[0] = i;
725 let mut j0 = 0usize;
726 let mut min_val = vec![inf; n + 1];
727 let mut used = vec![false; n + 1];
728 loop {
729 used[j0] = true;
730 let i0 = p[j0];
731 let mut delta = inf;
732 let mut j1 = 0usize;
733 for j in 1..=n {
734 if !used[j] {
735 let cur = cost[i0 - 1][j - 1] - u[i0] - v[j];
736 if cur < min_val[j] {
737 min_val[j] = cur;
738 way[j] = j0;
739 }
740 if min_val[j] < delta {
741 delta = min_val[j];
742 j1 = j;
743 }
744 }
745 }
746 for j in 0..=n {
747 if used[j] {
748 u[p[j]] += delta;
749 v[j] -= delta;
750 } else {
751 min_val[j] -= delta;
752 }
753 }
754 j0 = j1;
755 if p[j0] == 0 {
756 break;
757 }
758 }
759 loop {
760 let j1 = way[j0];
761 p[j0] = p[j1];
762 j0 = j1;
763 if j0 == 0 {
764 break;
765 }
766 }
767 }
768 let mut ans = vec![0usize; n];
769 for j in 1..=n {
770 if p[j] != 0 {
771 ans[p[j] - 1] = j - 1;
772 }
773 }
774 let total_cost = (1..=n).map(|i| cost[i - 1][ans[i - 1]]).sum();
775 (total_cost, ans)
776}
777pub fn uniform_matroid_greedy(weights: &[f64], k: usize) -> Vec<usize> {
780 let mut indexed: Vec<(usize, f64)> = weights.iter().enumerate().map(|(i, &w)| (i, w)).collect();
781 indexed.sort_by(|a, b| b.1.partial_cmp(&a.1).unwrap_or(std::cmp::Ordering::Equal));
782 indexed[..k.min(indexed.len())]
783 .iter()
784 .map(|&(i, _)| i)
785 .collect()
786}
787pub fn greedy_set_cover(n_elem: usize, sets: &[Vec<usize>]) -> Vec<usize> {
789 let mut covered = vec![false; n_elem];
790 let mut remaining = n_elem;
791 let mut selected = vec![];
792 while remaining > 0 {
793 let best = sets
794 .iter()
795 .enumerate()
796 .filter(|(i, _)| !selected.contains(i))
797 .max_by_key(|(_, s)| s.iter().filter(|&&e| !covered[e]).count());
798 match best {
799 Some((idx, set)) => {
800 let newly_covered = set.iter().filter(|&&e| !covered[e]).count();
801 if newly_covered == 0 {
802 break;
803 }
804 selected.push(idx);
805 for &e in set {
806 if !covered[e] {
807 covered[e] = true;
808 remaining -= 1;
809 }
810 }
811 }
812 None => break,
813 }
814 }
815 selected
816}
817pub fn knapsack_01(items: &[(i64, i64)], capacity: i64) -> (i64, Vec<usize>) {
821 let _n = items.len();
822 let mut best = 0i64;
823 let mut best_sol = vec![];
824 fn backtrack(
825 idx: usize,
826 cap: i64,
827 current_val: i64,
828 current_sol: &mut Vec<usize>,
829 best: &mut i64,
830 best_sol: &mut Vec<usize>,
831 items: &[(i64, i64)],
832 ) {
833 if current_val > *best {
834 *best = current_val;
835 *best_sol = current_sol.clone();
836 }
837 if idx == items.len() {
838 return;
839 }
840 let ub: i64 = current_val
841 + items[idx..]
842 .iter()
843 .scan(cap, |c, &(v, w)| {
844 if *c >= w {
845 *c -= w;
846 Some(v)
847 } else {
848 let frac = (*c * v) / w.max(1);
849 *c = 0;
850 Some(frac)
851 }
852 })
853 .sum::<i64>();
854 if ub <= *best {
855 return;
856 }
857 let (v, w) = items[idx];
858 if cap >= w {
859 current_sol.push(idx);
860 backtrack(
861 idx + 1,
862 cap - w,
863 current_val + v,
864 current_sol,
865 best,
866 best_sol,
867 items,
868 );
869 current_sol.pop();
870 }
871 backtrack(
872 idx + 1,
873 cap,
874 current_val,
875 current_sol,
876 best,
877 best_sol,
878 items,
879 );
880 }
881 let mut sol = vec![];
882 backtrack(0, capacity, 0, &mut sol, &mut best, &mut best_sol, items);
883 (best, best_sol)
884}
885pub fn tsp_nearest_neighbor(dist: &[Vec<f64>]) -> (Vec<usize>, f64) {
887 let n = dist.len();
888 if n == 0 {
889 return (vec![], 0.0);
890 }
891 let mut visited = vec![false; n];
892 let mut tour = vec![0usize];
893 visited[0] = true;
894 let mut cost = 0.0;
895 for _ in 1..n {
896 let last = *tour
897 .last()
898 .expect("tour is non-empty: initialized with element 0");
899 let next = (0..n).filter(|&j| !visited[j]).min_by(|&a, &b| {
900 dist[last][a]
901 .partial_cmp(&dist[last][b])
902 .unwrap_or(std::cmp::Ordering::Equal)
903 });
904 if let Some(next) = next {
905 cost += dist[last][next];
906 tour.push(next);
907 visited[next] = true;
908 }
909 }
910 cost += dist[*tour
911 .last()
912 .expect("tour is non-empty: initialized with element 0")][tour[0]];
913 (tour, cost)
914}
915pub fn tsp_held_karp(dist: &[Vec<f64>]) -> f64 {
917 let n = dist.len();
918 if n <= 1 {
919 return 0.0;
920 }
921 let full = 1usize << n;
922 let inf = f64::INFINITY;
923 let mut dp = vec![vec![inf; n]; full];
924 dp[1][0] = 0.0;
925 for mask in 1..full {
926 for u in 0..n {
927 if dp[mask][u] == inf {
928 continue;
929 }
930 if mask & (1 << u) == 0 {
931 continue;
932 }
933 for v in 0..n {
934 if mask & (1 << v) != 0 {
935 continue;
936 }
937 let next_mask = mask | (1 << v);
938 let new_cost = dp[mask][u] + dist[u][v];
939 if new_cost < dp[next_mask][v] {
940 dp[next_mask][v] = new_cost;
941 }
942 }
943 }
944 }
945 let last_mask = full - 1;
946 (1..n)
947 .filter_map(|u| {
948 let c = dp[last_mask][u] + dist[u][0];
949 if c < inf {
950 Some(c)
951 } else {
952 None
953 }
954 })
955 .fold(inf, f64::min)
956}
957#[cfg(test)]
958mod tests {
959 use super::*;
960 #[test]
961 fn test_bipartite_matching_perfect() {
962 let mut g = BipartiteMatchingGraph::new(3, 3);
963 for u in 0..3 {
964 for v in 0..3 {
965 g.add_edge(u, v);
966 }
967 }
968 let (size, _, _) = g.hopcroft_karp();
969 assert_eq!(size, 3);
970 }
971 #[test]
972 fn test_bipartite_matching_partial() {
973 let mut g = BipartiteMatchingGraph::new(3, 2);
974 g.add_edge(0, 0);
975 g.add_edge(1, 0);
976 g.add_edge(1, 1);
977 g.add_edge(2, 1);
978 let (size, _, _) = g.hopcroft_karp();
979 assert_eq!(size, 2);
980 }
981 #[test]
982 fn test_flow_network_max_flow() {
983 let mut net = FlowNetwork::new(4);
984 net.add_edge(0, 1, 3);
985 net.add_edge(0, 2, 2);
986 net.add_edge(1, 3, 3);
987 net.add_edge(2, 3, 2);
988 let flow = net.max_flow(0, 3);
989 assert_eq!(flow, 5);
990 }
991 #[test]
992 fn test_hungarian_algorithm() {
993 let cost = vec![vec![4, 1, 3], vec![2, 0, 5], vec![3, 2, 2]];
994 let (total, assignment) = hungarian(&cost);
995 let mut cols: Vec<usize> = assignment.clone();
996 cols.sort();
997 cols.dedup();
998 assert_eq!(cols.len(), 3, "Assignment must be a permutation");
999 let computed: i64 = assignment
1000 .iter()
1001 .enumerate()
1002 .map(|(i, &j)| cost[i][j])
1003 .sum();
1004 assert_eq!(total, computed);
1005 assert!(total <= 6, "Should find near-optimal solution");
1006 }
1007 #[test]
1008 fn test_uniform_matroid_greedy() {
1009 let weights = [0.5, 3.0, 1.2, 2.8, 0.9];
1010 let selected = uniform_matroid_greedy(&weights, 3);
1011 assert_eq!(selected.len(), 3);
1012 assert!(selected.contains(&1));
1013 assert!(selected.contains(&3));
1014 }
1015 #[test]
1016 fn test_greedy_set_cover() {
1017 let sets = vec![vec![0, 1, 2], vec![1, 3, 4], vec![2, 4, 5]];
1018 let cover = greedy_set_cover(6, &sets);
1019 let mut covered = [false; 6];
1020 for &idx in &cover {
1021 for &e in &sets[idx] {
1022 covered[e] = true;
1023 }
1024 }
1025 assert!(covered.iter().all(|&c| c), "All elements should be covered");
1026 }
1027 #[test]
1028 fn test_knapsack_01() {
1029 let items = vec![(4, 2), (5, 3), (3, 2), (7, 4)];
1030 let (val, sel) = knapsack_01(&items, 7);
1031 assert_eq!(val, 12);
1032 let total_weight: i64 = sel.iter().map(|&i| items[i].1).sum();
1033 assert!(total_weight <= 7, "selected items exceed capacity");
1034 let total_value: i64 = sel.iter().map(|&i| items[i].0).sum();
1035 assert_eq!(total_value, 12, "selected items should have value 12");
1036 }
1037 #[test]
1038 fn test_tsp_held_karp_small() {
1039 let d = vec![
1040 vec![0.0, 10.0, 15.0, 20.0],
1041 vec![10.0, 0.0, 35.0, 25.0],
1042 vec![15.0, 35.0, 0.0, 30.0],
1043 vec![20.0, 25.0, 30.0, 0.0],
1044 ];
1045 let opt = tsp_held_karp(&d);
1046 assert!(
1047 (opt - 80.0).abs() < 1e-9,
1048 "Held-Karp should find optimal tour of length 80, got {}",
1049 opt
1050 );
1051 }
1052 #[test]
1053 fn test_build_combinatorial_optimization_env() {
1054 let mut env = Environment::new();
1055 let result = build_combinatorial_optimization_env(&mut env);
1056 assert!(
1057 result.is_ok(),
1058 "build_combinatorial_optimization_env failed: {:?}",
1059 result.err()
1060 );
1061 }
1062}
1063#[cfg(test)]
1064mod spec_wrapper_tests {
1065 use super::*;
1066 #[test]
1067 fn test_flow_network_spec() {
1068 let mut net = FlowNetworkSpec::new(4);
1069 net.add_edge(0, 1, 3);
1070 net.add_edge(0, 2, 2);
1071 net.add_edge(1, 3, 3);
1072 net.add_edge(2, 3, 2);
1073 assert_eq!(net.max_flow_ford_fulkerson(0, 3), 5);
1074 assert_eq!(net.min_cut(0, 3), 5);
1075 assert!(net.augmenting_path(0, 3));
1076 }
1077 #[test]
1078 fn test_matching_problem_bipartite() {
1079 let mut mp = MatchingProblem::new(true, 3, 3);
1080 mp.add_edge(0, 0, 1.0);
1081 mp.add_edge(1, 1, 2.0);
1082 mp.add_edge(2, 2, 3.0);
1083 assert_eq!(mp.max_matching(), 3);
1084 }
1085 #[test]
1086 fn test_spanning_tree_kruskal() {
1087 let mut st = SpanningTree::new(4);
1088 st.add_edge(0, 1, 1.0);
1089 st.add_edge(1, 2, 2.0);
1090 st.add_edge(2, 3, 3.0);
1091 st.add_edge(0, 3, 10.0);
1092 let mst = st.kruskal();
1093 assert_eq!(mst.len(), 3);
1094 let w: f64 = mst.iter().map(|&(_, _, w)| w).sum();
1095 assert!((w - 6.0).abs() < 1e-9);
1096 }
1097 #[test]
1098 fn test_shortest_path_dijkstra() {
1099 let mut sp = ShortestPath::new(3);
1100 sp.add_edge(0, 1, 1.0);
1101 sp.add_edge(1, 2, 2.0);
1102 sp.add_edge(0, 2, 5.0);
1103 let dist = sp.dijkstra(0);
1104 assert!((dist[2] - 3.0).abs() < 1e-9);
1105 }
1106 #[test]
1107 fn test_bellman_ford() {
1108 let mut sp = ShortestPath::new(3);
1109 sp.add_edge(0, 1, 1.0);
1110 sp.add_edge(1, 2, 2.0);
1111 sp.add_edge(0, 2, 5.0);
1112 let dist = sp.bellman_ford(0);
1113 assert!((dist[2] - 3.0).abs() < 1e-9);
1114 }
1115 #[test]
1116 fn test_floyd_warshall() {
1117 let mut sp = ShortestPath::new(3);
1118 sp.add_edge(0, 1, 1.0);
1119 sp.add_edge(1, 2, 2.0);
1120 sp.add_edge(0, 2, 5.0);
1121 let d = sp.floyd_warshall();
1122 assert!((d[0][2] - 3.0).abs() < 1e-9);
1123 }
1124 #[test]
1125 fn test_knapsack_solver() {
1126 let solver = KnapsackSolver::new(vec![(4, 2), (5, 3), (3, 2), (7, 4)], 7);
1127 let (val, _) = solver.dynamic_programming();
1128 assert_eq!(val, 12);
1129 }
1130 #[test]
1131 fn test_graph_coloring() {
1132 let gc = GraphColoring::new(3, vec![(0, 1), (1, 2), (0, 2)]);
1133 let (k, _) = gc.greedy_color();
1134 assert!(k >= 3);
1135 assert!(gc.chromatic_number_bound() >= 3);
1136 let (dk, _) = gc.dsatur();
1137 assert!(dk >= 3);
1138 }
1139 #[test]
1140 fn test_traveling_salesman() {
1141 let d = vec![
1142 vec![0.0, 10.0, 15.0, 20.0],
1143 vec![10.0, 0.0, 35.0, 25.0],
1144 vec![15.0, 35.0, 0.0, 30.0],
1145 vec![20.0, 25.0, 30.0, 0.0],
1146 ];
1147 let tsp = TravelingSalesman::new(d);
1148 let hk = tsp.held_karp();
1149 assert!((hk - 80.0).abs() < 1e-9);
1150 let nn = tsp.nearest_neighbor();
1151 assert!(nn > 0.0);
1152 }
1153 #[test]
1154 fn test_steiner_tree() {
1155 let st = SteinerTree::new(5, vec![0, 2, 4]);
1156 assert!(st.approx_2() >= 0.0);
1157 assert!(st.dreyfus_wagner() >= 0.0);
1158 }
1159 #[test]
1160 fn test_matroid_intersection() {
1161 let mi = MatroidIntersection::new("graphic", "partition");
1162 assert!(mi.exchange_property());
1163 }
1164}
1165#[cfg(test)]
1166mod extended_comb_opt_tests {
1167 use super::*;
1168 #[test]
1169 fn test_branch_bound() {
1170 let mut bb = BranchBoundData::integer_program("most-fractional", "LP relaxation");
1171 bb.explore(100);
1172 assert_eq!(bb.nodes_explored, 100);
1173 assert!(bb.description().contains("B&B"));
1174 }
1175 #[test]
1176 fn test_cutting_plane() {
1177 let mut cp = CuttingPlane::gomory();
1178 cp.add_cut();
1179 cp.add_cut();
1180 assert_eq!(cp.num_cuts_added, 2);
1181 assert!(cp.description().contains("Gomory"));
1182 }
1183 #[test]
1184 fn test_set_cover() {
1185 let sc = SetCoverData::greedy(5, 10);
1186 assert!(sc.approximation_ratio > 1.0);
1187 assert!(sc.approx_description().contains("Greedy"));
1188 }
1189 #[test]
1190 fn test_facility_location() {
1191 let fl = FacilityLocation::new(3, 10, vec![1.0, 2.0, 3.0], 5.0);
1192 assert_eq!(fl.total_opening_cost(), 6.0);
1193 assert!((FacilityLocation::jms_approximation_ratio() - 1.861).abs() < 0.001);
1194 }
1195 #[test]
1196 fn test_vehicle_routing() {
1197 let vr = VehicleRouting::capacitated(3, 20, 100.0);
1198 assert_eq!(vr.num_vehicles, 3);
1199 assert!(vr.christofides_description().contains("Christofides"));
1200 }
1201}