"""
A Hybrid Tabu Search for Lot Sizing Problems
Copyright (C) Joao Pedro Pedroso, LIACC/UP, 2005

these routines are dependent on:
  * the model specified in "fri-jt.mod"
      - variables named x[*,*,*] and y[*,*,*] are read from glpk's set of variables
      
  * data/instance format specied in "fridata.py"
      - requires .Prod, .Compat[p], .Period
      - uses the data class to produce the (ampl) model's ".dat" in a temporary file


"""



from glpkpi import *
from glpk import *
import random

LOG = 0
EPS = 1.e-9

lpx_no_print()

import resource
def clock():
    """clock() -> floating point number

    Return the CPU time in seconds (user + system) since the start of the
    process.  This is done via a call to resource.getrusage, so it avoids
    the wraparound problems in time.clock()."""

    res = resource.getrusage(resource.RUSAGE_SELF)
    return res[0]+res[1]


class LotSizing(lpx):
    def __init__(self, mod_file, data, CUT = {}):
        self.data = data
        # "produce ampl data file from data information"
        tmpdat = 'test.dat'	# make a temporary ampl data file with this name
        self.data.mkampl(tmpdat, CUT)
        lpx.__init__(self, mod_file, tmpdat)


    def getsol(self):
        """copy lot sizing model's x, y, obj solution into dictionaries
        x[p,b,t] - (continuous) production
        y[p,b,t] - (0-1) setup
        z - objective
        *** NOTE ***
        supposing a (heuristic or exact) solution in available in self.sol
        """
        x = {}
        y = {}
        for p in self.data.Prod:
            for b in self.data.Compat[p]:
                for t in self.data.Period:
                    xname = "x[%s,%s,%d]" % (p,b,t)
                    yname = "y[%s,%s,%d]" % (p,b,t)
                    x[p,b,t] = self.sol[xname]
                    y[p,b,t] = self.sol[yname]
        z = self.x[0]
        return x,y,z
    

    def printsol(self):
        """print solution in a easy-to-read format"""
        print "%10s%10s%10s%16s%16s" % ("product", "block", "period", "y", "x")
        for p in self.data.Prod:
            for b in self.data.Compat[p]:
                for t in self.data.Period:
                    yname = "y[%s,%s,%d]" % (p,b,t)
                    xname = "x[%s,%s,%d]" % (p,b,t)
                    if self.sol[yname] != 0:
                        print "%10s%10s%10i%16g%16g" % (p, b, t, self.sol[yname], self.sol[xname])
        print "(remaining variables are 0)"
        print "objective:", self.x[0]


    def writesol(self,filename):
        """write x and y solution as text in 'filename'"""
        f = open(filename,'w')
        keys = self.sol.keys()
        keys.sort()
        for i in keys:
            print >>f, "%s = %r" % (i, self.sol[i])
        f.close()

        
    def printstat(self):
        "print the current status of all the integer variables"
        print "bnds,\tkind,\t[lb,val,ub]"
        for i in self.ivars:
            print lpx_get_col_name(self.lp,i),\
                  lpx_get_col_type(self.lp,i), "\t", lpx_get_col_kind(self.lp,i),\
                  "[", lpx_get_col_lb(self.lp,i), lpx_mip_col_val(self.lp,i),\
                  lpx_get_col_ub(self.lp,i),"]"


    def update_sol(self, y):
        """update the current lp to reflect the choices in the y variable
           y is a dictionary: y[
        """
        for p in self.data.Prod:
            for b in self.data.Compat[p]:
                for t in self.data.Period:
                    yname = "y[%s,%s,%d]" % (p,b,t)
                    i = self.indices[yname]
                    value = y[p,b,t]
                    lpx_set_col_bnds(self.lp, i, LPX_FX, value, 0.)
        

    def destruct(self, KIND):
        """partial destruct a solution
        KIND is the glpk kind of variable to set (LPX_CV or LPX_IV)
        """

        unfix = self.ivars[:]
        random.shuffle(unfix)
        alpha = random.random()
        nfix = int(len(self.ivars)*alpha)
        del unfix[nfix:]
        if LOG:
            print "alpha:", alpha
            print "unfix:", unfix
        
        for ivar in range(len(unfix)):
            i = unfix[ivar]
            lpx_set_col_bnds(self.lp, i, LPX_DB, self.lb[ivar], self.ub[ivar])
            lpx_set_col_kind(self.lp, i, KIND)
        if LOG:
            print "made", unfix, " variables", KIND, "..."
        return unfix


    def rnd_fix(self):
        """fix a y variable at a random value; return its index on the lp"""
        ### this is rather inefficient, but later might be required... !!!
        # t = random.choice(self.data.Period)
        # p = random.choice(self.data.Prod)
        # b = random.choice(self.data.Compat[p])
        # i = self.indices["y[%s,%s,%d]" % (p,b,t)]
        # if LOG:
        #     print "changing setup for y[%s,%s,%d]" % (p,b,t)
        i = random.choice(self.ivars)
        ival = 1 - self.x[i]
        lpx_set_col_bnds(self.lp, i, LPX_FX, ival, 0.)
        if LOG:
            print "changed setup for %s to %d" % (lpx_get_col_name(self.lp,i), ival)
        return i



    def limited_solve(self, tmlim):
        """solve with branch & bound, with a timeout limit
        (note: might return before finding a feasible solution!)
        """

        lpx_set_real_parm(self.lp, LPX_K_TMLIM, tmlim)
        lpx.solve(self)
        stat = lpx_mip_status(self.lp)

        if stat == LPX_I_OPT or stat == LPX_I_FEAS:
            # fix variables at the solution found
            for ivar in range(len(self.ivars)):
                i = self.ivars[ivar]
                value = self.x[i]
                lpx_set_col_bnds(self.lp, i, LPX_FX, value, 0.)
        else:
            print "failed, B&B could not find integer solution"
            raise Exception

        lpx_set_real_parm(self.lp, LPX_K_TMLIM, -1)

        
    def relax_and_fix(self):
        """relax-and-fix heuristic:
        make the variables:
          * for the current period: integer
          * for the previous periods: fixed
          * for the subsequent periods: continuous
        solve for all periods t=1...T
        """
        for ivar in range(len(self.ivars)):
            i = self.ivars[ivar]
            lpx_set_col_bnds(self.lp, i, LPX_DB, self.lb[ivar], self.ub[ivar])
            lpx_set_col_kind(self.lp, i, LPX_CV)
        if LOG:
            print "made all integer variables continuous..."
        # order vars
        for t in self.data.Period:
            var = []
            per = ",%d]" % t
            for i in self.ivars:
                nm = self.names[i]
                if nm.find(per) != -1:
                    var.append(i)
            for i in var:
                lpx_set_col_kind(self.lp, i, LPX_IV)
            if LOG:
                print 'made', var, 'integer'
            lpx.solve(self)
            for i in var:
                value = round(self.x[i])
                if LOG:
                    print "fixing %s (%f) to %f" % (self.names[i], self.x[i], value)
                lpx_set_col_bnds(self.lp, i, LPX_FX, value, 0.)
                ### lpx_set_col_kind(self.lp, i, LPX_CV) ???
        

    def relax_and_fix_one_prod(self):
        """a more restricted relax-and-fix:
        make integer only the current periods'variables
        that concern ONE product 
        (this seems to get many times the same sol...)
        """
        for ivar in range(len(self.ivars)):
            i = self.ivars[ivar]
            lpx_set_col_bnds(self.lp, i, LPX_DB, self.lb[ivar], self.ub[ivar])
            lpx_set_col_kind(self.lp, i, LPX_CV)
        if LOG:
            print "made all integer variables continuous..."

        # order vars
        for t in self.data.Period:
            rndprods = self.data.Prod[:]	# make a product list in random order
            random.shuffle(rndprods)
            for p in rndprods:
                if LOG:
                    print "period", t, "product", p,
                var = []
                per = ",%d]" % t
                prd = "[%s," % p
                for i in self.ivars:
                    nm = self.names[i]
                    if nm.find(per) != -1 and nm.find(p) != -1:
                        var.append(i)
                if LOG:
                    print '; made', var, 'integer'
                for i in var:
                    lpx_set_col_kind(self.lp, i, LPX_IV)
                    if LOG:
                        print self.names[i]
                lpx.solve(self)
                for i in var:
                    value = round(self.x[i])
                    if LOG:
                        print "fixing %s (%f) to %f" % (self.names[i], self.x[i], value)
                    lpx_set_col_bnds(self.lp, i, LPX_FX, value, 0.)
                    ### lpx_set_col_kind(self.lp, i, LPX_CV)	### // ??? required ??? 


    def reconstruct_one_prod(self, unfix):
        """reconstruct, starting from a given solutiona more restricted relax-and-fix:
        make integer only the current periods'variables
        that concern ONE product 
        (this seems to get many times the same sol...)
        """

        # order vars
        for t in self.data.Period:
            rndprods = self.data.Prod[:]	# make a product list in random order
            random.shuffle(rndprods)
            for p in rndprods:
                if LOG:
                    print "period", t, "product", p,
                var = []
                per = ",%d]" % t
                prd = "[%s," % p
                for i in unfix:
                    nm = self.names[i]
                    if nm.find(per) != -1 and nm.find(p) != -1:
                        var.append(i)
                if var == []:
                    continue
                if LOG:
                    print '; made', var, 'integer'
                for i in var:
                    lpx_set_col_kind(self.lp, i, LPX_IV)
                    if LOG:
                        print self.names[i]
                lpx.solve(self)
                for i in var:
                    value = round(self.x[i])
                    if LOG:
                        print "fixing %s (%f) to %f" % (self.names[i], self.x[i], value)
                    lpx_set_col_bnds(self.lp, i, LPX_FX, value, 0.)
                    ### lpx_set_col_kind(self.lp, i, LPX_CV)	### // ??? required ??? 


    def relax_and_fix_one_var(self):
        """an even more restricted relax-and-fix:
        for each period, select ONLY ONE variable at a time for fixing;
        variable selection is done in a random order
        """
        for ivar in range(len(self.ivars)):
            i = self.ivars[ivar]
            lpx_set_col_bnds(self.lp, i, LPX_DB, self.lb[ivar], self.ub[ivar])
            lpx_set_col_kind(self.lp, i, LPX_CV)
        if LOG:
            print "made all integer variables continuous..."

        # order vars
        for t in self.data.Period:
            var = []
            per = ",%d]" % t
            for i in self.ivars:
                nm = self.names[i]
                if nm.find(per) != -1:
                    var.append(i)

            random.shuffle(var)

            for i in var:
                lpx_set_col_kind(self.lp, i, LPX_IV)
                lpx.solve(self)
                value = round(self.x[i])
                if LOG:
                    print "fixing %s (%f) to %f" % (self.names[i], self.x[i], value)
                lpx_set_col_bnds(self.lp, i, LPX_FX, value, 0.)
                ### lpx_set_col_kind(self.lp, i, LPX_CV) ???


