module Cs = Cstr.Rat

module Debug = DebugTypes.Debug(struct let name = "PLPCore" end)
module Profile = Profile.Profile(struct let name = "PLPCore" end)

module Vec = Cs.Vec
module Minimization = Min.Classic(Vec)
module Explore = PSplxExec.Explore
module Coeff = Vec.Coeff

module MapV = Map.Make(struct type t = int let compare = Stdlib.compare end)

type region_t = Cone | NCone

module Boundary = struct
	type t = Cs.t * Vec.t

	let equal : t -> t -> bool
		= fun (c1,_) (c2,_) ->
		Cs.equalSyn c1 c2

	let get_cstr : t -> Cs.t
		= fun (c,_) ->
		c

	let get_cstrs : t list -> Cs.t list
		= fun bs ->
		List.map get_cstr bs

	let to_string : t -> string
		= fun (c,p) ->
		"Boundary : " ^ (Cs.to_string Var.to_string c)
		^ " - point_other_side : " ^ (Vec.to_string Var.to_string p)

    let canon : t -> t
        = fun (c,v) ->
        (Cs.canon c, v)
end

let (strict_comp : Cs.t -> Cs.t)
	= fun cstr ->
	let cstr_comp = Cs.compl cstr in
	(match Cs.get_typ cstr_comp with
		| Cstr_type.Le -> {cstr_comp with Cs.typ = Cstr_type.Lt}
		| _ -> cstr_comp)

let (cstr_eq : Cs.t -> Cs.t)
	= fun cstr ->
	{cstr with Cs.typ = Cstr_type.Eq}

let is_trivial : Cs.t -> bool
		= fun c ->
		(Cs.Vec.equal (Cs.get_v c) Cs.Vec.nil)

module ExplorationPoint = struct

	type t =
		| Direction of int * Boundary.t (** (id of origin region, boundary) *)
		| Point of Vec.t

	let to_string : t -> string
		= function
		| Direction (id,(c,x)) ->
			Printf.sprintf "(%s, %s, %s)"
			(string_of_int id)
			(Cs.to_string Var.to_string c)
			(Vec.to_string Var.to_string x)
		| Point v ->
			Vec.to_string Var.to_string v

	let equal : t -> t -> bool
		= fun p1 p2 ->
		match (p1,p2) with
		| Direction(id1, (cstr1, pointToExplore1)), Direction (id2, (cstr2, pointToExplore2)) ->
				id1 = id2
				&& Cs.equalSyn cstr1 cstr2
				&& Vec.equal pointToExplore1 pointToExplore2
		| Point v1, Point v2 -> Vec.equal v1 v2
		| _,_ -> false
end

