#! /usr/bin/python3

import sys
from math import sqrt
from pyx import *

verbose = False

def plot_side (c, side):
   xm = side[0]
   ym = side[1]
   xl = side[2]
   x0 = xm - xl/2
   x1 = xm + xl/2
   c.stroke(path.line(x0, ym, x1, ym), [style.linewidth(0.1 * unit.w_cm), color.rgb.blue])

def plot_air (c, p, q):
   c.stroke(path.line(p[0], p[1], q[0], q[1]), [style.linestyle.dashed, color.rgb.red])

def plot_stroke (c, Rk, rbitk):
   pk = Rk['p'][rbitk]
   qk = Rk['p'][1-rbitk]
   c.stroke(path.line(pk[0], pk[1], qk[0], qk[1]), [style.linewidth(0.8 * Rk['w'] * unit.w_cm), style.linecap.round, color.gray(0.8)])

def plot_strokes (R):

   c = canvas.canvas()
   style.linewidth(0.2 * unit.w_cm)
   qend = None

   for k in range(len(R)):
        Rk = R[k]
        rbitk = R[k]['rbit']
        pk = Rk['p'][rbitk]
        qk = Rk['p'][1-rbitk]
        plot_stroke (c, Rk, 0)
        for i in range(len(Rk['links'][0])):
           linki = Rk['links'][0][i]  
           sidei = linki[0]
           plot_side (c, sidei)
        if k > 0: 
           plot_air (c, qend, pk)
        qend = qk

   c.writeEPSfile("strokes.eps")

def create_stroke (sid, p, q, w):
    '''
    This function creates a stroke object with endpoints {p,q} with width {w}
    and no sides and identifier {sid}.
    '''
    R = {'sid': sid, 'p':[p,q], 'w':w, 'e':extrude_time(p,q), 'links':[[].copy(),[].copy()], 'Tini': None, 'Tend': None, 'rbit': 0}
    return R

def add_side (R1, R2):
    '''
    Creates a side object which is an output of {R1} and input of {R2}. The two strokes
    must be adjacent. For now we assume that the strokes are horizontal, which R1
    below R2, and have the same width and starting point to the left of the endpoint.
    '''
    p1 = R1['p'][0]
    q1 = R1['p'][1]
    y1 = p1[1]
    assert y1 == q1[1] #R1 must be horizontal
    assert p1[0] <= q1[0] #R1 is left to right 
    p2 = R2['p'][0]
    q2 = R2['p'][1]
    y2 = p2[1]
    assert y2 == q2[1] #R2 must be horizontal
    assert p2[0] <= q2[0] #R2 is left to right 

    assert R1['w'] == R2['w'] #R1 and R2 must have the same width
    assert y2 == y1 + R1['w'] #R1 and R2 must be adjacent with R1 below R2
    
    x0 = max(p1[0], p2[0])
    x1 = min(q1[0], q2[0])
    assert x0 < x1
    xm = (x0 + x1)/2 #x of midpoint of side
    ym = (q1[1] + q2[1])/2 #y of side
    xl = x1 - x0 #lenght of side
    tm1 = extrude_time_passage (p1,q1,(xm,y1))
    tm2 = extrude_time_passage (p2,q2,(xm,y2))
    j1 = len(R1['links'][1])
    i2 = len(R2['links'][0])
    side = (xm, ym, xl, R1, j1, R2, i2)
    link1 = (side, tm1)
    link2 = (side, tm2)
    R1['links'][1].append(link1) #Append side to outputs of {R1}
    R2['links'][0].append(link2) #Append side to inputs of {R2}
    return side

def extrude_time (p,q):
    '''
    Returns the times of extruding a stroke from point {p} to point {q},
    accounting for acceleration and deceleration.
    '''
    dx = p[0] - q[0]
    dy = p[1] - q[1]
    d = sqrt (dx*dx + dy*dy)
    acc_time = 0 #For now we are not taking into account acceleration delay
    acc_dist = 0 
    vel_max = 20 #Steady-state velocity (mm/s)
    if d >= 2 * acc_dist:
       t = 2 * acc_time + (d - 2 * acc_dist)/vel_max
    else:   
       assert False
    return t  

