Vous avez reçu un message "Your GitLab account has been locked ..." ? Pas d'inquiétude : lisez cet article https://docs.gricad-pages.univ-grenoble-alpes.fr/help/unlock/

store.ml 29.9 KB
Newer Older
1
(*-----------------------------------------------------------------------
2
** Copyright (C) 2002, 2003 - Verimag.
3
4
5
6
** This file may only be copied under the terms of the GNU Library General
** Public License 
**-----------------------------------------------------------------------
**
7
** File: store.ml
8
9
10
11
** Main author: jahier@imag.fr
*)

open Formula
12
open Value
13
open Constraint
14
open Polyhedron
15
16
open Util
open List
17
open Draw
18
        
19
let debug_store = false
20

21
22
(* exported *)
exception No_numeric_solution	
23

24

25

26
27
28
29
30
31

type vars_domain =
    Unsat
  | Range of (Formula.var_name, Polyhedron.range) Hashtbl.t
  | Poly of (Formula.vnt list * Poly.t * (int -> string)) list

32
type t = {
33
34
35
36
37
38
39
40
41
42
  (** 
    This field is used to store where each variable ranges.  It is
    set to [Unsat] when the system becomes unsatisfiable, namely,
    when the range for one of the variable becomes empty.

    Some variables are represented by Ranges (polyhedron of dimension
    one). Some others by plain Polyhedron. The idea is that, at bdd
    leaves, if it remains some delayed constraints, we switch to a
    polyhedron representation.  *)
  var : vars_domain ;
43
44
45
46
47
48
49
50
51
52
53
54
  
  (**
    This field is used to substitute a variable by an expression. This 
    is to deal with equalities: when an equality is encountered, 
    we can remove one dimension by putting the equality into a 
    such a substitution. 
    
    Then we apply it to all the other relations. the value of the
    substituted varaible is then obtained once the other var have
    been drawn. 

    We add an element to this list if 
55
56
57
    - an equality is encountered during the drawing/bdd-traversal
    - whenever a variable become bounded (1) after a constraint is 
      added to the store
58
59
60
61
62
63
64
65
66
67
68
69
70
71
    
    (1) i.e., when the interval is reduced to one single point
  *) 
  substl : Ne.subst list;
  
  (** 
    When the dimension of an atomic formula is greater than 1, we
    delay its addition to the store until an equality makes it a
    constraint of dimension 1 (i.e., it contains only 1 var). At bdd
    leaves, if this list is not empty, it means that the current
    formula cannot be solved with an interval based solver.
    In that case, we use a polyhedron solver.
  *)
  delay : Constraint.ineq list
72
}
73
74
75


(* exported *)
76
let (new_store : Formula.vnt list -> t) =
77
  fun vnt_l -> 
78
79
    let (add_one_var : (Formula.var_name, range) Hashtbl.t -> Formula.vnt 
	   -> unit) =
80
81
      fun tbl (vn, vt) -> 
	let range = (
82
83
84
85
	  match vt with
	      Formula.IntT(min, max) -> RangeI(min, max)
	    | Formula.FloatT(min, max) -> RangeF(min, max) 
	    | Formula.BoolT -> assert false
86
87
	)
	in
88
	  Hashtbl.replace tbl vn range
89
    in
90
91
92
    let tbl = (Hashtbl.create (List.length vnt_l)) in
      List.iter (add_one_var tbl) vnt_l;
      {
93
	var = Range(tbl) ;
94
95
96
	substl = [];
	delay = []
      }
97
98

(* exported *)
99
let (copy : t -> t) =
100
101
102
103
104
105
106
107
108
109
110
111
112
  fun store ->
    match store.var with
	Unsat -> { var = Unsat; substl = store.substl ; delay = store.delay }
      | Range tbl -> 
	  { 
	    var = Range (Hashtbl.copy tbl) ; 
	    substl = store.substl ; 
	    delay = store.delay
	  }
      | Poly p -> 
	  (* ougth to be useless *)
	  assert false

113
114

(* Normalised atomic constraints *)
115
116
117
118
119
type nac = 
  | NSupF   of float (** >  *)
  | NSupEqF of float (** >= *)
  | NInfF   of float (** <  *)
  | NInfEqF of float (** <= *)
120
  | NEqF    of float (** = *)
121
122
123

  | NSupEqI of int (** >=  *)
  | NInfEqI of int (** <=  *)
124
  | NEqI    of int (** = *)
125

126
let eps = epsilon_float    
127

128
129
130
(****************************************************************************)
(* Pretty printing   *)

131

132
133
134
135
136
137
138
139
140
let (range_to_string : range -> string) =
  fun range -> 
    match range with 
    RangeI(min, max) ->
      ("[" ^ (string_of_int min) ^ ", " ^ (string_of_int max) ^ "] ")
  | RangeF(min, max) ->
      ("[" ^ (string_of_float min) ^ ", " ^ (string_of_float max) ^ "] ")