module Region = struct

	type 'c t = {
		id : int;
		r : (Boundary.t * int option) list;
		point : Vec.t; (* Un point dans la région *)
		sx : 'c PSplx.t (* Tableau de simplexe dont l'objectif a donné cette région *)
	}

	let mk : int -> Boundary.t list -> Vec.t -> 'c PSplx.t -> 'c t
		= fun id bounds point sx ->
		{id = id;
		 r = List.map (fun b -> (b,None)) bounds;
		 point = point;
		 sx = sx
		}

    let canon : 'c t -> 'c t
        = fun reg ->
        {reg with
            r = List.map (fun (bnd, i) -> (Boundary.canon bnd, i)) reg.r;
        }

	let copy : 'c t -> 'c t
		= fun reg -> { reg with
			sx = PSplx.copy reg.sx;
		}

	let get_cstrs : 'c t -> Cs.t list
		= fun reg ->
		List.map (fun (b,_) -> Boundary.get_cstr b) reg.r

	let to_string : 'c t -> string
		= fun r ->
		let cstrs = get_cstrs r in
		let vars = Cs.getVars cstrs |> Var.Set.elements in
		(Printf.sprintf "Region %i : \n\tBoundaries : %s\n\t Point : %s\n\t\n\n\tPlot : %s"
			r.id
			(Misc.list_to_string
				(fun (b,i) -> Printf.sprintf "%s -> %s"
					(Boundary.to_string b)
					(match i with None -> "unknown" | Some i -> string_of_int i))
				r.r "\n\t\t")
			(Vec.to_string Var.to_string r.point)
			("P = Polyhedron(ieqs = " ^ (Misc.list_to_string
				(fun c -> Cs.plot vars c)
				cstrs ",") ^ ")"))
		^
		(if Debug.is_enabled DebugTypes.Detail
		then Printf.sprintf "\n\tTableau : \n%s" (PSplx.to_string r.sx)
		else "")

	let extract : 'c PSplx.t -> Cs.t list
	  	= fun sx ->
        sx.PSplx.obj.Objective.lin
	  	|> Misc.rem_dupl (fun c c' -> Cs.equal c c')
	  	|> List.filter (fun c -> is_trivial c |> not)

	(** Returns a point in the regions' interior. The chosen point is a 1 from each boundary if possible. *)
	let getPointInside : region_t -> Var.t -> Cs.t list -> Vec.t option
		= let getPointInside_cone : Var.t -> Cs.t list -> Vec.t option
			= fun horizon cstrs ->
			match Splx.getPointInside_cone horizon cstrs with
			| None -> None
			| Some pl -> Some (Rtree.map Vec.ofSymbolic pl)
		and getPointInside_not_cone : Var.t -> Cs.t list -> Vec.t option
			= fun horizon cstrs ->
			let cstrs' = List.mapi
				(fun i cstr -> i, cstr)
				cstrs
			in
			match Opt.getAsg horizon cstrs' with
			| None -> None
			| Some pl -> Some (Rtree.map Vec.ofSymbolic pl)
		in
		function
		| Cone -> getPointInside_cone
		| NCone -> getPointInside_not_cone

	(** Returns true if the given point lies in the given region. *)
	let (contains : 'c t -> Vec.t -> bool)
		= fun reg p ->
		List.for_all
			(fun ((c,_),_) -> Cs.satisfy (Vec.toRat p) c)
			reg.r

	(** Returns true if the given point lies in the given region. *)
	let (contains' : Cs.t list -> Vec.t -> bool)
		= fun cstrs p ->
		List.for_all (Cs.satisfy (Vec.toRat p)) cstrs

	(** Returns true if the given point lies in the given relaxed ({i i.e.} non-strict) region. *)
	let (contains_large : 'c t -> Vec.t -> bool)
		= fun reg p ->
		List.for_all
			(fun ((c,_),_) -> Cs.satisfy (Vec.toRat p) {c with Cs.typ = Cstr_type.Le})
			reg.r

	(** Returns true if the given point lies in the given relaxed ({i i.e.} non-strict) region. *)
	let (contains'_large : Cs.t list -> Vec.t -> bool)
		= fun cstrs p ->
		List.for_all
			(fun c -> Cs.satisfy (Vec.toRat p) {c with Cs.typ = Cstr_type.Le})
			cstrs

	let contains_ep : 'c t -> ExplorationPoint.t -> bool
		= fun reg ->
		function
		| ExplorationPoint.Direction (_,(_,point)) ->
			contains reg point
		| ExplorationPoint.Point point ->
			contains reg point

	let contains_large_ep : 'c t -> ExplorationPoint.t -> bool
		= fun reg ->
		function
		| ExplorationPoint.Direction (_,(_,point)) ->
			contains_large reg point
		| ExplorationPoint.Point point ->
			contains_large reg point

	let (equalSyn : 'c t -> 'c t -> bool)
		= fun reg1 reg2 ->
		List.for_all2 Cs.equalSyn (get_cstrs reg1) (get_cstrs reg2)

	let (equal : 'c t -> 'c t -> bool)
		= fun reg1 reg2 ->
		Misc.list_eq (fun c cstrs -> List.exists (fun c' -> Cs.equal c c') cstrs)
			(get_cstrs reg1) (get_cstrs reg2)

	let eval_obj : 'c t -> Cs.Vec.t -> Cs.Vec.Coeff.t option
		= fun reg p ->
		let cstr = PSplx.objValue reg.sx in
		Some (Cs.eval cstr p)
end

type 'c mapRegs_t = 'c Region.t MapV.t

(*let next_id : int ref = ref 0*)

type 'c t = {
	regs : 'c mapRegs_t;
	todo: ExplorationPoint.t list
}

let empty : 'c t
	= {regs = MapV.empty ; todo = []}

let to_string : 'c t -> string
	= fun plp ->
	let regs = MapV.bindings plp.regs
		|> List.split
		|> Stdlib.snd
	in
	Printf.sprintf "Regions = %s\nTodo list : \n%s"
		(Misc.list_to_string Region.to_string regs "\n")
		(Misc.list_to_string ExplorationPoint.to_string plp.todo " ; ")

(** Module that checks results on the fly.
	It is based on {!module:Debug.Check}, thus these verifications cost almost nothing if checking is disabled. *)
module Check = struct

	let satisfy : Vec.t -> Cs.t -> bool
			= fun v c ->
			let v' = Vec.toRat v in
			Cs.satisfy v' c

	let check_boundary : Boundary.t -> unit
		= fun (cstr,point) ->
		Debug.Check.check
			(lazy (not (satisfy point cstr)))
			(lazy(Printf.sprintf "Point %s satisfies its boundary %s"
				(Vec.to_string Var.to_string point)
				(Cs.to_string Var.to_string cstr)))

	let check_region : 'c Region.t -> unit
		= fun r ->
		begin
			Debug.Check.check (lazy(Region.contains r r.Region.point))
				(lazy(Printf.sprintf "point not in region %s"
				(Region.to_string r)));
			List.iter (fun (b,_) -> check_boundary b) r.Region.r
		end

	let check_result : 'c t -> unit
		= fun plp ->
		MapV.iter (fun _ -> check_region) plp.regs
end

module Adjacency = struct

	(** Returns true if [reg1] and [reg2] are adjacent through [cstr1] and [cstr2].
	@param cstr1 must belong to [reg1]
	@param cstr2 must belong to [reg2] *)
	let adjacent_cstr_lp : 'c Region.t -> 'c Region.t -> Cs.t -> Cs.t -> bool
		= fun reg1 reg2 cstr1 cstr2->

			let cstrs_reg1 = {cstr1 with Cs.typ = Cstr_type.Le}
				:: (Misc.popAll Cs.equalSyn (Region.get_cstrs reg1) cstr1)
			and cstrs_reg2 = {cstr2 with Cs.typ = Cstr_type.Le}
				:: (Misc.popAll Cs.equal (Region.get_cstrs reg2) cstr2)
			in
			let cstrs = cstrs_reg1 @ cstrs_reg2 in
			let horizon = Cs.getVars cstrs
				|> Var.horizon
			in
			let cstrs = List.mapi (fun i c -> (i,c)) cstrs in
			match Splx.mk horizon cstrs |> Splx.checkFromAdd with
			| Splx.IsUnsat _ -> false
			| Splx.IsOk _ -> true

	(** Returns true if [reg1] and [reg2] are adjacent through [cstr1] and [cstr2].
	@param cstr1 must belong to [reg1]
	@param cstr2 must belong to [reg2] *)
	let adjacent_cstr : 'c Region.t -> 'c Region.t -> Cs.t -> Cs.t -> bool
		= fun reg1 reg2 cstr1 _ ->
		let p = Cs.get_saturating_point cstr1 in
		let obj1 = PSplx.objValue reg1.Region.sx
		and obj2 = PSplx.objValue reg2.Region.sx
		in
		Cs.Vec.Coeff.equal (Cs.eval obj1 p) (Cs.eval obj2 p)

	let adjacent' : int -> Cs.t -> ((int * Cs.t) * 'a) list -> 'c mapRegs_t -> int option
		= fun id_init cstr l regMap ->
		let reg = MapV.find id_init regMap in
		try
			let ((id',_),_) = List.find
				(fun ((id', cstr'), _) ->
					id_init <> id'
				  &&
				  	(adjacent_cstr reg (MapV.find id' regMap) cstr cstr'))
				l
			in
			Some id'
		with Not_found -> None

	(** returns [Some i] with [i] the index of the adjacent region if it exists, [None] otherwise. *)
	let adjacent : int -> Cs.t -> ((int * Cs.t) * 'a) list -> 'c mapRegs_t -> int option
		= fun id_init cstr l regMap ->
		Profile.start "adjacent";
		let res = adjacent' id_init cstr l regMap in
		Profile.stop "adjacent";
		res
end


(** Several methods to choose the next point to explore. *)
module Next_Point = struct

	module Min = Min.Make(Vec)(Vec)(Cs)(MinLP.Splx(Cs))

	let conv
		= Min.({vec_CsVec = (fun x -> Vec.toRat x);
			csVec_Vec  = (fun x -> Vec.ofRat x) ;
			csCoeff_Coeff = (fun x -> Vec.Coeff.ofQ x) ;
			csInput_Cs = (fun x -> x) ;
			vecInput_Vec = (fun x -> x);
			vec_VecInput = (fun x -> x)
		})

	module Point_in_Region = struct
		(** Returns the index of the region in which the given point is minimal w.r.t objective functions. *)
		let eval_regions : ExplorationPoint.t -> 'c t -> (int * 'c Region.t) option
			= fun p plp ->
			try
				match p with
				| ExplorationPoint.Point point -> begin
					let point' = conv.Min.vec_CsVec point in
					MapV.fold
						(fun id reg value_opt ->
							match Region.eval_obj reg point', value_opt with
							| Some v, None -> Some ((id,reg),v)
							| Some v, Some (_,v') -> if Cs.Vec.Coeff.lt v' v
								then Some ((id,reg),v)
								else value_opt
							| None, _ -> value_opt)
						plp.regs
						None
					|> function
						| None -> None
						| Some (i,_) -> Some i
				end

				| ExplorationPoint.Direction (id, (_, point)) ->
					let point' = conv.Min.vec_CsVec point in
					MapV.fold
						(fun id' reg value_opt ->
							if id = id'
							then value_opt
							else match Region.eval_obj reg point', value_opt with
								| Some v, None -> Some ((id,reg),v)
								| Some v, Some (_, v') -> if Cs.Vec.Coeff.lt v' v
									then Some ((id,reg),v)
									else value_opt
								| None, _ -> value_opt)
						plp.regs
						None
					|> function
						| None -> None
						| Some (i,_) -> Some i
			with Invalid_argument _ -> None

		let in_region' : ExplorationPoint.t -> 'c t -> (int * 'c Region.t) option
			= fun p plp ->
			match eval_regions p plp with
			| Some (i,reg) ->
				if Region.contains_large_ep reg p
				then Some (i,reg)
				else None
			| None -> None
		(** Returns a region containing the given point, if it exists. *)

		let in_region : ExplorationPoint.t -> 'c t -> (int * 'c Region.t) option
			= fun p plp ->
			Profile.start "in_region";
			let res = in_region' p plp in
			Profile.stop "in_region";
			res
	end

	(** Module that chooses the next point to explore by raytracing. *)
	module RayTracing = struct

		let eval : (int * 'c Region.t) list -> Min.direction_type -> ((int * Cs.t) * Min.direction) list
			= fun regs dir_type ->
			List.fold_left
				(fun res (id,reg) ->
					List.fold_left
					(fun res ((cstr',_),_) ->
						let v = Min.Sort.value dir_type cstr' in
						((id, cstr'), (dir_type, v)) :: res)
					res reg.Region.r)
				[] regs

		let eval : (int * 'c Region.t) list -> Min.direction_type -> ((int * Cs.t) * Min.direction) list
			= fun regs dir_type ->
			Profile.start "raytracing";
			let res = eval regs dir_type in
			Profile.stop "raytracing";
			res

		(*
		(** Removes the initial region [reg] (identified by [id]) and those that are on the same side of [cstr] that [reg]. *)
		let filter_regions : int -> Cs.t -> mapRegs_t -> (int * 'c Region.t) list
			= fun id cstr regMap ->
			MapV.fold
				(fun id' reg res ->
					if id = id' || (List.exists (fun ((cstr',_),_) -> Cs.equal cstr' cstr) reg.Region.r)
					then res
					else (id',reg) :: res )
				regMap
				[]
		*)
		(** Removes the initial region [reg] (identified by [id]) and those that are on the same side of [cstr] that [reg].
			TODO: use LP.*)
		let filter_regions : int -> Cs.t -> 'c mapRegs_t -> (int * 'c Region.t) list
			= fun id cstr regMap ->
			MapV.fold
				(fun id' reg res ->
					if id = id' || (List.exists (fun ((cstr',_),_) -> Cs.equal cstr' cstr) reg.Region.r)
					then res
					else (id',reg) :: res )
				regMap
				[]

		let debug_evals
			= fun l ->
			Debug.exec l DebugTypes.Detail
				(lazy(Printf.sprintf "Evals before post-processing : \n%s"
						(Misc.list_to_string
							(fun ((id,c),(_,v)) -> Printf.sprintf "id %i, %s -> %s"
								id
								(Cs.to_string Var.to_string c)
								(Cs.Vec.Coeff.to_float v |> string_of_float))
							l "\n")))

		let compute_next_point : Cs.t -> ('a * Min.direction) list list ->  Cs.Vec.Coeff.t option
			= fun cstr evals ->
			let i = Misc.findi (List.exists (fun ((_,cstr'),_) -> Cs.equalSyn cstr cstr')) evals
			in
			let eval_i = List.filter
				(fun ((_,cstr'),_) -> not (Cs.equal cstr cstr')
					&& not (Cs.equal (strict_comp cstr) cstr'))
				(List.nth evals i)
			in
			match eval_i with
			| [] ->
				if List.length evals = (i+1) (* il n'y a pas d'élément après cstr dans la liste*)
				then None
				else let (_,(_,v2)) = List.nth evals (i+1) |> List.hd in
				Some v2
			| (_,(_,v2))  :: _ -> Some v2

		let eval_to_string : ('a * Min.direction) list list -> string
			= fun l ->
			Misc.list_to_string
				(fun l -> Misc.list_to_string
					(fun ((_,c),(_,v)) ->
						(Cs.to_string Var.to_string c) ^ ", " ^ (Cs.Vec.Coeff.to_float v|>string_of_float)) l " ; ")
				l "\n"

		exception Adjacent of int

		let adjust'' : region_t -> (int * Boundary.t) -> 'c mapRegs_t -> (int * 'c Region.t) list -> ExplorationPoint.t
			= fun reg_t (id, (cstr, pointOtherSide)) regMap regions_to_eval ->
			Min.conv := conv;
			let x0 = (MapV.find id regMap).Region.point in
			let dir_type = Min.TwoPoints (x0,pointOtherSide) in
			let v1 = Min.Sort.value dir_type cstr in
			let evals = ((id,cstr), (dir_type, v1))
				:: (eval regions_to_eval dir_type)
				|> Min.Sort.sort
				|> Min.Sort.stack
			in
			Debug.log DebugTypes.Detail (lazy(Printf.sprintf "evals = %s"
				(eval_to_string evals)));
			try
				(* XXX: que se passe t-il si i n'est pas 0? *)
				let i = Misc.findi (List.exists (fun ((_,cstr'),_) -> Cs.equal cstr cstr')) evals in
				let evals =
					if List.length (List.nth evals i) > 1
					then (* On a trouvé des candidat adjacents, il faut vérifier qu'ils sont bien adjacents *)
						match Adjacency.adjacent id cstr (List.nth evals i) regMap with
						| None -> begin match reg_t with
							| Cone -> (Misc.sublist evals 0 i) @ [[(id,cstr), (dir_type, v1)]] @ (Misc.sublist evals (i+1) (List.length evals))
							| NCone -> evals
						end
						| Some id_adj -> Stdlib.raise (Adjacent id_adj) (* un des candidats était bien adjacent *)
					else evals
				in
				match compute_next_point cstr evals with
				| None -> ExplorationPoint.Direction(id, (cstr, pointOtherSide))
				| Some v2 ->
					let v' = Cs.Vec.Coeff.div
						(Cs.Vec.Coeff.add v1 v2)
						(Cs.Vec.Coeff.of_float 2.)
					in
					Debug.log DebugTypes.Detail
						(lazy (Printf.sprintf "computing middle of %s and %s : %s"
						(Cs.Vec.Coeff.to_string v1)
						(Cs.Vec.Coeff.to_string v2)
						(Cs.Vec.Coeff.to_string v')));
					let x = Min.Sort.getPoint cstr (dir_type,v') in
						Debug.log DebugTypes.Detail
						(lazy (Printf.sprintf "The next exploration point will be %s"
							(Vec.to_string Var.to_string x)));
					ExplorationPoint.Direction(id, (cstr, x))
			with Not_found -> Stdlib.failwith "PLP.adjust"

		let adjust' : region_t -> (int * Boundary.t) -> 'c mapRegs_t -> ExplorationPoint.t
			= fun reg_t (id, (cstr, pointOtherSide)) regMap ->
			let regions_to_eval = filter_regions id cstr regMap in
			adjust'' reg_t (id, (cstr, pointOtherSide)) regMap regions_to_eval

		let adjust : region_t -> (int * Boundary.t) -> 'c mapRegs_t -> ExplorationPoint.t
			= fun reg_t (id, (cstr, pointOtherSide)) regMap ->
			Debug.log DebugTypes.Detail
				(lazy (Printf.sprintf "Adjusting exploration point (%s, %s, %s)"
					(Vec.to_string Var.to_string((MapV.find id regMap).Region.point))
					(Cs.to_string Var.to_string cstr)
					(Vec.to_string Var.to_string pointOtherSide)));
			adjust' reg_t (id, (cstr, pointOtherSide)) regMap

		let apply_adjacency : 'c t -> int -> Boundary.t -> int -> 'c t
			= fun plp id_init boundary id_adj ->
			(* Updating the region map*)
			let (cstr,_) = boundary in
			let cstr_adj =  strict_comp cstr in
			let reg = MapV.find id_init plp.regs
			and reg_adj = MapV.find id_adj plp.regs in
			let reg = {reg with
				Region.r = List.map (fun (b,i) ->
					if Boundary.equal b boundary
					then (b,Some id_adj)
					else (b,i))
				 	reg.Region.r
			}
			and reg_adj = {reg_adj with
				Region.r = List.map (fun ((cstr',v),i) ->
					if Cs.equal cstr' cstr_adj(* XXX: equal?*)
					then ((cstr',v), Some id_init)
					else ((cstr',v),i))
				 	reg_adj.Region.r
			} in
			let map' = MapV.add id_init reg plp.regs
				|> MapV.add id_adj reg_adj in
			(* Updating the exploration point list*)
			let todo' = List.fold_left
				(fun res -> function
					| ExplorationPoint.Point _ as x -> x :: res
					| ExplorationPoint.Direction (id,(cstr,_)) as x ->
						if id = id_adj && (Cs.equal cstr cstr_adj)(* XXX: equal?*)
						then res
						else x :: res)
				[] plp.todo in
			{regs = map' ; todo = todo'}


		(** The PLP state [t] is modified if an adjacency is found. *)
		let rec get_next_point_rec : region_t -> 'c t -> 'c t
			= fun reg_t plp ->
			match plp.todo with
			| [] -> plp
			| (ExplorationPoint.Direction (id,b)) :: tl -> begin
				try
					let p = adjust reg_t (id,b) plp.regs in
					{plp with todo = p :: tl}
				with Adjacent id_adjacent ->
					apply_adjacency {plp with todo = tl} id b id_adjacent
					|> get_next_point_rec reg_t
				end
			| (ExplorationPoint.Point vec) :: tl ->
				if MapV.exists (fun _ reg -> Region.contains reg vec) plp.regs
				then (Debug.log DebugTypes.Detail
						(lazy (Printf.sprintf "Exploration point %s lies within an already discovered region"
						(Vec.to_string Var.to_string vec)));
					get_next_point_rec reg_t {plp with todo = tl})
				else {plp with todo = (ExplorationPoint.Point vec) :: tl}

		(** Returns the next point to explore, or [None] if there is none.
		The returned point is the first element of [plp.todo]. *)
		let get_next_point : region_t -> 'c t -> ExplorationPoint.t option * 'c t
			= fun reg_t plp ->
			let plp = get_next_point_rec reg_t plp in
			if plp.todo = []
			then (None, plp)
			else (Some (List.hd plp.todo),plp)
	end

	(** Other algorithm to find new exploration points *)
	module RayTracing_min = struct


		let eval : (int * 'c Region.t) list -> Min.direction_type -> ((int * Cs.t) * Min.direction) list
			= fun regs dir_type ->
			List.fold_left
				(fun res (id,reg) ->
					List.fold_left
					(fun res ((cstr',_),_) ->
						let v = Min.Sort.value dir_type cstr' in
						((id, cstr'), (dir_type, v)) :: res)
					res reg.Region.r)
				[] regs

		let debug_evals
			= fun l ->
			Debug.exec l DebugTypes.Detail
				(lazy(Printf.sprintf "Evals before post-processing : \n%s"
						(Misc.list_to_string
							(fun ((id,c),(_,v)) -> Printf.sprintf "id %i, %s -> %s"
								id
								(Cs.to_string Var.to_string c)
								(Cs.Vec.Coeff.to_float v |> string_of_float))
							l "\n")))

		let compute_next_point : Cs.t -> ('a * Min.direction) list list ->  Cs.Vec.Coeff.t option
			= fun cstr evals ->
			let i = Misc.findi (List.exists (fun ((_,cstr'),_) -> Cs.equalSyn cstr cstr')) evals
			in
			let eval_i = List.filter
				(fun ((_,cstr'),_) -> not (Cs.equal cstr cstr')
					&& not (Cs.equal (strict_comp cstr) cstr'))
				(List.nth evals i)
			in
			match eval_i with
			| [] ->
				if List.length evals = (i+1) (* il n'y a pas d'élément après cstr dans la liste*)
				then None
				else let (_,(_,v2)) = List.nth evals (i+1) |> List.hd in
				Some v2
			| (_,(_,v2))  :: _ -> Some v2

		let eval_to_string : ('a * Min.direction) list list -> string
			= fun l ->
			Misc.list_to_string
				(fun l -> Misc.list_to_string
					(fun ((_,c),(_,v)) ->
						(Cs.to_string Var.to_string c) ^ ", " ^ (Cs.Vec.Coeff.to_float v|>string_of_float)) l " ; ")
				l "\n"

		exception Adjacent of int

		let adjust : region_t -> 'c mapRegs_t -> (int * Boundary.t) -> (int * 'c Region.t) -> ExplorationPoint.t
			= fun reg_t regMap (id, (cstr, pointOtherSide)) reg ->
			Min.conv := conv;
			let x0 = (MapV.find id regMap).Region.point in
			let dir_type = Min.TwoPoints (x0,pointOtherSide) in
			let v1 = Min.Sort.value dir_type cstr in
			let evals = ((id,cstr), (dir_type, v1))
				:: (eval [reg] dir_type)
				|> debug_evals
				|> Min.Sort.sort
				|> Min.Sort.stack
			in
			Debug.log DebugTypes.Detail (lazy(Printf.sprintf "evals = %s"
				(eval_to_string evals)));
			try
				(* XXX: que se passe t-il si i n'est pas 0? *)
				let i = Misc.findi (List.exists (fun ((_,cstr'),_) -> Cs.equal cstr cstr')) evals in
				let evals =
					if List.length (List.nth evals i) > 1
					then (* On a trouvé des candidat adjacents, il faut vérifier qu'ils sont bien adjacents *)
						match Adjacency.adjacent id cstr (List.nth evals i) regMap with
						| None -> begin match reg_t with
							| Cone -> (Misc.sublist evals 0 i) @ [[(id,cstr), (dir_type, v1)]] @ (Misc.sublist evals (i+1) (List.length evals))
							| NCone -> evals
						end
						| Some id_adj -> Stdlib.raise (Adjacent id_adj) (* un des candidats était bien adjacent *)
					else evals
				in
				match compute_next_point cstr evals with
				| None -> ExplorationPoint.Direction(id, (cstr, pointOtherSide))
				| Some v2 ->
					let v' = Cs.Vec.Coeff.div
						(Cs.Vec.Coeff.add v1 v2)
						(Cs.Vec.Coeff.of_float 2.)
					in
					Debug.log DebugTypes.Detail
						(lazy (Printf.sprintf "computing middle of %s and %s : %s"
						(Cs.Vec.Coeff.to_string v1)
						(Cs.Vec.Coeff.to_string v2)
						(Cs.Vec.Coeff.to_string v')));
					let x = Min.Sort.getPoint cstr (dir_type,v') in
						Debug.log DebugTypes.Detail
						(lazy (Printf.sprintf "The next exploration point will be %s"
							(Vec.to_string Var.to_string x)));
					ExplorationPoint.Direction(id, (cstr, x))
			with Not_found -> Stdlib.failwith "PLP.adjust"

		let apply_adjacency : 'c t -> int -> Boundary.t -> int -> 'c t
			= fun plp id_init boundary id_adj ->
			(* Updating the region map*)
			let (cstr,_) = boundary in
			let cstr_adj =  strict_comp cstr in
			let reg = MapV.find id_init plp.regs
			and reg_adj = MapV.find id_adj plp.regs in
			let reg = {reg with
				Region.r = List.map (fun (b,i) ->
					if Boundary.equal b boundary
					then (b,Some id_adj)
					else (b,i))
				 	reg.Region.r
			}
			and reg_adj = {reg_adj with
				Region.r = List.map (fun ((cstr',v),i) ->
					if Cs.equal cstr' cstr_adj(* XXX: equal?*)
					then ((cstr',v), Some id_init)
					else ((cstr',v),i))
				 	reg_adj.Region.r
			} in
			let map' = MapV.add id_init reg plp.regs
				|> MapV.add id_adj reg_adj in
			(* Updating the exploration point list*)
			let todo' = List.fold_left
				(fun res -> function
					| ExplorationPoint.Point _ as x -> x :: res
					| ExplorationPoint.Direction (id,(cstr,_)) as x ->
						if id = id_adj && (Cs.equal cstr cstr_adj)(* XXX: equal?*)
						then res
						else x :: res)
				[] plp.todo in
			{regs = map' ; todo = todo'}

		(** The PLP state [t] is modified if an adjacency is found. *)
		let rec get_next_point_rec : region_t -> 'c t -> 'c t
			= fun reg_t plp ->
			match plp.todo with
			| [] -> plp
			| ((ExplorationPoint.Direction (id,b)) as ep):: tl -> begin
				match Point_in_Region.in_region ep plp with
				| None -> {plp with todo = ep :: tl}
				| Some (i,reg) ->
				try
					let p = adjust reg_t plp.regs (id,b) (i,reg) in
					 get_next_point_rec reg_t {plp with todo = p :: tl}
				with Adjacent id_adjacent ->
					apply_adjacency {plp with todo = tl} id b id_adjacent
					|> get_next_point_rec reg_t
				end
			| ((ExplorationPoint.Point vec) as ep) :: tl ->
				match Point_in_Region.in_region ep plp with
				| None -> {plp with todo = ep :: tl}
				| Some _ -> begin
					(Debug.log DebugTypes.Detail
					(lazy (Printf.sprintf "Exploration point %s lies within an already discovered region"
					(Vec.to_string Var.to_string vec)));
					get_next_point_rec reg_t {plp with todo = tl})
				end

		(** Returns the next point to explore, or [None] if there is none.
			The returned point is the first element of [plp.todo]. *)
		let get_next_point : region_t -> 'c t -> ExplorationPoint.t option * 'c t
			= fun reg_t plp ->
			let plp = get_next_point_rec reg_t plp in
			if plp.todo = []
			then (None, plp)
			else (Some (List.hd plp.todo),plp)
	end

	module Greedy = struct

		let get_next : 'c t -> int * Boundary.t -> ExplorationPoint.t option
			= fun plp (reg_id, (cstr,point)) ->
			let ipoint = (MapV.find reg_id plp.regs).Region.point in
			let dir = Min.TwoPoints (ipoint, point) in
			let t = Min.Sort.value dir cstr in
			let vec = Min.Sort.getPoint' (dir, Cs.Vec.Coeff.add t (Cs.Vec.Coeff.of_float 0.1)) in
			match Point_in_Region.in_region (ExplorationPoint.Point vec) plp with
			| None -> Some (ExplorationPoint.Point vec)
			| _ ->  None


		let rec get_next_point_rec : region_t -> 'c t -> 'c t
			= fun reg_t plp ->
			match plp.todo with
			| [] -> plp
			| (ExplorationPoint.Direction (id,b)) :: tl -> begin
				match get_next plp (id,b) with
				| Some ep -> {plp with todo = ep :: tl}
				| None -> get_next_point_rec reg_t {plp with todo = tl}
				end
			| (ExplorationPoint.Point vec) :: tl ->
				if MapV.exists (fun _ reg -> Region.contains reg vec) plp.regs
				then (Debug.log DebugTypes.Detail
						(lazy (Printf.sprintf "Exploration point %s lies within an already discovered region"
						(Vec.to_string Var.to_string vec)));
					get_next_point_rec reg_t {plp with todo = tl})
				else {plp with todo = (ExplorationPoint.Point vec) :: tl}

		(** Returns the next point to explore, or [None] if there is none.
			The returned point is the first element of [plp.todo]. *)
		let get_next_point : region_t -> 'c t -> ExplorationPoint.t option * 'c t
			= fun reg_t plp ->
			Min.conv := conv;
			let plp = get_next_point_rec reg_t plp in
			if plp.todo = []
			then (None, plp)
			else (Some (List.hd plp.todo),plp)
	end
end

let get_next_point
	= fun reg_t plp ->
	Profile.start "get_next_point";
	let res = match !Flags.plp with
		| Flags.Greedy -> Next_Point.Greedy.get_next_point reg_t plp
		| Flags.Adj_Raytracing -> Next_Point.RayTracing.get_next_point reg_t plp
		| Flags.Adj_Raytracing_min -> Next_Point.RayTracing_min.get_next_point reg_t plp
	in
	Profile.stop "get_next_point";
	res

let reg_id : int ref = ref 0

(** Returns a fresh id for a new region. *)
let get_fresh_id : unit -> int
	= fun () ->
	let id = !reg_id in
	reg_id := id + 1;
	id

(** Module that runs one exploration *)
module Exec = struct

	(* XXX: faut-il remettre la version ci-dessus pour les floats? *)

	(** [correct_point cstrs point] returns a point that lies in [cstrs]'s interior.
		If no such point exists, [None] is returned.
		It uses {!val:Region.getPointInside}.
	*)
	let correct_point : region_t -> Cs.t list -> Vec.t option
		= fun reg_t cstrs ->
		Debug.log DebugTypes.Normal (lazy("Correcting region"));
        let horizon = Cs.getVars cstrs |> Var.horizon in
		match Region.getPointInside reg_t horizon cstrs with
		| None -> Debug.exec None DebugTypes.Normal
			(lazy (Printf.sprintf "Region has empty interior:\n%s"
				(Cs.list_to_string cstrs)))
		| Some p when Region.contains' cstrs p -> begin
			Debug.log DebugTypes.Normal (lazy(Printf.sprintf "Changing point to %s"
				(Vec.to_string Var.to_string p)));
			Some p
			end
		| Some p -> Stdlib.failwith
			(Printf.sprintf "Point correction returned %s, which is not in region %s"
				(Vec.to_string Var.to_string p)
				(Cs.list_to_string cstrs))

	(** [get_boundaries cstrs point] returns the list of boundaries by minimizing [cstrs].
		@param point must be a point in the interior of [cstrs]. *)
	let (get_boundaries : region_t -> Cs.t list -> Vec.t -> Boundary.t list)
		= fun reg_t cstrs point ->
		Debug.log DebugTypes.Detail (lazy("Looking for true boundaries in new region"));
		match reg_t with
			| Cone -> Minimization.minimize_cone point cstrs
			| NCone -> Minimization.minimize point cstrs

	let exec' : region_t -> 'c PSplx.t -> Vec.t -> 'c Region.t option
	  	= fun reg_t sx pointToExplore ->
	  	Debug.log DebugTypes.Normal
	  		(lazy("Exec on the point " ^ (Vec.to_string Var.to_string pointToExplore)));
	  	Profile.start "LP" ;
        (try Explore.push false pointToExplore sx with ex -> begin
            Profile.stop "LP" ;
            raise ex
        end);
	 	Profile.stop "LP" ;
	 	Debug.log DebugTypes.Normal
	  		(lazy("Found solution " ^ (PSplx.objValue sx |> Cs.to_string Var.to_string)));
	 	let cstrs = Region.extract sx in
	 	match correct_point reg_t cstrs with
	 	| None -> None
	 	| Some point ->
	 		let bounds = get_boundaries reg_t cstrs point in
	 		let id = get_fresh_id() in
	 		let reg = Region.mk id bounds point sx in
	 		Check.check_region reg;
	 		Some reg

	 (** [exec strat sx map exp] returns the region resulting from the exploration of [exp].
		If this region has empty interior, it returns [None]. *)
	let exec : region_t -> 'c PSplx.t -> Vec.t -> 'c Region.t option
	  	= fun reg_t sx pointToExplore ->
		let sx' = PSplx.copy sx in
	  	Profile.start "exec";
	  	let res = try exec' reg_t sx' pointToExplore with ex -> begin
            Profile.stop "exec" ;
            raise ex
            end
        in
	  	Profile.stop "exec";
	  	res
end

let extract_points : 'c Region.t -> int -> ExplorationPoint.t list
	= fun reg reg_id->
	List.map
		(fun (b,_) -> ExplorationPoint.Direction (reg_id, b))
		reg.Region.r

type 'c config = {
	add_region : 'c Region.t option -> 'c Region.t -> ExplorationPoint.t -> 'c t -> 'c t;
	reg_t : region_t;
	points : ExplorationPoint.t list;
	regions : 'c Region.t list;
	}

(* TODO : tirer parti des égalités! *)
module Add_Region = struct

	(** [should_explore_again reg1 cstr reg2] checks whether the point that has been explored should be explored again.
	@param reg1 is the provenance region.
	@param cstr is constraint that has been crossed in the previous region.
	@param reg2 is the new discovered region.*)
	let should_explore_again : 'c Region.t -> Cs.t -> 'c Region.t -> bool
		= fun reg1 cstr reg2 ->
		let cstr_comp = strict_comp cstr in
		Adjacency.adjacent_cstr reg1 reg2 cstr cstr_comp
		|> not

	(* TODO: better equality test? *)
	(** [region_already_known reg plp] checks whether [reg] is already present in [plp]. *)
	let region_already_known : 'c Region.t -> 'c t -> bool
		= fun reg plp ->
		MapV.exists
			(fun _ reg' -> Region.equal reg reg')
			plp.regs

	(** Return a witness if the region has no interior, [None] otherwise. *)
	let has_interior : 'c Region.t -> bool
		= fun reg ->
		let cstrs = Region.get_cstrs reg in
		let nvar = Cs.getVars cstrs |> Var.horizon in
		let cstrs' = List.mapi (fun i c -> (i,c)) cstrs in
		match Splx.checkFromAdd (Splx.mk nvar cstrs') with
		| Splx.IsOk _ -> true
		| Splx.IsUnsat _ -> false (*match w with
			| [i,_] -> Some (List.assoc i cstrs')
			| _ -> Stdlib.failwith "PLP.has_interior"*)

	let remove_point : ExplorationPoint.t -> ExplorationPoint.t list -> ExplorationPoint.t list
		= fun point points ->
		Misc.pop ExplorationPoint.equal points point

	let do_not_add : 'c Region.t -> ExplorationPoint.t -> 'c t -> 'c t
		= fun _ explorationPoint plp ->
		Debug.log DebugTypes.Normal (lazy "Region not added");
		{plp with todo = (remove_point explorationPoint plp.todo)}

	let add : 'c config -> 'c Region.t option -> 'c Region.t -> ExplorationPoint.t -> 'c t -> 'c t
		= fun config reg_origin new_reg explorationPoint plp ->
		if region_already_known new_reg plp
		then do_not_add new_reg explorationPoint plp
		else if has_interior new_reg
			then config.add_region reg_origin new_reg explorationPoint plp
			else do_not_add new_reg explorationPoint plp
end

(** Leads to a feasible parametric simplex tableau. *)
let init_sx : 'c Factory.t -> 'c PSplx.t -> Vec.t -> bool
	= fun factory sx point ->
	Debug.log DebugTypes.Normal
		(lazy(Printf.sprintf "PLP initialization with point = %s"
			(Vec.to_string Var.to_string point)));
  	if PSplxExec.Init.findFeasibleBasis factory sx point
  	then Debug.exec true DebugTypes.Normal (lazy "Problem Sat")
	else Debug.exec false DebugTypes.Normal (lazy "Problem Unsat")



(** This function is used when an execution (ran on point [p]) led to a region with empty interior.
It returns the middle between the initial region's point and [p]. *)
let find_new_point : Vec.t -> Vec.t -> Vec.t
	= fun init_point explored_point ->
	Vec.middle init_point explored_point

let rec exec : 'c config -> 'c PSplx.t -> 'c t -> 'c t
	= fun config sx_init plp ->
	Debug.log DebugTypes.Normal
		(lazy (Printf.sprintf "Current exploration list : \n%s"
			(Misc.list_to_string ExplorationPoint.to_string plp.todo " ; ")));
	match get_next_point config.reg_t plp with
	| (None, plp') -> plp'
	| (Some explorationPoint, plp') ->
		match explorationPoint with
		| ExplorationPoint.Direction (id, (cstr, pointToExplore)) -> begin
			Debug.log DebugTypes.Detail
				(lazy (Printf.sprintf "Exec on constraint %s and point %s"
					(Cs.to_string Var.to_string cstr)
					(Vec.to_string Var.to_string pointToExplore)));
			let reg = MapV.find id plp'.regs in
			match Exec.exec config.reg_t reg.Region.sx pointToExplore with
			| None ->
				let newPointToExplore = find_new_point reg.Region.point pointToExplore in
				let newExplorationPoint = ExplorationPoint.Direction(id, (cstr, newPointToExplore)) in
				(* on remplace l'ancien ExplorationPoint.t de todo par le nouveau *)
				exec config sx_init {plp' with todo = List.tl plp'.todo @ [newExplorationPoint]}
			| Some reg' -> exec config sx_init (Add_Region.add config (Some reg) reg' explorationPoint plp')
			end
		| ExplorationPoint.Point pointToExplore ->
			match Exec.exec config.reg_t sx_init pointToExplore with
			| None -> exec config sx_init {plp' with todo = (Add_Region.remove_point explorationPoint plp'.todo)}
			| Some reg' -> exec config sx_init (Add_Region.add config None reg' explorationPoint plp')

let result_to_string : ('c Region.t * 'c Cons.t) list option -> string
	= fun res ->
	match res with
	| None -> "None"
	| Some res ->
		Printf.sprintf "Regions = %s\n"
			(Misc.list_to_string
				(fun (reg, (c,_)) ->
					(Region.to_string reg)
					^ "\nSolution = "
					^ (Cs.to_string Var.to_string c))
				res "\n")

(** Careful: if a region has no simplex tableau (because it was given in input), it does not appear in the result. *)
let get_results : 'c Factory.t -> 'c t -> ('c Region.t * 'c Cons.t) list
	= fun factory plp ->
	let results = List.map (fun (_,reg) ->
        let reg' = Region.canon reg in
        (reg', PSplx.objValueCert factory reg'.Region.sx)
    ) (MapV.bindings plp.regs)
	in
	Debug.exec results DebugTypes.MOutput
		(lazy(result_to_string (Some results)))

(* also initializes the fresh region id *)
let init_regions : 'c config -> 'c PSplx.t -> 'c mapRegs_t
	= fun config init_sx ->
	reg_id := 0;
	List.fold_left
		(fun map reg ->
			if !reg_id <= reg.Region.id
			then reg_id := reg.Region.id + 1;
			MapV.add reg.Region.id {reg with Region.sx = init_sx} map)
		MapV.empty config.regions

let choose_init_point : 'c PSplx.t -> 'c config -> Vec.t * 'c config
	= fun sx config ->
	if config.points = []
	then
		let point = PSplx.getParams sx
			|> Var.Set.elements
			|> List.map (fun p -> (Vec.Coeff.of_int (Random.int 100),p))
			|> Vec.mk
		in
	  	(point, {config with points = [ExplorationPoint.Point point]})
	else
		let v = match List.hd config.points with
		| ExplorationPoint.Point v -> v
		| ExplorationPoint.Direction(_,(_,v)) -> v
		in
		(v,config)

let init_and_exec : 'c Factory.t -> 'c config -> 'c PSplx.t -> 'c t option
	= fun factory config sx ->
	let (point,config) = choose_init_point sx config in
	if init_sx factory sx point
	then let plp = {
		regs = init_regions config sx;
		todo = config.points
		} in
		Some (exec config sx plp)
	else None

let run : 'c Factory.t -> 'c config -> 'c PSplx.t -> ('c Region.t * 'c Cons.t) list option
	= fun factory config sx ->
	Random.init 0;
	Debug.log DebugTypes.Title
		(lazy("PLP solver with scalar_type = " ^ (Vec.Coeff.name)));
	Debug.log DebugTypes.MInput
		(lazy (Printf.sprintf "Exploration points provided : %s\n%s"
			(Misc.list_to_string ExplorationPoint.to_string config.points " ; ")
			(PSplx.to_string sx)));
	Profile.start "run";
	let res = match init_and_exec factory config sx with
	| None -> None
	| Some plp -> Some (get_results factory plp)
	in
	Profile.stop "run";
	res