def air_time (p,q):
    '''
    Returns the times of moving the nozzle from point {p} to point {q},
    without extruding, accounting for nozzle lifting and dropping, acceleration and deceleration.
    '''
    if p == None or q == None:
        return 0
    dx = p[0] - q[0]
    dy = p[1] - q[1]
    d = sqrt (dx*dx + dy*dy)
    lift_time = 0.5 # time for lifting and dropping in seconds.
    acc_time = 0 #For now we are not taking into account acceleration delay
    acc_dist = 0 
    vel_max = 20 #Steady-state velocity (mm/s)
    if d >= 2 * acc_dist:
       t = 2 * acc_time + (d - 2 * acc_dist)/vel_max
    else:   
       assert False
    return t + 2 * lift_time 


def extrude_time_passage (p,q,m):
    '''
    Assumes that a stroke is being extruded from point {p} to point {q};
    returns the time from the start until when nozzle passes the point {m}.
    '''
    dx = p[0] - q[0]
    dy = p[1] - q[1]
    d = sqrt (dx*dx + dy*dy)
    gx = p[0] - m[0]
    gy = p[1] - m[1]
    g = sqrt (gx*gx + gy*gy)

    acc_time = 0 #For now we are not taking into account acceleration delay
    acc_dist = 0 
    vel_max = 20 #Steady-state velocity (mm/s)
    if d >= 2 * acc_dist:
        if g >= acc_dist and g <= d - acc_dist:  
           t = acc_time + (g - acc_dist)/vel_max
        else:   
           assert False
    else:   
       assert False
    return t   

def make_example ():
    w = 1.00 #1 mm
    N = 15
    R = [None]*N
    y = 0
    R[0] = create_stroke (0, (2,y), (4,y), w)
    R[1] = create_stroke (1, (5,y), (12,y), w)
    y = y + w
    R[2] = create_stroke (2, (3,y), (15,y), w)
    y = y + w
    R[3] = create_stroke (3, (2,y), (8,y), w)
    R[4] = create_stroke (4, (13,y), (18,y), w)
    y = y + w
    R[5] = create_stroke (5, (0,y), (5,y), w)
    R[6] = create_stroke (6, (11,y), (14,y), w)
    R[7] = create_stroke (7, (16,y), (21,y), w)
    y = y + w
    R[8] = create_stroke (8, (2,y), (7,y), w)
    R[9] = create_stroke (9, (11,y), (14,y), w)
    R[10] = create_stroke (10, (16,y), (21,y), w)
    y = y + w
    R[11] = create_stroke (11, (3,y), (21,y), w)
    y = y + w
    R[12] = create_stroke (12, (5,y), (12,y), w)
    R[13] = create_stroke (13, (16,y), (19,y), w)
    y = y + w
    R[14] = create_stroke (14, (7,y), (10,y), w)

    add_side (R[0], R[2])
    add_side (R[1], R[2])
    add_side (R[2], R[3])
    add_side (R[2], R[4])
    add_side (R[3], R[5])
    add_side (R[4], R[6])
    add_side (R[4], R[7])
    add_side (R[5], R[8])
    add_side (R[6], R[9])
    add_side (R[7], R[10])
    add_side (R[8], R[11])
    add_side (R[9], R[11])
    add_side (R[10], R[11])
    add_side (R[11], R[12])
    add_side (R[11], R[13])
    add_side (R[12], R[14])

    return R