(* exported *)
141
let (to_string : t -> string) =
142
143
144
  fun s -> 
    let var_str = ("\n*** Variable ranges: \n" ^ 
      match s.var with
145
	  Unsat -> "Empty store"
146
	| Poly _ -> " XXX Store.to_string:  finish me ! " 
147
	| Range(tbl) -> 
148
149
150
151
152
153
154
155
156
157
158
	    (Hashtbl.fold
	       (fun vn range acc ->
		  ("   " ^ vn ^ " in " ^ (range_to_string range) ^ "\n" ^ acc)
	       )
	       tbl
	       "\n"
	    )
		  )
    and substl_str = ("\n*** Substitutions: \n" ^ Ne.substl_to_string s.substl)
    and delay_str = ("\n*** Delayed constraints: \n" ^ 
      List.fold_left 
159
	(fun acc d -> acc ^ "\n" ^ (Constraint.ineq_to_string d))
160
161
162
163
164
165
	""
	s.delay)
    in
      (var_str ^ substl_str ^ delay_str)

(****************************************************************************)
166

167
168
169
170
171
172
173
174
175
176
177
(* 
   Note that we check the satisfiability of constraints over
   polyhedra at bdd leaves, which, in some circumstances, migth be
   inefficient. The point is that, if we chose to check the formula
   satisfiability during the bdd traversal, we take the risk that a
   very polyhedron is created whereas it was not necessary (because
   of forthcoming equalities that would reduce its dimension). And
   creating polyhedron with too big (>15) dimensions simply runs
   forever, which is really bad.
*)
(* exported *)
178
let (switch_to_polyhedron_representation : t -> t * t) =
179
180
181
182
183
184
185
186
187
188
189
190
  fun store ->
    (* handle delayed constraints using polyhedron *)
    match store.var with
	Unsat -> raise No_numeric_solution (* this ougth to be dead code ... *)
      | Poly _ -> assert false
      | Range tbl -> 
	  if 
	    store.delay = [] 
	  then 
	    (store, {var = Poly []; substl = []; delay = [] })
	  else
	    let poly_l = 
191
192
	      (* This functions removes from [store] variables that are 
		 moved to the polyhedron store *)
193
194
195
196
	      Polyhedron.build_poly_list_from_delayed_cstr tbl store.delay 
	    in
	      List.iter 
		(fun (_, poly, _) -> 
197
198
199
200
201
202
		   if Poly.is_empty poly then (
		     if debug_store then (
		       print_string "\n The polyhedron is empty .\n"; 
		       flush stdout );
		     raise No_numeric_solution
		   )
203
204
205
206
207
208
209
		)
		poly_l;
	      (store, {var = Poly poly_l; substl = []; delay = []})
		

(****************************************************************************)

210
211
212
213
214
215
216
217
218
219
let (div : int -> int -> int) =
  fun x y -> 
    (* I define my own integer division as the division of Pervasives
       does not consistently rounds its result (ie, the result is
       round to the least integer if it is positive, and to the
       greatest integer if it is negative). *)
    let xf = float_of_int x
    and yf = float_of_int y
    in 
      int_of_float (floor (xf /.yf))
220

221

222
let (normalise : Constraint.ineq -> var_name * nac ) =
223
  fun cstr -> 
224
    (* Transform atomic formula into a data type that is easier to
225
       process.
226

227
228
       Fails if [cstr] contains more than one variable (in which
       case the constraint should have been delayed).
229
    *)
230
231
232
    let (get_vn_and_constant : Ne.t -> (
	     Value.num  (* The constant *)
	   * Value.num  (* The coefficient of the variable *)
233
234
235
	   * var_name   (* The name of the variable *)
	 ) 
	) =
236
      fun ne -> 
237
	let list = Ne.fold (fun vn num acc -> (vn,num)::acc) ne [] in
238
	  match list with
239
240
241
242
243
244
245
246
247
	      (* 0 var *)
	      [("", cst)] -> 
		( match cst with
		      I(_) -> (cst, I(0) , "")
		    | F(_) -> (cst, F(0.), "")
		)

	      (* 1 var *)
	    | [("", cst); (vn, coeff)] -> (cst, coeff, vn)
248
249
250
251
252
253
	    | [(vn, coeff); ("", cst)] -> (cst, coeff, vn)
	    | [(vn, coeff)] -> 
		( match coeff with
		      I(_) -> (I(0), coeff, vn)
		    | F(_) -> (F(0.), coeff, vn)
		)
254
255
256

	      (* more than 1 var *)
	    | _ -> 
257
		assert false
258
    in
259
      match cstr with
260
	  GZ(ne) ->
