model "HRP"
	uses "mmxprs"; !gain access to the Xpress-Optimizer solver
	
	parameters
		INFILE= "hrp.txt"
		Version = 2
	end-parameters
	
	!sample declarations section
	declarations
		WORKERS: set of integer
		TASKS: set of integer
		TASKSETS: set of integer
		
		MaxHours: real
		TaskDuration: array(TASKS) of real
		MaxNumWorkers: integer
	end-declarations
	
	initialisations from INFILE
		MaxHours
		TaskDuration
	end-initialisations

	finalize(TASKS)
	
	forward function ConstructFeasible: integer
	forward procedure CreateVariables
	forward function IsInTaskSets(MyTasks: set of integer): boolean
		
	MaxNumWorkers:= ConstructFeasible
	
	WORKERS:= 1..MaxNumWorkers
	
	finalize(WORKERS)
	
	declarations
		z: array(TASKS,WORKERS) of mpvar		!fraction of service to client j from site i
		w: array(WORKERS) of mpvar				!1 if facility is constructed at site i; 0 otherwise
		x: dynamic array(TASKSETS) of mpvar
		
		TaskSet: array(TASKSETS) of set of integer
		
		AssignTask: array(TASKS) of linctr
	end-declarations

	if Version = 1 then	
		!construct a feasible solution
		
		forall(i in TASKS) do
			AssignTask(i):=
			sum(j in WORKERS) z(i,j) = 1
		end-do
		
		forall(j in WORKERS) do
			sum(i in TASKS) TaskDuration(i) * z(i,j) <= MaxHours * w(j)
		end-do
		
		forall(j in WORKERS) do
			w(j) is_binary
			
			forall(i in TASKS) do
				z(i,j) is_binary	
			end-do
		end-do
		
(!		forall(j in WORKERS | j >= 2) do
			sum(i in TASKS) z(i,j-1) >= sum(i in TASKS) z(i,j)
		end-do
!)	
		
		Obj:= sum(j in WORKERS) w(j)
	elif Version = 2 then
		CreateVariables
		
		forall(i in TASKS) do
			AssignTask(i):=
			sum(s in TASKSETS | i in TaskSet(s)) x(s) = 1
		end-do
		
		Obj:= sum(s in TASKSETS) x(s)		
	end-if
		
	minimise(Obj)