def make_example_simple ():
    w = 1.00 #1 mm
    N = 12
    R = [None]*N
    y = 0
    R[0] = create_stroke (0, (0,y), (60,y), w)
    y = y + w
    R[1] = create_stroke (1, (0,y), (20,y), w)
    R[2] = create_stroke (2, (40,y), (60,y), w)
    y = y + w
    R[3] = create_stroke (3, (0,y), (20,y), w)
    R[4] = create_stroke (4, (40,y), (60,y), w)
    y = y + w
    R[5] = create_stroke (5, (0,y), (20,y), w)
    R[6] = create_stroke (6, (40,y), (60,y), w)
    y = y + w
    R[7] = create_stroke (7, (0,y), (20,y), w)
    R[8] = create_stroke (8, (40,y), (60,y), w)
    y = y + w
    R[9] = create_stroke (9, (0,y), (20,y), w)
    R[10] = create_stroke (10, (40,y), (60,y), w)
    y = y + w
    R[11] = create_stroke (11, (0,y), (60,y), w)

    add_side (R[0], R[1])
    add_side (R[0], R[2])

    add_side (R[1], R[3])
    add_side (R[2], R[4])
    
    add_side (R[3], R[5])
    add_side (R[4], R[6])

    add_side (R[5], R[7])
    add_side (R[6], R[8])

    add_side (R[7], R[9])
    add_side (R[8], R[10])

    add_side (R[9], R[11])
    add_side (R[10], R[11])

    return R


def print_links (edge, Rk, rbitk, ind):
    indx = ' ' * ind
    rbit_both = (Rk['rbit'] + rbitk) % 2
    L = Rk['links'][edge]
    for i in range(len(L)):
        linki = L[i]
        sys.stderr.write("%s  %3d" % (indx, i))
        side = linki[0]
        sys.stderr.write(" (%6.2f %6.2f) " % (side[0], side[1]))
        sys.stderr.write(" len = %6.2f" % (side[2]))
        if edge == 0:
            Rl = side[3]
            jl = side[4]
            sys.stderr.write("  R%03d out [%d] " % (Rl['sid'],jl))
            assert (side[5] == Rk)
            assert (side[6] == i)
        else:     
            Rl = side[5]
            jl = side[6]
            sys.stderr.write("  R%03d inp [%d] " % (Rl['sid'],jl))
            assert (side[3] == Rk)
            assert (side[4] == i)
        if rbit_both == 0:   
           trel = linki[1]
        else:   
           trel = Rk['e'] - linki[1]  
        sys.stderr.write("  trel = %7.3f\n" % (trel))

def print_stroke (Rk, rbitk, label, ind):
   '''
     Prints the stroke {Rk} inverted if {rbitk} = 1. If the stroke
     already has 'rbit' equal to 1, they cancel each other. {ind}
     is the number of blank to indent.
   '''
   indx = ' ' * ind
   rbit_both = (Rk['rbit'] + rbitk) % 2
   sys.stderr.write("%s Stroke R%03d:%d " % (indx, Rk['sid'],rbit_both))
   if label != None:
       sys.stderr.write(" [%s] " % label)
   sys.stderr.write(" (%6.2f %6.2f) " % (Rk['p'][rbit_both][0], Rk['p'][rbit_both][1]))
   sys.stderr.write(" (%6.2f %6.2f) " % (Rk['p'][1-rbit_both][0], Rk['p'][1-rbit_both][1]))
   if Rk['Tini'] != None: sys.stderr.write(" t = %7.3f " % (Rk['Tini']))
   if Rk['Tend'] != None: sys.stderr.write(" t = %7.3f " % (Rk['Tend']))
   sys.stderr.write(" e = %6.2f " % Rk['e'])
   sys.stderr.write(" w = %4.2f\n" % Rk['w'])
   sys.stderr.write("%s  Inputs: \n" % indx)
   print_links (0, Rk, rbitk, ind + 2)
   sys.stderr.write("%s  Outputs: \n" % indx)
   print_links (1, Rk, rbitk, ind + 2)

def print_strokes (R):
    for Rk in R:
        print_stroke (Rk, 0, None, 0)
        sys.stderr.write("\n")

def print_cands (Cands, ind):
    indx = ' ' * ind
    for k in range(len(Cands)):
        Ck = Cands[k]
        print_stroke (Ck[0], Ck[1], ('C%03d' % k), ind)
        sys.stderr.write("%s  Score: %.6f \n" % (indx, Ck[2]))
        sys.stderr.write("\n")