261
262
263
264
265
266
267
268
269
270
271
272
273
	    let (cst, coeff, vn) = get_vn_and_constant ne in
	      ( match (cst, coeff) with
		    (I(i_cst), I(i_coeff)) -> 
		      let i = div (- i_cst)  i_coeff in 
			if i_coeff > 0
			then 
			  if (i_cst mod i_coeff) = 0
			  then (vn, NSupEqI(i+1))
			  else (vn, NSupEqI(i))
			else
			  if (i_cst mod i_coeff) = 0
			  then (vn, NInfEqI(i-1))
			  else (vn, NInfEqI(i-1))
274

275
276
277
278
		  | (F(f_cst), F(f_coeff)) ->
		      if f_coeff > 0.
		      then (vn, NSupF(-.f_cst /. f_coeff))
		      else (vn, NInfF(-.f_cst /. f_coeff))
279

280
281
		  | (I(i), F(f)) -> 
		      failwith ("*** Error: " ^ (string_of_int i) 
282
				^ " is an integer and " 
283
284
285
				^ (string_of_float f) ^ " is a float.\n")
		  | (F(f), I(i)) ->  
		      failwith ("*** Error: " ^ (string_of_int i) 
286
				^ " is an integer and " 
287
				^ (string_of_float f) ^ " is a float.\n")
288
	      )
289

290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
	| GeqZ(ne) -> 
	    let (cst, coeff, vn) = get_vn_and_constant ne in
	      ( match (cst, coeff) with
		    (I(i_cst), I(i_coeff)) ->  
		      let i =  div (- i_cst)  i_coeff in 
			if i_coeff > 0
			then 
			  if (i_cst mod i_coeff) = 0
			  then (vn, NSupEqI(i))
			  else (vn, NSupEqI(i))
			else
			  if (i_cst mod i_coeff) = 0
			  then (vn, NInfEqI(i))
			  else (vn, NInfEqI(i-1))
		  | (F(f_cst), F(f_coeff)) -> 
		      if f_coeff > 0.
		      then (vn, NSupEqF(-.f_cst /. f_coeff))
		      else (vn, NInfEqF(-.f_cst /. f_coeff))
		  | (I(_), F(_)) -> assert false
		  | (F(_), I(_)) -> assert false
	      )
311

312
313
314



315
let (make_subst : vn -> Value.num -> Ne.subst) =
316
317
  fun vn value -> 
    (* returns a subst from [vn] to [value] *)
318
319
320
321
322
323
324
325
326
327
328
329
    match value with
	I _ -> ((vn, (I 1)), Ne.make "" value)
      | F _ -> ((vn, (F 1.)), Ne.make "" value)



(** if a constraint [cstr] = [GeqZ(ne)] is such that the store
  contains the constraint [eqZ(-ne)] among its delayed variables,
  then [cstr] turns out to be an equality. In that case, this
  function returns [ne] as well as the store with the delayed
  constraint [eqZ(-ne)] removed. *)

330
(* let (is_ineq_cstr_an_eq : Constraint.t -> t -> (t * Ne.t) option) = *)
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
(*   fun cstr store -> *)
(*     let ne = *)
(*       match cstr with *)
(* 	  GZ(ne) -> ne *)
(* 	| GeqZ(ne) -> ne *)
(* 	| _ -> assert false *)
(*     in *)
(*     let rec get_cstr ne d acc = *)
(*       match d with *)
(* 	  [] -> raise Not_found *)
(* 	| GZ(ne2)::dtail ->  *)
(* 	    if ne = ne2 *)
(* 	    then GZ(ne2), List.rev_append acc dtail *)
(* 	    else get_cstr ne dtail (GZ(ne2)::acc) *)
(* 	| GeqZ(ne2)::dtail -> *)
(* 	    if ne = ne2 *)
(* 	    then GeqZ(ne2), List.rev_append acc dtail *)
(* 	    else get_cstr ne dtail (GeqZ(ne2)::acc) *)
(* 	| _ -> assert false *)
(*     in *)
(*     let d2 = *)
(*       try  *)
(* 	match  cstr, (get_cstr ne d []) with *)
(* 	  | GeqZ(ne), (GZ(_),  d2) -> GZ(ne)::d2 *)
(* 	      (* Indeed, GeqZ(ne) => GZ(ne) *) *)
(* 	  | _ -> d *)
(* 	with  *)
(* 	    Not_found -> d *)
(* 	in *)
(* 	let new_d = *)
(* 	  try  *)
(* 	    if *)
(* 	      match cstr, (get_cstr (Ne.neg ne) d []) with *)
(* 		  GeqZ(ne), (GeqZ(_), d2) -> true *)
(* 		| _ -> false *)
(* 	    then *)
(* 	       *)
(* 	    else *)
(* 	       *)
(* 	  with  *)
(* 	      Not_found -> d *)
(* 	in *)
373

374

375
(* exported *)
376
let rec (add_constraint : t -> Constraint.t -> t) =
377
378
379
380
381
  fun store cstr0 ->
    let cstr = Constraint.apply_substl store.substl cstr0 in      
      add_constraint_do store cstr

and (add_constraint_do : t -> Constraint.t -> t) =
382
  fun store cstr ->
383
384
    let _ =
      if debug_store then (
385
386
	print_string "\n***********************************************\n " ;   
	print_string (to_string  store);   
387
388
389
390
391
392
	print_string "\n**** add_constraint: " ;  
	print_string (Constraint.to_string cstr);
(* 	print_string "\n"; *)
	flush stdout 
      );  
    in
