model "HRP"
	uses "mmxprs"; !gain access to the Xpress-Optimizer solver
	
	parameters
		INFILE= "hrp-hard.txt"
	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 procedure CreateVariablesCOLGEN
	forward procedure ColumnGeneration
	forward function IsInTaskSets(MyTasks: set of integer): boolean
	forward function knapsack(C:array(TASKS) of real,
	                   A:array(TASKS) of real, 
	                   B:real, 
	                   ybest:array(TASKS) of integer,
	                   pass: integer):real
	
	forward procedure print_set(dj:real, vx: array(TASKS) of integer)
	
	declarations
		x: dynamic array(TASKSETS) of mpvar
		y: array(TASKS) of mpvar            ! Knapsack variables
		
		TaskSet: array(TASKSETS) of set of integer
		
		AssignTask: array(TASKS) of linctr
		Obj: linctr
	end-declarations

	CreateVariablesCOLGEN
	
	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)	
	
	ColumnGeneration
			
	minimise(Obj)

	writeln("Best integer solution: ", getobjval, " workers")
	writeln("  Task sets used: ")
	forall(i in TASKSETS | getsol(x(i)) > 0) do
		writeln(TaskSet(i), "\t\tDuration:", sum(j in TaskSet(i)) TaskDuration(j))
	end-do 

!*************functions and procedures

	procedure CreateVariablesCOLGEN
		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
	end-procedure
	
	
	procedure ColumnGeneration
	  declarations
	   dualdem: array(TASKS) of real 
	   ybest: array(TASKS) of integer
	   ubnd, zbest, objval: real
	   bas: basis
	  end-declarations
	
	  setparam("XPRS_CUTSTRATEGY", 0)      ! Disable automatic cuts 
	  setparam("XPRS_PRESOLVE", 0)         ! Switch presolve off
	  npatt:=getsize(TASKSETS)
	  npass:=1
	
	  while(true) do
	    minimize(XPRS_LIN, Obj)       ! Solve the LP 
	
	    savebasis(bas)                     ! Save the current basis 
	    objval:= getobjval                 ! Get the objective value      
	
	                                       ! Get the solution values 
	    forall(j in 1..npatt)  soluse(j):=getsol(x(j))
	    forall(i in TASKS)  dualdem(i):=getdual(AssignTask(i))
	                                       ! Solve a knapsack problem
	    zbest:= 1 - knapsack(dualdem, TaskDuration, MaxHours, ybest, npass)
		
		writeln(" Objective value: ", objval)
	    write("Pass ", npass, ": ")   
	    if zbest = 0 then
	      writeln("no profitable column found.\n")
	      break
	    else 
	      print_set(zbest, ybest)          ! Print the new pattern
	      npatt+=1
	      create(x(npatt))               ! Create a new var. for this pattern
	      x(npatt) is_integer
	
	      Obj+= x(npatt)            ! Add new var. to the objective
	      
	      forall(i in TASKS | ybest(i) > 0) do
	      	TaskSet(npatt)+={i} 
	      end-do
	      
	      forall(i in TASKS | i in TaskSet(npatt)) do
	         AssignTask(i)+= x(npatt)  ! Add new var. to demand constr.s
		  end-do
		  
	      loadprob(Obj)               ! Reload the problem 
	      loadbasis(bas)                   ! Load the saved basis 
	    end-if
	    npass+=1
	   end-do
	end-procedure
	
	 function knapsack(C:array(TASKS) of real,
	                   A:array(TASKS) of real, 
	                   B:real, 
	                   ybest:array(TASKS) of integer,
	                   pass: integer):real
	
	! Hide the demand constraints
	  forall(j in TASKS) sethidden(AssignTask(j), true)
	
	! Define the knapsack problem
	  KnapCtr := sum(j in TASKS) A(j)*y(j) <= B
	  KnapObj := sum(j in TASKS) C(j)*y(j)
	
	! Integrality condition
	  if pass=1 then
	   forall(j in TASKS) y(j) is_binary
	  end-if
	
	  maximize(KnapObj)
	  returned:=getobjval
	  forall(j in TASKS) ybest(j):=round(getsol(y(j)))
	
	! Reset knapsack constraint and objective, unhide demand constraints 
	  KnapCtr := 0
	  KnapObj := 0  
	  forall(j in TASKS) sethidden(AssignTask(j), false)
	
	 end-function
	
	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

	 procedure print_set(dj:real, vx: array(TASKS) of integer)
	  declarations 
	   dw: real
	  end-declarations
	  
	  writeln("New task set found with marginal cost ", dj)
	  write("   Task set: ")
	  dw:=0
	  forall(i in TASKS | vx(i) > 0) do 
	    write(i, " ")
	    dw += TaskDuration(i)*vx(i)
	  end-do 
	  writeln("   Total hours: ", dw)
	 end-procedure
	
end-model