def heuristic1 (R):
    '''
    Returns the optimal order {S} of the strokes according to heurist 1,
    and also the total extrusion time {Tend} of {S} and also the maximum cooling
    time {Cmax} for any side, and the corresponding side {smax}. For now we assume that the input strokes are given 
    in the proper order: bottom to top rows, and left to right within each row.
    It does not extrude the raster links for now.
    '''

    Tend = 0
    Cmax = 0
    smax = None
    qend = None
    for k in range(len(R)):
        Rk = R[k]
        assert Rk['rbit'] == 0
        pk = Rk['p'][0]
        qk = Rk['p'][1]
        ek = Rk['e']
        if k > 0:
           Tini = Tend + air_time (qend,pk)
        else:   
           Tini = 0 
        Rk['Tini'] = Tini 
        Tend = Tini + ek
        Rk['Tend'] = Tend
        qend = qk
        for ik in range(len(Rk['links'][0])):
            linkik = Rk['links'][0][ik]
            side = linkik[0]
            Tik = Tini + linkik[1]
            Rl = side[3]
            jl = side[4]
            assert side[5] == Rk
            assert side[6] == ik
            assert Rl['rbit'] == 0
            linkjl = Rl['links'][1][jl]
            assert linkjl[0] == side
            Tjl = Rl['Tini'] + linkjl[1]
            Ck = Tik - Tjl #cooling time of segment
            if Ck > Cmax:
                Cmax = Ck
                smax = side
    return R, Tend, Cmax, smax

def get_initial_front_0 (R):
    '''
    Returns a list with the stroke R[0]
    '''
    L = [R[0]]
    return L

def get_initial_front_all (R):
    '''
    Returns a list with all strokes in R
    '''
    L = R.copy()
    return L

def get_initial_front_fringe (R):
    '''
    Returns all strokes in {R} that have no input or output sides.
    '''
    L = [].copy()
    for Rk in R:
       inputs = Rk['links'][0]
       outputs = Rk['links'][1]
       #if len(inputs) == 0 or len(outputs) == 0:
       if len(inputs) == 0:
          L.append (Rk)
    return L

def worst_cooling (qend, Tend, Rf, rbitf, ind):
   '''
     Assumes that the nozzle is currently at {qend} and
     the current time is {Tend}. Enumerates the sides of 
     the front stroke {Rf} that have been already half-covered
     and computes their cooling times, assuming that {Rf} will
     be extruded in the direction {rbitf}. Returns the largest 
     of those cooling times.
   '''
   indx = ' ' * ind
   if verbose: 
       sys.stderr.write ('%s### enter worst_cooling (Tend: %.6f) \n' % (indx, Tend))
       print_stroke (Rf, rbitf, 'Rf', ind)
   assert Rf['Tini'] == None  
   assert Rf['rbit'] == 0
   pf = Rf['p'][rbitf]
   qf = Rf['p'][1-rbitf]
   ef = Rf['e']
   Tmove = Tend + air_time (qend,pf)
   Coolmax = 0
   for edge in range(2):
      linksf = Rf['links'][edge]  
      for jf in range(len(linksf)):
          linkfj = linksf[jf] 
          side = linkfj[0]
          if edge == 0:
             Ra = side[3]  
             ja = side[4]
             assert side[5] == Rf
             assert side[6] == jf
          else: 
             Ra = side[5]  
             ja = side[6]
             assert side[3] == Rf
             assert side[4] == jf
          
          if Ra['Tini'] != None:
             linkaj = Ra['links'][1-edge][ja]
             rbita = Ra['rbit']
             ea = Ra['e']
             if verbose: print_stroke (Ra, 0, 'Ra', ind + 2)
             if rbitf == 0: 
                Tfj = Tmove + linkfj[1]
             else: 
                Tfj = Tmove + (ef - linkfj[1])

             if rbita == 0: 
                Taj = Ra['Tini'] + linkaj[1] 
             else: 
                Taj = Ra['Tini'] + (ea - linkaj[1])
             if verbose: sys.stderr.write ('Ra[Tini]: %.6f, Taj: %.6f, Tend: %.6f, Tfj: %.6f\n' % (Ra['Tini'], Taj, Tend, Tfj))   
             assert Taj < Tfj
             assert Taj < Tend
             assert Tend < Tfj
             Coolj = Tfj - Taj
             if Coolj > Coolmax:
                Coolmax = Coolj
   if verbose:
      sys.stderr.write ('%sCoolmax: %.6f\n' % (indx, Coolmax))   
      sys.stderr.write ('%s### exiting worst_cooling \n' % indx)
   return Coolmax          