393
    let (var, sl, d) = (store.var, store.substl, store.delay) in
394
    let dim = Constraint.dimension cstr in
395
396
397
      if 
	dim = 0
      then
398
	( match cstr with
399
400
401
	      EqZ(ne) -> add_eq_to_store store ne
	    | Bv _ -> assert false
	    | Ineq GZ(ne) ->
402
		if 
403
404
405
406
407
408
		  ( match (Ne.find "" ne) with 
			Some v -> Value.num_sup_zero v
		      | None -> false 
		  )
		then 
		  store
409
410
411
		else
		  raise No_numeric_solution
	    | Ineq GeqZ(ne) -> 
412
		if 
413
414
415
416
417
418
419
420
		  ( match (Ne.find "" ne) with 
			Some v -> Value.num_supeq_zero v
		      | None -> true 
		  )
		then 
		  store
		else 
		  raise No_numeric_solution
421
	)
422
423
424
      else if 
	dim > 1
      then
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
	(* 
	   We could also choose not to delay those constraints and
	   switch to the polyhedron representation for the concerned
	   variables there.

	   What is the better solution really ???  
	*)
	( match cstr with
	      EqZ(ne) -> add_eq_to_store store ne
	    | Bv _ -> assert false
	    | Ineq ineq ->
		if debug_store then (
		  print_string "\n ==> delay  " ;  
		  print_string (Constraint.to_string cstr);
		  flush stdout 
		);
		{ var=var ; substl=sl ; delay = ineq::d }
	)
443
444
445
      else
	(* dim = 1 *)
	match store.var with
446
	    Unsat ->
