Skip to main content

oxilean_std/combinatorial_optimization/
functions.rs

1//! Auto-generated module
2//!
3//! ๐Ÿค– Generated with [SplitRS](https://github.com/cool-japan/splitrs)
4
5use 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}
62/// `Graph : Type` โ€” a simple undirected graph (vertices: Nat, edges: pairs).
63pub fn graph_ty() -> Expr {
64    type0()
65}
66/// `BipartiteGraph : Type` โ€” bipartite graph with two vertex sets.
67pub fn bipartite_graph_ty() -> Expr {
68    type0()
69}
70/// `Matching : Graph โ†’ Type` โ€” a set of non-adjacent edges.
71pub fn matching_ty() -> Expr {
72    arrow(graph_ty(), type0())
73}
74/// `PerfectMatching : Graph โ†’ Prop` โ€” every vertex is matched.
75pub fn perfect_matching_ty() -> Expr {
76    arrow(graph_ty(), prop())
77}
78/// `MaxMatching : Graph โ†’ Matching โ†’ Prop` โ€” a maximum cardinality matching.
79pub fn max_matching_ty() -> Expr {
80    arrow(graph_ty(), arrow(matching_ty(), prop()))
81}
82/// `Alternating_path : Graph โ†’ Matching โ†’ List Nat โ†’ Prop` โ€” alternating path wrt matching.
83pub fn alternating_path_ty() -> Expr {
84    arrow(
85        graph_ty(),
86        arrow(matching_ty(), arrow(list_ty(nat_ty()), prop())),
87    )
88}
89/// `AugmentingPath : Graph โ†’ Matching โ†’ List Nat โ†’ Prop` โ€” augmenting path.
90pub fn augmenting_path_ty() -> Expr {
91    arrow(
92        graph_ty(),
93        arrow(matching_ty(), arrow(list_ty(nat_ty()), prop())),
94    )
95}
96/// `Blossom : Graph โ†’ Matching โ†’ List Nat โ†’ Prop` โ€” odd cycle contracted in Edmonds' algorithm.
97pub fn blossom_ty() -> Expr {
98    arrow(
99        graph_ty(),
100        arrow(matching_ty(), arrow(list_ty(nat_ty()), prop())),
101    )
102}
103/// `HallCondition : BipartiteGraph โ†’ Prop`
104/// For every subset S of left vertices, |N(S)| โ‰ฅ |S|.
105pub fn hall_condition_ty() -> Expr {
106    arrow(bipartite_graph_ty(), prop())
107}
108/// `HallTheorem : โˆ€ G : BipartiteGraph, HallCondition G โ†” PerfectMatching G`
109pub 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}
121/// `KonigTheorem : MaxMatchingSize = MinVertexCoverSize` for bipartite graphs.
122pub 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}
134/// `TutteCondition : Graph โ†’ Prop`
135/// For every S โІ V, the number of odd components of G-S is โ‰ค |S|.
136pub fn tutte_condition_ty() -> Expr {
137    arrow(graph_ty(), prop())
138}
139/// `TutteTheorem : โˆ€ G, PerfectMatching G โ†” TutteCondition G`
140pub 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}
152/// `BergeTheorem : max matching augmenting path characterization`
153pub 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}
173/// `FlowNetwork : Type` โ€” directed graph with capacities.
174pub fn flow_network_ty() -> Expr {
175    type0()
176}
177/// `Flow : FlowNetwork โ†’ Type` โ€” feasible flow satisfying capacity and conservation.
178pub fn flow_ty() -> Expr {
179    arrow(flow_network_ty(), type0())
180}
181/// `FlowValue : FlowNetwork โ†’ Flow โ†’ Real` โ€” the value of a flow (net flow out of source).
182pub fn flow_value_ty() -> Expr {
183    arrow(flow_network_ty(), arrow(flow_ty(), real_ty()))
184}
185/// `MaxFlow : FlowNetwork โ†’ Real` โ€” maximum flow value.
186pub fn max_flow_ty() -> Expr {
187    arrow(flow_network_ty(), real_ty())
188}
189/// `Cut : FlowNetwork โ†’ (Nat โ†’ Bool) โ†’ Prop` โ€” s-t cut (partition of vertices).
190pub fn cut_ty() -> Expr {
191    arrow(flow_network_ty(), arrow(arrow(nat_ty(), bool_ty()), prop()))
192}
193/// `CutCapacity : FlowNetwork โ†’ (Nat โ†’ Bool) โ†’ Real` โ€” total capacity of a cut.
194pub fn cut_capacity_ty() -> Expr {
195    arrow(
196        flow_network_ty(),
197        arrow(arrow(nat_ty(), bool_ty()), real_ty()),
198    )
199}
200/// `MinCut : FlowNetwork โ†’ Real` โ€” minimum cut capacity.
201pub fn min_cut_ty() -> Expr {
202    arrow(flow_network_ty(), real_ty())
203}
204/// `MaxFlowMinCut : โˆ€ N, MaxFlow N = MinCut N` โ€” the max-flow min-cut theorem.
205pub 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}
217/// `FordFulkersonTermination : Ford-Fulkerson terminates on integer capacities.`
218pub 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}
229/// `DinicsComplexity : Dinic's algorithm runs in O(Vยฒ E).`
230pub fn dinics_complexity_ty() -> Expr {
231    pi(
232        BinderInfo::Default,
233        "N",
234        flow_network_ty(),
235        app(cst("DinicsTimeBound"), bvar(0)),
236    )
237}
238/// `CostMatrix : Nat โ†’ Nat โ†’ Real` โ€” cost matrix for assignment.
239pub fn cost_matrix_ty() -> Expr {
240    arrow(nat_ty(), arrow(nat_ty(), real_ty()))
241}
242/// `Assignment : Nat โ†’ (Nat โ†’ Nat) โ†’ Prop` โ€” a permutation assignment for nร—n matrix.
243pub fn assignment_ty() -> Expr {
244    arrow(nat_ty(), arrow(arrow(nat_ty(), nat_ty()), prop()))
245}
246/// `OptimalAssignment : CostMatrix โ†’ (Nat โ†’ Nat) โ†’ Prop` โ€” assignment minimizing total cost.
247pub fn optimal_assignment_ty() -> Expr {
248    arrow(cost_matrix_ty(), arrow(arrow(nat_ty(), nat_ty()), prop()))
249}
250/// `HungarianAlgorithm : CostMatrix โ†’ (Nat โ†’ Nat)` โ€” solves assignment in O(nยณ).
251pub fn hungarian_algorithm_ty() -> Expr {
252    arrow(cost_matrix_ty(), arrow(nat_ty(), nat_ty()))
253}
254/// `HungarianCorrectness : HungarianAlgorithm gives OptimalAssignment.`
255pub 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}
267/// `TSPInstance : Type` โ€” complete graph with edge weights.
268pub fn tsp_instance_ty() -> Expr {
269    type0()
270}
271/// `TSPTour : TSPInstance โ†’ List Nat โ†’ Prop` โ€” a Hamiltonian cycle.
272pub fn tsp_tour_ty() -> Expr {
273    arrow(tsp_instance_ty(), arrow(list_ty(nat_ty()), prop()))
274}
275/// `TSPOptimal : TSPInstance โ†’ Real` โ€” optimal TSP tour length.
276pub fn tsp_optimal_ty() -> Expr {
277    arrow(tsp_instance_ty(), real_ty())
278}
279/// `HeldKarpBound : TSPInstance โ†’ Real` โ€” Held-Karp lower bound.
280pub fn held_karp_bound_ty() -> Expr {
281    arrow(tsp_instance_ty(), real_ty())
282}
283/// `HeldKarpLowerBound : HeldKarpBound โ‰ค TSPOptimal.`
284pub 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}
296/// `ChristofidesApproximation : 3/2-approximation for metric TSP.`
297pub 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}
316/// `VertexCoverApproximation : 2-approx for vertex cover.`
317pub 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}
333/// `SetCoverApproximation : H_n-approximation (log n) for set cover.`
334pub 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}
350/// `Matroid : Type` โ€” a matroid (ground set + independent sets satisfying axioms).
351pub fn matroid_ty() -> Expr {
352    type0()
353}
354/// `IndependentSet : Matroid โ†’ List Nat โ†’ Prop`
355pub fn independent_set_ty() -> Expr {
356    arrow(matroid_ty(), arrow(list_ty(nat_ty()), prop()))
357}
358/// `MatroidBase : Matroid โ†’ List Nat โ†’ Prop` โ€” maximal independent set.
359pub fn matroid_base_ty() -> Expr {
360    arrow(matroid_ty(), arrow(list_ty(nat_ty()), prop()))
361}
362/// `MatroidRank : Matroid โ†’ Nat` โ€” rank of a matroid.
363pub fn matroid_rank_ty() -> Expr {
364    arrow(matroid_ty(), nat_ty())
365}
366/// `GreedyOptimality : Greedy algorithm is optimal for matroid weight maximization.`
367pub 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}
380/// `MatroidIntersection : common independent set in two matroids.`
381pub fn matroid_intersection_ty() -> Expr {
382    arrow(matroid_ty(), arrow(matroid_ty(), list_ty(nat_ty())))
383}
384/// `MatroidIntersectionOptimality : the max weight common independent set.`
385pub 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}
403/// `GraphicMatroid : Graph โ†’ Matroid` โ€” cycle matroid of a graph.
404pub fn graphic_matroid_ty() -> Expr {
405    arrow(graph_ty(), matroid_ty())
406}
407/// `UniformMatroid : Nat โ†’ Nat โ†’ Matroid` โ€” U(k,n) uniform matroid.
408pub fn uniform_matroid_ty() -> Expr {
409    arrow(nat_ty(), arrow(nat_ty(), matroid_ty()))
410}
411/// `SubmodularFunction : (List Nat โ†’ Real) โ†’ Prop`
412/// f is submodular if f(A) + f(B) โ‰ฅ f(AโˆชB) + f(AโˆฉB).
413pub fn submodular_function_ty() -> Expr {
414    arrow(arrow(list_ty(nat_ty()), real_ty()), prop())
415}
416/// `SubmodularMaximization : (List Nat โ†’ Real) โ†’ List Nat` โ€” greedy maximization.
417pub fn submodular_maximization_ty() -> Expr {
418    arrow(arrow(list_ty(nat_ty()), real_ty()), list_ty(nat_ty()))
419}
420/// `SubmodularGreedy_1_2_approx : greedy gives 1/2-approximation for monotone submodular max.`
421pub 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}
440/// `SupermodularFunction : (List Nat โ†’ Real) โ†’ Prop` โ€” supermodular (negation of submodular).
441pub fn supermodular_function_ty() -> Expr {
442    arrow(arrow(list_ty(nat_ty()), real_ty()), prop())
443}
444/// `PolymatroidRankFunction : axioms for polymatroid rank.`
445pub fn polymatroid_rank_ty() -> Expr {
446    arrow(arrow(list_ty(nat_ty()), real_ty()), prop())
447}
448/// `ILPInstance : Type` โ€” integer linear program: min cยทx s.t. Ax โ‰ค b, x โˆˆ โ„ค^n.
449pub fn ilp_instance_ty() -> Expr {
450    type0()
451}
452/// `ILPSolution : ILPInstance โ†’ List Int โ†’ Prop`
453pub fn ilp_solution_ty() -> Expr {
454    arrow(ilp_instance_ty(), arrow(list_ty(int_ty()), prop()))
455}
456/// `ILPOptimal : ILPInstance โ†’ List Int โ†’ Prop`
457pub fn ilp_optimal_ty() -> Expr {
458    arrow(ilp_instance_ty(), arrow(list_ty(int_ty()), prop()))
459}
460/// `GomoryCut : ILPInstance โ†’ ILPInstance` โ€” adds a Gomory cutting plane.
461pub fn gomory_cut_ty() -> Expr {
462    arrow(ilp_instance_ty(), ilp_instance_ty())
463}
464/// `BranchAndBound : ILPInstance โ†’ List Int` โ€” branch-and-bound solver.
465pub fn branch_and_bound_ty() -> Expr {
466    arrow(ilp_instance_ty(), list_ty(int_ty()))
467}
468/// `LPRelaxation : ILPInstance โ†’ Real` โ€” optimal LP relaxation value.
469pub fn lp_relaxation_ty() -> Expr {
470    arrow(ilp_instance_ty(), real_ty())
471}
472/// `LPRelaxationLowerBound : LPRelaxation โ‰ค ILP optimal.`
473pub 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}
485/// `Polytope : Type` โ€” a convex polytope (intersection of halfspaces).
486pub fn polytope_ty() -> Expr {
487    type0()
488}
489/// `Vertex : Polytope โ†’ List Real โ†’ Prop` โ€” extreme point of polytope.
490pub fn polytope_vertex_ty() -> Expr {
491    arrow(polytope_ty(), arrow(list_ty(real_ty()), prop()))
492}
493/// `TotallyUnimodular : (Nat โ†’ Nat โ†’ Int) โ†’ Prop` โ€” every square submatrix has det โˆˆ {-1,0,1}.
494pub fn totally_unimodular_ty() -> Expr {
495    arrow(arrow(nat_ty(), arrow(nat_ty(), int_ty())), prop())
496}
497/// `TUIntegralPolyhedra : TU matrix โ†’ LP has integer optimal vertex.`
498pub 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}
509/// `BipartiteIncidenceTU : incidence matrix of bipartite graph is TU.`
510pub 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}
521/// `FacetDefiningInequality : Polytope โ†’ (List Real โ†’ Prop) โ†’ Prop`
522pub fn facet_defining_ty() -> Expr {
523    arrow(
524        polytope_ty(),
525        arrow(arrow(list_ty(real_ty()), prop()), prop()),
526    )
527}
528/// `WeakDuality : for LP, dual objective โ‰ฅ primal objective.`
529pub 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}
541/// `StrongDuality : for LP, dual objective = primal when both feasible.`
542pub 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}
557/// Register all combinatorial optimization axioms into the kernel environment.
558pub 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}
711/// Solve the minimum-weight assignment problem using the Hungarian algorithm.
712/// Returns (total cost, assignment vector where assignment[i] = j means row i โ†’ col j).
713pub 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}
777/// Greedy maximum weight independent set in uniform matroid U(k, n).
778/// Returns indices of k elements with highest weights.
779pub 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}
787/// Greedy set cover: returns selected set indices covering all elements 0..n_elem.
788pub 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}
817/// Branch-and-bound for 0-1 knapsack (maximize).
818/// items: (value, weight), capacity: max weight.
819/// Returns (max_value, selected_item_indices).
820pub 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}
885/// Nearest-neighbor heuristic for TSP. Returns tour (list of vertex indices).
886pub 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}
915/// Held-Karp DP lower bound (exact TSP for small n โ‰ค 20).
916pub 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}