def check_candidate (qend, Tend, Rk, rbitk, Front, Delta, ind):
   '''
     Checks whether choosing candidate {Rk}
     in the direction {rbitk} for the next stroke would
     not immediately imply excessive cooling time for any side in the current cutset.
   '''
   indx = ' ' * ind
   if verbose:
      sys.stderr.write ('%s### enter check_candidate: \n' % indx)
      print_stroke (Rk, rbitk, 'Candidate Rk', ind) 
   assert Rk['rbit'] == 0
   pk = Rk['p'][rbitk]
   qk = Rk['p'][1-rbitk]
   Tendk = Tend + air_time (qend,pk) + Rk['e']
   for Rf in Front:
      if verbose: print_stroke (Rf, 0, 'Front Rf') 
      if Rf == Rk:
         Coolf = worst_cooling (qend, Tend, Rk, rbitk, ind)
      else:
         Coolf0 = worst_cooling (qk, Tendk, Rf, 0, ind) 
         Coolf1 = worst_cooling (qk, Tendk, Rf, 1, ind) 
         Coolf = min (Coolf0, Coolf1)
      if Coolf > Delta:
         sys.stderr.write ('%scooling time exceeded for candidate R%03d:%d (Tend: %.6f) between R%03d and some neighbor: %.6f\n' % (indx, Rk['sid'], rbitk, Tendk, Rf['sid'], Coolf)) 
         if verbose: sys.stderr.write ('%s### exiting check_candidate (false) \n' % indx)
         return False
   if verbose: sys.stderr.write ('%s### exiting check_candidate (true) \n' % indx)
   return True  

def exclude_bad_candidates (qend, Tend, Cands, Front, Delta, ind):
    '''
    Excludes from the list {Cands} every candidate which, if it was
    stroked next, would cause some half-covered side to exceed the
    max allowed cooling time {Delta}. Assumes that all half-covered sides
    are inputs or outputs of strokes in the {Front}. Assumes also that the nozzle
    is currently at {qend} and the current time is {Tend}. {Cands} is a list of triples
    {(R, rbit, tmove)} where {R} is a stroke that has not been extruded yet and has {R['rbit']=0},
    and {rbit} is a 0 or 1, and {tmove} is the air time from {qend} to the start of the candidate.
    '''
    indx = ' ' * ind
    if verbose: sys.stderr.write ('%s### enter exclude_bad_candidates (Tend: %.6f): \n' % (indx, Tend))
    if qend == None:
       if verbose: sys.stderr.write ('%s### exiting exclude_bad_candidates (trivial) \n' % indx)
       return Cands.copy()

    NewCands = [].copy()
    
    for Ck in Cands:
       assert len(Ck) == 3 
       Rk = Ck[0]
       rbitk = Ck[1]
       if check_candidate (qend, Tend, Rk, rbitk, Front, Delta, ind):
          NewCands.append(Ck)
    if verbose: sys.stderr.write ('%s### exiting exclude_bad_candidates (eliminated %d) \n' % (indx, len(Cands) - len(NewCands)))
    return NewCands      
    