447
	      raise No_numeric_solution
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
	  | Poly pl -> 
	      (* 
		 Polka is currently only called at bdd leaves, not
		 during the bdd traversal. Of course, if this changes,
		 some more code ougth to be written here.
	      *)
	      assert false
	  | Range(tbl) -> 
	      ( match cstr with
		    EqZ(ne) -> add_eq_to_store store ne
		  | Bv _ -> assert false
		  | Ineq ineq ->
		      let (vn, nac) = normalise ineq in
			( match (nac, (Hashtbl.find tbl vn)) with 
			      (NSupEqI(i), RangeI(min, max)) -> 
				if 
				  i <= min
				then 
				  {var=Range(tbl) ; substl=sl ; delay=d }
				else if 
				  i > max
				then 
				  {var=Unsat ; substl=sl ; delay=d }
				else (* min < i <= max *)
				  let sl1, d1, da = 
				    if 
				      i = max
				    then 
				      (*
					Whenever, a variable becomes bounded, we:
					- add it to the list of substitutions
					- remove it from the store.
					- check if delayed cstr should be awaked
					once this new subst has been applied
				      *)
				      let s = make_subst vn (I max) in 
				      let (d_awake, d_delay) = 
					List.partition
486
487
488
489
490
					  (fun ineq -> 
					     Constraint.dimension_ineq ineq <= 1)
					  (List.map 
					     (Constraint.apply_subst_ineq s)
					     d)
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
				      in
					Hashtbl.remove tbl vn;
					(s::sl, d_delay, d_awake)
				    else 
				      (
					Hashtbl.replace tbl vn (RangeI(i, max)); sl, 
					d,
					[]
				      )
				  in 
				    (* Applying the waked constraints *)
				    List.fold_left
				      (fun acc cstr -> 
					 if debug_store then (
					   print_string "\n <== awake ";
					   print_string (Constraint.ineq_to_string cstr);
					   flush stdout
					 );
509
					 add_constraint_do acc (Ineq cstr))
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
				      { var=Range(tbl) ; substl=sl1 ; delay=d1 }
				      da
			    | 
			      (NInfEqI(i), RangeI(min, max)) -> 
				if 
				  i < min
				then 
				  { var=Unsat ; substl=sl ; delay=d }
				else if 
				  i >= max
				then
				  { var=Range(tbl) ; substl=sl ; delay=d }
				else (* min <= i < max *)
				  let sl1, d1, da =
				    if 
				      i = min 
				    then 
				      let s = make_subst vn (I min) in
				      let (d_awake, d_delay) = 
					List.partition
530
531
532
533
534
					  (fun ineq -> 
					     Constraint.dimension_ineq ineq <= 1)
					  (List.map 
					     (Constraint.apply_subst_ineq s) 
					     d)
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
				      in
					Hashtbl.remove tbl vn;
					(s::sl, d_delay, d_awake)
				    else 
				      (
					Hashtbl.replace tbl vn (RangeI(min, i)); sl, 
					d,
					[]
				      )
				  in
				    (* Applying the waked constraints *)
				    List.fold_left
				      (fun acc cstr -> 
					 if debug_store then (
					   print_string "\n <== awake ";
					   print_string (Constraint.ineq_to_string cstr);
					   flush stdout
					 );
553
					 add_constraint_do acc (Ineq cstr))
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
				      { var=Range(tbl) ; substl=sl1 ; delay=d1 }
				      da
				      
			    (** **)
			    |  (NSupEqF(f), RangeF(min, max)) -> 
				 if 
				   f <= min
				 then
				   {var=Range(tbl) ; substl=sl ; delay=d }
				 else if 
				   f > max
				 then
				   {var=Unsat ; substl=sl ; delay=d }
				 else (* min < f <= max *)
				   let sl1, d1, da =
				     if 
				       f = max
				     then  
				       let s = make_subst vn (F max) in
				       let (d_awake, d_delay) = 
					 List.partition
575
576
577
578
579
					   (fun ineq -> 
					      Constraint.dimension_ineq ineq <= 1)
					   (List.map 
					      (Constraint.apply_subst_ineq s) 
					      d)
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
				       in
					 Hashtbl.remove tbl vn;
					 (s::sl, d_delay, d_awake)
				     else  
				       (
					 Hashtbl.replace tbl vn (RangeF(f, max)); sl,
					 d,
					 []
				       )
				   in
				     (* Applying the waked constraints *)
				     List.fold_left
				       (fun acc cstr -> 
					 if debug_store then (
					   print_string "\n <== awake ";
					   print_string (Constraint.ineq_to_string cstr);
					   flush stdout
					 );
598
					  add_constraint_do acc (Ineq cstr))
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
				       { var=Range(tbl) ; substl=sl1 ; delay=d1 }
				       da
				       
			    | 
			      (NInfEqF(f), RangeF(min, max)) -> 
				if 
				  f < min
				then 
				  {var=Unsat ; substl=sl ; delay=d }
				else if 
				  f >= max
				then
				  {var=Range(tbl) ; substl=sl ; delay=d }
				else (* min <= f < max *)
				  let sl1, d1, da =
				    if 
				      f = min 
				    then  
				      let s = make_subst vn (F min) in
				      let (d_awake, d_delay) = 
					List.partition
620
621
622
623
624
					  (fun ineq -> 
					     Constraint.dimension_ineq ineq <= 1)
					  (List.map
					     (Constraint.apply_subst_ineq s) 
					     d)
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
				      in
					Hashtbl.remove tbl vn;
					(s::sl, d_delay, d_awake)
				    else   
				      (
					Hashtbl.replace tbl vn (RangeF(min, f)); sl,
					d,
					[]
				      )
				  in
				    (* Applying the waked constraints *)
				    List.fold_left
				      (fun acc cstr -> 
					 if debug_store then (
					   print_string "\n <== awake ";
					   print_string (Constraint.ineq_to_string cstr);
					   flush stdout
					 );
643
					 add_constraint_do acc (Ineq cstr))
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
				      { var=Range(tbl) ; substl=sl1 ; delay=d1 }
				      da
		  		      
			    |  (NSupF(f), RangeF(min, max)) -> 
				 if 
				   f < min
				 then
				   {var=Range(tbl) ; substl=sl ; delay=d }
				 else if 
				   f >= max
				 then 
				   {var=Unsat ; substl=sl ; delay=d }
				 else (* min <= f < max *)
				   let sl1, d1, da =
				     if 
				       (f +. eps) = max
				     then
				       let s = make_subst vn (F max) in
				       let (d_awake, d_delay) = 
					 List.partition
664
665
666
667
668
					   (fun ineq ->
					      Constraint.dimension_ineq ineq <= 1)
					   (List.map
					      (Constraint.apply_subst_ineq s) 
					      d)
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
				       in
					 Hashtbl.remove tbl vn;
					 (s::sl, d_delay, d_awake)
				     else   
				       (
					 Hashtbl.replace tbl vn (RangeF(f+.eps, max)); sl,
					 d,
					 []
				       )
				   in
				     (* Applying the waked constraints *)
				     List.fold_left
				       (fun acc cstr -> 
					 if debug_store then (
					   print_string "\n <== awake ";
					   print_string (Constraint.ineq_to_string cstr);
					   flush stdout
					 );
687
					  add_constraint_do acc (Ineq cstr))
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
				       { var=Range(tbl) ; substl=sl1 ; delay=d1 }
				       da
				       
			    | 
			      (NInfF(f), RangeF(min, max)) -> 
				if 
				  f <= min
				then 
				  {var=Unsat ; substl=sl ; delay=d }
				else if 
				  f > max
				then
				  {var=Range(tbl) ; substl=sl ; delay=d }
				else (* min < f <= max *)
				  let sl1, d1, da =
				    if 
				      (f -. eps) = min 
				    then
				      let s = make_subst vn (F min) in
				      let (d_awake, d_delay) = 
					List.partition
709
710
711
712
713
					  (fun ineq -> 
					     Constraint.dimension_ineq ineq <= 1)
					  (List.map
					     (Constraint.apply_subst_ineq s)
					     d)
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
				      in
					Hashtbl.remove tbl vn;
					(s::sl, d_delay, d_awake)
				    else   
				      (
					Hashtbl.replace tbl vn (RangeF(min, f-.eps)); sl,
					d,
					[]
				      )
				  in
				    (* Applying the waked constraints *)
				    List.fold_left
				      (fun acc cstr -> 
					 if debug_store then (
					   print_string "\n <== awake ";
					   print_string (Constraint.ineq_to_string cstr);
					   flush stdout
					 );
732
					 add_constraint_do acc (Ineq cstr))
733
734
735
736
737
738
				      { var=Range(tbl) ; substl=sl1 ; delay=d1 }
				      da
				      
			    | _ -> assert false
			)
	      )