class Stats:
    def __init__(self):
        self.niter = 0
        self.ndestr = 0
        self.nalltabu = 0
        self.nreconstr = 0

        self.niunset = 0
        self.niset = 0
        self.nidestr = 0

    def __del__(self):
        print "\n\n\nRUN STATISTTICS"
        print "number of iterations:", self.niter
        print "number of destructions:", self.ndestr
        print "number of reconstructions:", self.nreconstr
        print "number of blocked situations:", self.nalltabu

        print "number of un-setup improvements on best:", self.niunset
        print "number of setup improvements on best:", self.niset
        print "number of destruct/reconstruct improvements on best:", self.nidestr

stats = Stats()	# quick and dirty


def tabumove(prob, par, curr, zcurr, zbest, tabu, it):
    """check current neighbors for an improving solution
       arguments:
       prob -> a LotSizing instance (hence holding a lpx structure
       par -> ?????
       curr -> current solution (y[ in the MIP, 0-1 vector)
       zbest -> best found solution's objective value
       tabu -> list of variables which are tabu
       """

    stats.niter +=1
    cand = None

    # get the set of y=1 on the current solution
    # required info:
    #   * the product corresponding to each variable
    #   * the variables (machines) corresponding to each product (i.e., variable)
    #   * 



    if LOG:
        print "TC: starting candidate search; initial solution:"
        prob.printsol()

    move = None
    for t in prob.data.Period:
        if LOG:
            print "TC: PERIOD", t
        rndprods = prob.data.Prod[:]	# make a product list in random order
        random.shuffle(rndprods)
        for p in rndprods:
            if LOG:
                print "TC: PRODUCT", p
            tabulen = int(len(prob.ivars) * random.random())	# !!! tabu length !!!
            # if p is being produced in this period
            setup = []
            nosetup = []
            for b in prob.data.Compat[p]:
                index = prob.indices["y[%s,%s,%d]" % (p,b,t)]
                if curr[p,b,t] == 1:
                    setup.append(index)
                else:
                    nosetup.append(index)

            for i in setup:
                lpx_set_col_bnds(prob.lp, i, LPX_FX, 0., 0.)
                prob.solve()
                x,y,z = prob.getsol()
                if LOG:
                    print "unsetting up", prob.names[i], "obj:", z

                if z < zbest-EPS:	# aspiration criterion
                    stats.niunset += 1
                    if LOG:
                        print "returning (best updated)"
                    return y, z

                if (cand == None or z < zcand-EPS) and \
                       it - tabu[t,p] > tabulen: 	# p is not tabu
                    cand = y.copy()
                    zcand = z
                    move = (i,0,None,None)
                    tmove = t,p	# tabu information
                    if LOG:
                        print "updated candidate: unsetting", i, "(%s)" % prob.names[i], "z=", zcand

                    if zcand < zcurr:
                        if LOG:
                            print "returning (improved current solution)"
                        return cand, zcand

                for j in nosetup:
                    lpx_set_col_bnds(prob.lp, j, LPX_FX, 1., 0.)
                    prob.solve()
                    x,y,z = prob.getsol()
                    if LOG:
                        print "setting up", prob.names[j], "obj:", z

                    if z < zbest-EPS:	# aspiration criterion
                        stats.niset += 1
                        if LOG:
                            print "returning (-- best updated)"
                        return y, z

                    if (cand == None or z < zcand-EPS) and \
                           it - tabu[t,p] > tabulen: 	# p is not tabu
                        cand = y.copy()
                        zcand = z
                        move = (i,0,j,1)
                        tmove = t,p	# tabu information
                        if LOG:
                            print "updated candidate: setting", j, "(%s)" % prob.names[j], "z=", zcand
                        if zcand < zcurr:
                            if LOG:
                                print "returning (improved current solution)"
                            return cand, zcand


                    lpx_set_col_bnds(prob.lp, j, LPX_FX, 0., 0.)
                lpx_set_col_bnds(prob.lp, i, LPX_FX, 1., 0.)
                            
    if move == None:
        "all moves were tabu, make a random move"
        stats.nalltabu += 1
        index = prob.rnd_fix()
        p = random.choice(prob.data.Prod)
        b = random.choice(prob.data.Compat[p])
        t = random.choice(prob.data.Period)
        i = prob.indices["y[%s,%s,%d]" % (p,b,t)]
        print "blocked; changing setup for y[%s,%s,%d]" % (p,b,t)
        ival = 1 - curr[p,b,t]
        tabu[(t,p)] = it
        lpx_set_col_bnds(prob.lp, i, LPX_FX, ival, 0.)
        x,y,z = prob.getsol()
        print "returnning; random move:", prob.names[i], "->", ival
        return y,z

    else:
        stats.ndestr += 1
        i,ival,j,jval = move

        unfix = prob.destruct(LPX_CV)

        lpx_set_col_bnds(prob.lp, i, LPX_FX, ival, 0.)
        try:
            unfix.remove(i)
        except Exception:
            pass
        if j != None:
            lpx_set_col_bnds(prob.lp, j, LPX_FX, jval, 0.)
            try:
                unfix.remove(j)
            except Exception:
                pass

        # to do: remove all the tabu variables from unfix !!!
        prob.reconstruct_one_prod(unfix)
        stats.nreconstr += 1
        if LOG:
            prob.printsol()
            print "candidate objective:", zcand
        x,cand,zcand = prob.getsol()
        if zcand < zbest-EPS:
            stats.nidestr += 1
        if LOG:
            print "candidate objective after reconstruction:", zcand
            print "returning: reconstructed move"

        tabu[tmove] = it
        return cand,zcand