def choose_candidate_dist (Cands, qend, ind):
    indx = ' ' * ind
    if verbose: sys.stderr.write ('%s### enter choose_candidate: \n' % indx)
    tmin = 9999999
    Cmin = None
    for Ck in Cands:
       Rk = Ck[0]
       rbitk = Ck[1]
       pk = Rk['p'][rbitk]
       if qend != None:
          tmove = air_time (qend,pk)
       else:   
          tmove = 0 
       if tmove < tmin:
          tmin = tmove
          Cmin = Ck
    if verbose:      
       sys.stderr.write ('%s### exiting choose_candidate \n' % indx)
       print_stroke(Cmin[0],Cmin[1],'choosen')
    return Cmin 

def choose_candidate_raster (Cands, qend, ind):
    indx = ' ' * ind
    if verbose: sys.stderr.write ('%s### enter choose_candidate: \n' % indx)
    xymin = 9999999
    Cmin = None
    for Ck in Cands:
       Rk = Ck[0]
       rbitk = Ck[1]
       xy = raster_order (Rk) 
       if xy < xymin:
          xymin = xy
          Cmin = Ck
    if verbose: 
       sys.stderr.write ('%s### exiting choose_candidate \n' % indx)
       print_stroke(Cmin[0],Cmin[1],'choosen')
    return Cmin 

def update_front (Front, Rk):
    '''
    Removes the stroke {Rk} from the {Front} and adds
    all output strokes of {Rk} that are not yet on the {Front}. Returns
    the updated {Front}.
    '''
    NewFront = [].copy()
    for Fk in Front:
        if Fk['sid'] != Rk['sid']:
           NewFront.append (Fk)
    for edge in range(2):
       for linki in Rk['links'][edge]:
          sidei = linki[0]
          if edge == 0: 
              Rl = sidei[3]
              assert sidei[5] == Rk
          else: 
              Rl = sidei[5]
              assert sidei[3] == Rk
          if Rl not in NewFront and Rl['Tini'] == None:
              NewFront.append (Rl)        
    return NewFront

def heuristic2 (R, Delta):
    '''
    Returns the optimal order {S} of the strokes according to heuristic 2,
    with maximum allowed cooling {Delta}, and also the total extrusion time {Tend} of {S} and also the maximum cooling
    time {Cmax} for any side, and the corresponding side {smax}. 
    It does not extrude the raster links for now. Each element of {S} is a pair (stroke, reversebit).

    The solution will start with stroke R[0] and each new stroke will be adjacent at least one previously extruded
    stroke.
    '''
    Tend = 0
    Cmax = 0
    smax = None
    qend = None
    S = [].copy() 
    #Front = get_initial_front_0 (R) # current front
    Front = get_initial_front_fringe (R) # current front
    for k in range(len(R)):
        Cands = [].copy() #List of pairs (stroke, reversebit) 
        for i in range(len(Front)): 
            assert Front[i]['rbit'] == 0
            Cands.append((Front[i],0,air_time(qend,Front[i]['p'][0])))
            Cands.append((Front[i],1,air_time(qend,Front[i]['p'][1])))
        Cands = exclude_bad_candidates (qend, Tend, Cands, Front, Delta, 0)
        if len(Cands) == 0:
           sys.stderr.write("*************FAILED\n")    
           assert False  
        Ck = choose_candidate_raster (Cands, qend)
        #Ck = choose_candidate_dist (Cands, qend)
        Rk = Ck[0]
        rbitk = Ck[1]
        pk = Rk['p'][rbitk]
        qk = Rk['p'][1-rbitk]
        ek = Rk['e']
        assert Rk['Tini'] == None
        assert Rk['Tend'] == None
        assert Rk['rbit'] == 0
        if k > 0:
           Tini = Tend + air_time (qend,pk)
        else:   
           Tini = 0 
        Rk['Tini'] = Tini 
        Tend = Tini + ek
        Rk['Tend'] = Tend
        qend = qk
        Rk['rbit'] = rbitk
        S.append(Rk)
        sys.stderr.write("---------------------------------------------\n")    
        print_stroke (Rk, rbitk, 'appended')
        sys.stderr.write("---------------------------------------------\n")
        Cmax_k, smax_k = worst_cooling_one (Rk)
        if Cmax_k > Cmax:
           Cmax = Cmax_k
           smax = smax_k
        #Atualize the front and cutset.
        Front = update_front (Front, Rk)
    return S, Tend, Cmax, smax