739

740
and (add_eq_to_store : t -> Ne.t -> t) =
741
  fun store ne -> 
742
743
744
745
746
747
748
749
750
751
752
    (* 
       [add_eq_to_store s e] returns the store [s] with the numeric
       constraint [EqZ(e)] added.
    *)
    let dim = Ne.dimension ne in
      if 
	dim = 0
      then
	let unsat_store = { var = Unsat ; substl = [] ; delay = [] } in
	  match (Ne.find "" ne) with
	      Some x -> 
753
754
		 if Value.num_eq_zero x then store else unsat_store
	     | None -> store
755
756
757
      else 
	(* dim > 0 *)
	let (var, sl, d) = (store.var, store.substl, store.delay) in    
758
	let (vn_val_opt, ne_rest) = Ne.split ne in
759
760
761
762
763
	  match vn_val_opt with
	      None ->
		raise No_numeric_solution
	    | Some (vn, nval0) ->
		let nval = Value.neg nval0 in
764
		let s = ((vn, nval), ne_rest) in
765
		let d2 = List.map (Constraint.apply_subst_ineq s) d in
766
		  (* some delayed constraints may have been awaken *)
767
768
		let (waked, new_d) = 
		  List.partition 
769
		    (fun cstr -> Constraint.dimension_ineq cstr <= 1) d2
770
771
772
773
		in
		let new_sl = s::sl in
		let range_vn, new_var = 
		  match var with
774
775
776
777
778
779
		      Unsat -> assert false
		    | Poly pl -> assert false
		    | Range(tbl) ->
			let tbl2 = Hashtbl.copy tbl in
			  Hashtbl.remove tbl2 vn;
			  (Hashtbl.find tbl vn, Range(tbl2))
780
781
782
		in
		let new_store = 
		  List.fold_left
783
784
785
786
787
788
		    (fun acc cstr -> 
		       if debug_store then (
			 print_string "\n <== awake ";
			 print_string (Constraint.ineq_to_string cstr);
			 flush stdout
		       );
789
		       add_constraint_do acc (Ineq cstr))
790
791
792
		    { var = new_var; substl = new_sl; delay = new_d }
		    waked
		in
793
794
795
796
		  if 
		    Ne.is_a_constant ne_rest 
		  then
		    new_store
797
798
799
		  else
		    let c1, c2 =
		      match range_vn with
800
801
802
803
804
805
806
807
808
809
810
811
			  RangeI(min, max) ->
			    let min1, max1 =
			      Ne.diff (Ne.make "" (Value.quot_num (I (-min)) nval0)) ne_rest,
			      Ne.add  (Ne.make "" (Value.quot_num (I  max) nval0))   ne_rest
			    in
			      if 
				Value.num_sup_zero nval0
			      then
				GeqZ(min1), GeqZ(max1)
			      else
				GeqZ(Ne.neg_nexpr min1), GeqZ(Ne.neg_nexpr max1)

812
			| RangeF(min, max) ->
813
814
815
816
817
818
819
820
821
822
			    let min1, max1 =
			      Ne.diff (Ne.make "" (Value.quot_num (F (-.min)) nval0)) ne_rest,
			      Ne.add  (Ne.make "" (Value.quot_num (F  max)    nval0)) ne_rest
			    in
			      if 
				Value.num_sup_zero nval0
			      then
				GeqZ(min1), GeqZ(max1)
			      else
				GeqZ(Ne.neg_nexpr min1), GeqZ(Ne.neg_nexpr max1)
823
		    in
824
825
		      add_constraint_do
			(add_constraint_do new_store (Ineq c1)) (Ineq c2)
826
		      
827
    
828
(*******************************************************************************)
829
(* exported *)
830
let (is_store_satisfiable : t -> bool) =
831
  fun store -> 
832
    store.var <> Unsat
833
834
  

835
836
(*******************************************************************************)

837
838
839
840
841
(** [get_subst_from_solved_system sl] ougth to take in input a list
  of pair "var_name_and_coef, normal_expr" such that the normal_expr
  is reduced to a constant c. It returns the list of substitution
  (var_name -> c/coef).
*)
842
let rec (get_subst_from_solved_system : t -> ((string * Value.num) * Ne.t) list -> 
843
	   (string * Value.num) list) =
844
  fun store sl ->
845
846
847
    match sl with
	[] -> []
      | ((vn, v), ne)::tail -> 
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
	  let _ = 
	    assert (
	      if 
		(Ne.dimension ne) > 1 
	      then
		(
		  print_string (to_string  store);   
		  print_string "\n*** The expression <<";
		  print_string (Ne.to_string ne) ;
		  print_string ">> ougth to be a constant.\n";
		  flush stdout;
		  false
		)
	      else 
		true 
	    )
	  in