!procedures and functions
	function ConstructFeasible: integer
		declarations
			NumWorkers: integer
			WORKERS_TMP: set of integer
		end-declarations
		
		WORKERS_TMP:=1..getsize(TASKS)
		
		declarations
			WorkerLoad: array(WORKERS_TMP) of real
		end-declarations
		
		forall(i in TASKS) do
			forall(j in WORKERS_TMP) do
				if WorkerLoad(j) + TaskDuration(i) <= MaxHours then
					WorkerLoad(j)+=TaskDuration(i)
					break
				end-if
			end-do
		end-do
		
		NumWorkers:= integer(sum(j in WORKERS_TMP | WorkerLoad(j) > 0) 1)
		
		returned:= NumWorkers
	end-function
	
	procedure CreateVariables
		declarations
			NumWorkers: integer
			WORKERS_TMP: set of integer
			TaskSetIndex: integer
			INDICES: list of integer
			
			L1: list of integer
			L2: list of integer
			pos: integer
			FakeDuration: array(TASKS) of real
		end-declarations
		
		WORKERS_TMP:=1..getsize(TASKS)
		
		declarations
			WorkerLoad: array(WORKERS_TMP) of real
			WorkerTasks: array(WORKERS_TMP) of set of integer
		end-declarations
		
		!First job to first worker
		forall(i in TASKS) do
			INDICES:= INDICES+[i]
		end-do
		
		forall(i in INDICES) do
			forall(j in WORKERS_TMP) do
				if WorkerLoad(j) + TaskDuration(i) <= MaxHours then
					WorkerLoad(j)+=TaskDuration(i)
					WorkerTasks(j)+={i}
					break
				end-if
			end-do
		end-do
		
		forall(j in WORKERS_TMP | WorkerLoad(j) > 0 and IsInTaskSets(WorkerTasks(j)) = false) do
			TaskSet(getsize(TASKSETS) + 1):= WorkerTasks(j)
			create(x(getsize(TASKSETS)))
			x(getsize(TASKSETS)) is_binary
		end-do
		
		!Sort task in non-decreasing order
		INDICES:=[]
		forall(i in TASKS) do
			!find position
			pos:= 1
			forall(j in INDICES | TaskDuration(j) <= TaskDuration(i)) do
				pos+= 1
			end-do
			
			L1:= splithead(INDICES,pos-1)
			
			INDICES:= L1+[i]+INDICES
		end-do

		forall(j in WORKERS_TMP) do
			WorkerLoad(j):=0
			WorkerTasks(j):={}
		end-do
		
		forall(i in INDICES) do
			forall(j in WORKERS_TMP) do
				if WorkerLoad(j) + TaskDuration(i) <= MaxHours then
					WorkerLoad(j)+=TaskDuration(i)
					WorkerTasks(j)+={i}
					break
				end-if
			end-do
		end-do
		
		forall(j in WORKERS_TMP | WorkerLoad(j) > 0 and IsInTaskSets(WorkerTasks(j)) = false) do
			TaskSet(getsize(TASKSETS) + 1):= WorkerTasks(j)
			create(x(getsize(TASKSETS)))
			x(getsize(TASKSETS)) is_binary
		end-do
		
		!Sort task in non-increasing order
		INDICES:=[]
		forall(i in TASKS) do
			!find position
			pos:= 1
			forall(j in INDICES | TaskDuration(j) >= TaskDuration(i)) do
				pos+= 1
			end-do
			
			L1:= splithead(INDICES,pos-1)
			
			INDICES:= L1+[i]+INDICES
		end-do

		forall(j in WORKERS_TMP) do
			WorkerLoad(j):=0
			WorkerTasks(j):={}
		end-do
		
		forall(i in INDICES) do
			forall(j in WORKERS_TMP) do
				if WorkerLoad(j) + TaskDuration(i) <= MaxHours then
					WorkerLoad(j)+=TaskDuration(i)
					WorkerTasks(j)+={i}
					break
				end-if
			end-do
		end-do
		
		forall(j in WORKERS_TMP | WorkerLoad(j) > 0 and IsInTaskSets(WorkerTasks(j)) = false) do
			TaskSet(getsize(TASKSETS) + 1):= WorkerTasks(j)
			create(x(getsize(TASKSETS)))
			x(getsize(TASKSETS)) is_binary
		end-do
		
		
		!Sort task in random order (n rounds)
		forall(n in 1..100) do
			forall(i in TASKS) do
				FakeDuration(i):= random
			end-do
		
			INDICES:=[]
			forall(i in TASKS) do
				!find position
				pos:= 1
				forall(j in INDICES | FakeDuration(j) >= FakeDuration(i)) do
					pos+= 1
				end-do
				
				L1:= splithead(INDICES,pos-1)
				
				INDICES:= L1+[i]+INDICES
			end-do
	
			forall(j in WORKERS_TMP) do
				WorkerLoad(j):=0
				WorkerTasks(j):={}
			end-do
			
			forall(i in INDICES) do
				forall(j in WORKERS_TMP) do
					if WorkerLoad(j) + TaskDuration(i) <= MaxHours then
						WorkerLoad(j)+=TaskDuration(i)
						WorkerTasks(j)+={i}
						break
					end-if
				end-do
			end-do
			
			forall(j in WORKERS_TMP | WorkerLoad(j) > 0 and IsInTaskSets(WorkerTasks(j)) = false) do
				TaskSet(getsize(TASKSETS) + 1):= WorkerTasks(j)
				create(x(getsize(TASKSETS)))
				x(getsize(TASKSETS)) is_binary
			end-do		
		end-do	
		
	end-procedure
	
	function IsInTaskSets(MyTasks: set of integer): boolean
		declarations
			result: boolean
		end-declarations
		
		result:= false
		forall(i in TASKSETS | MyTasks = TaskSet(i)) do
			result:= true
			break
		end-do
		
		returned:= result
	end-function
	

	
end-model