####        i,ival,j,jval = move
####        print "candidate is", prob.names[i], "->", ival
####        lpx_set_col_bnds(prob.lp, i, LPX_FX, ival, 0.)
####        if j != None:
####            print "candidate is", prob.names[j], "->", jval
####            lpx_set_col_bnds(prob.lp, j, LPX_FX, jval, 0.)
####
####        tabu[tmove] = it
####        print "returning: normal tabu move"
####        return cand,zcand
        
    

def tabusearch(prob, sol, zsol, zbest, tmlim):
    "tabu search main procedure"

    # initialize tabu list:
    tabu = {}
    for t in prob.data.Period:
        for p in prob.data.Prod:
            tabu[t,p] = -len(prob.ivars)

    print "initial solution:", zsol, "/", zbest
    best = sol.copy()
    it = 0
    while clock() < tmlim:
        sol, zsol = tabumove(prob, None, sol, zsol, zbest, tabu, it)
        it += 1
        print "iter", it, "current solution:", zsol, "/", zbest, "CPU", clock()
        if zsol < zbest:
            "improved, updated best solution"
            best = sol.copy()
            zbest = zsol

    print "best found solution:", best, zbest
    return best, zbest




if __name__ == "__main__":
    import sys
    args = sys.argv
    args = ["", 1, 9, 600]
    if len(args) != 4:
        print "usage:", args[0], "seed instance-number tmlim(s)"
        sys.exit(-1)
    seed = int(args[1])
    num = int(args[2])
    tmlim = float(args[3])
    print seed, num, tmlim
    
    random.seed(seed)

    import fridata
    # data = fridata.toy
    # data = fridata.jt
    # data = fridata.rnd
    data = fridata.inst[num]
    prob = LotSizing("fri-jt.mod", data)

    # prob.limited_solve(10)
    ### prob.printstat()
    ### raise Exception

    ### print "EXACT SOLUTION:"
    ### prob.solve()
    ### prob.printsol()
    ### 
    ### print "relax_and_fix SOLUTION:"
    ### prob.relax_and_fix()
    ### prob.printsol()
    
    print "relax_and_fix_one_prod SOLUTION:"
    prob.relax_and_fix_one_prod()
    x,y,z = prob.getsol()
    ### prob.printstat()
    prob.printsol()
    print clock()
    ### sys.exit(0)
    ### raise Exception

    print "\n\n\n\n\ntabu search + reconstruct_one_prod:"
    sol,zsol = tabusearch(prob, y, z, z, tmlim)
    print "final solution:", sol
    print "final objective:", zsol
    sys.exit(0)
    
    ### LOG=0
    ### print prob.indices
    ### for i in range(100):
    ###     print "\n\n\n\n\ndestroy and B&B:"
    ###     prob.destruct(LPX_IV)
    ###     prob.rnd_fix()
    ###     prob.limited_solve(20)
    ###     x,y,z = prob.getsol()
    ###     prob.printsol()
    ### raise Exception

    LOG=0
    for i in range(1000):
        unfix = prob.destruct(LPX_CV)
        fixed = prob.rnd_fix()
        try:
            unfix.remove(fixed)
        except Exception:
            pass
        prob.reconstruct_one_prod(unfix)
        ### x,y,z = prob.getsol()
        prob.printsol()
    raise Exception

    print "reconstruct_one_prod SOLUTION:"
    prob.reconstruct_one_prod()
    x,y,z = prob.getsol()
    print "x=", x
    print "y=", y
    prob.printsol()

    print "reconstruct_one_prod SOLUTION:"
    prob.reconstruct_one_prod()
    x,y,z = prob.getsol()
    print "x=", x
    print "y=", y
    prob.printsol()

###    tabusearch(prob, y, 10000, z)