865
866
867
868
869
	  let s = 
	    match (Ne.find "" ne) with
		Some(cst) -> (vn, Value.quot_num cst v) 
	      | None      -> (vn, Value.diff_num v v)
	  in
870
	    s::(get_subst_from_solved_system store tail)
871

872
873

let (deduce_remaining_vars : (string * Value.num) list -> 
874
       t -> (string * Value.num) list) =
875
  fun fsl store ->
876
877
    (* 
       Once the num var have been drawn, we can deduce the values of
878
879
880
       variables that have substituted when equalities were
       encountered. 
    *)
881
    let sl0 = store.substl in
882
883
    let sl1 = 
      List.map 
884
	(fun (x, ne) -> (x, List.fold_left (Ne.apply_simple_subst) ne fsl))
885
886
	sl0
    in
887
	get_subst_from_solved_system store sl1 
888

889

890
891
892
893
894
895
896
897
898
899
900
901
902
903
let print_points pl =
  List.iter
    (* Used for debugging... *)
    (fun fl ->   
       print_string "(";
       List.iter
	 (fun f -> print_float f ; print_string " ")
	 fl;  
       print_string ") "
    )	   
    pl;
  print_string "\n";
  flush stdout

904

905
906
907
908
909
910
911
912
let rec (draw_n_points : int -> point list -> point list) =
  fun n pl -> 
    (* draw n points among pl (ant point may be drawn several times) *)
    let _ = assert (pl <> []) in
      if n = 0 then [] else
	let i = Random.int (List.length pl) in
	let p = List.nth pl i in
	  (p::(draw_n_points (n-1) pl))
913

914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
let rec (draw_n_distinc_points : int -> point list -> point list) =
  fun n pl -> 
    (* draw n distinct points among pl *)
    let _ = assert (pl <> []) in
      if n = 0 then [] else
	let i = Random.int (List.length pl) in
	let p = List.nth pl i in
	let new_pl = List.filter (fun p1 -> p1 <> p) pl in
	  (p::(draw_n_distinc_points (n-1) new_pl))
		
let rec (draw_point : point list -> point) =
  fun pl -> 
    (* Draw values in the convex hull of points in pl *)
    match pl with
	[] -> assert false
      | [p] -> p
      | p1::p2::tail ->
	  let alpha = Random.float 1. in
	  let p = 
	    List.map2
	      (fun x y -> 
		 let xy = alpha *. x +. (1. -. alpha) *. y in
		   xy
	      )
	      p1
	      p2
	  in
	    draw_point (p::tail)

(* a generic draw_inside parametrized by the way how points are chosen
 in the list of generators.*)
let  (draw_inside_gen : (int -> point list -> point list) -> 
       t -> Formula.num_subst list -> Formula.num_subst list) =
  fun draw_n_pts store sl_prev -> 
    match store.var with 
	Unsat -> assert false
      | Poly vntl_poly_rtn_list -> 
	  (List.fold_left 
	     (fun sl_acc (vntl, poly, rank_to_name) -> 
953
		let pl = Draw.get_vertices poly rank_to_name in
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
		let lgt = List.length pl in
		let dim = Poly.dim poly in
		let pl2 = draw_n_pts (min lgt (dim+1)) pl in
		let point = draw_point pl2 in
		let sl =
(*		  print_points pl; *)
		  snd
		    (List.fold_left
		       (fun (j, sl) x ->
			  let vn = rank_to_name j in
			  let vt = List.assoc vn vntl in
			  let s =
			    match vt with
				BoolT -> assert false 
			      | FloatT(_,_) -> (vn, F(x)) 
			      | IntT(_,_)  ->  
				  (* XXX brrrr, this is wrong !! *)  
				  (vn, I(truncate x))
			  in
			    (j+1, s::sl)
		       )
		       (0, [])
		       point
		    )
		in
		  rev_append sl sl_acc
	     )
	     []
	     vntl_poly_rtn_list
	  )
      | Range(tbl) -> 
	  let sl0 = 
	    Hashtbl.fold
	      (fun vn range acc -> 
		 ( match range with
		     | RangeI(min, max) -> 
			 let ran = Random.int (max - min + 1) 
			 in
			   ((vn, I(min + ran))::acc)

		     | RangeF(min, max)  -> 
			 let n = max -. min in
			 let ran = Random.float n in
			   ((vn, F(min +. ran))::acc)
		 )
	      )
	      tbl
	      []
	  in
	  let sl = rev_append sl0 sl_prev in
	    List.rev_append sl (deduce_remaining_vars sl store)
1005
  
1006

1007
(* exported *)
1008
1009
1010
let (draw_inside : t -> Formula.num_subst list -> Formula.num_subst list) =
  draw_inside_gen (draw_n_distinc_points)

1011
let (draw_edges : t -> Formula.num_subst list -> Formula.num_subst list) =
1012
1013
1014
1015
  draw_inside_gen (draw_n_points)


