Source file simplex.ml
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
open NumCompat
open Q.Notations
open Polynomial
open Mutils
type ('a, 'b) sum = Inl of 'a | Inr of 'b
let debug = false
(** Exploiting profiling information *)
let profile_info = ref []
let nb_pivot = ref 0
type profile_info =
{ number_of_successes : int
; number_of_failures : int
; success_pivots : int
; failure_pivots : int
; average_pivots : int
; maximum_pivots : int }
let init_profile =
{ number_of_successes = 0
; number_of_failures = 0
; success_pivots = 0
; failure_pivots = 0
; average_pivots = 0
; maximum_pivots = 0 }
let get_profile_info () =
let update_profile
{ number_of_successes
; number_of_failures
; success_pivots
; failure_pivots
; average_pivots
; maximum_pivots } (b, i) =
{ number_of_successes = (number_of_successes + if b then 1 else 0)
; number_of_failures = (number_of_failures + if b then 0 else 1)
; success_pivots = (success_pivots + if b then i else 0)
; failure_pivots = (failure_pivots + if b then 0 else i)
; average_pivots = average_pivots + 1
; maximum_pivots = max maximum_pivots i }
in
let p = List.fold_left update_profile init_profile !profile_info in
profile_info := [];
{ p with
average_pivots =
( try (p.success_pivots + p.failure_pivots) / p.average_pivots
with Division_by_zero -> 0 ) }
type iset = unit IMap.t
(** Mapping basic variables to their equation.
All variables >= than a threshold rst are restricted.*)
type tableau = Vect.t IMap.t
module Restricted = struct
type t =
{ base : int (** All variables above [base] are restricted *)
; exc : int option (** Except [exc] which is currently optimised *) }
let pp o {base; exc} =
Printf.fprintf o ">= %a " LinPoly.pp_var base;
match exc with
| None -> Printf.fprintf o "-"
| Some x -> Printf.fprintf o "-%a" LinPoly.pp_var base
let is_exception (x : var) (r : t) =
match r.exc with None -> false | Some x' -> x = x'
let restrict x rst =
if is_exception x rst then {base = rst.base; exc = None}
else failwith (Printf.sprintf "Cannot restrict %i" x)
let is_restricted x r0 = x >= r0.base && not (is_exception x r0)
let make x = {base = x; exc = None}
let set_exc x rst = {base = rst.base; exc = Some x}
let fold rst f m acc =
IMap.fold
(fun k v acc -> if is_exception k rst then acc else f k v acc)
(IMap.from rst.base m) acc
end
let pp_row o v = LinPoly.pp o v
let output_tableau o t =
IMap.iter
(fun k v -> Printf.fprintf o "%a = %a\n" LinPoly.pp_var k pp_row v)
t
let output_env o t =
IMap.iter (fun k v -> Printf.fprintf o "%i : %a\n" k WithProof.output v) t
let output_vars o m =
IMap.iter (fun k _ -> Printf.fprintf o "%a " LinPoly.pp_var k) m
(** A tableau is feasible iff for every basic restricted variable xi,
we have ci>=0.
When all the non-basic variables are set to 0, the value of a basic
variable xi is necessarily ci. If xi is restricted, it is feasible
if ci>=0.
*)
let unfeasible (rst : Restricted.t) tbl =
Restricted.fold rst
(fun k v m -> if Vect.get_cst v >=/ Q.zero then m else IMap.add k () m)
tbl IMap.empty
let is_feasible rst tb = IMap.is_empty (unfeasible rst tb)
(** Let a1.x1+...+an.xn be a vector of non-basic variables.
It is maximised if all the xi are restricted
and the ai are negative.
If xi>= 0 (restricted) and ai is negative,
the maximum for ai.xi is obtained for xi = 0
Otherwise, it is possible to make ai.xi arbitrarily big:
- if xi is not restricted, take +/- oo depending on the sign of ai
- if ai is positive, take +oo
*)
let is_maximised_vect rst v =
Vect.for_all
(fun xi ai ->
if ai >/ Q.zero then false else Restricted.is_restricted xi rst)
v
(** [is_maximised rst v]
@return None if the variable is not maximised
@return Some v where v is the maximal value
*)
let is_maximised rst v =
try
let vl, v = Vect.decomp_cst v in
if is_maximised_vect rst v then Some vl else None
with Not_found -> None
(** A variable xi is unbounded if for every
equation xj= ...ai.xi ...
if ai < 0 then xj is not restricted.
As a result, even if we
increase the value of xi, it is always
possible to adjust the value of xj without
violating a restriction.
*)
type result =
| Max of Q.t (** Maximum is reached *)
| Ubnd of var (** Problem is unbounded *)
| Feas (** Problem is feasible *)
type pivot = Done of result | Pivot of int * int * Q.t
type simplex = Opt of tableau * result
(** For a row, x = ao.xo+...+ai.xi
a valid pivot variable is such that it can improve the value of xi.
it is the case, if xi is unrestricted (increase if ai> 0, decrease if ai < 0)
xi is restricted but ai > 0
This is the entering variable.
*)
let rec find_pivot_column (rst : Restricted.t) (r : Vect.t) =
match Vect.choose r with
| None -> failwith "find_pivot_column"
| Some (xi, ai, r') ->
if ai </ Q.zero then
if Restricted.is_restricted xi rst then find_pivot_column rst r'
else (xi, -1)
else
(xi, 1)
(** Finding the variable leaving the basis is more subtle because we need to:
- increase the objective function
- make sure that the entering variable has a feasible value
- but also that after pivoting all the other basic variables are still feasible.
This explains why we choose the pivot with the smallest score
*)
let min_score s (i1, sc1) =
match s with
| None -> Some (i1, sc1)
| Some (i0, sc0) ->
if sc0 </ sc1 then s
else if sc1 </ sc0 then Some (i1, sc1)
else if i0 < i1 then s
else Some (i1, sc1)
let find_pivot_row rst tbl j sgn =
Restricted.fold rst
(fun i' v res ->
let aij = Vect.get j v in
if Q.of_int sgn */ aij </ Q.zero then
let score' = Q.abs (Vect.get_cst v // aij) in
min_score res (i', score')
else res)
tbl None
let safe_find err x t =
try IMap.find x t
with Not_found ->
if debug then
Printf.fprintf stdout "safe_find %s x%i %a\n" err x output_tableau t;
failwith err
(** [find_pivot vr t] aims at improving the objective function of the basic variable vr *)
let find_pivot vr (rst : Restricted.t) tbl =
let v = safe_find "find_pivot" vr tbl in
match is_maximised rst v with
| Some mx -> Done (Max mx)
| None -> (
let _, v = Vect.decomp_cst v in
let j', sgn = find_pivot_column rst v in
match find_pivot_row rst (IMap.remove vr tbl) j' sgn with
| None -> Done (Ubnd j')
| Some (i', sc) -> Pivot (i', j', sc) )
(** [solve_column c r e]
@param c is a non-basic variable
@param r is a basic variable
@param e is a vector such that r = e
and e is of the form ai.c+e'
@return the vector (-r + e').-1/ai i.e
c = (r - e')/ai
*)
let solve_column (c : var) (r : var) (e : Vect.t) : Vect.t =
let a = Vect.get c e in
if a =/ Q.zero then failwith "Cannot solve column"
else
let a' = Q.minus_one // a in
Vect.mul a' (Vect.set r Q.minus_one (Vect.set c Q.zero e))
(** [pivot_row r c e]
@param c is such that c = e
@param r is a vector r = g.c + r'
@return g.e+r' *)
let pivot_row (row : Vect.t) (c : var) (e : Vect.t) : Vect.t =
let g = Vect.get c row in
if g =/ Q.zero then row else Vect.mul_add g e Q.one (Vect.set c Q.zero row)
let pivot_with (m : tableau) (v : var) (p : Vect.t) =
IMap.map (fun (r : Vect.t) -> pivot_row r v p) m
let pivot (m : tableau) (r : var) (c : var) =
incr nb_pivot;
let row = safe_find "pivot" r m in
let piv = solve_column c r row in
IMap.add c piv (pivot_with (IMap.remove r m) c piv)
let adapt_unbounded vr x rst tbl =
if Vect.get_cst (IMap.find vr tbl) >=/ Q.zero then tbl else pivot tbl vr x
module BaseSet = Set.Make (struct
type t = iset
let compare = IMap.compare (fun x y -> 0)
end)
let get_base tbl = IMap.mapi (fun k _ -> ()) tbl
let simplex opt vr rst tbl =
let b = ref BaseSet.empty in
let rec simplex opt vr rst tbl =
( if debug then
let base = get_base tbl in
if BaseSet.mem base !b then Printf.fprintf stdout "Cycling detected\n"
else b := BaseSet.add base !b );
if debug && not (is_feasible rst tbl) then begin
let m = unfeasible rst tbl in
Printf.fprintf stdout "Simplex error\n";
Printf.fprintf stdout "The current tableau is not feasible\n";
Printf.fprintf stdout "Restricted >= %a\n" Restricted.pp rst;
output_tableau stdout tbl;
Printf.fprintf stdout "Error for variables %a\n" output_vars m
end;
if (not opt) && Vect.get_cst (IMap.find vr tbl) >=/ Q.zero then
Opt (tbl, Feas)
else
match find_pivot vr rst tbl with
| Done r -> (
match r with
| Max _ -> Opt (tbl, r)
| Ubnd x ->
let t' = adapt_unbounded vr x rst tbl in
Opt (t', r)
| Feas -> raise (Invalid_argument "find_pivot") )
| Pivot (i, j, s) ->
if debug then begin
Printf.fprintf stdout "Find pivot for x%i(%s)\n" vr (Q.to_string s);
Printf.fprintf stdout "Leaving variable x%i\n" i;
Printf.fprintf stdout "Entering variable x%i\n" j
end;
let m' = pivot tbl i j in
simplex opt vr rst m'
in
simplex opt vr rst tbl
type certificate = Unsat of Vect.t | Sat of tableau * var option
(** [normalise_row t v]
@return a row obtained by pivoting the basic variables of the vector v
*)
let normalise_row (t : tableau) (v : Vect.t) =
Vect.fold
(fun acc vr ai ->
try
let e = IMap.find vr t in
Vect.add (Vect.mul ai e) acc
with Not_found -> Vect.add (Vect.set vr ai Vect.null) acc)
Vect.null v
let normalise_row (t : tableau) (v : Vect.t) =
let v' = normalise_row t v in
if debug then Printf.fprintf stdout "Normalised Optimising %a\n" LinPoly.pp v';
v'
let add_row (nw : var) (t : tableau) (v : Vect.t) : tableau =
IMap.add nw (normalise_row t v) t
(** [push_real] performs reasoning over the rationals *)
let push_real (opt : bool) (nw : var) (v : Vect.t) (rst : Restricted.t)
(t : tableau) : certificate =
if debug then begin
Printf.fprintf stdout "Current Tableau\n%a\n" output_tableau t;
Printf.fprintf stdout "Optimising %a=%a\n" LinPoly.pp_var nw LinPoly.pp v
end;
match simplex opt nw rst (add_row nw t v) with
| Opt (t', r) -> (
match r with
| Ubnd x ->
if debug then
Printf.printf "The objective is unbounded (variable %a)\n"
LinPoly.pp_var x;
Sat (t', Some x)
| Feas -> Sat (t', None)
| Max n ->
if debug then begin
Printf.printf "The objective is maximised %s\n" (Q.to_string n);
Printf.printf "%a = %a\n" LinPoly.pp_var nw pp_row (IMap.find nw t')
end;
if n >=/ Q.zero then Sat (t', None)
else
let v' = safe_find "push_real" nw t' in
Unsat (Vect.set nw Q.one (Vect.set 0 Q.zero (Vect.mul Q.minus_one v')))
)
(** One complication is that equalities needs some pre-processing.
*)
open Mutils
open Polynomial
let find_solution rst tbl =
IMap.fold
(fun vr v res ->
if Restricted.is_restricted vr rst then res
else Vect.set vr (Vect.get_cst v) res)
tbl Vect.null
let find_full_solution rst tbl =
IMap.fold (fun vr v res -> Vect.set vr (Vect.get_cst v) res) tbl Vect.null
let choose_conflict (sol : Vect.t) (l : (var * Vect.t) list) =
let esol = Vect.set 0 Q.one sol in
let rec most_violating l e (x, v) rst =
match l with
| [] -> Some ((x, v), rst)
| (x', v') :: l ->
let e' = Vect.dotproduct esol v' in
if e' <=/ e then most_violating l e' (x', v') ((x, v) :: rst)
else most_violating l e (x, v) ((x', v') :: rst)
in
match l with
| [] -> None
| (x, v) :: l ->
let e = Vect.dotproduct esol v in
most_violating l e (x, v) []
let rec solve opt l (rst : Restricted.t) (t : tableau) =
let sol = find_solution rst t in
match choose_conflict sol l with
| None -> Inl (rst, t, None)
| Some ((vr, v), l) -> (
match push_real opt vr v (Restricted.set_exc vr rst) t with
| Sat (t', x) -> (
match l with [] -> Inl (rst, t', x) | _ -> solve opt l rst t' )
| Unsat c -> Inr c )
let fresh_var l =
1
+
try
ISet.max_elt
(List.fold_left
(fun acc c -> ISet.union acc (Vect.variables c.coeffs))
ISet.empty l)
with Not_found -> 0
module PrfEnv = struct
type t = WithProof.t IMap.t
let empty = IMap.empty
let register prf env =
let fr = LinPoly.MonT.get_fresh () in
(fr, IMap.add fr prf env)
let set_prf i prf env = IMap.add i prf env
let find idx env = IMap.find idx env
let rec of_list acc env l =
match l with
| [] -> (acc, env)
| (((lp, op), prf) as wp) :: l -> (
match op with
| Gt -> raise Strict
| Ge ->
let f, env' = register wp env in
of_list ((f, lp) :: acc) env' l
| Eq ->
let f1, env = register wp env in
let wp' = WithProof.neg wp in
let f2, env = register wp' env in
of_list ((f1, lp) :: (f2, fst (fst wp')) :: acc) env l )
let map f env = IMap.map f env
end
let make_env (l : Polynomial.cstr list) =
PrfEnv.of_list [] PrfEnv.empty
(List.rev_map WithProof.of_cstr
(List.mapi (fun i x -> (x, ProofFormat.Hyp i)) l))
let find_point (l : Polynomial.cstr list) =
let vr = fresh_var l in
LinPoly.MonT.safe_reserve vr;
let l', _ = make_env l in
match solve false l' (Restricted.make vr) IMap.empty with
| Inl (rst, t, _) -> Some (find_solution rst t)
| _ -> None
let optimise obj l =
let vr = fresh_var l in
LinPoly.MonT.safe_reserve vr;
let l', _ = make_env l in
let bound pos res =
match res with
| Opt (_, Max n) -> Some (if pos then n else Q.neg n)
| Opt (_, Ubnd _) -> None
| Opt (_, Feas) -> None
in
match solve false l' (Restricted.make vr) IMap.empty with
| Inl (rst, t, _) ->
Some
( bound false (simplex true vr rst (add_row vr t (Vect.uminus obj)))
, bound true (simplex true vr rst (add_row vr t obj)) )
| _ -> None
(** [make_certificate env l] makes very strong assumptions
about the form of the environment.
Each proof is assumed to be either:
- an hypothesis Hyp i
- or, the negation of an hypothesis (MulC(-1,Hyp i))
*)
let make_certificate env l =
Vect.normalise
(Vect.fold
(fun acc x n ->
let _, prf = PrfEnv.find x env in
ProofFormat.(
match prf with
| Hyp i -> Vect.set i n acc
| MulC (_, Hyp i) -> Vect.set i (Q.neg n) acc
| _ -> failwith "make_certificate: invalid proof"))
Vect.null l)
let find_unsat_certificate (l : Polynomial.cstr list) =
let l', env = make_env l in
let vr = fresh_var l in
match solve false l' (Restricted.make vr) IMap.empty with
| Inr c -> Some (make_certificate env c)
| Inl _ -> None
open Polynomial
open ProofFormat
let make_farkas_certificate (env : PrfEnv.t) v =
Vect.fold
(fun acc x n -> add_proof acc (mul_cst_proof n (snd (PrfEnv.find x env))))
Zero v
let make_farkas_proof (env : PrfEnv.t) v =
Vect.fold
(fun wp x n ->
WithProof.addition wp (WithProof.mult (Vect.cst n) (PrfEnv.find x env)))
WithProof.zero v
let frac_num n = n -/ Q.floor n
type ('a, 'b) hitkind =
| Forget
| Hit of 'a
| Keep of 'b
let violation sol vect =
let sol = Vect.set 0 Q.one sol in
let c = Vect.get 0 vect in
if Q.zero =/ c then Vect.dotproduct sol vect
else Q.abs (Vect.dotproduct sol vect // c)
let cut env rmin sol (rst : Restricted.t) tbl (x, v) =
let n, r = Vect.decomp_cst v in
let fn = frac_num (Q.abs n) in
if fn =/ Q.zero then Forget
else
let n0_5 = Q.one // Q.two in
let ccoeff_prop2 =
let tmin =
if fn </ n0_5 then
Q.ceiling (n0_5 // fn)
else
Q.neg (Q.ceiling (n0_5 // (Q.one -/ fn)))
in
("Prop2", fun v -> frac_num (v */ tmin))
in
let ccoeff = ccoeff_prop2 in
let cut_vector ccoeff =
Vect.fold
(fun acc x n ->
if Restricted.is_restricted x rst then Vect.set x (ccoeff n) acc
else acc)
Vect.null r
in
let lcut =
( fst ccoeff
, make_farkas_proof env (Vect.normalise (cut_vector (snd ccoeff))) )
in
let check_cutting_plane (p, c) =
match WithProof.cutting_plane c with
| None ->
if debug then
Printf.printf "%s: This is not a cutting plane for %a\n%a:" p
LinPoly.pp_var x WithProof.output c;
None
| Some (v, prf) ->
if debug then (
Printf.printf "%s: This is a cutting plane for %a:" p LinPoly.pp_var x;
Printf.printf "(viol %f) %a\n"
(Q.to_float (violation sol (fst v)))
WithProof.output (v, prf) );
Some (x, (v, prf))
in
match check_cutting_plane lcut with
| Some r -> Hit r
| None ->
let has_unrestricted =
Vect.fold
(fun acc v vl -> acc || not (Restricted.is_restricted v rst))
false r
in
if has_unrestricted then Keep (x, v) else Forget
let merge_result_old oldr f x =
match oldr with
| Hit v -> Hit v
| Forget -> (
match f x with Forget -> Forget | Hit v -> Hit v | Keep v -> Keep v )
| Keep v -> (
match f x with Forget -> Keep v | Keep v' -> Keep v | Hit v -> Hit v )
let merge_best lt oldr newr =
match (oldr, newr) with
| x, Forget -> x
| Hit v, Hit v' -> if lt v v' then Hit v else Hit v'
| _, Hit v | Hit v, _ -> Hit v
| Forget, Keep v -> Keep v
| Keep v, Keep v' -> Keep v'
let find_cut nb env u sol rst tbl =
if nb = 0 then
IMap.fold
(fun x v acc -> merge_result_old acc (cut env u sol rst tbl) (x, v))
tbl Forget
else
let lt (_, ((v1, _), p1)) (_, ((v2, _), p2)) =
ProofFormat.pr_size p1 </ ProofFormat.pr_size p2
in
IMap.fold
(fun x v acc -> merge_best lt acc (cut env u sol rst tbl (x, v)))
tbl Forget
let find_split env tbl rst =
let is_split x v =
let v, n =
let n, _ = Vect.decomp_cst v in
if Restricted.is_restricted x rst then
let n', v = Vect.decomp_cst (fst (fst (PrfEnv.find x env))) in
(v, n -/ n')
else (Vect.set x Q.one Vect.null, n)
in
if Restricted.is_restricted x rst then None
else
let fn = frac_num n in
if fn =/ Q.zero then None
else
let fn = Q.abs fn in
let score = Q.min fn (Q.one -/ fn) in
let vect = Vect.add (Vect.cst (Q.neg n)) v in
Some (Vect.normalise vect, score)
in
IMap.fold
(fun x v acc ->
match is_split x v with
| None -> acc
| Some (v, s) -> (
match acc with
| None -> Some (v, s)
| Some (v', s') -> if s' >/ s then acc else Some (v, s) ))
tbl None
let var_of_vect v = fst (fst (Vect.decomp_fst v))
let eliminate_variable (bounded, env, tbl) x =
if debug then
Printf.printf "Eliminating variable %a from tableau\n%a\n" LinPoly.pp_var x
output_tableau tbl;
let vr = LinPoly.MonT.get_fresh () in
let vr1 = LinPoly.MonT.get_fresh () in
let vr2 = LinPoly.MonT.get_fresh () in
let z = LinPoly.var vr1 in
let zv = var_of_vect z in
let t = LinPoly.var vr2 in
let tv = var_of_vect t in
let xdef = Vect.add z (Vect.uminus t) in
let xp = ((Vect.set x Q.one (Vect.uminus xdef), Eq), Def vr) in
let zp = ((z, Ge), Def zv) in
let tp = ((t, Ge), Def tv) in
let tbl = IMap.map (fun v -> Vect.subst x xdef v) tbl in
let env =
PrfEnv.map
(fun lp ->
let (v, o), p = lp in
let ai = Vect.get x v in
if ai =/ Q.zero then lp
else WithProof.addition (WithProof.mult (Vect.cst (Q.neg ai)) xp) lp)
env
in
let env =
PrfEnv.set_prf vr xp (PrfEnv.set_prf zv zp (PrfEnv.set_prf tv tp env))
in
let bounded = IMap.add x (vr, zv, tv) bounded in
if debug then (
Printf.printf "Tableau without\n %a\n" output_tableau tbl;
Printf.printf "Environment\n %a\n" output_env env );
(bounded, env, tbl)
let integer_solver lp =
let insert_row vr v rst tbl =
match push_real true vr v rst tbl with
| Sat (t', x) ->
Inl (Restricted.restrict vr rst, t', x)
| Unsat c -> Inr c
in
let vr0 = LinPoly.MonT.get_fresh () in
let l', env =
PrfEnv.of_list [] PrfEnv.empty (List.rev_map WithProof.of_cstr lp)
in
let nb = ref 0 in
let rec isolve env cr res =
incr nb;
match res with
| Inr c ->
Some
(Step
( LinPoly.MonT.get_fresh ()
, make_farkas_certificate env (Vect.normalise c)
, Done ))
| Inl (rst, tbl, x) -> (
if debug then begin
Printf.fprintf stdout "Looking for a cut\n";
Printf.fprintf stdout "Restricted %a ...\n" Restricted.pp rst;
Printf.fprintf stdout "Current Tableau\n%a\n" output_tableau tbl;
flush stdout
end;
if !nb mod 3 = 0 then
match find_split env tbl rst with
| None ->
None
| Some (v, s) -> (
let vr = LinPoly.MonT.get_fresh () in
let wp1 = ((v, Ge), Def vr) in
let wp2 = ((Vect.mul Q.minus_one v, Ge), Def vr) in
match (WithProof.cutting_plane wp1, WithProof.cutting_plane wp2) with
| None, _ | _, None ->
failwith "Error: splitting over an integer variable"
| Some wp1, Some wp2 -> (
if debug then
Printf.fprintf stdout "Splitting over (%s) %a:%a or %a \n"
(Q.to_string s) LinPoly.pp_var vr WithProof.output wp1
WithProof.output wp2;
let v1', v2' = (fst (fst wp1), fst (fst wp2)) in
if debug then
Printf.fprintf stdout "Solving with %a\n" LinPoly.pp v1';
let res1 = insert_row vr v1' (Restricted.set_exc vr rst) tbl in
let prf1 = isolve (IMap.add vr ((v1', Ge), Def vr) env) cr res1 in
match prf1 with
| None -> None
| Some prf1 ->
let prf', hyps = ProofFormat.simplify_proof prf1 in
if not (ISet.mem vr hyps) then Some prf'
else (
if debug then
Printf.fprintf stdout "Solving with %a\n" Vect.pp v2';
let res2 = insert_row vr v2' (Restricted.set_exc vr rst) tbl in
let prf2 =
isolve (IMap.add vr ((v2', Ge), Def vr) env) cr res2
in
match prf2 with
| None -> None
| Some prf2 -> Some (Split (vr, v, prf1, prf2)) ) ) )
else
let sol = find_full_solution rst tbl in
match find_cut (!nb mod 2) env cr sol rst tbl with
| Forget ->
None
| Hit (cr, ((v, op), cut)) -> (
let vr = LinPoly.MonT.get_fresh () in
if op = Eq then
Some (Step (vr, CutPrf cut, Done))
else
let res = insert_row vr v (Restricted.set_exc vr rst) tbl in
let prf =
isolve (IMap.add vr ((v, op), Def vr) env) (Some cr) res
in
match prf with
| None -> None
| Some p -> Some (Step (vr, CutPrf cut, p)) )
| Keep (x, v) -> (
if debug then
Printf.fprintf stdout "Remove %a from Tableau\n" LinPoly.pp_var x;
let bounded, env, tbl =
Vect.fold
(fun acc x n ->
if x <> 0 && not (Restricted.is_restricted x rst) then
eliminate_variable acc x
else acc)
(IMap.empty, env, tbl) v
in
let prf = isolve env cr (Inl (rst, tbl, None)) in
match prf with
| None -> None
| Some pf ->
Some
(IMap.fold
(fun x (vr, zv, tv) acc ->
ExProof (vr, zv, tv, x, zv, tv, acc))
bounded pf) ) )
in
let res = solve true l' (Restricted.make vr0) IMap.empty in
isolve env None res
let integer_solver lp =
nb_pivot := 0;
if debug then
Printf.printf "Input integer solver\n%a\n" WithProof.output_sys
(List.map WithProof.of_cstr lp);
match integer_solver lp with
| None ->
profile_info := (false, !nb_pivot) :: !profile_info;
None
| Some prf ->
profile_info := (true, !nb_pivot) :: !profile_info;
if debug then
Printf.fprintf stdout "Proof %a\n" ProofFormat.output_proof prf;
Some prf