def sort_candidates_none (qend, Cands):
   '''
     Sorts the candidate list based on air time from {qend} to the starting point.
     Returns the sorted list.
   '''
   return Cands.copy()    

def sort_candidates_dist (qend, Cands):
   '''
     Sorts the candidate list based on air time from {qend} to the starting point.
     Returns the sorted list.
   '''
   return(sorted(Cands, key = lambda x: x[2]))    

def raster_order (Rk):
   '''
     Returns combination of {Y} and {X} coordinates of stroke {Rk} adequate for
     raster order sorting.
   '''
   assert Rk['rbit'] == 0
   y = Rk['p'][0][1]
   x = Rk['p'][0][0]
   assert y == Rk['p'][1][1] #Stroke must be horizontal
   assert x < Rk['p'][1][0] #Stroke must be left to right
   return 10000*y + x

def sort_candidates_raster (qend, Cands):
   '''
     Sorts the candidate list based on Y with ties broken by X.
     Returns the sorted list.
   '''
   return(sorted(Cands, key = lambda Ck: raster_order(Ck[0])+1000000*Ck[1]))    


def heuristic2R (R, Delta):
    '''
    Returns the optimal order {S} of the strokes according to heuristic 2,
    with maximum allowed cooling {Delta}, and also the total extrusion time {Tend} of {S} and also the maximum cooling
    time {Cmax} for any side, and the corresponding side {smax}. 
    It does not extrude the raster links for now. Each element of {S} is a pair (stroke, reversebit).

    The solution will start with stroke R[0] and each new stroke will be adjacent at least one previously extruded
    stroke. Uses backtracking.
    '''
    Front = get_initial_front_fringe (R) # current front
    S, Cmax, smax = heuristic2R_aux ([].copy(), Front, Delta)
    if S == None:
        Tend = None
    else:    
        Tend = S[len(S)-1]['Tend']
    return S, Tend, Cmax, smax


def worst_cooling_one (Rk):
    '''
    Computes the worst cooling time for the sides of a stroke {Rk} that
    has just been added to the solution list. Returns the worst cooling time
    {Cmax} and the corresponding side {smax}.
    '''
    assert Rk['Tini'] != None
    Cmax = 0
    smax = None
    ek = Rk['e']
    for edge in range(2): 
      linksk = Rk['links'][edge]
      for ik in range(len(linksk)):
        linkik = linksk[ik]
        side = linkik[0]
        if Rk['rbit'] == 0:
           Tik = Rk['Tini'] + linkik[1]
        else:  
           Tik = Rk['Tini'] + (ek - linkik[1])
        if edge == 0:
          Rl = side[3]
          jl = side[4]
          assert side[5] == Rk
          assert side[6] == ik
        else:
          Rl = side[5]
          jl = side[6]
          assert side[3] == Rk
          assert side[4] == ik
        linkjl = Rl['links'][1-edge][jl]
        assert linkjl[0] == side
        if Rl['Tini'] != None:
           if Rl['rbit'] == 0: 
               Tjl = Rl['Tini'] + linkjl[1]
           else:    
               Tjl = Rl['Tini'] + (Rl['e'] - linkjl[1])
           Cooljk = Tik - Tjl #cooling time of segment 
           if Cooljk > Cmax:
              Cmax = Cooljk
              smax = side
    return Cmax, smax          