(* exported *)
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
let (draw_vertices : t -> Formula.num_subst list -> Formula.num_subst list) =
  fun store sl_prev ->
    match store.var with 
	Unsat -> assert false
      | Poly vntl_poly_rtn_list -> 
	  (List.fold_left 
	     (fun sl_acc (vntl, poly, rank_to_name) -> 
		let pl = Draw.get_vertices poly rank_to_name in
		let i = Random.int (List.length pl) in
		let p = List.nth pl i in
		let sl =     
(* 		  print_points pl; *)
		  snd
		    (List.fold_left
		       (fun (j, sl) x -> 
			  let vn = rank_to_name j in
			  let vt = List.assoc vn vntl in
			  let s =
			    match vt with
				BoolT -> assert false 
			      | FloatT(_,_) -> (vn, F(x)) 
			      | IntT(_,_)  ->  
				  (* XXX brrrr, this is wrong !! *)  
				  (vn, I(truncate x))
			  in
			    (j+1, s::sl)
		       )
		       (0, [])
		       p
		    )
		in
		  rev_append sl sl_acc
	     )
	     []
	     vntl_poly_rtn_list
	  )
      | Range(tbl) -> 
	  let sl0 =
	    Hashtbl.fold
	      (fun vn range acc -> 
		 ( match range with
		     | RangeI(min, max) -> 
			 let ran = Random.int 2
			 in
			   if ran = 0 
			   then ((vn, I(min))::acc)
			   else ((vn, I(max))::acc)
		     | RangeF(min, max)  -> 
			 let  ran = Random.int 2 in
			   if ran = 0 
			   then ((vn, F(min))::acc)
			   else ((vn, F(max))::acc)
		 ))
	      tbl
	      []
	  in
	  let sl = rev_append sl0 sl_prev in
	    List.rev_append sl (deduce_remaining_vars sl store)

1075
1076
1077
1078
1079
1080
1081
1082
(* (* exported *) *)
(* let (draw_edges : t -> Formula.num_subst list -> Formula.num_subst list) = *)
(*   fun store sl_prev -> *)
(*     match store.var with  *)
(* 	Unsat -> assert false *)
(*       | Poly vntl_poly_rtn_list ->  *)
(* 	  (List.fold_left  *)
(* 	     (fun sl_acc (vntl, poly, rank_to_name) ->  *)
1083
(* 		let pl = Draw.get_vertices poly rank_to_name in *)
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
(* 		let i = Random.int (List.length pl) in *)
(* 		let k = Random.int (List.length pl) in *)
(* 		let p1 = List.nth pl i in *)
(* 		let p2 = List.nth pl k in *)
(* 		let alpha = Random.float 1. in *)
(* 		let sl =      *)
(* (* 		  print_points pl; *) *)
(* 		  snd *)
(* 		    (List.fold_left2 *)
(* 		       (fun (j, sl) x y -> *)
(* 			  let xy = alpha *. x +. (1. -. alpha) *. y in *)
(* 			  let vn = rank_to_name j in *)
(* 			  let vt = List.assoc vn vntl in *)
(* 			  let s = *)
(* 			    match vt with *)
(* 				BoolT -> assert false  *)
(* 			      | FloatT(_,_) -> (vn,F(xy))  *)
(* 			      | IntT(_,_)  ->   *)
(* 				  (* XXX brrrr, this is wrong !! *)   *)
(* 				  (vn, I(truncate xy)) *)
(* 			  in *)
(* 			    (j+1, s::sl) *)
(* 		       ) *)
(* 		       (0, []) *)
(* 		       p1 p2 *)
(* 		    ) *)
(* 		in *)
(* 		  rev_append sl sl_acc *)
(* 	     ) *)
(* 	     [] *)
(* 	     vntl_poly_rtn_list *)
(* 	  ) *)
(*       | Range(tbl) ->  *)
(* 	  let sl0 = *)
(* 	    Hashtbl.fold *)
(* 	      (fun vn range acc ->  *)
(* 		 ( match range with *)
(* 		     | RangeI(min, max) ->  *)
(* 			 let ran = Random.int (max - min + 1)  *)
(* 			 in *)
(* 			   ((vn, I(min + ran))::acc) *)
(* 		     | RangeF(min, max)  ->  *)
(* 			 let n = max -. min in *)
(* 			 let ran = Random.float n in *)
(* 			   ((vn, F(min +. ran))::acc) *)
(* 		 )) *)
(* 	      tbl *)
(* 	      [] *)
(* 	  in *)
(* 	  let sl = rev_append sl0 sl_prev in *)
(* 	    List.rev_append sl (deduce_remaining_vars sl store) *)
1135

1136
(* exported *)
1137
(* let (get_all_vertices : t -> Formula.num_subst list list) =
1138
1139
1140
  fun  ->  *)
  
(* exported *)
1141
(* let (get_center_of_gravity : t -> Formula.num_subst list) =
1142
  fun  ->  *)
1143