def heuristic2R_aux (S, Front, Delta):
    '''
      Given a partial solution {S} for the stroke order optimization problem, and
      the corresponding list of {Front} strokes tries to complete it into a full
      solution {Sfull}. If it succeds returns the {Sfull}, the maximum cooling time
      {Cmax} among the new strokes added, and the corresponding side {smax}. If it fails returns None in place of
      {Sfull}. Assumes that the {Front} has at least one stroke in each component of the connectivity graph that
      has not been yet completey extruded.
    '''
    ind = 2 * len(S)
    indx = ' ' * ind
    if len(Front) == 0:
        return S, 0, None
    N = len(S)
    if N == 0:
        Tend = 0
        qend = None
    else:    
        Tend = S[len(S)-1]['Tend']
        rbitend = S[len(S)-1]['rbit']
        qend = S[len(S)-1]['p'][1-rbitend]
    Cands = [].copy() #List of pairs (stroke, reversebit) 
    for i in range(len(Front)): 
        assert Front[i]['rbit'] == 0
        Cands.append((Front[i],0,air_time(qend,Front[i]['p'][0])))
        Cands.append((Front[i],1,air_time(qend,Front[i]['p'][1])))
    Cands = exclude_bad_candidates (qend, Tend, Cands, Front, Delta, ind)
    Cands = sort_candidates_raster (qend, Cands)
    #Cands = sort_candidates_dist (qend, Cands)
    #Cands = sort_candidates_none (qend, Cands)
    if len(S) == 999:
        print_cands (Cands, ind)
        sys.exit()
    for k in range(len(Cands)):
        Ck = Cands[k]
        Rk = Ck[0]
        assert Rk['rbit'] == 0
        rbitk = Ck[1]
        pk = Rk['p'][rbitk]
        qk = Rk['p'][1-rbitk]
        ek = Rk['e']
        assert Rk['Tini'] == None
        assert Rk['Tend'] == None
        if qend != None:
            Tinik = Tend + air_time (qend,pk)
        else:   
            Tinik = 0 
        Rk['Tini'] = Tinik 
        Tendk = Tinik + ek
        Rk['Tend'] = Tendk
        qend = qk
        Rk['rbit'] = rbitk
        S.append(Rk)
        if verbose:
           sys.stderr.write("%s---------------------------------------------\n" % indx)    
           print_stroke (Rk, 0, 'appended', ind)
           sys.stderr.write("%s---------------------------------------------\n" % indx)
        else:   
           sys.stderr.write("%sappended R%03d:%d, Tini = %.6f, Tend = %.6f \n" % (indx, Rk['sid'], Rk['rbit'], Rk['Tini'], Rk['Tend']))    
        Cmax_k, smax_k = worst_cooling_one (Rk)
        NewFront = update_front (Front, Rk)
        S_rec, Cmax_rec, smax_rec = heuristic2R_aux (S, NewFront, Delta)
        if S_rec != None:
            if Cmax_rec > Cmax_k:
                return S_rec, Cmax_rec, smax_rec
            else:
                return S_rec, Cmax_k, smax_k
        else:
            Rk['rbit'] = 0
            Rk['Tini'] = None
            Rk['Tend'] = None
            S.pop()
            if verbose:
               sys.stderr.write("%s---------------------------------------------\n" % indx)    
               print_stroke (Rk, rbitk, 'deleted', ind)
               sys.stderr.write("%s---------------------------------------------\n" % indx)
            else:   
               sys.stderr.write("%sdeleted R%03d:%d \n" % (indx, Rk['sid'], rbitk))    

    sys.stderr.write("%s**FAILED\n" % indx)    
    return None, None, None

R = make_example ()
#R = make_example_simple ()

#print_strokes (R)

S, Tend, Cmax, smax = heuristic1 (R)

#Delta = 8 
#S, Tend, Cmax, smax = heuristic2 (R, Delta)

#Delta = 5
#S, Tend, Cmax, smax = heuristic2R (R, Delta)

if S != None:
   sys.stderr.write("*************SUCCEDED\n")
   print_strokes (S)
   sys.stderr.write ('Total extrusion time: %7.3f\n' % Tend)
   sys.stderr.write ('Maximum cooling time: %7.3f ' % Cmax)
   sys.stderr.write ('from stroke R%03d output [%d] ' % (smax[3]['sid'],smax[4]))
   sys.stderr.write ('to stroke R%03d output [%d]\n' % (smax[5]['sid'],smax[6]))
   plot_strokes (S)
else:
   sys.stderr.write("*************FAILED\n")