#!/usr/bin/env python
# batchruntomo - run one or more data sets in batch mode
#
# Author: David Mastronarde
#
# $Id: batchruntomo,v 4d0fa707e305 2025/04/01 21:24:27 mast $
#

progname = 'batchruntomo'
prefix = 'ERROR: ' + progname + ' - '

# Index of batch directives in allDirectives array, and of first template
batDictInd = 4
tmplDictInd = 1

# Variables that need to be global, probably superfluous but here as an FYI
global datasetDir, dualAxis, setName, scanHeader, ifMontage, defocus, pixelSize
global fidSizeNm, fidSizePix, montFrameData, numSurfaces
global rawXsize, rawYsize, zsize, fiducialless, coarseBinning, centerOnGold, posSampleType
global xtiltNeeded, fidThickness, fidIncShift, reconThickness, didLocalAlign, madeZfactors
global aliBinning, aliXunbinned, aliYunbinned, patchTrack, totalDelTilt, latestMessages
global suppressAbort, absDirectiveFile, userTemplateDir, summaryMessage, logFile
global finalRetVal, correctCTF, eraseGold, axisLet, axisNum, useVolMatch, nxRecA, nyRecA
global finishSetAndQuit, finalAlignResid, axisRotation, transposeForAli, edfLines
global edfChanged, excludedViewsA, excludedViewsB, noXaxisTilt, stackExtension
global fromExtension, originalStackExt, autofitCTF, ctf3dSlabThick, doBothRecons
global rawFor3dCTF, skipAlignStackMods, filterIn2D, myParallelGPU, myGPUlist
global needGPUrelease, alignLines, makeSymLinks, lastTranslateInd, dsetDirTranslateInd
global expandFactor

# Print a string to both standard out and a log file if one is defined
def prnLog(string, end  = '\n', flush = False):
   prnstr(string, end = end, flush = flush)
   if logFile:
      prnstr(string, file = logFile, end = end, flush = flush)


# Close the log file, Backup and save the edf file if it is modified
def closeLogFileWriteEdf():
   global logFile, edfLines, edfChanged
   if logFile:

      # Etomo needs the comma
      prnstr('Batchruntomo finished with data set, ' + datetime.datetime.now().ctime() +\
                '  [brt12]')
      logFile.close()
      logFile = None

   if edfChanged:
      edfFile = setName + '.edf'
      makeBackupFile(edfFile + '~')
      makeBackupFile(edfFile)
      if writeTextFile(edfFile, edfLines, True):
         prnstr('Error writing modified ' + edfFile)
      edfChanged = False

   edfLines = []


# Print a warning with prefix WARNING: and blank line at the end, to log and output
# by default, and flushes
def warning(strings, toLog = True):
   if isinstance(strings, str):
      strings = [strings]
   strings.append(' ')
   strings[0] = 'WARNING: ' + strings[0]
   for line in strings:
      if toLog:
         prnLog(line, flush = True)
      else:
         prnstr(line, flush = True)
         
      
# Send an email if there is an address
def sendEmail(subject, message):
   if not emailAddress:
      return
   try:
      msg = MIMEText(message)
      msg['Subject'] = subject
      msg['From'] = 'batchruntomo'
      msg['To'] = emailAddress
      s = smtplib.SMTP(SMTPserver)
      s.sendmail(emailAddress, [emailAddress], msg.as_string())
      s.quit()
   except Exception:
      if SMTPserver == 'localhost':
         warning('Failed to send email notification; you probably need to specify ' + \
                 'an SMTP server')
      else:
         warning('Failed to send email notification')
   

# Abort either an axis or a data set with the error string
def abortSet(errString):
   global summaryMessage, finalRetVal
   if suppressAbort:
      return
   if renamingOnly:
      exitError('renaming only - ' + errString)
      
   abortStr = ['ABORT SET: ', 'ABORT AXIS: ', 'ABORT AXIS: ']
   prnLog(abortStr[axisNum] + errString)
   prnLog('', flush = True)
   if dualAxis and axisNum:
      message = 'Batchruntomo aborted axis ' + axisLet.upper() + ' of'
   elif dualAxis:
      message = 'Batchruntomo aborted combine of'
   else:
      message = 'Batchruntomo aborted'
   message += ' dataset ' + setName + ' after error:\n' + errString + '\n'
   sendEmail('Batchruntomo error on ' + setName, message)
   summaryMessage += message
   if not axisNum:
      finalRetVal += 1
   if exitOnError:
      sys.exit(1)


# Report the error from running a program and then abort
def reportImodError(abortText = None):
   errStrings =  getErrStrings()
   num = len(errStrings)
   for ind in range(num):
      l = errStrings[ind]
      if ind == num - 1:
         l = 'ERROR: ' + l
      prnLog(l, end = '')
   if abortText:
      abortSet(abortText)


# Renames a file, aborts set if failure with formatted message.  formatStr must have two 
# {} for from and to name
def renameAndAbort(fromName, toName, formatStr = 'Renaming {} to {}'):
   try:
      os.rename(fromName, toName)
   except OSError:
      abortSet(fmtst(formatStr, fromName, toName))
      return 1
   return 0


# Test whether a value returned by lookupDirective is a string and abort with message
def testDirectiveValue(val, directive, dtype):
   if isinstance(val, str):
      abortSet('An error occurred converting the value of the directive ' + directive + \
                  ' to a ' + dtype)
      return 1
   return 0


# Print lines starting with tags in a log file.  Tags is an array of duples, the first
# element of each is the tag itself, the second is the sum of 1 if it is used for
# multiline messages and 2 if it should be stripped and 4 if it is not at start of a line.
# Optionally, a tag can have a third element, a suffix to put at the end of the line
# when it is printed.  This suffix should include all desired leading spaces.  The tag
# itself may have ONE '.*' in order to match lines containing elements on both sides
def printTaggedMessages(logfile, tags):
   global latestMessages
   blankNeeded = ('ERROR', 'INFO', 'WARNING')
   needFinalBlank = False
   needBlank = ''
   latestMessages = []
   anyError = False
   if isinstance(logfile, str):
      loglines = readTextFile(logfile, None, True)
      if isinstance(loglines, str):
         warning('Error ' + loglines)
         return
   else:
      loglines = logfile
      
   readingMulti = False
   for l in loglines:
      if readingMulti:
         prnLog(l)
         needFinalBlank = True
         latestMessages.append(l)
         if l.strip() == '':
            readingMulti = False
            needBlank = ''
            needFinalBlank = False
            
      else:
         for tag in tags:
            if '.*' in tag[0]:
               index = tag[0].find('.*')
               if tag[1] & 4:
                  match = tag[0][:index] in l and tag[0][index+2:] in l
               else:
                  match = l.startswith(tag[0][:index]) and tag[0][index+2:] in l
            elif tag[1] & 4:
               match = tag[0] in l
            else:

               # Workaround to badly fixed bug 2503 in Etomo: suppress No GUI
               match = l.startswith(tag[0]) and l.find('No GUI') < 0
               
            if match:
               # Suppress uninformative vmstopy output if there is already a message
               if not anyError or not (l.startswith('ERROR:') and \
                                          l.endswith('exited with status 1')):
                  if tag[0].startswith('ERROR'):
                     anyError = True

                  # If we needed a blank for a particular kind of message, and this is
                  # not another of that type, put the blank out
                  if needBlank and not tag[0].startswith(needBlank):
                     prnLog('')
                     needBlank = ''
                     needFinalBlank = False
                     
                  # Print message one way or the other
                  latestMessages.append(l)
                  lOut = l
                  if tag[1] & 2:
                     lOut = l[len(tag[0]):].lstrip()
                  if len(tag) > 2:
                     lOut += tag[2]
                  prnLog(lOut)
                  if len(tag) > 2 and tag[2] == logSuffixTag:
                     prnLog('')
                  needFinalBlank = True

                  # See if a blank is needed after contiguous lines of this type
                  if not needBlank:
                     for need in blankNeeded:
                        if tag[0].startswith(need):
                           needBlank = need
                           break
                  
               if tag[1] & 1:
                  readingMulti = True
               break

   if needFinalBlank:
      prnLog(' ')
            

# Look for a tag in a set of lines, and get the value after the separator
def findTaggedValue(lines, tag, separator, valType):
   for l in lines:
      ind = l.find(separator)
      if ind > 0 and tag in l and ind < len(l) - 1:
         valAll = l[ind + 1:].strip()
         if valType == STRING_VALUE:
            return valAll
         vsplit = valAll.split()
         if not len(vsplit):
            return None
         try:
            if valType == INT_VALUE:
               value = int(vsplit[0])
            else:
               value = float(vsplit[0])
            return value
         except:
            return None
   return None


# Write a text file and report a returned error
def writeTextFileReportErr(filename, lines):
   err = writeTextFile(filename, lines, True)
   if err:
      abortSet('Error ' + err)
   return err


# Read a text file and report a returned error and return [] on error
def readTextFileReportErr(filename, message = None):
   lines = readTextFile(filename, message, True)
   if isinstance(lines, str):
      abortSet('Error ' + lines)
      return []
   if len(lines) == 0:
      abortSet('File ' + filename + ' is empty')
   return lines
   

# When doing parallel batch, translate the path put out originally to one needed here
def translateParallelPath(dfile):
   global lastTranslateInd
   fileIn = dfile
   lastTranslateInd = -1
   if fromParallelPaths and toParallelPaths:
      for ind in range(numTranslations):
         fromPath = fromParallelPaths[ind]
         if dfile.startswith(fromPath):
            dfile = dfile.replace(fromPath, toParallelPaths[ind])
            lastTranslateInd = ind
            break
   
   return dfile


# Reverse translate a path using the given translation index or skip if negative
# This gets translated paths back to the form on the Etomo computer for messages it wants
def reverseTranslatePath(dfile, transInd):
   if not fromParallelPaths or transInd < 0:
      return dfile
   toPath = toParallelPaths[transInd]
   if dfile.startswith(toPath):
      dfile = dfile.replace(toPath, fromParallelPaths[transInd])

   return dfile


# Looks for stacks with the given name (including directory prefix and dot) and all
# possible extensions, giving priority to standard ones.  Returns duple with index of
# first and second extension found, either among the standard ones if any is found, or
# among alternate ones, and their Z sizes
def findPossibleStacks(stackRoot, alreadyStackExt):
   numPossible = len(possibleStackExts)
   exists1 = -1
   exists2 = -1
   size1 = 0
   size2 = 0
      
   for ind in range(numPossible):
      possible = possibleStackExts[ind]
      if fromExtension and fromExtension != possible and \
         (not alreadyStackExt or alreadyStackExt != possible):
         continue
      name = stackRoot + possible
      if os.path.exists(name):
         try:
            (nx, ny, nz) = getmrcsize(name)
            if exists1 < 0:
               exists1 = ind
               size1 = nz

            # Make sure a file with expected extension is included in the two
            elif exists2 < 0 or (typeExtension and stackExtension and \
                                 stackExtension == possible):
               exists2 = ind
               size2 = nz

         except ImodpyError:
               warning('Could not read header of possible stack file ' + name)
            

      # At end of standard extensions, simply stop if there is one
      if ind == len(standardTypeExts) -1:
         if exists1 >= 0:
            break

   # Ignore two files if either one has a size of 1
   if exists2 >= 0 and size1 == 1:
      exist1 = exists2
      size1 = size2
      exists2 = -1
   if exists2 >= 0 and size2 == 1:
      exists2 = -1

   return (exists1, exists2, size1)


# Check existence of stack, rename if necessary.  Stack includes dot at end
def checkRenameStack(stack):
   global stackExtension, fromExtension, originalStackExt
   (existInd1, existInd2, size1) = findPossibleStacks(stack, extIfAlreadySet)
   posExt = possibleStackExts[max(0, existInd1)]

   if existInd1 >= 0 and existInd2 >= 0:
      abortSet('There are two possible stack files in the dataset directory : ' + 
               stack + posExt + ' and ' + stack + possibleStackExts[existInd2])
      return 1

   if existInd1 < 0:
      if fromExtension:
         abortSet('Stack file does not exist: ' + stack + fromExtension)
      else:
         abortSet('Stack file does not exist with any allowed extension: ' + stack)
      return 1

   # Keep track of original stack extension for putting in edf file
   if not originalStackExt:
      originalStackExt = posExt

   # If there is a type extension, keep the file as is and record the stack extension
   fromExtension = ''
   if not stackExtension and (typeExtension or not existInd1):
      stackExtension = '.' + posExt
      return 0

   if (existInd1 > 0 and not typeExtension) or (typeExtension and stackExtension and \
                                                stackExtension[1:] != posExt):

      # Otherwise rename the stack to st or to the extension established by first axis
      oldName = stack + posExt
      newName = stack + 'st'
      if stackExtension:
         newName = stack + stackExtension[1:]
      if os.path.exists(newName):
         abortSet('Cannot rename stack from ' + oldName + ' to: ' + newName + \
                  ' because a single-image file with that name already exists')
         return 1

      try:

         # Etomo is looking for "to:" to find new name
         os.rename(oldName, newName)
         prnLog('[brt9]  Renamed stack from ' + oldName + ' to: ' + newName)
         if not stackExtension:
            stackExtension = '.st'
      except OSError:
         abortSet('Error renaming stack from ' + oldName + ' to ' + newName)
         return 1

   return 0


# Test for whether symbolic links should be made instead of renaming/moving
def testForSymLink():
   global makeSymLinks
   if makeSymLinks == None:
      makeSymLinks = False
      if lookupDirective(setupPrefix, 'makeSymbolicLinks', 0, BOOL_VALUE):
         if 'win32' in sys.platform or 'cygwin' in sys.platform:
            prnLog('Directive to make symbolic links instead of delivering is ' + \
                   'ignored on Windows')
         else:
            makeSymLinks = True
            prnLog('Making symbolic links instead of delivering files to data ' + \
                   'set directory')
   return makeSymLinks


# Delivery some associated file with the given source and destination names
def deliverAncillary(source, dest, typeName, fullExt):
   if os.path.exists(source) and not os.path.exists(dest):
      try:
         if testForSymLink():
            os.symlink(source, dest)
         else:
            os.rename(source, dest)
      except OSError:
         abortSet(fmtstr('Error moving {} file {}.{} from {} to {}', setName, fullExt,
                         deliverFromDir, datasetDir))
         return 1
   return 0
   
      
# Move a stack to the dataset directory
def deliverStack(axis):
   global stackExtension, fromExtension
   testForSymLink()
   stackRoot = setName + axis
   stackRootDot = stackRoot + '.'
   sourceRoot = os.path.join(deliverFromDir, stackRootDot)
   destRoot = os.path.join(datasetDir, stackRootDot)
   
   # Etomo may change the current dir after delivery, so if they match, skip
   if deliverFromDir == datasetDir:
      return 0

   # Get existence of files and their Z sizes from both locations
   (srcInd1, srcInd2, srcSize1) = findPossibleStacks(sourceRoot, '')
   (destInd1, destInd2, destSize1) = findPossibleStacks(destRoot, extIfAlreadySet)
   srcExt = possibleStackExts[max(0, srcInd1)]
   destExt = possibleStackExts[max(0, destInd1)]

   # Handle clear conflict cases
   if destInd1 >= 0 and destInd2 >= 0:
      abortSet('There are two possible stack files in the dataset directory : ' + 
               stackRootDot + destExt + ' and ' + 
               stackRootDot + possibleStackExts[destInd2])
      return 1

   if srcInd1 >= 0 and srcInd2 >= 0:
      abortSet('There are two possible stack files in ' + deliverFromDir + ' : ' +
               stackRootDot + srcExt + ' and ' + 
               stackRootDot + possibleStackExts[srcInd2])
      return 1

   # Nothing already there and nothing available
   if srcInd1 < 0 and destInd1 < 0:
      abortSet('Stack file ' + stackRoot + ' does not exist in ' + \
                deliverFromDir + ' with any allowed extension')
      return 1

   # Something already there and something available, and not only one is single-image
   if srcInd1 >= 0 and destInd1 >= 0 and ((srcSize1 > 1 and destSize1 > 1) or \
                                          (srcSize1 == 1 and destSize1 == 1)):
      
      mess = 'A stack file ' + stackRootDot + destExt + ' already exists in ' + \
               datasetDir + ' but there is a possible stack file ' + stackRootDot +\
               srcExt + ' in ' + deliverFromDir
      if makeSymLinks:
         prnLog('WARNING: ' + mess)
         prnLog('Assuming this is the result of making a symbolic link previously')
      else:
         abortSet(mess)
         return 1
      
   # Nothing available or single-image available but stack already there
   renaming = False
   delivering = False
   fromExtension = ''
   if srcInd1 < 0 or (destInd1 >= 0 and (srcSize1 == 1 or makeSymLinks)):
      prnLog('Stack file is already delivered to dataset directory; assuming ' + \
             'associated files are too')

      # The existing extension is OK if: no type extension and it is st, or
      # type extension and either no stack extension yet, or it matches established one
      if (not typeExtension and not destInd1) or \
         (typeExtension and (not stackExtension or stackExtension[1:] == destExt)):
         if not stackExtension:
            stackExtension = '.' + destExt
         return 0
      renaming = True;
      if not stackExtension:
         stackExtension = '.' + destExt
      source = destRoot + destExt

   # True delivery
   else:
      delivering = True
      source = sourceRoot + srcExt
      if not stackExtension:
         if typeExtension:
            stackExtension = '.' + srcExt
         else:
            stackExtension = '.st'
      renaming = srcExt != stackExtension[1:]

   dest = destRoot + stackExtension[1:]

   # Test if there is a file there: which should be single-image given tests above
   if os.path.exists(dest):
      abortSet('A single-image file named ' + os.path.basename(dest) + \
               ' already exists in ' + datasetDir)
      return 1
   try:
      if delivering and makeSymLinks:
         os.symlink(source, dest)
      else:
         os.rename(source, dest)

      # Etomo is looking for "to:" to find new name/location
      if delivering:
         prnLog('[brt8]  Delivered stack ' +
                reverseTranslatePath(source, delFromTranslateInd) + ' to: ' +
                reverseTranslatePath(dest, dsetDirTranslateInd))
      if renaming:
         prnLog('[brt9]  Renamed stack from ' +
                reverseTranslatePath(os.path.basename(source), delFromTranslateInd) + \
                ' to: ' + reverseTranslatePath(stackRoot, dsetDirTranslateInd) + \
                stackExtension)
   except OSError:
      abortSet('Error moving/renaming stack file from ' + source + ' to ' + dest)
      return 1

   # Move .mdoc, .log, .rawtlt files also
   if deliverAncillary(source + '.mdoc', dest + '.mdoc', 'metadata', 
                       stackExtension + '.mdoc'):
      return 1
   if deliverAncillary(sourceRoot + 'log', destRoot + 'log', 'tilt series log', 'log'):
      return 1
   if deliverAncillary(sourceRoot + 'rawtlt', destRoot + 'rawtlt', 'raw tilt angle',
                       'rawtlt'):
      return 1

   return 0
   

# Process option and environment variable for queue command and max jobs
def getQueueOptions(queueOption, queueEnvVar, maxOption, maxEnvVar, gpuText):
   runningOnIt = False
   maxJobs = 0
   if os.getenv(queueEnvVar):
      command = os.environ[queueEnvVar]
      if command == 'None':
         command = None
      else:
         runningOnIt = True
   else:
      command = PipGetString(queueOption, '')

   if command:
      if cpuList and not coresPerClusterJob:
         exitError('You cannot enter a CPU machine list with a queue command') 
      if not parallelRoot:
         exitError('Use of a cluster queue is allowed only when running batch in ' + \
                   'parallel')
      
      maxJobs = PipGetInteger(maxOption, 0)
      maxEntered = 1 - PipGetErrNo()
      if os.getenv(maxEnvVar):
         try:
            maxJobs = int(os.environ[maxEnvVar])
            maxEntered = 1
         except ValueError:
            exitError('Converting environment variable ' + maxEnvVar + ' to integer')

      if maxJobs <= 0 or not maxEntered:
         exitError(fmtstr('Maximum number of {0}jobs must be entered if a {0}cluster' + \
                          ' command is entered', gpuText))
         
   return (command, runningOnIt, maxJobs)


# Process option and environment variable for cores or GPUs per job
def getClusterJobOption(option, envVar):
   entered = 0
   if os.getenv(envVar) != None:
      envar = os.environ[envVar]
      if envar == 'None':
         perJob = 0
      else:
         try:
            perJob = int(envar)
            entered = 1
         except ValueError:
            exitError('Environment variable ' + envVar + ' must be an integer')

   else:
      perJob = PipGetInteger(option, 0)
      entered = 1 - PipGetErrNo()

   if entered:
      if perJob <= 0:
         exitError('The value for ' + option + ' or ' + envVar + ' must be positive')
      if not parallelRoot:
         exitError('Cluster node options can be entered only when running batch in ' + \
                   'parallel')

   return perJob


# Generate two sed commands for adding an edf file entry, deleting if necessary
# An alternate delimiter can be used as long as it used for all commands passed on one
# call to modifyEdfLines
def edfDelAndAdd(option, value, delim = '/'):
   option = 'batchruntomo.' + option
   return [fmtstr('{0}{1}={0}d', delim, option),
           fmtstr('{0}{1}{0}a{0}{2}={3}{0}', delim, 'Setup.DatasetName=', option, value)]


# Return a 'true' or 'false' string for placing in the edf file
def boolStringForEdf(value):
   if value:
      return 'true'
   return 'false'


# Read the edf file if necessary
def readEdfFileIfNeeded():
   global edfLines
   if edfLines:
      return 0
   edfLines = readTextFileReportErr(setName + '.edf')
   if not edfLines:
      if bypassEtomo:
         prnLog('Ignore that error!  Starting a blank edf file.')
         edfLines = ['Setup.DatasetName=' + setName]
         return 0
      return 1
   return 0


# Modify the edf lines
def modifyEdfLines(sedcom):
   global edfLines, edfChanged
   if readEdfFileIfNeeded():
      return 1
   newLines =  pysed(sedcom, edfLines, retErr = True)
   if isinstance(newLines, str):
      abortSet('Error modifying lines from .edf file')
      return 1
   edfLines = newLines
   edfChanged = True


# See if cutviews.info exists and has views matching current exclude list
# Returns 0 if they do not match, 1 if they do, and -1 for read error
def checkExcludedViewsRemoved(setroot, excludeList):

   # Get all possible files and read the last one
   cutList1 = glob.glob(setroot + '_cutviews[0-9].info')
   cutList2 = glob.glob(setroot + '_cutviews[1-9][0-9].info')
   cutList1.sort()
   cutList2.sort()
   cutList1 += cutList2
   if not cutList1:
      return 0
   cutLines = readTextFileReportErr(cutList1[-1])
   if not cutLines:
      return -1

   # Get the list from the file and see if it matches the exclude list
   cutSplit = cutLines[-1].replace(',', ' ').split()
   excludeSplit = excludeList.replace(',', ' ').split()
   if len(cutSplit) != len(excludeSplit):
      return 0
   for cut in cutSplit:
      if cut not in excludeSplit:
         return 0
   return 1
   

# Print out errors from directive file reading or validation
def printDirectiveErrors(errors):
   if errors:
      prnLog('ERROR: Incorrect directive(s) as listed below:')
      for l in errors:
         prnLog(l)
      prnLog('')
      abortSet('Bad directives')
   return len(errors)
   

# Read a directive or template file and convert to a dictionary, skipping comment
# lines and lines without an =
def readDirectiveOrTemplate(filename, index):
   global absDirectiveFile, excludedViewsA, excludedViewsB, userTemplateDir
   global dsetDirTranslateInd
   errRet = (1, [], 0, False, '')
   directLines = readTextFile(filename, 'directive/template file', True)
   if isinstance(directLines, str):
      abortSet('Error ' + directLines)
      return errRet

   # For the batch file, find existing lines of various things
   if index == batDictInd:
      rewriteBatch = False
      nodi = (-1, '')
      skipViewsA = copyPrefix + 'skip'
      skipViewsB = copyPrefix + 'bskip'
      lineNumDict = {rootNameText : nodi, dataDirText : nodi, scopeTmplText : nodi,
                     sysTmplText : nodi, userTmplText : nodi, skipViewsA : nodi,
                     skipViewsB : nodi, cpTomoExtText : nodi}

      # Add raw boundary model variations if doing delivery
      if doDelivery:
         for operation in (patchTrackText, autoSeedText):
            for axlet in ('.a.', '.b.', '.any.'):
               direc = operation + axlet + 'rawBoundaryModel'
               lineNumDict[direc] = nodi

      # Find lines
      for ind in range(len(directLines)):
         line = directLines[ind].lstrip()
         lsplit = line.split('=')
         if len(lsplit) < 2:
            continue
         lineDirec = lsplit[0].strip()
         if lineDirec in lineNumDict:
            lineNumDict[lineDirec] = (ind, lsplit[1].strip())
         

      rootname = lineNumDict[rootNameText][1]
      datadir = lineNumDict[dataDirText][1]
      
      # Get true root name and fix or add line
      if numRootOpts:
         rootname = rootByOption[dfileInd]
         rootDirect = copyPrefix + 'name = ' + rootname
         rewriteBatch = True
         if lineNumDict[rootNameText][0] >= 0:
            directLines[lineNumDict[rootNameText][0]] = rootDirect
         else:
            directLines.append(rootDirect)

      # Take care of the data directory now if there are current directory entries
      if numCurrent and rootname:
         if not doDelivery:
            datadir = imodAbsPath(currentDirs[dfileInd])
            dsetDirTranslateInd = curDirTranslateInds[dfileInd]

         # If delivery, now we can make the directory
         # But we can't deliver file(s) until we know about dual/single
         else:
            if deliverDirs:
               delInd = min(dfileInd, len(deliverDirs) - 1)
               if not os.path.isdir(deliverDirs[delInd]):
                  abortSet('Cannot make directory for dataset ' + rootname + ' because '+\
                              deliverDirs[delInd] + ' is not an existing directory')
               topDir = deliverDirs[delInd]
               dsetDirTranslateInd = deliverTranslateInds[delInd]
            else:
               topDir = currentDirs[min(dfileInd, len(currentDirs) - 1)]
               dsetDirTranslateInd = curDirTranslateInds[min(dfileInd,
                                                            len(currentDirs) - 1)]

            datadir = imodAbsPath(os.path.join(topDir, rootname))

            # Detect Etomo error when it loses the original stack location and adjust
            # for it
            if startingStep > 0. and topDir.endswith(rootname) and \
               not (os.path.exists(datadir) and os.path.isdir(datadir)) and \
               os.path.exists(os.path.join(topDir, rootname + '.edf')):
               datadir = topDir
               warning(['CurrentLocation entry for ' + rootname + ' appears to ' + \
                       'incorrectly be the dataset directory', 
                        'instead of the original location, so it will ' +\
                        'be used as the current dataset location'])

            if not (os.path.exists(datadir) and os.path.isdir(datadir)):
               if os.path.exists(datadir):
                  abortSet(datadir + ' already exists and is not a directory')
                  return errRet
               try:
                  os.mkdir(datadir)
               except OSError:
                  abortSet('Error making directory for dataset, ' + datadir)
                  return errRet

               prnLog('Created dataset directory ' + datadir)

         dataDirect = dataDirText + ' = ' + datadir
         rewriteBatch = True
         if lineNumDict[dataDirText][0] >= 0:
            directLines[lineNumDict[dataDirText][0]] = dataDirect
         else:
            directLines.append(dataDirect)

      # Now check and resolve pathless templates
      for (tmplText, ind) in ((scopeTmplText, tmplDictInd), (sysTmplText, tmplDictInd+ 1),
                              (userTmplText, tmplDictInd + 2)):
         tmplName = lineNumDict[tmplText][1]
         if tmplName:
            (absTmplName, err, userTemplateDir, 
             errMess) = absTemplatePath(tmplName, ind - tmplDictInd, userTemplateDir, 
                                        direcFileErrorNames[ind])
            if err < 0: 
               abortSet(message)
               return errRet
            if err > 0:
               rewriteBatch = True
               directLines[lineNumDict[tmplText][0]] = tmplText + ' = ' + absTmplName

      # Deliver raw boundary model(s) and adjust directive to local file
      if doDelivery:
         for operation in (patchTrackText, autoSeedText):
            for axlet in ('.a.', '.b.', '.any.'):
               direc = operation + axlet + 'rawBoundaryModel'
               if lineNumDict[direc][0] >= 0:
                  srcPath = lineNumDict[direc][1]
                  rawFile = os.path.basename(srcPath)
                  destPath = os.path.join(datadir, rawFile)
                  srcExists = os.path.exists(srcPath)
                  destExists = os.path.exists(destPath)
                  if srcExists and not destExists:
                     try:
                        os.rename(srcPath, destPath)
                        prnLog('Moved boundary model ' + srcPath + ' to ' + datadir)
                     except OSError:
                        abortSet('Error renaming/moving raw boundary model ' + srcPath + \
                                    ' to ' + destPath)
                        return errRet

                  elif not (destExists and not srcExists):
                     abortSet('Raw boundary model ' + srcPath + ' does not exist')
                     return errRet

                  rewriteBatch = True
                  directLines[lineNumDict[direc][0]] = direc + ' = ' + rawFile

      # Check for excluded views to be removed
      excludedViewsA = ''
      excludedViewsB = ''
      checkRoot = os.path.join(datadir, rootname)
      if lineNumDict[skipViewsA][0] >= 0 and lineNumDict[skipViewsA][1]:

         # For first axis, check for no axis letter or a
         matchSingle = checkExcludedViewsRemoved(checkRoot, lineNumDict[skipViewsA][1])
         if matchSingle < 0:
            return errRet
         matchA = checkExcludedViewsRemoved(checkRoot + 'a', lineNumDict[skipViewsA][1])
         if matchA < 0:
            return errRet

         # If either matches, remove the skip entry
         if matchA or matchSingle:
            rewriteBatch = True
            directLines[lineNumDict[skipViewsA][0]] = skipViewsA + ' ='
            prnLog('Excluded views ' + lineNumDict[skipViewsA][1] +
                   ' were already removed from A/only axis')
         else:
            excludedViewsA = lineNumDict[skipViewsA][1]

      # Do the same for B
      if lineNumDict[skipViewsB][0] >= 0 and lineNumDict[skipViewsB][1]:
         matchB = checkExcludedViewsRemoved(checkRoot + 'b', lineNumDict[skipViewsB][1])
         if matchB < 0:
            return errRet
         if matchB:
            rewriteBatch = True
            directLines[lineNumDict[skipViewsB][0]] = skipViewsB + ' ='
            prnLog('Excluded views ' + lineNumDict[skipViewsB][1] +
                   ' were already removed from B axis')
         else:
            excludedViewsB = lineNumDict[skipViewsB][1]
               
   # Back to general processing of all kinds of files
   validErr = []
   for ind in range(len(directLines)):
      line = directLines[ind].lstrip()
      if line.startswith('#') or line == '':
         continue
      lsplit = line.split('=')
      if len(lsplit) < 2:
         validErr.append('Directive from ' + filename + ' lacks an = separator: ' + line)
         continue
      if lsplit[0].strip() == '':
         validErr.append('Directive from ' + filename + ' lacks a key: ' + line)
         continue
      allDirectives[index][lsplit[0].strip()] = (lsplit[1].strip(), ind)

   err = printDirectiveErrors(validErr)
   if not err and index == batDictInd:
      return (0, directLines, lineNumDict[cpTomoExtText], rewriteBatch, datadir)
   return (err, [], 0, False, '')
   

# Look for a directive with the given prefix and "option" (which might include a process)
# starting in the given directive dictionary, and converting by the type
# Returns None if the directive is not present, except for booleans where it returns 0
def lookupDirective(prefix, option, startDct, valType):
   bestDict = -1
   if prefix.startswith(comPrefix):
      keys = [prefix + '.' + option, prefix + 'a.' + option, prefix + 'b.' + option]
   elif prefix == setupPrefix:
      keys = [prefix + option, prefix + option, prefix + option]
   else:
      keys = [prefix + '.any.' + option, prefix + '.a.' + option, prefix + '.b.' + option]
   keyCheck = [(0, 1), (0, 1), (0, 2)]
   for dct in range(startDct, batDictInd + 1):
      for keyInd in keyCheck[axisNum]:
         if keys[keyInd] in allDirectives[dct]:
            better = True
            newInd = allDirectives[dct][keys[keyInd]][1]

            # If two entries are equivalent with regard to axis preference, new one is
            # better if it comes from later dictionary or was later in file
            equivBetter = dct > bestDict or newInd > bestInd
            if bestDict >= 0:

               # For dual axis, new one is better if it matches the current axis and
               # previous one did not; or if they are for same axis and this one is later
               if dualAxis:
                  better = (keyInd == axisNum and bestKeyInd != axisNum) or \
                      (keyInd == bestKeyInd and equivBetter)
               else:

                  # For single axis, all "any" and "a" entries are equivalent
                  better = equivBetter

            if better:
               bestDict = dct
               bestInd = newInd
               bestKeyInd = keyInd

   # If nothing was found, return None, or 0 for a boolean
   if bestDict < 0:
      if valType == BOOL_VALUE:
         return 0
      return None

   # Otherwise return 1 for a boolean only if it is specifically 1, or return the
   # converted value, or the string
   value = allDirectives[bestDict][keys[bestKeyInd]][0]
   if valType == BOOL_VALUE:
      if value == '1':
         return 1
      return 0
   elif valType == INT_VALUE or valType == FLOAT_VALUE:
      if value == '':
         return None
      try:
         if valType == INT_VALUE:
            numval = int(value)
         else:
            numval = float(value)
         return numval
      except ValueError:
         return 'ERROR'
   else:
      return value


# Get all the directives for making a com file after the fact
def laterComDirectives(startInd):
   lines = []
   if startInd < 1 and len(allDirectives[0]):
      lines.append('ChangeParametersFile ' + defaultsFile)
   if startInd <= tmplDictInd and scopeTmplText in allDirectives[batDictInd]:
      lines.append('ChangeParametersFile ' + allDirectives[batDictInd][scopeTmplText][0])
   if startInd <= tmplDictInd + 1 and sysTmplText in allDirectives[batDictInd]:
      lines.append('ChangeParametersFile ' + allDirectives[batDictInd][sysTmplText][0])
   if startInd <= tmplDictInd + 2 and userTmplText in allDirectives[batDictInd]:
      lines.append('ChangeParametersFile ' + allDirectives[batDictInd][userTmplText][0])
   lines.append('ChangeParametersFile ' + absDirectiveFile)
   return lines
   

# "Use" a file to replace one in sequence, with options to save previous one as _orig
# or to make a backup file out of it
def useFileAsReplacement(useFile, oldFile, saveOrig, makeBackup):
   (base, ext) = os.path.splitext(oldFile)
   origname = base + '_orig' + ext
   try:
      if saveOrig and not os.path.exists(origname):
            err = oldFile + ' to ' + origname
            os.rename(oldFile, origname)
      elif makeBackup:
         makeBackupFile(oldFile)
      else:
         cleanupFiles([oldFile])
      err = useFile + ' to ' + oldFile
      os.rename(useFile, oldFile)
   except OSError:
      abortSet('Error renaming ' + err + ' : ' + str(sys.exc_info()[1]))
      return 1


# Run imodtrans to get a boundary model onto (prealigned) stack
# In case there was any screwup in pixel size change of boundary model, make sure it
# fits the current stack before transforming
def transformRawBoundaryModel(modelIn, modelOut):
   try:
      imfile = datasetFilename('.preali')
      (panx, pany, panz) = getmrcsize(imfile)
      comstr = fmtstr('imodtrans -I "{}{}" -i "{}" -2 "{}.prexg" -S {} -tx {}' +\
                         ' -ty {} "{}" "{}"', dataName, stackExtension, imfile, dataName,
                      1. / coarseBinning, (panx - rawXsize // coarseBinning) / 2.,
                      (pany - rawYsize // coarseBinning) / 2., modelIn, modelOut)
      prnLog('Transforming ' + modelIn + ' to ' + modelOut + ' with:\n' + comstr)
      runcmd(comstr)
      return 0
   except ImodpyError:
      reportImodError('Could not transform boundary model to match stack')
      return 1


# Responds to the result of checking the check file; issues message and exits or sets
# finish flag as appropriate
def processQuitAction(action, message):
   global finishSetAndQuit
   if not action:
      return 
   if action == 'Q':
      if not message:
         message = 'RECEIVED SIGNAL TO QUIT, JUST EXITING'
      prnLog(message + '   [brt5]')
      if parallelRoot:
         sys.exit(1)
      else:
         sys.exit(0)
   if action == 'F':
      finishSetAndQuit = True

      
# Check for Q and F in the check file and process the result
def checkForQuit():
   action = checkForProChunksQuit(checkFile, proChunkCheckFile)
   processQuitAction(action, '')
   

# Runs a com file or com chunks using processchunks
def runOneProcess(comfile, single = True, usingGPU = False, message = '',
                  useMostCPUs = False):
   startTime = time.time()
   if (single or ('sirt' not in comfile and 'ctf3d' not in comfile)) and not \
      os.path.exists(comfile):
      abortSet('Command file ' + comfile + ' does not exist')
      return 1
   
   # Check for quitting then compose the rest of the command array
   checkForQuit()
   comArray = ['-n', str(niceness)]
   outfile = 'processchunks' + axisLet + '.out'
   if remoteDataDir:
      comArray += ['-w', remoteDataDir]
   machines = cpuList
   (comroot, ext) = os.path.splitext(comfile)
   if message:
      mess = message + ' (running ' + comfile
   else:
      mess = 'Running ' + comfile
      
   # single case no GPU:  regular run if no queue
   #                      queue run if not running on queue already
   #                      regular run if queue but running on queue already
   # Single case GPU:     simple run if no queue with gpu list
   #                      use gpuQueueCommand instead, no G option needed
   # chunk case, no GPU:  cpuList if no queue
   #                      queue run if queue
   # chunk case, GPU:     regular run if no queue with -G and current GPU list
   #                      if gpusPerClusterJob and no gpuQueueCommand, regular run with -G 
   #                      and GPU list made of localhost entries
   #                      if gpuQueueCommand, run on that

   if single:

      # Figure out the machine to use for single and its thread limit if possible
      threads = firstCpuLimit
      if useMostCPUs:
         machines = topCPUmachine
         threads = topCpuLimit
      elif not useFirstCPUforSingle:
         machines = '1'
         threads = localCpuLimit
      comArray += ['-s', '-e', '1']
      if threads > 0:
         comArray += ['-O', str(threads)]
      comuse = comfile
   else:
      mess += ' in multiple chunks'
      comuse = comroot

   if queueCommand and not usingGPU and not (single and runningOnQueue):
      machines = queueCommand
      comArray += ['-q', str(maxQueueJobs)]

   if usingGPU:
      if gpuQueueCommand:
         machines = gpuQueueCommand
         comArray += ['-q', str(maxGpuQueueJobs)]
      else:
         machines = myGPUlist
         if machines != '1':
            comArray.append('-G')
      mess += ' using GPU'

   if message:
      mess += ')'
   mess += '   [brt2]'
   if not machines:
      machines = '1'
   comArray += [machines, comuse]

   # Run the process detached
   prnLog(mess, flush=True)
   os.environ['PIP_PRINT_ENTRIES'] = '1'
   (error, finished, topQuit, numDone, mess) = \
       runProcesschunks(comArray, outfile, checkFile, proChunkCheckFile)
   os.environ['PIP_PRINT_ENTRIES'] = '0'
   if finished == -1:
      prnLog(mess)

   processQuitAction(topQuit, mess)
   if error < 0:
      prnLog('ERROR: ' + mess)
      abortSet('Cannot start processchunks to run ' + comfile)
      return 1
   if error:
      abortSet(mess)
      return 1

   doPrint = True
   if finished == 1:
      for comskip in handlingMessages:
         if comskip in comroot:
            doPrint = False
            break
      
   # Get error and warnings from logs
   if single and doPrint:
      printTaggedMessages(comroot + '.log', standardTags)
   else:
      printTaggedMessages(outfile, [('WARNING:', 0)])
   
   # After loop, one last check for quit and some more set aborts
   checkForQuit()
   if finished == -1:
      abortSet('An error occurred running ' + comfile)
   elif finished == -2:
      abortSet('Strangely, processchunks quit but Q was not detected in the check file')
   elif finished == 1:
      (minutes, seconds, frac) = elapsedTimeComponents(startTime);
      prnLog(fmtstr('Successfully finished {}   in {:02d}:{:02d}.{}   [brt3]\n', comfile,
                    minutes, seconds, frac), flush=True)
   return finished < 0


# If multiple parallel BRT jobs are running with GPUs, try to get an allocation to use
# in the coming job
def manageGPUallocation():
   global myParallelGPU, myGPUlist, needGPUrelease
   checkInterval = 5.
   warnInterval = 60. * 5.
   myParallelGPU = parallelGPU
   myGPUlist = gpuList
   if not useGPU or not maxParallelGPUs or queueCommand or coresPerClusterJob:
      return 0

   # Loop forever!
   startTime = datetime.datetime.now()
   lastWarn = startTime
   while True:
      try:
         allocLines = runcmd(fmtstr('gpuallocator -root {} -full "{}" -comm "{}" ' + \
                                    '-max {} -con {} -pid {}', parallelRoot, gpuList, 
                                    startingDir, maxParallelGPUs, hostname, myPID))
         myParallelGPU = 0
         myGPUlist = ''
         if not allocLines:
            abortSet('Gpuallocator gave no output, cannot proceed')
            return 1

         # Got an allocation: set variables for using them appropriately including flag
         # that they need to be released
         needGPUrelease = True
         for line in allocLines:
            myParallelGPU += 1
            if myGPUlist:
               myGPUlist += ','
            myGPUlist += line.rstrip('\r\n')
         
         prnLog(fmtstr('Received an allocation of {} GPUs in {} seconds', len(allocLines),
                       int(round((datetime.datetime.now() - startTime).seconds))))
         return 0

      except ImodpyError:

         # Got an error: See if it No GPUs and wait, or report error and abort
         errStrings = getErrStrings()
         if not errStrings:
            abortSet('Gpuallocator gave an error with no message, cannot proceed')
            return 1
         for line in errStrings:
            if '[GPA1]' in line:
               nowt = datetime.datetime.now()
               if (nowt - lastWarn).seconds > warnInterval:
                  warning(fmtstr('Have not been able to get a GPU allocation for {} ' + \
                                 'minutes', int(round((nowt - startTime).seconds / 60.))))
                  lastWarn = nowt
               break
         else:   # ELSE ON FOR
            for line in errStrings:
               prnLog(line)
            abortSet('Gpuallocator exited with an error, cannot proceed')
            return 1

         time.sleep(checkInterval)
         continue
               

# Release an allocation of GPUs if one was obtaied
def releaseGPUallocation():
   global needGPUrelease
   if not needGPUrelease:
      return
   try:
      runcmd(fmtstr('gpuallocator -root {} -comm "{}" -con {} -pid {}', parallelRoot, 
                    startingDir, hostname, myPID))
   except ImodpyError:
      errStrings = getErrStrings()
      mess = ': '
      if errStrings:
         mess += errStrings[-1]
      warning('Error trying to release GPU allocation' + mess)
   needGPUrelease = False


# Get the number of cores or GPUs to be used for a process that can use either
def numberOfProcessingUnits():
   if manageGPUallocation():
      return -1
   numProc = 1
   if myParallelGPU > 1:
      numProc = myParallelGPU
   elif queueCommand and not gpuQueueCommand:
      numProc = maxQueueJobs;
   elif parallelCPU and not useGPU:
      numProc = parallelCPU
   return numProc


# Give lines for running makecomfile, add the output file entry, get  the com file, and
# run it
def makeAndRunOneCom(comlines, comfile, message = '', useGPU = False, skipRun = False):
   comlines.insert(0, 'OutputFile	' + comfile)
   comlines.append('NamingStyle	' + str(nameStyle))
   comlines.append('StackExtension	' + stackExtension[1:])
   try:
      runcmd('makecomfile -StandardInput', comlines)
   except ImodpyError:
      reportImodError('Error making ' + comfile)
      return 1

   if skipRun:
      return 0;

   if runOneProcess(comfile, usingGPU = useGPU, message = message):
      return 1
   return 0


# Use the sed commands to modify a command file, write it, and run it
# Reads the command file first unless lines are supplied in inLines
def modifyWriteAndRunCom(comfile, sedcom, inLines = None, message = '', skipRun = False):
   if not inLines:
      inLines = readTextFileReportErr(comfile)
      if not inLines:
         return 1
               
   if pysed(sedcom, inLines, comfile, retErr=True):
      abortSet('Error modifying ' + comfile)
      return 1
   if skipRun:
      return 0
   if runOneProcess(comfile, message = message):
      return 1
   return 0
   

# Finds and converts a single value after a token like '=' on the given line
# Throws ValueError if token doesn't exist or value has conversion error
def getOneValueAfterToken(line, token, valType):
   ind = line.find(token) + 1
   if ind < 1:
      junk = int('=')
   if valType == INT_VALUE:
      return int(line[ind:])
   else:
      return float(line[ind:])


# Analyze tiltalign log file for angle, shift, and thickness
def analyzeAlignLog(twoSurf, angleArr, noLogOK = False):

   # Fraction of shift to subtract from fiducial-based thickness to get reconstruction
   # thickness, and minimum on each side before doing that
   thickShiftFrac = 1.
   minFidAdjustThick = 4
   angleArr[0:5] = 5 * [0.]
   if noLogOK and not os.path.exists('align' + axisLet + '.log'):
      return 0
   
   try:
      loglines = runcmd('alignlog -a align' + axisLet + '.log')
   except ImodpyError:
      reportImodError('Extracting angle analysis from align log')
      return 1
   
   tags = ['Total tilt angle change', 'X axis tilt needed', 'Unbinned thickness',
           'Incremental unbinned shift', 'Total unbinned shift']
   numBot, numTop = -1, -1
   try:
      for line in loglines:
         if '# of points' in line:

            # There can be one or three of these entries, so the second one replaces
            # numBot and the third one gives numTop
            num = getOneValueAfterToken(line, '=', INT_VALUE)
            if numBot < 0:
               numBot = num
            elif numTop < 0:
               numBot = num
               numTop = 0
            else:
               numTop = num
         for ind in range(len(tags)):
            if tags[ind] in line:
               angleArr[ind] = getOneValueAfterToken(line, '=', FLOAT_VALUE)
   except ValueError:
      abortSet('Error extracting information from align log')
      return 1

   # Adjust the thickness by the X-tilt if no tilt is to be used: the spacing itself
   # is increased by cosine of the angle, and the pitch adds tangent times extent in Y
   if noXaxisTilt:
      angleArr[2] = angleArr[2] / math.cos(math.radians(angleArr[1])) + \
                    0.9 * rawYsize * math.tan(math.fabs(math.radians(angleArr[1])))

   # For two surfaces, also compute a reconstruction thickness (thickness are floats here)
   if twoSurf:
      angleArr[5] = angleArr[2]
      if numBot >= minFidAdjustThick and numTop >= minFidAdjustThick and not centerOnGold:
         angleArr[5] -= thickShiftFrac * math.fabs(angleArr[4])
   angleArr[6] = numBot
   angleArr[7] = max(0, numTop)

   return 0


# Determine parameter modifications for a montage
# Inputs are the raw X and Y size of the montage and the desired aligned stack size,
# and whether sizes are being transposed for a rotation close to 90 degrees
# Values put into frameArr are:
# 0,1: actual aligned stack X and Y size
# 2,3: starting X and Y coordinates for blend
# 4,5: ending X and Y coordinates for blend
# 6,7: actual blended raw/prealigned stack size = FULLIMAGE
# 8,9: SUBSETSTART X and Y
# This function can be moved to imodpy so etomo can access it from a script
def montageFrameValues(nxmont, nymont, transpose, nxali, nyali, frameArr):

   (nxOut, nyOut) = (nxmont, nymont)
   if transpose:
      (nxOut, nyOut) = (nymont, nxmont)
   
   # aligned stacksize
   (frameArr[0], frameArr[1]) = runGoodframe(nxali, nyali)
   (frameArr[6], frameArr[7]) = runGoodframe(nxOut, nyOut)
   if frameArr[0] < 0 or frameArr[6] < 0:
      return min(frameArr[0], frameArr[6])

   # starting coordinates X and Y for blend, centered on center in old orientation
   frameArr[2] = -((frameArr[0] - nxmont) // 2)
   frameArr[3] = -((frameArr[1] - nymont) // 2)
   
   # ending for blend, and SUBSETSTART
   for ind in (0, 1):
      frameArr[ind + 4] = frameArr[ind + 2] + frameArr[ind] - 1
      frameArr[ind + 8] = -((frameArr[ind] - frameArr[ind + 6]) // 2)

   return 0


# Look up directives for both a raw and prealign boundary model with the given
# prefix, if there is a raw one, transform it; of there is a prealigned one use it
# directly for A; for B axis look for a specific one then try to transform the A
# axis one with eitherthe transferfid output or a run of transferfid when patch tracking
def getOrXformBoundaryModel(prefixText, dfltModelSuffix, doTransfer):

   # Boundary model has to be transformed if indicated
   # A model on A raw stack has to be transformed to preali
   aRawBound = lookupDirective(prefixText, 'rawBoundaryModel', 0, STRING_VALUE)
   aBoundaryMod = lookupDirective(prefixText, 'prealiBoundaryModel', 0, STRING_VALUE)
   if not axisInd:
      rawBound = aRawBound
      boundaryMod = aBoundaryMod
   else:

      # If there is an A raw model and not an A preali model, set default name
      # and then set the variable that it actually exists
      boundaryMod = ''
      rawBound = ''
      if aBoundaryMod == None and aRawBound != None:
         aBoundaryMod = setName + 'a' + dfltModelSuffix
      haveAaxisBound = aBoundaryMod != None and os.path.exists(aBoundaryMod)

      # A raw model for B has to be specific to the B axis
      key = prefixText + '.b.rawBoundaryModel'
      if key in allDirectives[batDictInd]:
         rawBound = allDirectives[batDictInd][key][0]
      key = prefixText + '.b.prealiBoundaryModel'
      if key in allDirectives[batDictInd]:
         boundaryMod = allDirectives[batDictInd][key][0]
      if boundaryMod:
         prnLog('GOT B MODEL')

      # But if there is not a raw model for B, see if there was one for A
      # that can be transferred with the transform
      if not rawBound and not boundaryMod and haveAaxisBound and doTransfer:
         boundaryMod = dataName + dfltModelSuffix
         if patchTrack:

            # Patch tracking: run transferfid to find best matching views and transform
            comlines = ['RootNameOfDataFiles       ' + setName]
            comlines += laterComDirectives(0)
            comlines += [fmtstr('OneParameterChange {}transferfid.transferfid.' +\
                                   'BoundaryModel={}', comPrefix, aBoundaryMod),
                         fmtstr('OneParameterChange {}transferfid.transferfid.' +\
                                   'SeedModel={}', comPrefix, boundaryMod),
                         fmtstr('OneParameterChange {}transferfid.transferfid.' +\
                                'CorrespondingCoordFile=', comPrefix)]
            if makeAndRunOneCom(comlines, 'transferfid' + comExt,
                                'Transferring boundary model from A to B axis'):
               return (1, None)
                         
         else:

            # Fiducials: try to get the change in Z from the transfer log
            deltaZ = 0
            loglines = readTextFile('transferfid.log', None, True)
            if not isinstance(loglines, str):
               fromView = 0
               toView = 0
               for line in loglines:
                  if 'from view' in line and 'to view' in line:
                     lsplit = line.split()
                     try:
                        for ind in range(len(lsplit) - 2):
                           if lsplit[ind] == 'from' and lsplit[ind + 1] == 'view':
                              fromView = int(lsplit[ind + 2])
                           if lsplit[ind] == 'to' and lsplit[ind + 1] == 'view':
                              toView = int(lsplit[ind + 2])
                        deltaZ = toView - fromView
                     except Exception:
                        pass

            # Then run imodtrans to transform model
            try:
               (panx, pany, panz) = getmrcsize(datasetFilename('.preali'))
               comstr = fmtstr('imodtrans -2 "{}_AtoB.xf" -l 0 -n {},{},{} -tz {} ' +
                               '"{}" "{}"', setName, panx, pany, panz, deltaZ,
                               aBoundaryMod, boundaryMod)
               runcmd(comstr)
            except ImodpyError:
               reportImodError('Failed to transform boundary model from A to B')
               return (1, None)

   # Transform a raw model
   if rawBound and not boundaryMod:
      boundaryMod = dataName + dfltModelSuffix
      if transformRawBoundaryModel(rawBound, boundaryMod): 
         return (1, None)
   return (0, boundaryMod)


# Common place to fetch the name of com file and process for newstack or blendmont
def comAndProcessForAlignedStack(prefix):
   if ifMontage:
      comRoot = prefix + 'blend'
      process = 'blendmont'
   else:
      comRoot = prefix + 'newst'
      process = 'newstack'
   return (comRoot, process)


# Run a tilt comfile, splitting it if appropriate
def splitAndRunTilt(comfile, message, numProc = 0):
   if not numProc:
      numProc = numberOfProcessingUnits()
      if numProc < 1:
         return 1

   if numProc > 1:
      try:
         runcmd('splittilt -n ' + str(numProc) + ' ' + comfile)
      except ImodpyError:
         reportImodError('Error running splittilt on ' + comfile)
         releaseGPUallocation()
         return 1

   err = runOneProcess(comfile, numProc < 2, useGPU, message = message)
   releaseGPUallocation()
   return err


# Run findsection on a file, produce a tomopitch model or fill array with Z limits
# Returns -1 for an error from findsection, 1 for failures to write or read log or
# convert numbers
def runFindSection(filename, numScales, boxSize, pitchMod = None, topBots = None,
                   block = None, failureOK = False, comRoot = 'findsection_tmp'):
   global suppressAbort
   limitKeys = ['Median Z values', 'autopatchfit combine', 'Absolute limits']
   comLines = ['$findsection -StandardInput',
               'TomogramFile ' + filename,
               'NumberOfDefaultScales ' + str(numScales),
               fmtstr('SizeOfBoxesInXYZ {0},1,{0}', boxSize)]
   if pitchMod:
      comLines.append('TomoPitchModel ' + pitchMod)
      comLines.append('NumberOfSamples 5')
      if block:
         comLines.append('BlockSize ' + str(block))
   if writeTextFileReportErr(comRoot + comExt, comLines):
      return 1
   mess = 'Finding Z limits of the material in the tomogram'
   if pitchMod:
      mess = 'Getting a model for tomogram positioning'
   suppressAbort = failureOK
   err = runOneProcess(comRoot + comExt, message = mess)
   suppressAbort = False
   if err:
      return -1
   if topBots:
      loglines = readTextFileReportErr(comRoot + '.log')
      if not loglines:
         return 1
      for line in loglines:
         for ind in (0, 1, 2):
            if limitKeys[ind].lower() in line.lower():
               lsplit = line.split()
               try:
                  topBots[2 * ind] = int(lsplit[-2])
                  topBots[2 * ind + 1] = int(lsplit[-1])
               except Exception:
                  abortSet('Error converting Z limits in: ' + line)
                  return 1

   return 0


# Parse a log from tomopitch to get the thickness, angles, and Z shift
def parseTomopitchLog(angleArr):
   pitchFile = 'tomopitch' + axisLet + '.log'
   tags = ['x-tilted  lines', 'X axis tilt -', 'Angle offset -', 'Z shift -']
   original = [0, 0., 0., 0.]
   untilted = [0, 0., 0., 0.]
   gotAllLines = 0

   if posSampleType <= 0:
      return 1
   if not os.path.exists(pitchFile):
      return 1
   pitchLines = readTextFileReportErr(pitchFile)
   if not pitchLines:
      return -1

   for line in pitchLines:
      if 'ERROR: ' in line:
         return 2

      # Look for the standard tagged lines but also get original values
      for ind in range(4):
         if tags[ind] in line:
            lsplit = line.split()
            try:
               if ind:
                  angleArr[ind] = float(lsplit[-1])
                  if noXaxisTilt:
                     for word in range(len(lsplit)):
                        if 'Original:' in lsplit[word] and word + 1 < len(lsplit):
                           original[ind] = float(lsplit[word + 1])
               else:
                  angleArr[ind] = int(lsplit[-1])
            except ValueError:
               abortSet('Error converting tomopitch output of ' + \
                           tags[ind].replace(' -', '') + ' to a number')
               return -1

      # But if want no Xtilt, look for the line that starts that output,
      # Then process the specific lines after it with the untilted values
      if noXaxisTilt and 'all line pairs' in line:
         gotAllLines = 1
      elif gotAllLines:
         try:
            if gotAllLines == 2:
               ind = line.find('add')
               toInd = line.find('to', ind + 3)
               untilted[2] = float(line[ind + 3:toInd])
            if gotAllLines == 3:
               lsplit = line.split()
               untilted[0] = int(lsplit[-1])
               ind = line.find('shift of')
               toInd = line.find(';', ind)
               untilted[3] = float(line[ind + 8:toInd])
            gotAllLines += 1
         except ValueError:
            abortSet('Error converting tomopitch output for non x-tilted lines ' + \
                           ' to a number')
            return -1
               
   # Add original and no X tilt results (returning original X tilt is consistent with
   # the other returned values)
   if noXaxisTilt and gotAllLines:
      for ind in range(4):
         angleArr[ind] = original[ind] + untilted[ind]

   return 0


# Look up whether SIRT or fake SIRT or one of those plus regular BP are to be done
def getReconTypes():
   doSIRT = lookupDirective(runtimePrefix + 'Reconstruction', 'useSirt', 0, BOOL_VALUE)
   doBPalso = lookupDirective(runtimePrefix + 'Reconstruction', 'doBackprojAlso', 0,
                              BOOL_VALUE)
   doRegBPalso = lookupDirective(runtimePrefix + 'Reconstruction',
                                 'doRegularBPalsoIfFakeSIRT', 0, BOOL_VALUE)
   doFakeSIRT = lookupDirective(comPrefix + 'tilt', 'tilt.FakeSIRTiterations', 0,
                                STRING_VALUE)
   doBoth = (doSIRT and doBPalso) or (doFakeSIRT and (doRegBPalso or doBPalso))
   return (doSIRT, doFakeSIRT, doBoth)


# Return flags for SIRT, fake SIRT, and doing both that and regular BP
# Compose the name of the SIRT rec from the leave iterations directive.  Return a name 
# of None if SIRT is not being done or an integer 1 if it is but the name can't be made
# Return a name for the fake SIRT recon, or None if not being done
def getSIRTrecName(recRoot):
   (doSIRT, doFakeSIRT, doBoth) = getReconTypes()
   
   recName = None
   fakeName = None
   if doSIRT:
      leaveList = lookupDirective(comPrefix + 'sirtsetup', 'sirtsetup.LeaveIterations', 0,
                                  STRING_VALUE)
      if leaveList:
         lsplit = leaveList.replace('-',',').split(',')
         lastOne = lsplit[-1:][0]
         if len(lastOne) < 2:
            lastOne = '0' + lastOne
         recName =  datasetFilename('.srec' + lastOne, root = recRoot)
      else:
         recName = 1

   if doFakeSIRT:
      fakeName = datasetFilename('.rec', root = recRoot)
      if doBoth:
         fakeName = datasetFilename('.slfrec', root = recRoot)

   return (doSIRT, doFakeSIRT, doBoth, recName, fakeName)


# Run clip stats on raw or fixed stack
def runClipStats(command, stackSuffix, descrip):
   comfile = dataName + stackSuffix + '.st.stats' + comExt
   if writeTextFileReportErr(comfile, [command]):
      return None
   err = runOneProcess(comfile, message = 'Getting statistics for ' + descrip + ' stack')
   cleanupFiles([comfile])
   if err:
      return None

   # Extract the summary lines
   clipLines = readTextFileReportErr(dataName + stackSuffix + '.st.stats.log')
   return clipLines


# Test whether a step should be run
def needStep(step):
   return step >= startingStep - 0.005 and step <= endingStep + 0.005


# Report a step was reached if it is in the range to run and in the list to report
def reportReachedStep(step):
   reachPoints = (0, 5, 6, 7, detect3DstepNum, 13, 14)
   if step in reachPoints and needStep(step):
      prnLog('Reached step ' + str(step))
      prnLog('', flush = True)


# Run a command entered as a directive either instead of or after a step
# Returns 1 for error, 0 if no action or -1 for success when replacing a step,
# or 0 for no action OR success when running after a step
def replaceOrRunAfterStep(step, runAfter):
   if not needStep(step):
      return 0
   direc = 'eplaceStep'
   goodRet = -1
   if runAfter:
      direc = 'unAfterStep'
      goodRet = 0
   comfile = 'r' + direc + str(step) + axisCom
   command = lookupDirective(runtimePrefix + 'R' + direc, str(step), 0, STRING_VALUE)
   if not command:
      return 0
   command = command.replace('%{setname}', dataName)
   if writeTextFileReportErr(comfile, ['$' + command]):
      return 1
   if runOneProcess(comfile):
      return 1
   return goodRet


# SINGLE-CALL FUNCTIONS FOR INITIAL STEPS

# Read in a validation file as a csv and store the directives in dictionaries
def processValidationFile():
   comPref = comPrefix[:-1]
   runPref = runtimePrefix[:-1]
   try:
      message = 'Opening '
      csvfile = open(validateFile, 'r')
      message = 'Using csv reader on '
      reader = csv.reader(csvfile)
      validLines = []
      for row in reader:
         validLines.append(row)
   except Exception:
      exitError(message + validateFile)

   for line in validLines:
      if len(line) > 1 and len(line[1]) > 0:
         lsplit = line[0].split('.')
         if len(lsplit) < 2:
            continue
         batchOK = len(line) > 3 and line[3].strip() == 'Y'
         templateOK = len(line) > 4 and line[4].strip() == 'Y'
         isBool = len(line) > 2 and line[2].lower().strip() == 'bool'
         if lsplit[0] == comPref:
            if len(lsplit) < 4:
               continue
            combase = lsplit[1]
            if combase not in baseComDict:
               baseComDict[combase] = combase
               baseComDict[combase + 'a'] = combase
               baseComDict[combase + 'b'] = combase
            key = lsplit[1] + '.' + lsplit[2] + '.' + lsplit[3]
            validComDict[key.lower()] = (key, templateOK, batchOK, isBool)

         elif lsplit[0] == runPref:
            if len(lsplit) < 4:
               continue
            key = lsplit[1] + '.' + lsplit[3]
            validRunDict[key.lower()] = (key, templateOK, batchOK, isBool)

         else:
            validOtherDict[line[0].lower()] = (line[0], templateOK, batchOK, isBool)

            
# Look for every directive in the lists of valid ones from master file
def checkAllDirectives(mainFile = 'directive'):
   errors = []
   source = direcFileErrorNames[0:4] + [mainFile]
   types = ['template', 'batch directive']
   comPref = comPrefix[:-1]
   runPref = runtimePrefix[:-1]
   for ind in range(batDictInd + 1):
      dirType = ind // batDictInd
      if mainFile == 'template':
         dirType = 0
      for direc in allDirectives[ind]:
         messfrom = 'irective from ' + source[ind] + ' file'
         dsplit = direc.split('.')
         badCase = False
         OKforFile = True
         emptyBool = False
         valNotBool = allDirectives[ind][direc][0] != '0' and \
             allDirectives[ind][direc][0] != '1'
         if len(dsplit) < 2 or ((dsplit[0] == comPref or dsplit[0] == runPref) \
                                   and len(dsplit) < 4):
            errors.append('D' + messfrom + ' too short: ' + direc)

         # comparam first, start by checking the com file is in list
         elif dsplit[0] == comPref:
            comlow = dsplit[1].lower()
            if comlow not in baseComDict:
               errors.append('D' + messfrom + ' does not include a known com file: '
                             + direc)
            else:

               # Check for a match other than incorrect case
               combase = baseComDict[comlow]
               key = dsplit[1] + '.' + dsplit[2] + '.' + dsplit[3]
               if key.lower() in validComDict:
                  badCase = key != validComDict[key.lower()][0]
                  OKforFile = validComDict[key.lower()][dirType + 1]
                  emptyBool = validComDict[key.lower()][3] and valNotBool

               else:

                  # Look for a process match next
                  proclow = dsplit[2].lower()
                  for vkey in validComDict:
                     vsplit = vkey.split('.')
                     if vsplit[0] == combase and vsplit[1].lower() == proclow:
                        badCase =  dsplit[1] not in baseComDict or vsplit[1] != dsplit[2]
                        
                        # If the process matches, now try to find and read autodoc
                        optFile = PipOpenInstalledAdoc(vsplit[1])
                        errmess = ''
                        if optFile:
                           adocLines = readTextFile(optFile, None, True)
                           if isinstance(adocLines, str):
                              errmess = 'error reading it'
                        else:
                           errmess = 'error opening it at standard location'

                        # If no autodoc, issue a warning because we just can't tell
                        if errmess:
                           warning(['Unknown d' + messfrom + ': ' + direc,
                                    '  Could not check ' + vsplit[1] +
                                    '.adoc because of '+ errmess])
                        else:

                           # Otherwise look for a case-insensitive match and report a
                           # bad case if any, otherwise report error if no match
                           target = fmtstr(r'\[ *Field *= *{} *\]', dsplit[3])
                           optMatch = re.compile(target)
                           optLC = re.compile(target, re.IGNORECASE)
                           for line in adocLines:
                              if re.match(optLC, line):
                                 if not re.match(optMatch, line):
                                    badCase = True
                                 break
                           else:  # ELSE ON FOR
                              errors.append('Unknown d' + messfrom + ': ' + direc)
                           
                        break
                  else:   # ELSE ON FOR
                     errors.append('D' + messfrom + ' does not include a known process: '
                                   + direc)

         # runtime next, check for a/b/any; all else must match
         elif dsplit[0] == runPref:
            if dsplit[2] not in ['any', 'a', 'b']:
               errors.append('D' + messfrom + ' does not include a/b/any: ' + direc)
            else:
               key = dsplit[1] + '.' + dsplit[3]
               if key.lower() in validRunDict:
                  badCase = key != validRunDict[key.lower()][0]
                  OKforFile = validRunDict[key.lower()][dirType + 1]
                  emptyBool = validRunDict[key.lower()][3] and valNotBool
               else:
                  errors.append('Unknown d' + messfrom + ': ' + direc)

         # Any other directives must match entirely
         else:
            if direc.lower() in validOtherDict:
               badCase = direc != validOtherDict[direc.lower()][0]
               OKforFile = validOtherDict[direc.lower()][dirType + 1]
               emptyBool = validOtherDict[direc.lower()][3] and valNotBool
            else:
               errors.append('Unknown d' + messfrom + ': ' + direc)
            
         if badCase:
            errors.append('D' + messfrom + ' has incorrect case: ' + direc)
         if not OKforFile:
            errors.append('D' + messfrom + ' is not intended for use in a ' + \
                             types[dirType] + ' file: ' + direc)
         if emptyBool:
            errors.append('D' + messfrom + ' is a boolean and must be 0 or 1: ' + direc)

   return printDirectiveErrors(errors)


# Look through directives for templates and basic setup parameters
def scanSetupDirectives():
   global datasetDir, dualAxis, setName, scanHeader, ifMontage, defocus, pixelSize
   global fidSizeNm, fidSizePix
   source = direcFileErrorNames

   batDirec = allDirectives[batDictInd]
   ifMontage = False
   dualAxis = False
   scanHeader = False
   pixelSize, fidSizeNm, fidSizePix = 0., None, 0.
   datasetDir, setName = '', ''
   defocus = -1000000.
   if scopeTmplText in batDirec:
      err = readDirectiveOrTemplate(batDirec[scopeTmplText][0], tmplDictInd)
      if err[0]:
         return 1
   if sysTmplText in batDirec:
      err = readDirectiveOrTemplate(batDirec[sysTmplText][0], tmplDictInd + 1)
      if err[0]:
         return 1
   if userTmplText in batDirec:
      err = readDirectiveOrTemplate(batDirec[userTmplText][0], tmplDictInd + 2)
      if err[0]:
         return 1

   if copyPrefix + 'name' in batDirec:
      setName = batDirec[copyPrefix + 'name'][0]

   # This was placed into the dictionary in readDirectiveOrTemplate based on entries
   # for currentDir, deliveryDir, and makeSubDir
   if dataDirText in batDirec:
      datasetDir = batDirec[dataDirText][0]
   for indDir in range(batDictInd + 1):
      direc = allDirectives[indDir]
      try:
         if scanHeadText in direc:
            scanHeader = direc[scanHeadText][0] == '1'
         if copyPrefix + 'montage' in direc:
            ifMontage = direc[copyPrefix + 'montage'][0] != '0'
         if copyPrefix + 'dual' in direc:
            dualAxis = direc[copyPrefix + 'dual'][0] != '0'
         dirKey = copyPrefix + 'pixel'
         if dirKey in direc:
            ftext = direc[dirKey][0]
            if ftext:
               pixelSize = float(ftext)
         dirKey = copyPrefix + 'gold'
         if dirKey in direc:
            ftext = direc[dirKey][0]
            if ftext:
               fidSizeNm = float(ftext)
         dirKey = copyPrefix + 'defocus'
         if dirKey in direc:
            ftext = direc[dirKey][0]
            if ftext:
               defocus = float(ftext)

      except ValueError:
         abortSet(fmtstr('Error converting the string "{}" for directive {} to float ' + \
                            'in {} file', ftext, dirKey, source[indDir]))
         return 1

   if validation > 0:
      return 0
   
   if (not pixelSize and not scanHeader):
      abortSet('Pixel size missing from directives, and header is not being scanned')
      return 1
   if not setName or not datasetDir:
      abortSet('Set name or dataset directory missing from directives')
      return 1

   return 0


# Look for probably old default values in adoc that are masking a template value
def checkDefaultsInBatchFile():
   numDuplicate = 0
   direcMasked = []
   origDefaults = (('setupset.scanHeader', '1'),
                   ('runtime.Preprocessing.any.archiveOriginal', '1'),
                   ('comparam.eraser.ccderaser.LineObjects', '2'),
                   ('comparam.eraser.ccderaser.BoundaryObjects', '3'),
                   ('comparam.eraser.ccderaser.AllSectionObjects', '1-3'),
                   ('comparam.track.beadtrack.RoundsOfTracking', '4'),
                   ('runtime.BeadTracking.any.numberOfRuns', '2'),
                   ('comparam.align.tiltalign.RobustFitting', '1'),
                   ('setupset.copyarg.gold', '0'))
   
   for (key, value) in origDefaults:
      if key in allDirectives[batDictInd] and value == allDirectives[batDictInd][key][0]:
         numDuplicate += 1
         for ind in range(tmplDictInd, batDictInd):
            if key in allDirectives[ind] and value != allDirectives[ind][key][0]:
               direcMasked.append('    ' + key + ' = ' + value)
               
   if numDuplicate >= 4 and len(direcMasked):
      mess = ['Incorrect directives in batch file may be masking template values.',
              str(numDuplicate) + ' directives in the batch file match the value of a ' +\
              'batch default.', 
              'These came from an old starting batch file that contained batch defaults',
              '(from before IMOD 4.10.15).  Directives overriding a template value are:']
      mess += direcMasked
      mess += ['Check these directives and remove any other unintended directives that',
               'occur in ' + defaultsFile + ' from your starting batch files.']
      warning(mess)


# Get some basic processing flags, prealign binning, and dataset set
def getAxisInitialParameters():
   global rawXsize, rawYsize, zsize, fiducialless, coarseBinning, patchTrack, eraseGold
   global correctCTF, ctf3dSlabThick, doBothRecons, rawFor3dCTF, skipAlignStackMods
   global filterIn2D, autofitCTF
   
   fiducialless = lookupDirective(runtimePrefix + 'Fiducials', 'fiducialless', 0,
                                  BOOL_VALUE)
   trackMethod = lookupDirective(runtimePrefix + 'Fiducials', 'trackingMethod', 0,
                                INT_VALUE)
   patchTrack = isinstance(trackMethod, int) and trackMethod == 1

   eraseGold = lookupDirective(runtimePrefix + 'AlignedStack', 'eraseGold', 0, INT_VALUE)
   if isinstance(eraseGold, str):
      abortSet('Error converting eraseGold entry to integer')
      return 1

   filterIn2D = lookupDirective(runtimePrefix + 'AlignedStack', 'filterStack', 0,
                                BOOL_VALUE)
   
   if eraseGold == 1 and (fiducialless or patchTrack):
      abortSet('Cannot erase gold with fiducials after fiducialless processing or ' +\
                  'patch tracking')
      return 1

   (doSIRT, doFakeSIRT, doBothRecons) = getReconTypes()

   # Check CTF parameters
   correctCTF = lookupDirective(runtimePrefix + 'AlignedStack', 'correctCTF', 0,
                                BOOL_VALUE)
   if correctCTF:
      autofitCTF = lookupDirective(runtimePrefix + 'CTFplotting', 'autoFitRangeAndStep', 
                                   0, STRING_VALUE)
   if correctCTF and defocus < -999999:
      if autofitCTF:
         if not lookupDirective(comPrefix + 'ctfplotter', 'ctfplotter.ScanDefocusRange',
                                0, STRING_VALUE):
            abortSet('Either defocus or a range to scan must be entered to use ' + \
                     'Ctfplotter')
            return 1
            
      else:
         abortSet('Defocus must be entered to correct CTF without using Ctfplotter')
         return 1

   ctf3dSlabThick = lookupDirective(comPrefix + 'ctf3dsetup', 
                                    'ctf3dsetup.SlabThicknessInNm', 0, INT_VALUE)
   if isinstance(ctf3dSlabThick, str):
      abortSet('Error converting ctf3dsetup.SlabThicknessInNm entry to integer')
      return 1
   if ctf3dSlabThick and not correctCTF:
      abortSet('Directive for correcting CTF must also be present to do 3-D ' + \
               'CTF-corrected reconstruction')
      return 1
   if ctf3dSlabThick and ((doSIRT or doFakeSIRT) and not doBothRecons):
      abortSet('You cannot do just a SIRT or SIRT-like reconstruction with 3-D CTF ' + \
               'correction')
      return 1

   rawFor3dCTF = lookupDirective(comPrefix + 'ctf3dsetup', 
                                 'ctf3dsetup.UseUnalignedImages', 0, BOOL_VALUE)
   skipAlignStackMods = ctf3dSlabThick and not doBothRecons
   superSample = lookupDirective(comPrefix + 'tilt', 'tilt.SuperSampleFactor', 0,
                                 INT_VALUE)
   if not isinstance(superSample, str) and superSample and superSample > 1 and \
      lookupDirective(comPrefix + 'tilt', 'tilt.ExpandInputLines', 0, BOOL_VALUE):
      abortSet('You cannot expand input lines when doing 3D CTF correction with ' + \
               'unaligned images')
      return 1
      
   (comRoot, process) = comAndProcessForAlignedStack('pre')
   coarseBinning = lookupDirective(comPrefix + comRoot, process + '.BinByFactor', 0,
                                   INT_VALUE)
   if isinstance(coarseBinning, str):
      abortSet('Error converting binning value in directive to integer')
      return 1
   if not coarseBinning:
      coarseBinning = 1
   
   try:
      if ifMontage:
         (rawXsize, rawYsize, zsize) = getMontageSize(dataName + stackExtension, 
                                                      dataName + '.pl')
      else:
         (rawXsize, rawYsize, zsize) = getmrcsize(dataName + stackExtension)
   except ImodpyError:
      reportImodError('Error getting size of ' + dataName + stackExtension)
      return 1

   return 0


# Run imodtrans on a model to make it match existing raw stack,
# Then replace original file.  A failure is not fatal unless it leaves no raw model file
def fixModelHeaderForImage(model, image):
   if not os.path.exists(image) or not os.path.exists(model):
      return 0
   tmpName = model + '.transtemp'
   try:
      runcmd(fmtstr('imodtrans -I "{}" "{}" "{}"', image, model, tmpName))
   except ImodpyError:
      warning('Could not fix coordinate information in model ' + model)
      reportImodError(None)
      return 0
   makeBackupFile(model)
   try:
      os.rename(tmpName, model)
   except OSError:
      prnLog('Error renaming fixed model ' + tmpName + ' to ' + model)
      return 1
   prnLog('Modified model ' + model + ' to adjust for change in pixel size of raw stack')
   return 0
            

# If pixel size was changed for an FEI file, make sure all raw boundary and X-ray models
# are changed so that they can still load correctly on the raw stacks
def fixAllSuppliedModelHeaders():
   modSet = set()
   bModSet = set()
   axList = ['a', 'any']
   if dualAxis:
      image = setName + 'a' + stackExtension
      modSet.add(setName + 'a.erase')
      axList.append('b')
      bModSet.add(setName + 'b.erase')
   else:
      image = setName + stackExtension
      modSet.add(setName + '.erase')

   for operation in ('SeedFinding', 'PatchTracking'):
      for axLet in axList:
         key = runtimePrefix + operation + '.' + axLet + '.rawBoundaryModel'
         if key in allDirectives[batDictInd]:
            rawBound = allDirectives[batDictInd][key][0]
            if axLet == 'b':
               bModSet.add(rawBound)
            else:
               modSet.add(rawBound)

   for ind in range(len(modSet)):
      rawBound = modSet.pop()
      if fixModelHeaderForImage(rawBound, image):
         return 1

   for ind in range(len(bModSet)):
      rawBound = bModSet.pop()
      if fixModelHeaderForImage(rawBound, setName + 'b' + stackExtension):
         return 1


# SINGLE-CALL FUNCTIONS FOR PROCESSING STEPS AND A FEW HELPERS FOR THEM

# Look at SDs of images and exclude views (mean not used because base unknown)
def analyzeSDsAdjustExcludes(clipLines):
   numAverage = 5
   maxExclude = 3
   crit = lookupDirective(runtimePrefix + 'Preprocessing', 'endExcludeCriterion', 0,
                          FLOAT_VALUE)
   if testDirectiveValue(crit, 'Preprocessing.endExcludeCriterion', 'float'):
      return 1
   if crit == None or crit <= 0.:
      return 0

   prnLog('Analyzing image SDs to see if views should be excluded at ends of series')

   # If there are no incoming lines from fixed stack, run on raw stack
   if not clipLines:
      clipLines = runClipStats('$clip stat ' + dataName + stackExtension, '', 'raw')
      if not clipLines:
         return 1

   # Extract the sds as second to last number
   sds = []
   for line in clipLines:
      if 'mean' in line or '-----' in line:
         continue
      if 'all' in line:
         break;
      lsplit = line.split()
      try:
         sds.append(float(lsplit[-1]))
      except Exception:
         abortSet('Error converting standard deviation to float in clip stats output')
         return 1

   if len(sds) < numAverage + 2 + 2 * maxExclude:
      return 0

   # Find out if there should be dark region analysis too
   darkRatio = lookupDirective(runtimePrefix + 'Preprocessing', 'darkExcludeRatio',
                               0, FLOAT_VALUE)
   darkFraction = lookupDirective(runtimePrefix + 'Preprocessing', 'darkExcludeFraction',
                                  0, FLOAT_VALUE)
   if testDirectiveValue(darkRatio, 'Preprocessing.darkExcludeRatio', 'float') or \
          testDirectiveValue(darkFraction, 'Preprocessing.darkExcludeFraction', 'float'):
      return 1
   if darkRatio == None:
      darkRatio = dfltDarkExcludeRatio
   if not darkFraction:
      darkFraction = dfltDarkExcludeFraction

   # Run histogram analysis if ratio is positive
   if darkRatio > 0.:
      numHist = min(zsize, 4 * maxExclude)
      clipcom = fmtstr('$clip hist -2d -s -iz 0-{},{}-{} {}{}', numHist / 2 - 1,
                       zsize - numHist / 2, zsize - 1, dataName, stackExtension)
      comfile = 'cliphist' + axisLet + comExt
      if writeTextFileReportErr(comfile, [clipcom]):
         return 1
      if runOneProcess(comfile, message = 'Analyzing histograms for large dark areas'):
         return 1
      cleanupFiles([comfile])
      histLines = readTextFileReportErr('cliphist' + axisLet + '.log')
      if not histLines:
         return 1

      # Find which ones pass the criteria and add to list
      darkExcludes = []
      for line in histLines:
         if line.startswith('Slice ') and 'no histogram dip' not in line:
            
            lsplit = line.split()
            try:
               zval = int(lsplit[1].replace(':', ''))
               lowPeak = float(lsplit[4])
               highPeak = float(lsplit[6])
               fraction = float(lsplit[-1])
            except Exception:
               abortSet('Error converting output from clip hist')
               return 1
            if lowPeak < darkRatio * highPeak and fraction > darkFraction:
               darkExcludes.append(zval)

   # Look for low images at ends of series
   excludes = []
   baseInd = 0
   baseView = 1
   for direc in (1, -1):

      # use varying indentations before computing the mean; get the mean
      for indent in range(2, maxExclude + 2):
         goodMean = 0.
         for ind in range(numAverage):
            goodMean += sds[baseInd + direc * (ind + indent)] / numAverage

         # See how many consecutive views are below the criterion; stop at first above
         numLow = 0
         for ind in range(indent):
            if sds[baseInd + direc * ind] < crit * goodMean:
               numLow += 1
            else:
               break

         # Done if found any and there is one not below criterion (to be conservative)
         if numLow and numLow < indent:
            for ind in range(numLow):
               excludes.append(baseView + direc * ind)
            break

      # If looking for dark regions too, now look from the first non-excluded one for
      # consecutive ones that would be excluded on that basis
      if darkRatio > 0.:
         for ind in range(numLow, 2 * maxExclude):
            view = baseView + direc * ind
            if view - 1 in darkExcludes:
               excludes.append(view)
            else:
               break

      # Set up to repeat on ending views
      baseInd = -1
      baseView = len(sds)

   if not len(excludes):
      prnLog('No views have low SDs or dark regions near end of series')
      return 0
   
   # Get existing skip list
   skipLet = ''
   if axisLet == 'b':
      skipLet = 'b'
   directive = copyPrefix + skipLet + 'skip'
   excludes.sort()
   if directive in allDirectives[batDictInd]:
      userExcludes = allDirectives[batDictInd][directive][0]
      message = 'Added these views to existing list of views to skip: '
      numAdd = 0
      userList = parselist(userExcludes)
      for view in excludes:
         if view not in userList:
            userExcludes += ',' + str(view)
            message += ' ' + str(view)
            numAdd += 1
      if not numAdd:
         prnLog('Views with low SDs were already in the list of views to skip')
         return 0

   else:
      userExcludes = ''
      message = 'These views have low SDs or dark regions and will be skipped: '
      comma = ''
      for view in excludes:
         userExcludes += comma + str(view)
         comma = ','
         message += ' ' + str(view)

   prnLog(message)

   # Modify all four com files
   for (comRoot, option, afterOpt) in [('xcorr', 'SkipViews', 'RotationAngle'),
                                       ('track', 'SkipViews', 'RotationAngle'),
                                       ('align', 'ExcludeList', 'RotationAngle'),
                                       ('tilt', 'EXCLUDELIST2', 'XTILTFILE')]:
      comfile = comRoot + axisCom
      sedcom = sedDelAndAdd(option, userExcludes, afterOpt)
      inLines = readTextFileReportErr(comfile)
      if not inLines:
         return 1
      if pysed(sedcom, inLines, comfile, retErr=True):
         abortSet('Error modifying ' + comfile)
         return 1
      
   return 0


# Operations needed when using fiducialless alignment for final alignment
def fidlessFileOperations():
   global didLocalAlign, madeZfactors
   didLocalAlign = False
   madeZfactors = False
   
   rotfile = 'rotation' + axisLet + '.xf'
   sinrot = math.sin(math.radians(axisRotation))
   if writeTextFileReportErr(rotfile, [fmtstr('{0:.6f} {1:.6f} {2:.6f} {0:.6f} 0. 0.',
                                              math.cos(math.radians(axisRotation)), 
                                              sinrot, -sinrot)]):
      return 1

   # Take care of xtilt file
   preLines = readTextFileReportErr(dataName + '.prexf')
   if not preLines:
      return 1
   xtilts = []
   for ind in range(len(preLines)):
      xtilts.append('0.')
   if writeTextFileReportErr(dataName + '.xtilt', xtilts):
      return 1

   try:
      runcmd('xftoxg -nfit 0 ' + dataName + '.prexf')
      runcmd(fmtstr('xfproduct {0}.prexg {1} {0}_nonfid.xf', dataName, rotfile))
      shutil.copyfile(dataName + '_nonfid.xf', dataName + '.xf')
      shutil.copyfile(dataName + '.rawtlt', dataName + '.tlt')
   except ImodpyError:
      reportImodError('Cannot prepare transformations')
      return 1
   except Exception:
      abortSet('Error copying xf or tlt file')
      return 1
   return 0


# Run tiltxcorr to do patch tracking
def runPatchTracking():
   contourPieces = lookupDirective(patchTrackText, 'contourPieces', 0, INT_VALUE)
   if isinstance(contourPieces, str):
      abortSet('Error converting contour pieces in directive to integer')
      return 1

   if not lookupDirective(comPrefix + 'xcorr_pt',
                          'tiltxcorr.SizeOfPatchesXandY', 0, STRING_VALUE):
      abortSet('Size of patches must be specified to use patch tracking')
      return 1

   if lookupDirective(comPrefix + 'xcorr_pt',
                      'tiltxcorr.NumberOfPatchesXandY', 0, STRING_VALUE) and \
      lookupDirective(comPrefix + 'xcorr_pt',
                      'tiltxcorr.OverlapOfPatchesXandY', 0, STRING_VALUE):
      abortSet('You cannot enter both the number of patches to track and the ' + \
               'fractional overlap of patches (xcorr_pt.tiltxcorr.NumberOfPatchesXandY' +\
               ' and xcorr_pt.tiltxcorr.OverlapOfPatchesXandY)')
      return 1

   # Start command list for running makecomfile
   comlines = ['InputFile xcorr' + axisCom,
               'BinningOfImages ' + str(coarseBinning),
               'RootNameOfDataFiles ' + dataName]
   comlines += laterComDirectives(0)

   # Look for boundary models to transform
   (error, boundaryMod) = getOrXformBoundaryModel(patchTrackText, '_ptbound.mod', True)
   if error:
      return 1
   if boundaryMod:
      comlines.append(fmtstr('OneParameterChange {}xcorr_pt{}.tiltxcorr.' +
                             'BoundaryModel={}', comPrefix, axisLet, boundaryMod))

   # Convert runtime directive for contour pieces into com directive for imodchopconts
   if contourPieces and contourPieces > 1:
      pieceLength = lookupDirective(comPrefix + 'xcorr_pt', 
                                    'imodchopconts.LengthOfPieces', 0, INT_VALUE)

      # Check if length is entered, allow the default to be overridden
      if isinstance(pieceLength, str):
         abortSet('Error converting length of contour pieces in directive to integer')
         return 1
      if pieceLength and pieceLength > 0:
         abortSet('Both xcorr_pt.imodchopconts.LengthOfPieces > 0 and ' +\
                     'PatchTracking.contourPieces were entered; only one is allowed')
         return 1
      comlines.append(fmtstr('OneParameterChange {}xcorr_pt{}.imodchopconts.' + \
                                'NumberOfPieces={}', comPrefix, axisLet, contourPieces))

   if makeAndRunOneCom(comlines, 'xcorr_pt' + axisCom,
                       'Tracking patches to make alignment model'):
      return 1


# Test whether alignment after patch tracking gave a tilt angle adjustment above limits
def OKtoAdjustPatchTrack():
   maxTiltAdjust = lookupDirective(patchTrackText, 'maxTiltAdjustment', 0, 
                                   FLOAT_VALUE)
   if isinstance(maxTiltAdjust, str) or not maxTiltAdjust:
      maxTiltAdjust = dfltPTmaxTiltAdjust
   maxAdjustedAngle = lookupDirective(patchTrackText, 'maxAdjustedAngle', 0, 
                                      FLOAT_VALUE)
   if isinstance(maxAdjustedAngle, str) or not maxAdjustedAngle:
      maxAdjustedAngle = dfltPTmaxAdjustedAngle

   if math.fabs(totalDelTilt) > maxTiltAdjust:
      warning(fmtstr('Patch tracking not being iterated: tilt angle adjustment ' + \
                     'by {:.1f} is above the limit of {:.1f}', totalDelTilt,
                     maxTiltAdjust))
      return False
   
   # Forgive any problems reading the tilt angles
   lines = readTextFile(setName + axisLet + '.rawtlt', returnOnErr = True)
   if isinstance(lines, str):
      return True
   maxAngle = 0.;
   for line in lines:
      lsplit = line.split()
      try:
         maxAngle = max(maxAngle, math.fabs(float(lsplit[0]) + totalDelTilt))
      except Exception:
         return True

   if maxAngle > maxAdjustedAngle:
      warning(fmtstr('Patch tracking not being iterated: tilt angle adjustment ' + \
                     'by {:.1f} would make the highest angle be {:.1f}, above the ' + \
                     'limit of {:.1f}', totalDelTilt, maxAngle, maxAdjustedAngle))
      return False

   return True

            
# Get fiducial model by seeding and tracking or by RAPTOR
def makeSeedAndTrack():
   trackingMethod = lookupDirective(runtimePrefix + 'Fiducials', 'trackingMethod',
                                    0, INT_VALUE)
   numBTRuns = lookupDirective(runtimePrefix + 'BeadTracking', 'numberOfRuns',
                               0, INT_VALUE)
   seedingMethod = lookupDirective(runtimePrefix + 'Fiducials', 'seedingMethod',
                                   0, INT_VALUE)
   
   if isinstance(trackingMethod, str) or isinstance(numBTRuns, str) or \
          isinstance(seedingMethod, str):
      abortSet('Error converting tracking method or number of runs to integer')
      return 1
   if not trackingMethod:
      trackingMethod = 0
   if trackingMethod < 0 or trackingMethod > 2:
      abortSet('trackingMethod must be between 0 and 2')
   if (seedingMethod == None or seedingMethod < 1 or seedingMethod > 3 or \
          (axisInd == 0 and seedingMethod == 2)) and trackingMethod == 0:
      abortSet('seedingMethod must be between 1 and 3 and not be 2 for first/only axis')
      return 1

   # If runs is not there, assume 0 for RAPTOR and 1 for other
   if numBTRuns == None:
      if trackingMethod == 2:
         numBTRuns = 0
      else:
         numBTRuns = 1
   if numBTRuns <= 0 and trackingMethod != 2:
       abortSet('Number of beadtrack runs must be > 0 unless using RAPTOR')
       return 1

   # Modify track.com with ImageBinned entry
   btlines = readTextFileReportErr('track' + axisCom)
   if not btlines:
      return 1
   sedcom = sedDelAndAdd('ImagesAreBinned', coarseBinning, 'OutputModel')
   if pysed(sedcom, btlines, 'track' + axisCom, retErr=True):
      abortSet('Error modifying track' + axisCom)
      return 1
   lightBeads = optionValue(btlines, 'LightBeads', BOOL_VALUE)
   lightInt = 0
   if lightBeads:
      lightInt = 1
   edfcom = edfDelAndAdd('Track.' + axisUpperLet + '.LightBeads', lightInt)
   if modifyEdfLines(edfcom):
      return 1

   # RAPTOR
   if needStep(4) and trackingMethod == 2:
      markers = lookupDirective(runtimePrefix + 'RAPTOR', 'numberOfMarkers', 0,
                                INT_VALUE)
      if markers == None or isinstance(markers, str) or markers <= 0:
         abortSet('numberOfMarkers for RAPTOR missing, not positive, or gave ' +\
                     'conversion error')
         return 1

      if lookupDirective(runtimePrefix + 'RAPTOR', 'useAlignedStack', 0,
                         BOOL_VALUE):
         rapStack = datasetFilename('.preali')
         beadint = int(round(fidSizePix / coarseBinning))
      else:
         rapStack = dataName + stackExtension
         beadint = int(round(fidSizePix))

      line = fmtstr('$runraptor -mark {} -diam {} {}', markers, beadint, rapStack)
      if writeTextFileReportErr('runraptor' + axisCom,
                          ['# Command file to run raptor', line]):
         return 1
      if runOneProcess('runraptor' + axisCom, message = 'Tracking fiducials with RAPTOR'):
         return 1
      if useFileAsReplacement(dataName + '_raptor.fid', dataName + '.fid', False, True):
         return 1

   # Seed and track:
   elif needStep(4):

      # Transferfid for axis b
      skipAuto = False
      if axisInd > 0 and (seedingMethod & 2) != 0:
         comlines = ['RootNameOfDataFiles	' + setName]
         comlines += laterComDirectives(0)
         comlines.append(fmtstr('OneParameterChange {}transferfid.transferfid.' +\
                                   'LowestTiltTransformFile={}_AtoB.xf', comPrefix,
                                setName))
         if makeAndRunOneCom(comlines, 'transferfid' + comExt,
                             'Transferring fiducials from A to B axis'):
            return 1
         printTaggedMessages('transferfid.log',
                             standardTags + [('fiducials that failed', 4)])
         numFailed = findTaggedValue(latestMessages, 'fiducials that failed', ':',
                                     INT_VALUE)
         if numFailed != None and numFailed == 0:
            skipAuto = True

      # Autoseed
      if (seedingMethod & 1) != 0 and not skipAuto:
         twoSurf = 0
         if numSurfaces > 1:
            twoSurf = 1

         # If there is not a specific directive for two surfaces, set it based on
         # tiltalign entry
         comlines = laterComDirectives(0)
         if axisInd > 0 and (seedingMethod & 2) != 0:
            comlines.append(fmtstr('OneParameterChange {}autofidseedb.autofidseed.' +
                                   'AppendToSeedModel=1', comPrefix))
         if lookupDirective(comPrefix + 'autofidseed', 'autofidseed.TwoSurfaces', 0,
                            STRING_VALUE) == None:
            comlines.append(fmtstr('OneParameterChange {}autofidseed{}.autofidseed.' +
                                   'TwoSurfaces={}', comPrefix, axisLet, twoSurf))

         # If transferfid was not done, do not append to dummy model
         if axisInd > 0 and (seedingMethod & 2) == 0:
            comlines.append(fmtstr('OneParameterChange {}autofidseed{}.autofidseed.' +
                                   'AppendToSeedModel=0', comPrefix, axisLet))
            
         (error, boundaryMod) = getOrXformBoundaryModel \
                                (autoSeedText, '_afsbound.mod', (seedingMethod & 2) != 0)
         if error:
            return 1
                                                        
         # Add boundary model to com
         if boundaryMod:
            comlines.append(fmtstr('OneParameterChange {}autofidseed{}.autofidseed.' +
                                   'BoundaryModel={}', comPrefix, axisLet, boundaryMod))
            
         if makeAndRunOneCom(comlines, 'autofidseed' + axisCom,
                             'Finding seed points for fiducial model'):
            return 1
         printTaggedMessages('autofidseed' + axisLet + '.log', standardTags +
                             [('candidate points, including', 4), ('Final:   total', 0),
                              ('[AFS1]', 4), ('[AFS3]', 4)])

         # If track file was adjusted for either reason, use the adjusted file
         for line in latestMessages:
            if '[AFS1]' in line or '[AFS3]' in line:
               realTrack = 'track' + axisCom
               adjFile = 'track' + axisLet + '_adjusted' + comExt
               if not os.path.exists(adjFile):
                  abortSet('Autofidseed adjusted tracking parameters but ' + adjFile + \
                           ' does not exist')
                  return 1

               origTrack = 'track' + axisLet + '_orig' + comExt
               if os.path.exists(origTrack):
                  makeBackupFile(realTrack)
               elif renameAndAbort(realTrack, origTrack, 
                                   'Renaming {} to {} before adopting adjusted file'):
                  return 1
               if renameAndAbort(adjFile, realTrack, 
                              'Renaming adjusted track command file {} to {}'):
                  return 1

               break
                  
      edfcom = edfDelAndAdd(axisEdfLet + '.SeedingDone', 'true')
      if modifyEdfLines(edfcom):
         return 1

   # Run bead tracking indicated number of times
   if not needStep(5):
      numBTRuns = 0
   for trackInd in range(max(0, numBTRuns)):

      # After first time, save seed as _orig or back it up
      if (trackInd or trackingMethod == 2) and \
             useFileAsReplacement(dataName + '.fid', dataName + '.seed', True, True):
         return 1
      if runOneProcess('track' + axisCom, message =
                       'Tracking beads with Beadtrack, run # ' + str(trackInd + 1)):
         return 1
      printTaggedMessages('track' + axisLet + '.log',
                          standardTags + [('Total points missing =', 4)])
      missing = findTaggedValue(latestMessages, 'Total points missing', '=', INT_VALUE)
      if missing != None and missing == 0:
         break


# Add items from most recent alignment to the edf file
def postProcessTiltalign():
   zfacs = optionValue(alignLines, 'XStretchOption', INT_VALUE, numVal = 1)
   localAli = optionValue(alignLines, 'LocalAlignments', BOOL_VALUE)
   zShift = optionValue(alignLines, 'AxisZShift', FLOAT_VALUE, numVal = 1)
   angle = optionValue(alignLines, 'AngleOffset', FLOAT_VALUE, numVal = 1)
   edfcom = edfDelAndAdd('MadeZFactors' + axisUpperLet,
                         boolStringForEdf(zfacs != None and zfacs > 0))
   edfcom += edfDelAndAdd('UsedLocalAlignments' + axisUpperLet,
                          boolStringForEdf(localAli))
   edfcom += edfDelAndAdd(axisEdfLet + '.align.AxisZShift', round(zShift, 2))
   edfcom += edfDelAndAdd(axisEdfLet + '.align.AngleOffset', round(angle, 3))
   return modifyEdfLines(edfcom)
   
      
# Set up tiltalign command file from alignLines with the changes in sedcom, then run
# restrictalign unless skipping it and report results, then run tiltalign
#
def modifyRestrictAndRunAlign(comfile, sedcom, localAli, message, skipRestrict):
   global suppressAbort, alignLines
   
   if pysed(sedcom, alignLines, comfile, retErr=True):
      abortSet('Error modifying ' + comfile)
      return (-1, 0)
   alignLines =  pysed(sedcom, alignLines, retErr=True)
   if isinstance(alignLines, str):
      abortSet('Error modifying alignment lines: ' + str)
      return (-1, 0)

   retval = 0
   noRobust = 0
   ranAlready = False
   if not skipRestrict:

      # Get directives to modify defaults if any when running restrictalign
      minRatio = lookupDirective(runtimePrefix + 'RestrictAlign', 'minMeasurementRatio',
                                 0, FLOAT_VALUE)
      if not minRatio:
         minRatio = 0.
      targetRatio = lookupDirective(runtimePrefix + 'RestrictAlign',
                                    'targetMeasurementRatio', 0, FLOAT_VALUE)
      if not targetRatio:
         targetRatio = 0.
      resOrder = lookupDirective(runtimePrefix + 'RestrictAlign', 'orderOfRestrictions',
                                 0, STRING_VALUE)
      skipBT = lookupDirective(runtimePrefix + 'RestrictAlign', 'skipBeamTiltWithOneRot',
                               0, BOOL_VALUE)
      comlines = ['InputFile ' + comfile,
                  'OutputFile ' + 'restrictalign' + axisCom]
      
      locText = 'global '
      if localAli:
         comlines.append('LocalAlignValidation 3')
         locText = 'local '
      if minRatio or targetRatio:
         comlines.append(fmtstr('TargetAndMinRatios {},{}', targetRatio, minRatio))
      if resOrder:
         comlines.append(fmtstr('OneParameterChange {}restrictalign{}.restrictalign.' + \
                                'OrderOfRestrictions={}', comPrefix, axisLet, resOrder))
      if skipBT:
         comlines.append('SkipBeamTiltWithOneRot 1')

      comlines += laterComDirectives(0)

      if makeAndRunOneCom(comlines, 'restrictalign' + axisCom, 'Running restrictalign' + \
                          ' to optimize ' + locText + 'alignment parameters'):
         suppressAbort = False
         return (-1, 0)

      fromLog = 'restrictalign' + axisLet + '.log'
      resLines = readTextFileReportErr(fromLog)
      if not resLines:
         suppressAbort = False
         return (-1, 0)

      # Copy log to _global if there isn't one, or to _local
      if not localAli:
         toLog = 'restrictalign' + axisLet + '_global.log'
         copyGlobal = not isFileNewer(toLog, 'origcoms/align' + axisCom)
      if localAli or copyGlobal:
         if localAli:
            toLog = 'restrictalign' + axisLet + '_local.log'
         makeBackupFile(toLog)
         try:
            shutil.copyfile(fromLog, toLog)
         except Exception:
            prnLog('WARNING: Failed to copy ' + fromLog + ' to ' + toLog)

      # Echo all the summary output
      doOutput = False
      retval = 1
      for line in resLines:
         if line.startswith('restrictalign:'):
            doOutput = True
         if 'CHUNK DONE' in line or 'COMPLETED' in line or 'ERROR:' in line:
            break
         if 'No restriction' in line:
            retval = 0
         if '[rsa2]' in line:
            noRobust = 1
            if 'fail' in line:
               noRobust = 2
         if '[rsa1]' in line:
            ranAlready = True
            doOutput = False
         if doOutput:
            prnLog(line.strip('\r\n'))
      
      prnLog(' ')
      if noRobust:
         message = message.replace('with robust', 'without robust')
      if retval:
         alignLines = readTextFileReportErr(comfile)
         if not alignLines:
            suppressAbort = False
            return (-1, 0)
         
   if ranAlready:
      prnLog(message.replace('Doing', 'Did') + '     [brt3]')
   elif runOneProcess(comfile, message = message):
      return (-2, noRobust)
   if postProcessTiltalign():
      return (-2, noRobust)
   return (retval, noRobust)


# Run tiltalign in 2 or 3 stages
def runTiltalign():
   global xtiltNeeded, fidThickness, fidIncShift, reconThickness, didLocalAlign
   global madeZfactors, totalDelTilt, suppressAbort, finalAlignResid, alignLines

   minTotGlblStretch = 12
   minSurfGlblStretch = 4
   minRatioGlblStretch = 0.125
   minSurfLocalStretch = 1.0
   minTotSkipRestrict = 20
   minTotSkipCrossVal = 80
   minTotToSetAngle = 4
   
   # The command file should be configured by various inputs, so let's find out some
   patchSizeArr = optionValue(alignLines, 'TargetPatchSizeXandY', INT_VALUE, numVal = 2)
   minFidsArr = optionValue(alignLines, 'MinFidsTotalAndEachSurface', INT_VALUE,
                            numVal = 2)
   if not patchSizeArr or not minFidsArr:
      abortSet('Problem finding some options in align' + comExt)
   (nxpatch, nypatch) = patchSizeArr
   (minFidsTot, minFidsSurf) = minFidsArr
   alignCom = 'align' + axisCom
   alignLog = 'align' + axisLet + '.log'
   doRobust = lookupDirective(comPrefix + 'align', 'tiltalign.RobustFitting', 0,
                              BOOL_VALUE)
   robustOrig = doRobust

   enableStretch = lookupDirective(runtimePrefix + 'TiltAlignment', 'enableStretching',
                                   0, BOOL_VALUE)
   didLocalAlign = lookupDirective(comPrefix + 'align', 'tiltalign.LocalAlignments',
                                   0, BOOL_VALUE)
   
   crossVal = lookupDirective(comPrefix + 'restrictalign',
                              'restrictalign.UseCrossValidation', 0, INT_VALUE)
   if crossVal == None:
      crossVal = 1
                                  
   if patchTrack:
      enableStretch = False

   # If skipping align, get a few more values from com file, analyze log
   if skipTiltalign:
      glbStretch = optionValue(alignLines, 'XStretchOption', INT_VALUE, numVal = 1)
      madeZfactors = glbStretch != None and glbStretch > 0
      angleArr = [0., 0., 0., 0., 0., 0., 0, 0]
      if analyzeAlignLog(numSurfaces > 1, angleArr):
         return 1
      (totalDelTilt, xtiltNeeded, fidThickness, fidIncShift, fidTotalShift, \
          reconThickness, numBot, numTop) = angleArr
      return 0
   
   # first time, turn off local align and make sure stretch/skew off
   # But loop twice in case have to turn off robust fitting
   for loop in (0, 1):
      robText = 'without'
      if doRobust:
         robText = 'with'
      mess = 'Doing fine alignment with no distortion or local alignments, ' + \
          robText + ' robust fitting'
      suppressAbort = not loop and doRobust
      sedcom = [sedModify('SurfacesToAnalyze', numSurfaces),
                sedModify('LocalAlignments', 0),
                sedModify('XStretchOption', 0),
                sedModify('SkewOption', 0)] + \
                sedDelAndAdd('ImagesAreBinned', coarseBinning, 'OutputTransformFile') + \
                sedDelAndAdd('ScaleShifts', '1,' + str(coarseBinning), 'InputFile2')+ \
                sedDelAndAdd('RobustFitting', doRobust, 'OutputTransformFile')

      # Do restrictalign only the first time
      (restricted, noRobust) =  modifyRestrictAndRunAlign(alignCom, sedcom, 0, mess,
                                                          loop > 0)
      globalNoRobust = noRobust
      if noRobust:
         doRobust = 0
      if restricted < 0:
         if not suppressAbort:
            return 1
         suppressAbort = False
         
         # Look for robust failure: this shouldn't happen with cross-validation
         for line in latestMessages:
            if 'TOO FEW DATA POINTS TO DO ROBUST' in line.upper():
               prnLog('Trying again without robust fitting')
               doRobust = 0
               break
         else:    # ELSE ON FOR
            return 1

         continue

      messageTags = standardTags + [('Residual error', 4), ('leave-out error', 4),
                                    ('Benefit from', 4)]
      printTaggedMessages(alignLog, messageTags)
      break

   suppressAbort = False
   angleArr = [0., 0., 0., 0., 0., 0., 0, 0]
   if analyzeAlignLog(numSurfaces > 1, angleArr):
      return 1
   (totalDelTilt, xtiltNeeded, fidThickness, fidIncShift, fidTotalShift, \
       reconThickness, numBot, numTop) = angleArr
   if numSurfaces > 1:
      totalFid = numBot + numTop
      minOnSurf = min(numBot, numTop)
   else:
      totalFid = numBot

   localAlign = 0
   glbStretch = 0
   locStretch = 0
   useMinFids = 0
   numruns = 1
   madeZfactors = False
   if didLocalAlign and (not restricted or crossVal):
      localAlign = 1
      messageTags = standardTags + [('(Global)', 4), ('error local mean:', 4),
                                    ('Local.*leave-out error', 4), ('Benefit from', 4)]
   else:
      didLocalAlign = 0
   if (not restricted or crossVal) and (didLocalAlign or enableStretch):
      if didLocalAlign:
         numruns = 2
         
         # Turn on robust again for locals if they were on before and were turned off
         # because of no benefit - failure gives a 2
         if crossVal and noRobust == 1:
            doRobust = robustOrig

      # are there enough fids for stretch?
      if enableStretch:
         if totalFid > minTotGlblStretch and \
                (numSurfaces == 1 or \
                    (minOnSurf > minSurfGlblStretch and \
                        min(numBot, numTop) / float(totalFid) > minRatioGlblStretch)):
            glbStretch = 3
            madeZfactors = True
            numruns = 2
            if didLocalAlign and numSurfaces > 1:
               mindens = minOnSurf / float(rawXsize * rawYsize)
               minInArea = mindens * nxpatch * nypatch
               if minInArea > minSurfLocalStretch:
                  locStretch = 3
                  useMinFids = minFidsSurf
               else:
                  prnLog('Too few fiducials on minority surface to enable local ' +\
                            'stretching solution')

         else:
            prnLog('Too few fiducials on minority surface to enable ' +\
                            'stretching solution')
            
   
   # Do restrictalign if there are enough fiducials for old ratio-based restricting,
   # or if doing either locals or stretch, as long as there are enough fiducials if
   # adding just stretch.  But skip it if run before on this data set
   skipRestrict = (not crossVal and totalFid >= minTotSkipRestrict) or \
      (crossVal and not didLocalAlign and \
       (not glbStretch or totalFid >= minTotSkipCrossVal))
   if isFileNewer('restrictalign' + axisLet + '_local.log',
                  'origcoms/align' + axisCom):
      skipRestrict = True

   reloaded = False
   
   
   # Just run alignment once or twice more.  Again, if robust gives a problem, try it
   # without; and this time if local alignments failed, drop that out too
   for run in range(numruns):
      for loop in (0, 2):
         robText = 'without'
         if doRobust:
            robText = 'with'
         mess = 'Doing fine alignment with '
         if not glbStretch:
            mess += 'no'
         mess += ' distortion, with'
         if (restricted and not crossVal) or not didLocalAlign:
            mess += 'out'
         mess += ' local alignments, ' + robText + ' robust fitting'
         suppressAbort = not loop and doRobust
         delTilt = totalDelTilt
         if totalFid < minTotToSetAngle or \
            lookupDirective(runtimePrefix + 'TiltAlignment', 'noAngleOffset', 0,
                            BOOL_VALUE):
            delTilt = 0.
         sedcom = [sedModify('SurfacesToAnalyze', numSurfaces),
                   sedModify('LocalAlignments', localAlign),
                   sedModify('AngleOffset', delTilt)] + \
                   sedDelAndAdd('ImagesAreBinned', coarseBinning, 'OutputTransformFile')+\
                   sedDelAndAdd('ScaleShifts', '1,' + str(coarseBinning), 'InputFile2')+ \
                   sedDelAndAdd('RobustFitting', doRobust, 'OutputTransformFile')
         if not reloaded:
            sedcom += [sedModify('XStretchOption', glbStretch),
                       sedModify('SkewOption', glbStretch),
                       sedModify('LocalStretchOption', locStretch),
                       sedModify('LocalSkewOption', locStretch),
                       sedModify('MinFidsTotalAndEachSurface', fmtstr('{},{}', minFidsTot,
                                                                      useMinFids))]

         if madeZfactors:
            sedcom += sedDelAndAdd('OutputZFactorFile', dataName + '.zfac',
                                   'OutputTransformFile')
                
         (restricted, noRobust) =  modifyRestrictAndRunAlign \
             (alignCom, sedcom, localAlign, mess, skipRestrict or run > 0 or loop > 0)
         if restricted > 0 and crossVal and not run and not loop:
            reloaded = True
         if noRobust:
            doRobust = 0

         if restricted < 0:
            if not suppressAbort:
               return 1
            suppressAbort = False
            for line in latestMessages:
               if 'TOO FEW DATA POINTS TO DO ROBUST' in line.upper():
                  prnLog('Trying again without robust fitting')
                  doRobust = 0
                  break
               if 'Minimum numbers of fiducials are too high' in line:
                  prnLog('Trying again without local alignments')
                  didLocalAlign = 0
                  localAlign = 0
                  break
            else:     # ELSE ON FOR
               return 1

            continue
            
         printTaggedMessages(alignLog, messageTags)
         break

      suppressAbort = False
      if analyzeAlignLog(numSurfaces > 1, angleArr):
         return 1
      (totalDelTilt, xtiltNeeded, fidThickness, fidIncShift, fidTotalShift, \
          reconThickness, numBot, numTop) = angleArr

   # Find the last residual output message and echo it as the final result, save for eval
   for ind in range(len(latestMessages) - 1, -1, -1):
      line = latestMessages[ind]
      colonInd = line.find(':')
      if colonInd > 0 and 'error' in line and 'mean' in line:
         prnLog('Final align - ' + line.strip() + logSuffixTag)
         prnLog('')
         lsplit = line[colonInd + 1:].split()
         finalAlignResid = 0.
         try:
            finalAlignResid = float(lsplit[0])
         except Exception:
            pass
         break

   # Get the ta...log files
   optsNames = (('m', 'Mappings'), ('e', 'Error'), ('s', 'Solution'), ('l', 'Locals'),
                ('c', 'Coordinates'), ('a', 'Angles'), ('b', 'Beamtilt'), ('r', 'Robust'))
   for opt in optsNames:
      name = 'ta' + opt[1] + axisLet + '.log'
      try:
         partLines = runcmd(fmtstr('alignlog -{} {}', opt[0], alignLog))
         if partLines:
            for ind in range(len(partLines)):
               partLines[ind] = partLines[ind].rstrip('\r\n')
         err = writeTextFile(name, partLines, True)
         if err:
            warning('Error ' + err)
      except ImodpyError:
         warning('Error running alignlog -' + opt[0])

   return 0
     

# Make the aligned stack: also sets some globals before deciding whether it is needed
def makeAlignedStack(positionBinning):
   global aliBinning, aliXunbinned, aliYunbinned, correctCTF, transposeForAli
   global expandFactor

   alipre = runtimePrefix + 'AlignedStack'
   if positionBinning > 0:
      aliBinning = positionBinning
      outsizeText = ''
      mess = 'Making binned aligned stack for whole tomogram positioning'

   else:
      aliBinning = lookupDirective(alipre, 'binByFactor', 0, INT_VALUE)
      if testDirectiveValue(aliBinning, 'AlignedStack.binByFactor', 'integer'):
         return 1
      outsizeText = lookupDirective(alipre, 'sizeInXandY', 0, STRING_VALUE)
      mess = 'Making final aligned stack'
      
   linear = lookupDirective(alipre, 'linearInterpolation', 0, BOOL_VALUE)
   if not linear:
      linear = 0
   aliXunbinned = rawXsize
   aliYunbinned = rawYsize
   if not aliBinning:
      aliBinning = 1

   edfcom = ['/batchruntomo.Stack.' + axisUpperLet + '.SizeToOutputInXandY/d']
   if outsizeText:
      splits = outsizeText.replace(',', ' ').split()
      try:
         aliXunbinned = int(splits[0])
         aliYunbinned = int(splits[1])
         edfcom = edfDelAndAdd('Stack.' + axisUpperLet + '.SizeToOutputInXandY',
                               fmtstr('{},{}', aliXunbinned, aliYunbinned))
      except:
         abortSet('Error converting aligned stack output size')
         return 1

   if transposeForAli:
      (aliXunbinned, aliYunbinned) = (aliYunbinned, aliXunbinned)

   if positionBinning < 0 and skipAlignStackMods and rawFor3dCTF:
      return 0

   (comRoot, process) = comAndProcessForAlignedStack('')
   if ifMontage:
      err = montageFrameValues(rawXsize, rawYsize, transposeForAli, aliXunbinned,
                               aliYunbinned, montFrameData)
      if err == 1:
         reportImodError("Error running goodframe on montage sizes")
         return 1
      if err:
         abortSet("Error converting output of goodframe to integers")
         return 1
      sedcom = [sedModify('StartingAndEndingX', fmtstr('{},{}', montFrameData[2],
                                                       montFrameData[4])),
                sedModify('StartingAndEndingY', fmtstr('{},{}', montFrameData[3],
                                                       montFrameData[5])),
                '/^InterpolationOrder/d']
      if linear:
         sedcom.append('/^TransformFile/a/InterpolationOrder    1/')
                
      aliXunbinned = montFrameData[0]
      aliYunbinned = montFrameData[1]

   else:
      expandFactor = lookupDirective(comPrefix + 'newst', 'newstack.ExpandByFactor', 0,
                                     FLOAT_VALUE)
      if not expandFactor:
         expandFactor = 1.
      sedcom =  [sedModify('SizeToOutputInXandY',
                           fmtstr('{},{}',
                                  int(expandFactor * (aliXunbinned // aliBinning)),
                                  int(expandFactor * (aliYunbinned // aliBinning))))] + \
                       sedDelAndAdd('LinearInterpolation', linear, 'TransformFile')
      
   sedcom += sedDelAndAdd('BinByFactor', aliBinning, 'TransformFile')

   # Evaluate whether the aligned stack needs to be rebuilt, if it has already been
   # corrected and correction is to be done
   needRebuild = False
   if not (positionBinning > 0 or needStep(8)) and needStep(ctfCorrStepNum) and \
          correctCTF and os.path.exists(datasetFilename('.ali')):
      try:
         headLines = runcmd('header ' + datasetFilename('.ali'))
      except ImodpyError:
         reportImodError('Error reading header of existing aligned stack')
         return 1
      for line in headLines:
         if 'ctfPhaseFlip' in line or 'CTF correct' in line:
            needRebuild = True
            prnLog('Remaking aligned stack because it has already been CTF corrected')
      
   if positionBinning > 0 or needStep(8) or needRebuild:
      if modifyWriteAndRunCom(comRoot + axisCom, sedcom, message = mess):
         return 1
      if positionBinning <= 0:
         edfcom += edfDelAndAdd('ImageRotationForAliStack', axisRotation)
         edfcom += edfDelAndAdd('Stack.' + axisUpperLet + '.UseLinearInterpolation',
                                boolStringForEdf(linear))
         if modifyEdfLines(edfcom):
            return 1

   return 0
   

# Do optional CTF correction on aligned stack
def CTFPlotAlignedStack():
   replaced = replaceOrRunAfterStep(ctfPlotStepNum, False)
   if replaced:
      return max(0, replaced)
   
   # Do autofitting with ctfplotter if enabled
   if needStep(ctfPlotStepNum) and correctCTF and autofitCTF:
      ctfLines = readTextFileReportErr('ctfplotter' + axisCom)
      if not ctfLines:
         return 1
      sedcom = ['/^ExpectedDefocus/a/SaveAndExit	1/',
                '/^ExpectedDefocus/a/AutoFitRangeAndStep	' + autofitCTF + '/',
                '/^AngleRange/d']

      # Make existing defocus file a backup to avoid error messages
      makeBackupFile(dataName + '.defocus')
      if modifyWriteAndRunCom('ctfplotter_auto' + axisCom, sedcom, ctfLines,
                              'Finding defocus for CTF correction with Ctfplotter'):
         return 1

      return replaceOrRunAfterStep(ctfPlotStepNum, True)

   return 0


# Do correction with ctfphaseflip
def CTFCorrectAlignedStack():

   # Write the simple defocus file and use it if there is no defocus file; but if there
   # is, make sure it is in the com file
   focFileName = dataName + '.defocus'
   simpleName = dataName + '_simple.defocus'
   if not os.path.exists(focFileName):
      if writeTextFileReportErr(simpleName, [fmtstr('{} {} 0. 0. {}', zsize // 2,
                                                    zsize // 2, defocus)]):
         return 1
      sedcom = [sedModify('DefocusFile', simpleName)]
   else:
      sedcom = [sedModify('DefocusFile', focFileName)]

   sedcom.append(sedModify('PixelSize', aliBinning * pixelSize))
   if useGPU:
      sedcom += sedDelAndAdd('UseGPU', 0, 'DefocusFile')

   # See if X axis tilt is to be corrected; if so try to get it and add option
   correctXtilt = lookupDirective(runtimePrefix + 'CTFcorrection', 'correctForXAxisTilt',
                                  0, BOOL_VALUE)
   if correctXtilt and not noXaxisTilt:
      angleArr = [0., 0., 0., 0., 0., 0., 0, 0]
      if analyzeAlignLog(numSurfaces > 1, angleArr, noLogOK = True):
         return 1
      xtiltNeeded = angleArr[1]
      err = parseTomopitchLog(angleArr)
      if err < 0:
         return 1
      if not err:
         xtiltNeeded = angleArr[1]

      if xtiltNeeded:
         sedcom += sedDelAndAdd('XAxisTilt', xtiltNeeded, 'DefocusFile')
   
   comfile = 'ctfcorrection' + axisCom
   ctfLines = readTextFileReportErr(comfile)
   if not ctfLines:
      return 1
   if pysed(sedcom, ctfLines, comfile, retErr=True):
      abortSet('Error modifying ' + comfile)
      return 1

   if skipAlignStackMods:
      return 0

   numProc = numberOfProcessingUnits()
   if numProc < 1:
      return 1;

   if numProc > 1:
      target = 2 * numProc
      maxSlices = (zsize + target - 1) // target
      try:
         runcmd(fmtstr('splitcorrection -m {} {}', maxSlices, comfile))
      except ImodpyError:
         reportImodError('Error trying to run ctfcorrection in parallel')
         releaseGPUallocation()
         return 1
   
   err = runOneProcess(comfile, numProc < 2, useGPU,
                       message = 'Correcting for CTF with Ctfphaseflip')
   releaseGPUallocation()
   if err:
      return 1

   if useFileAsReplacement(datasetFilename('_ctfcorr.ali'), datasetFilename('.ali'),
                           False, True):
      return 1

   return 0
      

# Optionally detect gold for erasing with 3D method
def detectGoldIn3D():

   # Get the binning or use the binning that etomo would assign, bead size / 5 rounded
   # to an integer with a minimum binned size of 4
   if eraseGold and eraseGold > 1 and needStep(detect3DstepNum):
      binning = lookupDirective(runtimePrefix + 'GoldErasing', 'binning', 0, INT_VALUE)
      if not binning:

         # desired reduction is expFac * fidSize / optimal
         # actual reduction will be binning / expFac so binning is expFac * desired
         binning = int(round(expandFactor**2 * fidSizePix / fb3dOptimalBinnedSize))
         if binning > 1 and expandFactor**2 * fidSizePix / binning < fb3dMinBinnedSize:
            binning -= 1

      # Get the aligned stack if binning differs
      if binning != aliBinning:
         (comRoot, process) = comAndProcessForAlignedStack('')
         comlines = ['InputFile ' + comRoot + axisCom,
                     'BinningOfImages ' + str(binning),
                     'RootNameOfDataFiles ' + dataName]
         if not ifMontage:
            comlines.append(fmtstr('OneParameterChange {}newst_3dfind{}.newstack.' + \
                                      'SizeToOutputInXandY={},{}', comPrefix, axisLet,
                                   int(aliXunbinned * expandFactor) // binning,
                                   int(aliYunbinned * expandFactor) // binning))
         comlines += laterComDirectives(0)
         if makeAndRunOneCom(comlines, comRoot + '_3dfind' + axisCom,
                             'Making aligned stack for erasing gold'):
            return 1

      # Set up to get the reconstruction; get a thickness if any entered
      thickness = lookupDirective(runtimePrefix + 'GoldErasing', 'thickness', 0,
                                  INT_VALUE)

      # use fid alignment if two surfaces, otherwise there needs to be a directive
      if not thickness:
         if numSurfaces > 1:
            useSize = max(15., fidSizeNm) / (pixelSize / expandFactor)
            thickness = 2 * ((int(round(1.1 * fidThickness + 6. * useSize)) + 1) // 2)
         else:
            abortSet('A GoldErasing.thickness directive must be supplied for 3D gold ' + \
                     'finding')
            return 1

      # Get the com file
      comfile = 'tilt_3dfind' + axisCom
      comlines = ['InputFile tilt' + axisCom,
                  'OutputFile ' + comfile,
                  'NamingStyle	' + str(nameStyle),
                  'StackExtension	' + stackExtension[1:],
                  'BinningOfImages ' + str(binning),
                  'RootNameOfDataFiles ' + dataName,
                  'ThicknessToMake ' + str(thickness),
                  'ShiftInY ' + str(fidIncShift)]
      if binning != aliBinning:
         comlines.append('Use3dfindAliInput 1')
      comlines += laterComDirectives(0)

      try:
         runcmd('makecomfile -StandardInput', comlines)
      except ImodpyError:
         reportImodError('Error making ' + comfile)
         return 1

      # Make the reconstruction
      if splitAndRunTilt(comfile, 'Making tomogram for finding beads'):
         return 1

      # Find the beads
      comlines = ['RootNameOfDataFiles ' + dataName,
                  'BinningOfImages ' + str(binning),
                  'BeadSize ' + str(fidSizePix),
                  fmtstr('OneParameterChange {}findbeads3d{}.findbeads3d.' + \
                            'StorageThreshold=-1', comPrefix, axisLet)]

      btlines = readTextFileReportErr('track' + axisCom)
      if not btlines:
         return 1
      if optionValue(btlines, 'LightBeads', BOOL_VALUE):
         comlines.append(fmtstr('OneParameterChange {}findbeads3d{}.findbeads3d.' + \
                                   'LightBeads=1', comPrefix, axisLet))
      comlines += laterComDirectives(0)
      if makeAndRunOneCom(comlines, 'findbeads3d' + axisCom, 'Finding beads in tomogram'):
         return 1

      sedcom = edfDelAndAdd('AlignedStack.eraseGold', 2)
      modifyEdfLines(sedcom)

   return 0


# Erase gold in aligned stack one way or another
def eraseGoldInAlignedStack():
   if eraseGold > 1:

      # Get the reprojection if using 3D method
      comlines = ['InputFile tilt_3dfind' + axisCom,
                  'RootNameOfDataFiles ' + dataName]
      if makeAndRunOneCom(comlines, 'tilt_3dfind_reproject' + axisCom,
                          'Reprojecting model of beads in tomogram onto aligned stack'):
         return 1

   # Or transform fiducial model
   else:
      try:
         filledIn = ''
         if os.path.exists(dataName + '_nogaps.fid'):
            filledIn = '_nogaps'
         runcmd(fmtstr('xfmodel -xf {0}.tltxf {0}{1}.fid {0}_erase.fid', dataName,
                       filledIn))
      except ImodpyError:
         reportImodError('Could not transform fiducials for erasing gold')
         return 1

   # Erase the beads
   extraDiam = lookupDirective(runtimePrefix + 'GoldErasing', 'extraDiameter', 0,
                               FLOAT_VALUE)
   if not extraDiam:
      extraDiam = 0.
   comlines = ['RootNameOfDataFiles ' + dataName,
                  'BeadSize ' + str(expandFactor * fidSizePix / aliBinning + extraDiam)]
   comlines += laterComDirectives(0)
   if lookupDirective(comPrefix + 'golderaser', 'ccderaser.ExpandCircleIterations', 0,
                      STRING_VALUE) == None:
      comlines.append(fmtstr('OneParameterChange {}golderaser{}.ccderaser.' + \
                                'ExpandCircleIterations=2', comPrefix, axisLet))
   if makeAndRunOneCom(comlines, 'golderaser' + axisCom,
                       'Erasing beads from aligned stack', skipRun = skipAlignStackMods):
      return 1
   if skipAlignStackMods:
      return 0
   
   if useFileAsReplacement(datasetFilename('_erase.ali'), datasetFilename('.ali'), False,
                           True):
      return 1

   sedcom = edfDelAndAdd('AlignedStack.eraseGold', min(2, eraseGold))
   modifyEdfLines(sedcom)
   return 0


# 2D filtering, very simple operation
def filterAlignedStack():
   thisStep = 13
   replaced = replaceOrRunAfterStep(thisStep, False)
   if replaced:
      return max(0, replaced)
   if filterIn2D and needStep(thisStep):
      filtLines = readTextFileReportErr('mtffilter' + axisCom)
      if not filtLines:
         return 1
      sedcom = [sedModify('PixelSize', aliBinning * pixelSize)]
      if modifyWriteAndRunCom('mtffilter' + axisCom, sedcom, filtLines,
                              '2D filtering the aligned stack with Mtffilter', 
                              skipRun = skipAlignStackMods):
         return 1
      if skipAlignStackMods:
         return 0
      if useFileAsReplacement(datasetFilename('_filt.ali'), datasetFilename('.ali'), 
                              False, True):
         return 1
      return replaceOrRunAfterStep(thisStep, True)

   return 0


# Change many entries in tilt.com for the current situation
def modifyTiltComFile(sampleThickness):
   global xtiltNeeded
   
   comfile = 'tilt' + axisCom
   comlines = readTextFileReportErr(comfile)
   if not comlines:
      return 1

   if sampleThickness:
      thickness = sampleThickness

   else:
      
      # Get different variants on thickness entry and insist on only one
      thickness = lookupDirective(comPrefix + 'tilt', 'tilt.THICKNESS', 0, INT_VALUE)
      binnedThick = lookupDirective(runtimePrefix + 'Reconstruction', 'binnedThickness',
                                    0, INT_VALUE)
      fallbackThick = lookupDirective(runtimePrefix + 'Reconstruction',
                                      'fallbackThickness', 0, INT_VALUE)
      if testDirectiveValue(thickness, 'tilt.THICKNESS', 'integer') or \
             testDirectiveValue(binnedThick, 'Reconstruction.binnedThickness',
                                'integer'):
         return 1
      if testDirectiveValue(fallbackThick,  'Reconstruction.fallbackThickness',
                            'integer'):
         return 1
             
      if thickness and binnedThick:
         abortSet('Both tilt.THICKNESS and Reconstruction.binnedThickness were entered')
         return 1
      if thickness and fallbackThick:
         abortSet('Both tilt.THICKNESS and Reconstruction.fallbackThickness were entered')
         return 1
      if fallbackThick and binnedThick:
         abortSet('Both Reconstruction.fallbackThickness and ' +\
                     'Reconstruction.binnedThickness were entered')
         return 1

      if binnedThick:
         thickness = aliBinning * binnedThick

      # If that did not provide a thickness, set it to fallback if there is nothing from
      # the alignment; otherwise take it from the alignment, with possible extra amount
      # but use the fallback if this is too thin
      if not thickness:

         # First look for a thickness from tomopitch
         # Parse this log afresh; results from align log are already in globals
         angleArr = [0, 0., 0., 0.]
         err = parseTomopitchLog(angleArr)
         if err < 0:
            return 1
         if not err:
            thickness = angleArr[0]
            xtiltNeeded = angleArr[1]

         # Or use the align thickness if no positioning available
         elif reconThickness:
            thickness = 2 * ((int(round(reconThickness)) + 1) // 2)

         # In either case, add the extra thickness if defined
         if thickness:
            extra = lookupDirective(runtimePrefix + 'Reconstruction', 'extraThickness', 0,
                                    INT_VALUE)
            if testDirectiveValue(extra, 'Reconstruction.extraThickness', 'integer'):
               return 1
            if extra:
               thickness += extra
               
         # In any of these cases, fallback rules if thickness is too low or not set
         if fallbackThick and (not thickness or
                               thickness < fallbackThick * useFallbackRatio):
            warning(fmtstr('Using fallback thickness of {} because computed '
                           'thickness, {}, is less than {} of fallback', fallbackThick,
                           thickness, useFallbackRatio))
            thickness = fallbackThick

         if not thickness:
            abortSet('No thickness was specified and neither fiducial alignment nor ' + \
                        'positioning gave a thickness to use')
            return 1

                  
   if ifMontage:
      fullx = montFrameData[6]
      fully = montFrameData[7]
      sssx = montFrameData[8]
      sssy = montFrameData[9]
   else:
      fullx = rawXsize
      fully = rawYsize
      if transposeForAli:
         (fullx, fully) = (fully, fullx)
      sssx = -((aliXunbinned - fullx) // 2)
      sssy = -((aliYunbinned - fully) // 2)

   gpuVal = -1
   if useGPU:
      gpuVal = 0

   # Output 0 X axis tilt if that is desired
   if noXaxisTilt:
      xtiltNeeded = 0.
   sedcom = [sedModify('IMAGEBINNED', aliBinning),
             sedModify('XAXISTILT', xtiltNeeded),
             sedModify('FULLIMAGE', fmtstr('{} {}', fullx, fully)),
             sedModify('SUBSETSTART', fmtstr('{} {}', sssx, sssy)),
             sedModify('THICKNESS', thickness)] + \
             sedDelAndAdd('UseGPU', gpuVal, 'XTILTFILE')
   if not axisNum:
      sedcom.append(sedModify('OutputFile', datasetFilename('_full.rec')))
   if useGPU:
      sedcom += sedDelAndAdd('ActionIfGPUFails', '1,2', 'XTILTFILE');
   else:
      sedcom.append('/^ActionIfGPUFails/d')
   if didLocalAlign:
      sedcom += sedDelAndAdd('LOCALFILE', dataName + 'local.xf', 'XTILTFILE')
   else:
      sedcom.append('/^LOCALFILE/d')
   if madeZfactors:
      sedcom += sedDelAndAdd('ZFACTORFILE', dataName + '.zfac', 'XTILTFILE')
   else:
      sedcom.append('/^ZFACTORFILE/d')

   # Handle change in scaling if they turned off log and did not supply a scale
   logbase = optionValue(comlines, 'LOG', FLOAT_VALUE)
   if not logbase and not lookupDirective(comPrefix + 'tilt', 'tilt.SCALE', 0,
                                          STRING_VALUE):
      scaleArr = optionValue(comlines, 'SCALE', FLOAT_VALUE)
      if not scaleArr or len(scaleArr) < 2:
         abortSet('Cannot modify SCALE value in tilt' + axisLet + comExt + \
                  ' for linear ' + 'scaling')
         return 1

      # Copytomocoms produces scales from 1000 to 40 depending on X size.
      # Linear requires scale to be reduced by 5000, giving numbers from 0.2 to 0.008
      # 3 is a safe dividing point for deciding whether this has already happened
      if scaleArr[1] > 3.:
         sedcom.append(sedModify('SCALE', fmtstr('{} {:.3f}',
                                                 scaleArr[0], scaleArr[1] / 5000.)))

   if pysed(sedcom, comlines, comfile, retErr=True):
      abortSet('Error modifying ' + comfile)
      return 1
   return 0
   

# Make 3D CTF corrected tomogram
def make3dCtfCorrectedTomogram(tiltCom, numProc, recName):
   if doBothRecons and not rawFor3dCTF and makeAlignedStack(0):
      return 1

   comlines = ['InputFile ' + tiltCom,
               'ThicknessToMake ' + str(ctf3dSlabThick)]
   if not lookupDirective(comPrefix + 'ctf3dsetup', 'ctf3dsetup.RunSlabsInParallel', 0,
                          BOOL_VALUE):
      comlines.append('NumberOfProcessors ' + str(numProc))
   comlines += laterComDirectives(0)
   if eraseGold:
      comlines.append(fmtstr('OneParameterChange {}ctf3dsetup{}.ctf3dsetup.' + \
                             'EraseFiducials=1', comPrefix, axisLet))
   if lookupDirective(runtimePrefix + 'AlignedStack', 'filterStack', 0, BOOL_VALUE):
      comlines.append(fmtstr('OneParameterChange {}ctf3dsetup{}.ctf3dsetup.' + \
                             'FilterIn2D=1', comPrefix, axisLet))
   if (makeAndRunOneCom(comlines, 'ctf3dsetup' + axisCom,
                    'Making command files to make 3D CTF-corrected tomogram')):
      return 1
   if runOneProcess('ctf3d' + comExt, False, useGPU, message =
                    'Making 3D CTF corrected tomogram'):
      return 1

   ctfName = datasetFilename('_3dctf.rec')
   if renameAndAbort(ctfName, recName,
                     'Renaming {} to {} after make CTF corrected reconstruction'):
      return 1


# Make tomogram by backprojection and/or SIRT
def generateTomogram():
   recRoot = setName + axisLet
   if not dualAxis:
      recRoot = setName + '_full'
   recName = datasetFilename('.rec', root = recRoot)
   (doSIRT, doFakeSIRT, doBoth, SIRTname, fakeName) = getSIRTrecName(recRoot)
   if doSIRT and doFakeSIRT:
      abortSet('You cannot do both SIRT and a SIRT-like filter')
      return 1

   tiltCom = 'tilt' + axisCom

   # This is needed for both SIRT and 3D CTF
   numProc = numberOfProcessingUnits()
   if numProc < 1:
      return 1

   if doBoth or (not doSIRT and not ctf3dSlabThick):
      if splitAndRunTilt(tiltCom, 'Making tomogram by back-projection', numProc):
         releaseGPUallocation()
         return 1
   if doFakeSIRT and doBoth:
      tempComName = 'tilt' + axisLet + '_slf' + axisCom
      if renameAndAbort(recName, fakeName, 
                        'Renaming {} to {} before regular back-projection'):
         releaseGPUallocation()
         return 1
      if renameAndAbort(tiltCom, tempComName):
         releaseGPUallocation()
         return 1
      if pysed(['/FakeSIRTiterations/d'], tempComName, tiltCom, retErr=True):
         abortSet('Error modifying ' + tempComName)
         releaseGPUallocation()
         return 1

      if ctf3dSlabThick:
         if make3dCtfCorrectedTomogram(tiltCom, numProc, recName):
            releaseGPUallocation()
            return 1

      elif splitAndRunTilt(tiltCom, 'Making tomogram by regular back-projection', 
                           numProc):
         releaseGPUallocation()
         return 1

      makeBackupFile(tiltCom)
      if renameAndAbort(tempComName, tiltCom):
         releaseGPUallocation()
         return 1

   if doSIRT:
      comlines = laterComDirectives(0)
      comlines.append(fmtstr('OneParameterChange {}sirtsetup{}.sirtsetup.' + \
                                'NumberOfProcessors={}', comPrefix, axisLet, numProc))
      if makeAndRunOneCom(comlines, 'sirtsetup' + axisCom,
                          'Making command files to run SIRT'):
         releaseGPUallocation()
         return 1
      if runOneProcess('tilt' + axisLet + '_sirt' + comExt, False, useGPU, message =
                       'Making tomogram with SIRT'):
         releaseGPUallocation()
         return 1
         
   if not (doFakeSIRT and doBoth) and ctf3dSlabThick and \
      make3dCtfCorrectedTomogram(tiltCom, numProc, recName):
      releaseGPUallocation()
      return 1

   # Put the correct size of the reconstruction in the edf file when finished
   releaseGPUallocation()
   if doSIRT and not doBoth:
      recName = SIRTname
   try:
      (recx, recy, recz) = getmrcsize(recName)
   except ImodpyError:
      reportImodError('Could not get size of reconstruction ' + recName)
      return 1

   prefix = axisUpperLet + '.tomogramSize.'
   sedcom = edfDelAndAdd(prefix + 'Columns', recx) + \
            edfDelAndAdd(prefix + 'Rows', recy) + \
            edfDelAndAdd(prefix + 'Sections', recz)
   if modifyEdfLines(sedcom):
      return 1
   return 0


# Run the final tiltalign and maintain alignLines
def runFinalTiltalign(sedcom, mess):
   global alignLines
   hasSkip = os.getenv('TILTALIGN_SKIP_CROSS_VAL')
   if not hasSkip:
      os.environ['TILTALIGN_SKIP_CROSS_VAL'] = '1'
   err = modifyWriteAndRunCom('align' + axisCom, sedcom, alignLines,
                           'Doing final alignment ' + mess)
   if not hasSkip:
      del os.environ['TILTALIGN_SKIP_CROSS_VAL']
   if err:
      return 1
   alignLines =  pysed(sedcom, alignLines, retErr=True)
   if isinstance(alignLines, str):
      abortSet('Error modifying alignment lines: ' + str)
      return 1
   if postProcessTiltalign():
      return 1


# Run positioning (with whole tomogram for now)
def positionTomogram():
   global suppressAbort, alignLines
   numScales = [4, 2, 2, 1]
   boxSizes = [32, 20, 16, 12]
   blockSizes = [100, 60, 48, 36]
   minCenterFid = 4

   # If no positioning and conditions are satisfied to center on gold, get align results
   # and run align if there are enough on each surface
   if posSampleType <= 0:
      if not centerOnGold or numSurfaces < 2:
         return 0
      angleArr = [0., 0., 0., 0., 0., 0., 0, 0]
      if analyzeAlignLog(True, angleArr, noLogOK = True):
         return 1
      if angleArr[6] < minCenterFid or angleArr[7] < minCenterFid:
         return 0
      sedcom = [sedModify('AxisZShift', angleArr[4])]
      if runFinalTiltalign(sedcom, 'to center on gold'):
         return 1

   # Get thickness
   thickness = lookupDirective(runtimePrefix + 'Positioning', 'thickness', 0, INT_VALUE)
   if testDirectiveValue(thickness, 'Positioning.thickness', 'integer'):
      return 1

   # Cryosample
   if posSampleType > 1:
      if not thickness:
         abortSet('Thickness must be entered for cryopositioning')
         return 1
      if manageGPUallocation():
         return 1;

      findGold = 2
      if patchTrack or fiducialless:
         findGold = 0
         if lookupDirective(runtimePrefix + 'Positioning', 'hasGoldBeads', 0, BOOL_VALUE):
            if fidSizePix <= 0.:
               abortSet('Gold bead size must be entered to find beads in cryopositioning')
            findGold = 2

      comlines = ['RootNameOfDataFiles ' + dataName,
                  'ThicknessToMake ' + str(thickness),
                  'FindBeadsInVolume ' + str(findGold)]
      if gpuList:
         comlines.append('UseGPU 0')

      suppressAbort = True
      err = makeAndRunOneCom(comlines, 'cryoposition' + axisCom, 
                             'Positioning tomogram with Cryoposition',
                             useGPU = gpuList != '')
      suppressAbort = False
      cplines = readTextFileReportErr('cryoposition' + axisLet + '.log')
      if not cplines:
         return 1
      if err:
         warning('Tomogram positioning did not work; proceeding if possible')
         return 0

      aliBinning = findTaggedValue(cplines, 'stack with binning', '=', INT_VALUE)
      if not aliBinning:
         abortSet('Cannot find binning in cryoposition' + axisLet + '.log')
         return 1
      
   else:

      # Plastic section: Get binning and make aligned stack
      size = max(rawXsize, rawYsize)
      aliBinning = lookupDirective(runtimePrefix + 'Positioning', 'binByFactor', 0,
                                   INT_VALUE)
      if testDirectiveValue(aliBinning, 'Positioning.binByFactor', 'integer'):
         return 1
      if not aliBinning:
         for bin in (1, 2, 3, 4):
            if size / bin < positionBinningTarget or bin == 4:
               aliBinning = bin
               prnLog('Binning for positioning tomogram set to ' + str(aliBinning))
               break
         
      if makeAlignedStack(aliBinning):
         return 1

      # Set default thickness if needed and make tomogram
      if not thickness:
         for ind in range(len(dfltPosThicknesses) - 1, -1, -1):
            if size > sizesForPosThicknesses[ind]:
               thickness = dfltPosThicknesses[ind]
               prnLog('Thickness for positioning tomogram set to ' + str(thickness))
               break
            
      if modifyTiltComFile(thickness) or \
             splitAndRunTilt('tilt' + axisCom, 'Making tomogram for positioning'):
         return 1

      # Run findsection to get the model
      filename = datasetFilename('.rec')
      if not axisNum:
         filename = datasetFilename('_full.rec')
      binInd = min(aliBinning, len(boxSizes)) - 1
      err = runFindSection(filename, numScales[binInd], boxSizes[binInd],
                           pitchMod = 'tomopitch' + axisLet + '.mod',
                           comRoot = 'findsection_pos' + axisLet,
                           block = blockSizes[binInd], failureOK = True)

      # Distinguish an OK findsection failure from other errors that issued ABORT SET
      if err < 0:
         warning('Finding the section for positioning did not work; proceeding if ' +\
                    'possible')
         return 0
      if err > 0:
         return 1

   # Get existing Z shift, angle offset, and X tilt
   tiltLines = readTextFileReportErr('tilt' + axisCom)
   if not tiltLines:
      return 1
   zShiftOrig = 0.
   angleOffsetOrig = 0.
   if fiducialless:
      zShiftArr = optionValue(tiltLines, 'SHIFT', FLOAT_VALUE, numVal = 2)
      if zShiftArr and len(zShiftArr) > 1:
         zShiftOrig = zShiftArr[1]
      offset = optionValue(tiltLines, 'OFFSET', FLOAT_VALUE, numVal = 1)

   else:
      zShift = optionValue(alignLines, 'AxisZShift', FLOAT_VALUE, numVal = 1)
      if zShift:
         zShiftOrig = zShift
      offset = optionValue(alignLines, 'AngleOffset', FLOAT_VALUE, numVal = 1)

   if offset:
      angleOffsetOrig = offset
   xtiltOrig = 0.
   xtilt = optionValue(tiltLines, 'XAXISTILT', FLOAT_VALUE, numVal = 1)
   if noXaxisTilt:
      xtilt = 0
   if xtilt and posSampleType == 1:
      xtiltOrig = xtilt

   edfcom = edfDelAndAdd(axisEdfLet + '.sample.AngleOffset', round(angleOffsetOrig, 3))
   edfcom += edfDelAndAdd(axisEdfLet + '.sample.AxisZShift', round(zShiftOrig, 2))
   edfcom += edfDelAndAdd(axisEdfLet + '.sample.XAXISTILT', round(xtiltOrig, 3))
   edfcom += edfDelAndAdd(axisEdfLet + '.positioning.PosSampleType', posSampleType)
   if modifyEdfLines(edfcom):
      return 1

   # Modify tomopitch and run it
   sedcom = sedDelAndAdd('ScaleFactor', aliBinning, 'ModelFile') + \
       sedDelAndAdd('AngleOffsetOld', angleOffsetOrig, 'ModelFile') + \
       sedDelAndAdd('ZShiftOld', zShiftOrig, 'ModelFile') + \
       sedDelAndAdd('XAxisTiltOld', xtiltOrig, 'ModelFile')
   if posSampleType > 1 and not lookupDirective(comPrefix + 'tomopitch',
                                             'tomopitch.ExtraThickness', 0, INT_VALUE):
      sedcom.append(sedModify('ExtraThickness', cryoPosExtraThick))

   suppressAbort = True
   err = modifyWriteAndRunCom('tomopitch' + axisCom, sedcom,
                              message = 'Finding angles and thickness')
   suppressAbort = False
   if err:
      warning('Tomopitch failed; proceeding if possible')
      return 1;
   
   # Parse this log
   angleArr = [0, 0., 0., 0.]
   err = parseTomopitchLog(angleArr)
   if err < 0:
      return 1
   (thickness, xtilt, offset, zShift) = angleArr
   if err or not thickness:
      prnLog('Tomopitch did not produce a good result; proceeding if possible')
      return 0
   
   # Write tilt.com with the X-axis tilt and shift/angle offset for fidless
   prnLog(fmtstr('Positioning gave thickness {}, X-axis tilt {:.2f}, angle offset ' + \
                    '{:.2f}, Z shift {:.1f}', thickness, xtilt, offset, zShift))
   sedcom = [sedModify('XAXISTILT', xtilt), 
             sedModify('THICKNESS', thickness)]
   if fiducialless:
      sedcom += sedDelAndAdd('OFFSET', offset, 'OutputFile') + \
          sedDelAndAdd('SHIFT', '0. ' + str(zShift), 'OutputFile')
   if pysed(sedcom, tiltLines, 'tilt' + axisCom, retErr=True):
      abortSet('Error modifying tilt' + axisCom)
      return 1

   # Write align.com with Z shift and offset and run it
   if not fiducialless:
      sedcom = [sedModify('AxisZShift', zShift),
                sedModify('AngleOffset', offset)]
      if runFinalTiltalign(sedcom, 'with new position'):
         return 1

   return 0
      

# COMBINE FUNCTIONS
#
# Get the final patch size and extra targets for autopatchfid
def getAutoPatchfitParams():
   extraTargets = lookupDirective(combinePrefix, 'extraTargets', 0, STRING_VALUE)
   finalPatchSize = lookupDirective(combinePrefix, 'finalPatchSize', 0, STRING_VALUE)
   if not extraTargets:
      extraTargets = dfltExtraWarpTargets
   if not finalPatchSize:
      finalPatchSize = dfltFinalPatchSize
   return (finalPatchSize, extraTargets)


# Analyze the tomogram for Z limits and run setupcombine
def setupCombine():
   global axisLet, axisNum, useVolMatch, nxRecA, nyRecA
   topBotInd = 2
   patchSize = 'M'
   useVolMatch = 0
   edfVol = 'false'
   topBotA = 6 * [0]
   topBotB = 6 * [0]
   if fiducialless or patchTrack or not os.path.exists('transferfid.coord'):
      useVolMatch = 1
      edfVol = 'true'

   # Get the right files in place if SIRT is involved
   useSIRTifBoth = lookupDirective(combinePrefix, 'doSIRTifBoth', 0, INT_VALUE)
   if testDirectiveValue(useSIRTifBoth, 'Combine.doSIRTifBoth', 'integer'):
      return 1
   recRoot = setName + 'a'
   for axisNum in (1, 2):
      recName = datasetFilename('.rec', root = recRoot)
      (doSIRT, doFakeSIRT, doBoth, SIRTname, fakeName) = getSIRTrecName(recRoot)
      useSIRT = (doSIRT and (useSIRTifBoth or not doBoth)) or \
                (doFakeSIRT and doBoth and useSIRTifBoth)
      if useSIRT:
         useSIRTname = SIRTname
         likeText = ''
         if doFakeSIRT:
            useSIRTname = fakeName
            likeText = '-like'
         SIRTexists = os.path.exists(useSIRTname)
         recExists = os.path.exists(recName)
         bpName = datasetFilename('_BP.rec', root = recRoot)
         if recExists and not SIRTexists:
            prnLog('The file ' + recName + ' exists and the SIRT' + likeText + \
                   ' reconstruction ' + useSIRTname + \
                   ' does not; assuming it was already renamed')
         elif not SIRTexists:
            abortSet('Neither the SIRT' + likeText + ' reconstruction ' + useSIRTname + \
                     ' nor the file ' + recName + ' exist; cannot proceed')
            return 1
         else:
            if recExists and doBoth:
               if useFileAsReplacement(recName, bpName, False, True):
                  return 1
               prnLog('Renamed back-projection reconstruction ' + recName + ' to: ' + \
                         bpName)
            if useFileAsReplacement(useSIRTname, recName, False, True):
               return 1
            prnLog('Renamed SIRT' + likeText + ' reconstruction ' + useSIRTname + \
                   ' to: ' + recName)

      recRoot = setName + 'b'

   axisNum = 0

   # Get default top and bottom Z limits in case of failure
   try:
      (nx, zRangeA, nz) = getmrcsize(datasetFilename('a.rec'))
      topBotA[topBotInd] = 1
      topBotA[topBotInd + 1] = zRangeA
      (nx, zRangeB, nz) = getmrcsize(datasetFilename('b.rec'))
      topBotB[topBotInd] = 1
      topBotB[topBotInd + 1] = zRangeB
   except ImodpyError:
      reportImodError('Error running header on reconstruction files')
      return 1

   # Get the directives and use the defaults if none
   matchAtoBratio = lookupDirective(combinePrefix, 'matchAtoBThickRatio', 0, FLOAT_VALUE)
   if testDirectiveValue(matchAtoBratio, 'Combine.matchAtoBThickRatio', 'float'):
      return 1
   if not matchAtoBratio:
      matchAtoBratio = dfltMatchAtoBratio
   numScales = lookupDirective(combinePrefix, 'findSecNumScales', 0, INT_VALUE)
   if testDirectiveValue(numScales, 'Combine.findSecNumScales', 'integer'):
      return 1
   if not numScales:
      numScales = dfltCombineNumScales
   boxSize = lookupDirective(combinePrefix, 'findSecBoxSize', 0, INT_VALUE)
   if testDirectiveValue(boxSize, 'Combine.findSecBoxSize', 'integer'):
      return 1
   if not boxSize:
      boxSize = dfltCombineBoxSize

   # Need to find limits in A unless match ratio is very high
   needFindA = matchAtoBratio < 4.
   if needFindA:
      errA = runFindSection(datasetFilename('a.rec'), numScales, boxSize,
                            topBots = topBotA, comRoot = 'findsection_zlima',
                            failureOK = True)
      if errA > 0:
         return 1
      if not errA:
         zRangeA = topBotA[topBotInd + 1] - topBotA[topBotInd]
         if zRangeA < 4:
            abortSet('Findsection failed to find a useful Z range in A tomogram')

   # Need to find limits in B unless match ratio is very low
   needFindB = matchAtoBratio > 0.25
   if needFindB:
      errB = runFindSection(datasetFilename('b.rec'), numScales, boxSize, 
                            topBots = topBotB, comRoot = 'findsection_zlimb', 
                            failureOK = True)
      if errB > 0:
         return 1
      if not errB:
         zRangeB = topBotB[topBotInd + 1] - topBotB[topBotInd]
         if zRangeB < 4:
            abortSet('Findsection failed to find a useful Z range in B tomogram')

   # Match A to B if ratio was very high or B is sufficiently smaller
   matchAtoB = not needFindA or (needFindB and not errA and not errB and \
                                    zRangeB <= matchAtoBratio * zRangeA)
   edfMatch = 'B_to_A'
   if matchAtoB:
      topBotA = topBotB
      edfMatch = 'A_to_B'
      errA = errB

   # Give warning and proceed if the matched to tomogram had error in Z limits
   if errA:
      warning('Findsection failed to find Z limits for combine, proceeding anyway')
      
   # Determine number of surfaces, make sure 2 is still OK
   # If in doubt, just set to 1 and see what happens
   numSurf = numSurfaces
   if useVolMatch:
      numSurf = 2
   else:
   
      # Do the appropriate axis of align log
      if matchAtoB:
         axisLet = 'b'
      angleArr = [0., 0., 0., 0., 0., 0., 0, 0]
      if analyzeAlignLog(numSurfaces > 1, angleArr, noLogOK = True):
         return 1
      axisLet = 'a'
      (totalDelTilt, xtiltNeeded, fidThickness, fidIncShift, fidTotalShift, \
          reconThickness, numBot, numTop) = angleArr
      if numBot + numTop < solvematchMinFids:
         useVolMatch = 1
         prnLog('Using Dualvolmatch instead of Solvematch because there are too few ' +\
                   'fiducials')
      elif numSurf > 1:
         if fidThickness < (topBotA[topBotInd + 1] - topBotA[topBotInd]) / 2:
            numSurf = 1
            mess = 'because fiducial extent is much smaller than Z range'
         else:
            fallbackThick = lookupDirective(runtimePrefix + 'Reconstruction',
                                            'fallbackThickness', 0, INT_VALUE)
            if isinstance(fallbackThick, int) and \
                   reconThickness < fallbackThick * useFallbackRatio:
               numSurf = 1
               mess = 'because fiducial extent is much smaller than fallback thickness'

         if numSurf > 1 and (numBot < solvematchMinEachSide or \
                                numTop < solvematchMinEachSide or \
                                min(float(numBot) / numTop, float(numTop) / numTop) < \
                                solvematchMinSideRatio):
            numSurf = 1
            mess = 'because there are too few fiducials on one surface'

         if numSurf == 1:
            prnLog('Setting number of surfaces for Solvematch to 1 ' + mess)

   edfSurf = 'BothSides'
   if numSurf == 1:
      edfSurf = 'OneSide'
      
   # Get the rest of the combine setup directives
   # TODO: what if there is a blank directive to cancel extra targets?
   patchIn = lookupDirective(combinePrefix, 'patchSize', 0, STRING_VALUE)
   if patchIn:
      patchSize = patchIn
   (finalPatchSize, extraTargets) = getAutoPatchfitParams()
   
   wedgeFrac = lookupDirective(combinePrefix, 'wedgeReduction', 0, FLOAT_VALUE)
   lowRadius = lookupDirective(combinePrefix, 'lowFromBothRadius', 0, FLOAT_VALUE)
   if testDirectiveValue(wedgeFrac, 'Combine.wedgeReduction', 'float') or \
          testDirectiveValue(lowRadius, 'Combine.lowFromBothRadius', 'float'):
      return 1
   comLines = ['RootName ' + setName,
               'WarningsToStandardOut 1',
               'NamingStyle ' + str(nameStyle),
               'StackExtension ' + stackExtension[1:],
               'TransferPointFile transferfid.coord',
               'SurfaceModelType ' + str(numSurf),
               fmtstr('ZLowerAndUpper {},{}', topBotA[topBotInd], topBotA[topBotInd + 1]),
               'PatchTypeOrXYZ ' + patchSize,
               'AutoPatchFinalSize ' + finalPatchSize]
   if extraTargets:
      comLines.append('ExtraResidualTargets ' + extraTargets)
   if useVolMatch:
      comLines.append('InitialVolumeMatching 1')
   comLines += laterComDirectives(0)

   if matchAtoB:
      comLines.append('MatchAtoB 1')
   if wedgeFrac:
      comLines.append('WedgeReductionFraction ' + str(wedgeFrac))
   if lowRadius:
      comLines.append('LowFromBothRadius ' + str(lowRadius))

   try:
      runcmd('setupcombine -StandardInput', comLines, 'stdout')
   except ImodpyError:
      reportImodError('Error running setupcombine')
      return 1

   # Now try to modify the edf file.  Z limits are fatal.  First word -> batchruntomo
   sedcom = edfDelAndAdd('Combine.MatchMode', edfMatch) + \
       edfDelAndAdd('Combine.FiducialMatch', edfSurf) + \
       edfDelAndAdd('Combine.InitialVolumeMatching', edfVol) + \
       edfDelAndAdd('Combine.PatchSize', patchSize) + \
       edfDelAndAdd('Combine.FinalPatchSize', finalPatchSize) + \
       edfDelAndAdd('Combine.ExtraResidualTargets', extraTargets)

   if modifyEdfLines(sedcom):
      return 1

   # Finally, modify matchvol1.com to preserve a larger B volume size
   # The transposed ny,nz match setupcombine usage for better or worse
   try:
      (nxRecA, nza, nyRecA) = getmrcsize(datasetFilename('a.rec'))
      (nxb, nzb, nyb) = getmrcsize(datasetFilename('b.rec'))

   except ImodpyError:
      reportImodError('Error getting sizes of axis reconstruction files')
      return 1

   sedcom = [sedModify('OutputSizeXYZ', fmtstr('{} {} {}', nxRecA, max(nza, nzb),
                                               nyRecA))]
   comFile = 'matchvol1' + comExt
   if pysed(sedcom, comFile, comFile, retErr = True):
      abortSet('Error modifying ' + comFile)
      return 1
      
   return 0


# Run solvematch and make needed modifications if there is a center shift
def initialCombineMatch():
   global suppressAbort, useVolMatch

   # If we are coming in at this step, need to reconstruct whether volume matching
   # is supposed to be used
   if useVolMatch < 0:
      edfLines = readTextFileReportErr(setName + '.edf')
      if not edfLines:
         return 1
      for line in edfLines:
         if 'InitialVolumeMatching' in line:
            if 'true' in line:
               useVolMatch = 1
            elif 'false' in line:
               useVolMatch = 0
            break
      if useVolMatch < 0:
         abortSet('Cannot determine if initial volume matching was set up to be used')
         return 1

   # Set up for type of process and whether may need to run it twice
   if useVolMatch:
      comRoot = 'dualvolmatch'
      numLoop = 1
      tags = [('unbinned mean residual', 4), ('implies a center', 4), ('Falling back', 4)]
   else:
      comRoot = 'solvematch'
      numLoop = 2
      tags = [('Mean residual', 4), ('Scaling along', 0), ('Local fits', 0),
              ('Average mean', 4)]

   for loop in range(numLoop):

      # Run the process and process and read log regardless
      suppressAbort = True
      err = runOneProcess(comRoot + comExt, message = 'Getting initial alignment ' +\
                             'between volumes')
      suppressAbort = False

      # Warnings got printed regardless, errors got printed if it failed, so set tags to
      # avoid duplicate messages
      useTags = [('ERROR:', 0)] + tags
      if err:
         useTags = tags
      printTaggedMessages(comRoot + '.log', useTags);
      logLines = readTextFileReportErr(comRoot + '.log')
      if not logLines:
         return 1

      # For vol match, see if need to change the matchvol thickness
      if useVolMatch and not err:
         for line in logLines:
            if 'may need to set thickness' in line:
               lsplit = line.split()
               try:
                  newThick = int(lsplit[-1])
               except Exception:
                  abortSet('Error trying to get suggested thickness for matchvol')

               sedcom = [sedModify('OutputSizeXYZ', fmtstr('{} {} {}', nxRecA,
                                                              newThick, nyRecA))]
               comFile = 'matchvol1' + comExt
               if pysed(sedcom, comFile, comFile, retErr = True):
                  abortSet('Error modifying' + comFile)
                  return 1
               prnLog(fmtstr('Thickness for matched volume set to {} as suggested',
                             newThick))
               break

      # For solvematch look first to see if it suggests using one surface, modify the
      # com and loop for another run
      if not useVolMatch:
         redo = False
         for line in logLines:
            if 'Try specifying' in line and 'on one surface' in line:
               comFile = 'solvematch' + comExt
               if pysed([sedModify('SurfacesOrUseModels', 1)], comFile,
                        comFile, retErr = True):
                  abortSet('Error modifying' + comFile)
                  return 1
               prnLog('Rerunning with one surface as suggested')
               redo = True
               break

         if redo:
            continue
            
      # Otherwise done if no error
      if not err:
         return 0
      initShift = []
      shiftLimit = 0

      # If there was an error, try to find initial shift as well as center shift
      for line in logLines:

         if 'InitialShiftXYZ' in line and 'needs' in line:
            lsplit = line.split()
            try:
               for ind in (0, 1, 2):
                  initShift.append(int(lsplit[-3 + ind]))
            except Exception:
               abortSet('Error trying to get initial shift from :' + line)
         if 'CenterShiftLimit' in line and 'avoid stopping' in line:
            lsplit = line.split()
            try:
               shiftLimit = int(lsplit[-1])
            except ValueError:
               pass

         # If error line is found that says it's OK, modify patchcorr.com
         if 'ERROR:' in line and (('INITIAL SHIFT' in line.upper() and \
                                      'SOLUTION IS OK' in line.upper()) or \
                                     'Initial shift needs' in line):
            if initShift:
               prnLog(fmtstr('Setting initial shift for Corrsearch3d to {} {} {}',
                             initShift[0], initShift[2], initShift[1]))
               sedcm = sedDelAndAdd('InitialShiftXYZ', fmtstr('{},{},{}', initShift[0],
                                                              initShift[1], initShift[2]),
                                    'FlipYZMessages')
               comFile = 'patchcorr' + comExt
               if pysed(sedcm, comFile, comFile, retErr = True):
                  abortSet('Error modifying ' + comFile)
                  return 1

               # If that succeeded, modify com file with shift limit
               if shiftLimit:
                  if pysed([sedModify('CenterShiftLimit', shiftLimit)], comRoot + comExt,
                           comRoot + comExt, retErr = True):
                     abortSet('Error modifying solvematch.com with new shift limit')
                     return 1
               return 0

            abortSet('Initial shift is required to go on, but was not found in ' +\
                        comRoot + '.log')
            return 1

   abortSet('Error getting initial alignment between volumes')
   return 1
   

# Run matchvol, autopatchfit, and volcombine
def alignAndCombineAxes():
   if needStep(17) and runOneProcess('matchvol1' + comExt, useMostCPUs = True,
                                     message = 'Making initial matching volume'):
      return 1

   if needStep(18):
      (finalPatchSize, extraTargets) = getAutoPatchfitParams()
      comLines = ['$autopatchfit -StandardInput',
                  'FinalPatchTypeOrXYZ ' + finalPatchSize]
      if extraTargets:
         comLines.append('ExtraResidualTargets ' + extraTargets)
      if writeTextFileReportErr('autopatchfit' + comExt, comLines):
         return 1
      err = runOneProcess('autopatchfit' + comExt, message = 'Doing patch correlation ' +\
                             'and fitting to local patches', useMostCPUs = True)
      printTaggedMessages('autopatchfit.log', standardTags +
                          [('Using ', 4), ('Changing ', 4), ('Adding ', 4),
                           ('FINDWARP -', 0), ('found a good', 4, logSuffixTag)])
      if err:
         abortSet('Cannot align the two tomograms')
         return 1

   if needStep(19):
      if parallelCPU > 1:
         try:
            runcmd('splitcombine')
         except ImodpyError:
            reportImodError('Error trying to run splitcombine to combine in parallel')
            return 1
            
      if runOneProcess('volcombine' + comExt, parallelCPU < 2,
                       message = 'Combining the two volumes'):
         return 1
      

# Run trimvol if doTrimvol is true; do not do it if it is present and false;
# if it is not present do it if ANY of the original directives are present
def trimVolume():
   trimPrefix = runtimePrefix + 'Trimvol'
   sizeKeys = ('sizeInX', 'sizeInY', 'thickness', 'scaleFromX', 'scaleFromY', 
               'scaleFromZ', 'findSecAddThickness')
   replaced = replaceOrRunAfterStep(trimStepNum, False)
   if replaced:
      return max(0, replaced)

   if not lookupDirective(runtimePrefix + 'Postprocess', 'doTrimvol', 0, BOOL_VALUE):
      if lookupDirective(runtimePrefix + 'Postprocess', 'doTrimvol', 0, STRING_VALUE):
         return 0
      if lookupDirective(trimPrefix, 'reorient', 0, INT_VALUE) == None:
         for key in sizeKeys:
            if lookupDirective(trimPrefix, key, 0, FLOAT_VALUE) != None:
               break
         else:    # ELSE ON FOR
            return 0

   if axisNum > 0 and not lookupDirective(trimPrefix, 'doAorBofDualAxis', 0, BOOL_VALUE):
      return 0

   useFindSec = lookupDirective(trimPrefix, 'findSecAddThickness', 0, FLOAT_VALUE) != None
   if lookupDirective(trimPrefix, 'thickness', 0, FLOAT_VALUE) != None and useFindSec:
      abortSet('Cannot use both Trimvol.findSecAddThickness and Trimvol.thickness' + \
                  ' directives')
      return 1
      
   reorient = lookupDirective(trimPrefix, 'reorient', 0, INT_VALUE)
   if reorient == None:
      reorient = 2
   if isinstance(reorient, str) or reorient < 0 or reorient > 2:
      abortSet('The value for "reorient" must be 0, 1, or 2')
      return 1
   edfcom = edfDelAndAdd('TrimvolFlipped', boolStringForEdf(reorient != 0))
   edfcom += edfDelAndAdd('Trimvol.SwapYZ', boolStringForEdf(reorient == 1))
   edfcom += edfDelAndAdd('Trimvol.RotateX', boolStringForEdf(reorient == 2))
   
   # Figure out what volume to trim and get its size
   recRoot = dataName
   trimNames = [datasetFilename('.rec')]
   fsName = 'findsection_trim'
   if dualAxis and axisNum > 0:
      trimNames = [datasetFilename('_trim.rec')]
      fsName += axisLet
   elif dualAxis:
      recRoot = 'sum'
   else:
      recRoot = dataName + '_full'
   recNames = [datasetFilename('.rec', root = recRoot)]

   if not dualAxis or axisNum > 0:
      (doSIRT, doFakeSIRT, doBoth, SIRTname, fakeName) = getSIRTrecName(recRoot)
      if doSIRT or doFakeSIRT:
         useSIRT = True
         if doBoth:
            trimBoth = lookupDirective(trimPrefix, 'doSIRTifBoth', 0, INT_VALUE)
            if testDirectiveValue(trimBoth, 'Trimvol.doSIRTifBoth', 'integer'):
               return 1
            if trimBoth == None:
               trimBoth = 0
            useSIRT = trimBoth > 0

         if useSIRT:
            if doSIRT and not SIRTname:
               abortSet('Cannot trim SIRT output, cannot find LeaveIterations directive')
               return 1

            useSIRTname = SIRTname
            if doFakeSIRT:
               useSIRTname = fakeName

            recNames = [useSIRTname]
            if doBoth and trimBoth > 1:
               (bpRoot, ext) = os.path.splitext(trimNames[0])
               recNames.append(datasetFilename('.rec', root = recRoot))
               trimNames.append(datasetFilename('_BP.rec', 
                                                root = bpRoot))

   for name in recNames:
      if not os.path.exists(name):
         abortSet('The file to be trimmed, ' + name + ', does not exist')
         return 1
   try:
      (nxrec, nzrec, nyrec) = getmrcsize(recNames[0])
   except ImodpyError:
      reportImodError('Could not get size of volume to trim')
      return 1

   edfcom += edfDelAndAdd('TrimVol.Input.NColumns', nxrec)
   edfcom += edfDelAndAdd('TrimVol.Input.NRows', nzrec)
   edfcom += edfDelAndAdd('TrimVol.Input.NSections', nyrec)
   
   # Determine size parameters, convert fractions, and check them
   sizesScales = []
   baseVals = (nxrec, nyrec, nzrec, nxrec, nyrec, nzrec, nzrec)
   for key in sizeKeys:
      val = lookupDirective(trimPrefix, key, 0, FLOAT_VALUE)
      if testDirectiveValue(val, 'Trimvol.' + key, 'float'):
         return 1
      size = None
      if val != None:
         base = baseVals[len(sizesScales)]
         lowLim = 40
         if len(sizesScales) % 3 == 2:
            lowLim = 4
         size = int(round(val))
         if size <= 1.:
            size = int(round(val * base))
         if (len(sizesScales) < 6 and (val < 0.02 or size < lowLim or size > base)) or \
                (len(sizesScales) == 6 and (size < 0 or size > base)):
            abortSet('The size specified by the "Trimvol.' + key + \
                     '" directive is out of the allowed range')
            return 1
      
      sizesScales.append(size)

   # Run findsection on the BP file if any
   if useFindSec:
      zLimits = 6 * [0]
      if runFindSection(recNames[-1], 4, 32, topBots = zLimits, comRoot = fsName):
         return 1
      if zLimits[5] - zLimits[4] <= 0:
         abortSet('Findsection on ' + recNames[-1] + \
                     ' did not find limits for the section')
         return 1

   # Set up size options in command
   trimcom = '$trimvol -f'
   if reorient == 1:
      trimcom += ' -yz'
   elif reorient == 2:
      trimcom += ' -rx'

   optNames = (' -nx ', ' -ny ', ' -nz ')
   edfNames = ('X', 'Y', 'Z')
   for ind in range(3):
      if sizesScales[ind]:
         trimcom += optNames[ind] + str(sizesScales[ind])
         halfTrim = (baseVals[ind] - sizesScales[ind]) // 2
         edfcom += edfDelAndAdd('Trimvol.' + edfNames[ind] + 'Min', halfTrim + 1)
         edfcom += edfDelAndAdd('Trimvol.' + edfNames[ind] + 'Max',
                                baseVals[ind] - halfTrim)

   if useFindSec:
      tzmin = max(1, zLimits[4] - sizesScales[6])
      tzmax = min(nzrec, zLimits[5] + sizesScales[6])
      trimcom += fmtstr(' -z {},{}', tzmin, tzmax)
      edfcom += edfDelAndAdd('Trimvol.ZMin', tzmin)
      edfcom += edfDelAndAdd('Trimvol.ZMax', tzmax)
      
   # If any scaling is specified, add options.  Set default for Z, leave others at
   # trimvol defaults
   scaleMeanSD = lookupDirective(trimPrefix, 'scaleToMeanSD', 0, STRING_VALUE)
   if scaleMeanSD:
      trimcom += fmtstr(' -meansd "{}"', scaleMeanSD)
   doScaling = sizesScales[3] or sizesScales[4] or sizesScales[5] or scaleMeanSD
   edfcom += edfDelAndAdd('Trimvol.ConvertToBytes', boolStringForEdf(doScaling))
   if doScaling:
      if not sizesScales[5]:
         sizesScales[5] = min(nzrec, max(4, nzrec // 3))
      for ind in range(3, 6):
         edfNames = ('ScaleXM', 'ScaleYM', 'ScaleSectionM')
         if sizesScales[ind]:
            start = 1 + (baseVals[ind] - sizesScales[ind]) // 2
            end = start + sizesScales[ind] - 1
            trimcom += fmtstr(' -s{} {},{}', chr(ord('x') + ind - 3), start, end)
            edfcom += edfDelAndAdd('Trimvol.' + edfNames[ind - 3] + 'in', start)
            edfcom += edfDelAndAdd('Trimvol.' + edfNames[ind - 3] + 'ax', end)

   # Do one or two volumes
   for ind in range(len(recNames)):
      mess = 'Running trimvol on ' + recNames[ind] + ' to create ' + trimNames[ind]
      if writeTextFileReportErr('trimvol' + axisCom,
                                [trimcom +
                                 fmtstr(' "{}" "{}"', recNames[ind], trimNames[ind])]):
         return 1
      if runOneProcess('trimvol' + axisCom, message = mess):
         return 1

   if modifyEdfLines(edfcom):
      return 1

   return replaceOrRunAfterStep(trimStepNum, True)


# Run anisotropic diffusion in chunks
def runNAD():
   recFile = datasetFilename('.rec')
   nadFile = datasetFilename('.nad', root = os.path.splitext(recFile)[0])
   comFile = 'autoNAD' + comExt

   # Get the directives and make sure both are there, and make sure file exists
   iterations = lookupDirective(runtimePrefix + 'NAD', 'iterations', 0, INT_VALUE)
   Kvalue = lookupDirective(runtimePrefix + 'NAD', 'Kvalue', 0, FLOAT_VALUE)
   if testDirectiveValue(iterations, 'NAD.iterations', 'integer') or \
          testDirectiveValue(Kvalue, 'NAD.Kvalue', 'float'):
      return 1
   if not iterations and not Kvalue:
      return 0
   if (iterations and not Kvalue) or (Kvalue and not iterations):
      abortSet('Both "iterations" and "Kvalue" directives need to be present to run NAD')
      return 1
   if not os.path.exists(recFile):
      abortSet('You must post-process with Trimvol in order to run NAD')
      return 1

   # Get optional memory entry, compute the chunksize, write the dummy com file
   memory = lookupDirective(runtimePrefix + 'NAD', 'chunkMemoryMB', 0, INT_VALUE)
   if testDirectiveValue(iterations, 'NAD.chunkMemoryMB', 'integer'):
      return 1
   if not memory:
      memory = 512
   chunkSize = max(5, memory // 36)
   nadcom = fmtstr('$nad_eed_3d -n {} -k {} INPUTFILE OUTPUTFILE', iterations, Kvalue)
   if writeTextFileReportErr(comFile, [nadcom]):
      return 1

   # Make the chunk coms
   padding = max(8, iterations)
   chunkCom = fmtstr('chunksetup -m {} -p {} -no {} {} {}', chunkSize, padding, comFile,
                     recFile, nadFile)
   try:
      runcmd(chunkCom)
   except ImodpyError:
     reportImodError('Error running chunksetup on ' + comFile)
     return 1

   return runOneProcess(comFile, False, message = 'Filtering with anisotropic diffusion')


# Run Reducefiltvol for reduction and filtering
def runReduceFiltVol():
   if not lookupDirective(runtimePrefix + 'Postprocess', 'doReduceFilt', 0, BOOL_VALUE):
      return 0
   binning = lookupDirective(comPrefix + 'reducefiltvol', 'reducefiltvol.ReductionFactor',
                             0, FLOAT_VALUE)
   if testDirectiveValue(binning, 'reducefiltvol.ReductionFactor', 'float'):
      return 1
   if binning == None:
      binning = 1
   comFile = 'reducefiltvol' + comExt
   recFile = datasetFilename('.rec')
   if not os.path.exists(recFile):
      abortSet('You must post-process with Trimvol in order to run Reducefiltvol')
      return 1

   comlines = ['RootNameOfDataFiles ' + dataName,
               'InputFile ' + recFile,
               'BinningOfImages ' + str(binning)]
   comlines.append('OneParameterChange ' + comPrefix + \
                   'reducefiltvol.reducefiltvol.SetupChunksIfMemoryError=1')

   comlines += laterComDirectives(0)
      
   if makeAndRunOneCom(comlines, comFile, 'Reducing and/or filtering final tomogram'):
      return 1

   rfvLines = readTextFileReportErr('reducefiltvol.log')
   if not rfvLines:
      return 1
   
   itRan = True
   comName = 'rfvfilter' + comExt
   for line in rfvLines:
      if '[MTF1]' in line:
         itRan = False
      if 'processchunks' in line:
         lsplit = line.split()
         comName = lsplit[-1]
         
   if itRan:
      return 0

   return runOneProcess(comName, False,
                        message = 'Filtering volume in chunks to overcome memory limit')


# Run tomocleanup on the data set
def runCleanup():
   keep = lookupDirective(runtimePrefix + 'Cleanup', 'doCleanup', 0, BOOL_VALUE)
   if not keep:
      return 0

   cmdstr = 'tomocleanup'
   if lookupDirective(runtimePrefix + 'Cleanup', 'keepAligned', 0, BOOL_VALUE):
      cmdstr += ' -aligned'
   if lookupDirective(runtimePrefix + 'Cleanup', 'keepUntrimmed', 0, BOOL_VALUE):
      cmdstr += ' -untrimmed'
   if lookupDirective(runtimePrefix + 'Cleanup', 'keepAxis', 0, BOOL_VALUE):
      cmdstr += ' -axis'
   if lookupDirective(runtimePrefix + 'Cleanup', 'keepSIRT', 0, BOOL_VALUE):
      cmdstr += ' -sirt'

   try:
      prnLog('Running ' + cmdstr + ' on data set')
      cleanLines = runcmd(cmdstr + ' .')
   except ImodpyError:
      reportImodError('Error running tomocleanup on data set')
      return 1

   printTaggedMessages(cleanLines, [('WARNING:', 0)])
   return 0
   

# Do all the operations on an axis
def runOneAxis():
   global xtiltNeeded, fidThickness, fidIncShift, reconThickness, transposeForAli
   global numSurfaces, axisRotation, centerOnGold, posSampleType, noXaxisTilt, alignLines
   xtiltNeeded, fidThickness, fidIncShift, reconThickness = 0., 0., 0., 0.
   clipLines = []
   
   if getAxisInitialParameters():
      return 1
   
   # Xray removal
   if lookupDirective(runtimePrefix + 'Preprocessing', 'removeXrays', 0, BOOL_VALUE) and \
          needStep(1):

      xrayMess = 'Removing X-rays with Ccderaser'
      sedcom = []

      # Copy the model file to the standard local name and switch to that unless the
      # file already exists.  But do it with imodtrans to make sure pixel size fits the
      # current stack
      model = lookupDirective(comPrefix + 'eraser', 'ccderaser.ModelFile', 0,STRING_VALUE)
      if model:
         eraseMod = setName + axisLet + '.erase'
         if os.path.basename(model) != eraseMod and not os.path.exists(eraseMod):
            try:
               runcmd(fmtstr('imodtrans -I "{}" "{}" "{}"', 
                             setName + axisLet + stackExtension, model, eraseMod))
            except ImodpyError:
               reportImodError('Error transforming manual erasing model ' + model + \
                                  ' to ' + eraseMod)
               return 1
            sedcom = [sedModify('ModelFile', eraseMod)]

      if sedcom:
         if modifyWriteAndRunCom('eraser' + axisCom, sedcom, message = xrayMess):
            return 1
      else :
         if runOneProcess('eraser' + axisCom, message = xrayMess):
            return 1

      # Get clip stats output for fixed stack
      try:
         (tmpX, tmpY, rawZ) = getmrcsize(dataName + stackExtension)
      except ImodpyError:
         reportImodError('Could not get size of raw stack')
         return 1
      length = min(30, max(15, rawZ // 4))
      montOpt = ''
      if ifMontage:
         montOpt = '-P ' + dataName + '.pl -O -10,-10 '
      command = fmtstr('$clip stats {} -10,-10 -n 2.5 -l {} {}', montOpt, length,
                       dataName + '_fixed' + stackExtension)

      # Run clip and extract the summary lines
      clipLines = runClipStats(command, '_fixed', 'fixed')
      if not clipLines:
         return 1
      prnLog(' ')
      for line in clipLines:
         if 'all' in line or 'extreme' in line:
            prnLog(line)
      prnLog(' ')

      # Use fixed stack
      if useFileAsReplacement(dataName + '_fixed' + stackExtension, 
                              dataName + stackExtension, True, False):
         return 1

      # See if there is already an archive file to rename and find next number for that
      compName = dataName + '_xray' + stackExtension + '.gz'
      if os.path.exists(compName):
            compList1 = glob.glob(compName + '.[0-9]')
            compList2 = glob.glob(compName + '.[0-9][0-9]')
            compList1.sort()
            compList2.sort()
            compList = compList1 + compList2
            nextNum = 1
            if len(compList):
               nextNum = int(os.path.splitext(compList[-1])[1][1:]) + 1
            try:
               os.rename(compName, compName + '.' + str(nextNum))
            except OSError:
               abortSet('Error renaming ' + compName + ' before archiving again')
               return 1

      # Run archiveorig
      if lookupDirective(runtimePrefix + 'Preprocessing', 'archiveOriginal', 0, 
                         BOOL_VALUE):
         if writeTextFileReportErr('archiveorig' + axisCom,
                                   ['$archiveorig -d ' + dataName + stackExtension]):
            return 1
         err = runOneProcess('archiveorig' + axisCom, message =
                             'Archiving original stack as compressed difference file')
         cleanupFiles(['archiveorig' + axisCom])
         if err:
            return 1

   # Possibly a stopgap: look for low-count views at end of series and exclude them
   if not ifMontage and needStep(1) and analyzeSDsAdjustExcludes(clipLines):
      return 1
         
   # coarse alignment
   if needStep(2):
      if runOneProcess('xcorr' + axisCom, message = 'Finding coarse ' + \
                          'alignment by cross-correlation with Tiltxcorr'):
         return 1
      if ifMontage:
         edfcom = edfDelAndAdd('xcorr.blendmont.' + axisEdfLet + '.WasRun', 'true')
         if modifyEdfLines(edfcom):
            return 1

   # Get the tiltalign file regardless, since it will be operated on
   alignLines = readTextFileReportErr('align' + axisCom)
   if not alignLines:
      return 1
   rotarr = optionValue(alignLines, 'RotationAngle', FLOAT_VALUE)
   if not rotarr:
      abortSet('Cannot find RotationAngle in align' + axisCom)
      return 1
   axisRotation = rotarr[0]
   while axisRotation > 180.:
      axisRotation -= 360.
   while axisRotation < -180.:
      axisRotation += 360.
   transposeForAli = math.fabs(axisRotation) > 45. and math.fabs(axisRotation) < 135.

   # Get the number of surfaces from directives in case it changed, fall back to align.com
   numSurfaces = lookupDirective(comPrefix + 'align', 'tiltalign.SurfacesToAnalyze', 0,
                                 INT_VALUE)
   if testDirectiveValue(numSurfaces, 'tiltalign.SurfacesToAnalyze', 'integer'):
      return 1
   
   if numSurfaces == None:
      numSurfaces = optionValue(alignLines, 'SurfacesToAnalyze', INT_VALUE, numVal = 1)
      if numSurfaces == None:
         abortSet('Failed to find SurfacesToAnalyze in align' + axisCom)
   if numSurfaces <= 0:
      abortSet('SurfacesToAnalyze is 0 or negative in directives or in align' + axisCom)
      return 1

   if fiducialless or patchTrack:
      numSurfaces = 1

   # Get whether X axis tilt should be kept at 0
   noXaxisTilt = lookupDirective(runtimePrefix + 'Reconstruction', 'noXAxisTilt', 0, 
                                 BOOL_VALUE)

   # Get positioning and set up for centerOnGold if no positioning and 2 surfaces
   posSampleType = lookupDirective(runtimePrefix + 'Positioning', 'sampleType', 0,
                                   INT_VALUE)
   if testDirectiveValue(posSampleType, 'Positioning.sampleType', 'integer'):
      return 1

   if posSampleType == None:
      posSampleType = 0
   centerOnGold = False
   if posSampleType <= 0 and numSurfaces > 1 and \
          lookupDirective(runtimePrefix + 'Positioning', 'centerOnGold', 0, BOOL_VALUE):
      centerOnGold = True
      
   # For fiducialless, make up the final xf
   if fiducialless:
      if fidlessFileOperations():
         return 1
      reportReachedStep(5)

   # Otherwise do many things
   else:

      # Make prealigned stack
      (comRoot, process) = comAndProcessForAlignedStack('pre')
      if needStep(3):
         err = runOneProcess(comRoot + axisCom, message = 'Making coarse aligned stack')
         if ifMontage:
            edfcom = edfDelAndAdd('InvalidEdgeFunctions' + axisUpperLet,
                                  boolStringForEdf(err))
            if modifyEdfLines(edfcom):
               return 1
         if err:
            return 1

      # Patch tracking
      if needStep(4) and patchTrack:
         if runPatchTracking():
            return 1

      # Otherwise various ways of getting a fiducial model
      else:
         if makeSeedAndTrack():
            return 1

      # Iterate patch tracking with angle offset if desired, afer running align
      iteratePatchTrack = lookupDirective(patchTrackText, 'adjustTiltAngles', 0,
                                          BOOL_VALUE) and needStep(4) and patchTrack
      if iteratePatchTrack:
         if (runTiltalign()):
            return 1
         if OKtoAdjustPatchTrack():
            sedcom = sedDelAndAdd('AngleOffset', totalDelTilt, 'SizeOfPatchesXandY')
            if modifyWriteAndRunCom('xcorr_pt' + axisCom, sedcom, message = fmtstr(
                  'Running patch tracking again with angle offset of {:.1f}', 
                  totalDelTilt)):
               return 1

      # Tilt alignment results needed if aligned stack or reconstruction being made
      reportReachedStep(5)
      if endingStep >= 6 and startingStep <= 14 and runTiltalign():
         return 1

   # Positioning
   reportReachedStep(6)
   if needStep(7) and positionTomogram():
      return 1
   
   # Make aligned stack (tests need, sets correctCTF global, remakes if already corrected)
   reportReachedStep(7)
   if makeAlignedStack(-1):
      return 1
         
   # Do CTF plotter step (tests need)
   if CTFPlotAlignedStack():
      return 1

   # Set up the tilt com file now that output size is known and montage frame data set
   if endingStep >= detect3DstepNum and startingStep <= 14 and modifyTiltComFile(0):
      return 1

   # Do 3D gold detection (tests need, sets eraseGold global)
   if detectGoldIn3D():
      return 1
   
   # Optionally CTF correct the stack
   reportReachedStep(detect3DstepNum)
   if needStep(ctfCorrStepNum) and correctCTF and CTFCorrectAlignedStack():
      return 1

   # Optionally erase gold
   if needStep(eraseStepNum) and eraseGold and eraseGoldInAlignedStack():
      return 1

   # Optionally filter
   if filterAlignedStack():
      return 1

   # Make the reconstruction
   reportReachedStep(13)
   if needStep(14) and generateTomogram():
      return 1

   # Run trimvol on full data set or at axis level if directive for it
   reportReachedStep(14)
   if ((not dualAxis and needStep(trimStepNum)) or (dualAxis and needStep(14.5))) and \
          trimVolume():
      return 1

   # Run NAD on full data set only
   if not dualAxis and needStep(nadStepNum) and runNAD():
      return 1
   
   # Run Reducefiltvol on full data set only
   if not dualAxis and needStep(redFiltStepNum) and runReduceFiltVol():
      return 1
   
   # Run Cleanup on full data set only
   if not dualAxis and needStep(cleanStepNum) and runCleanup():
      return 1

   return 0


# Or run the steps for combine, return on error or when done
def runCombine():
   if needStep(15) and setupCombine():
      return 1

   if needStep(16) and initialCombineMatch():
      return 1

   if alignAndCombineAxes():
      return 1

   if needStep(trimStepNum) and trimVolume():
      return 1

   if needStep(nadStepNum) and runNAD():
      return 1

   if needStep(redFiltStepNum) and runReduceFiltVol():
      return 1

   if needStep(cleanStepNum):
      return runCleanup()

   return 0
   

#### MAIN PROGRAM  ####
#
# load System Libraries
import os, sys, shutil, math, csv, re, time, datetime, platform, socket

#
# Setup runtime environment
if os.getenv('IMOD_DIR') != None:
   IMOD_DIR = os.environ['IMOD_DIR']
   if sys.platform == 'cygwin' and sys.version_info[0] > 2:
      IMOD_DIR = IMOD_DIR.replace('\\', '/')
      if IMOD_DIR[1] == ':' and IMOD_DIR[2] == '/':
         IMOD_DIR = '/cygdrive/' + IMOD_DIR[0].lower() + IMOD_DIR[2:]
   sys.path.insert(0, os.path.join(IMOD_DIR, 'pylib'))
   from imodpy import *
   addIMODbinIgnoreSIGHUP()
else:
   sys.stdout.write(prefix + " IMOD_DIR is not defined!\n")
   sys.exit(1)

#
# load IMOD Libraries
from pip import *
from pysed import *
from comchanger import *
from prochunks import *

# Fallbacks from ../manpages/autodoc2man 3 1 batchruntomo
options = ["directive:DirectiveFile:FNM:", "root:RootName:FNM:",
           "current:CurrentLocation:FNM:", "deliver:DeliverToDirectory:FNM:",
           "make:MakeSubDirectory:B:", "one:ProcessOneAxis:I:",
           "cpus:CPUMachineList:CH:", "single:SingleOnFirstCPU:B:",
           "gpus:GPUMachineList:CH:", "parallel:ParallelBatchRootName:CH:",
           "maxGPUs:MaxGPUsInParallelBatch:I:", "queue:QueueCommand:CH:",
           "jobs:MaxJobsOnQueue:I:", "jcores:CoresPerClusterJob:I:",
           "jgpus:GPUsPerClusterJob:I:", "gqueue:GPUQueueCommand:CH:",
           "gjobs:MaxGPUJobsOnQueue:I:", "frompath:TranslatePathsFrom:CHM:",
           "topath:TranslatePathsTo:CHM:", "remote:RemoteDirectory:FN:",
           "nice:NiceValue:I:", "limit:LimitLocalThreads:I:", "check:CheckFile:FN:",
           "email:EmailAddress:CH:", "SMTP:SMTPserver:CH:",
           "validation:ValidationType:I:", "end:EndingStep:F:", "start:StartingStep:F:",
           "first:StartForFirstSetOnly:B:", "use:UseExistingAlignment:B:",
           "exit:ExitOnError:B:", "style:NamingStyle:I:", "pcm:MakeComExtensionPcm:I:",
           "ext:StackExtension:CH:", "axis:AxisOfExtension:I:", "etomo:EtomoDebug:I:",
           "bypass:BypassEtomo:B:", ":PID:B:", "help:usage:B:"]

copyPrefix = 'setupset.copyarg.'
setupPrefix = 'setupset.'
runtimePrefix = 'runtime.'
comPrefix = 'comparam.'
scopeTmplText = setupPrefix + 'scopeTemplate'
userTmplText = setupPrefix + 'userTemplate'
sysTmplText = setupPrefix + 'systemTemplate'
dataDirText = setupPrefix + 'datasetDirectory'
rootNameText = copyPrefix + 'name'
cpTomoExtText = copyPrefix + 'stackext'
scanHeadText = setupPrefix + 'scanHeader'
patchTrackText = runtimePrefix + 'PatchTracking'
autoSeedText = runtimePrefix + 'SeedFinding'
combinePrefix = runtimePrefix + 'Combine'
localBatchCopy = 'batchDirective.adoc'
myPID = os.getpid()
proChunkCheckFile = 'processchunks.cmds'
logSuffixTag = '    [:LOG]'
montFrameData = 10 * [0]
suppressAbort = False
altExtension = 'mrc'
userTemplateDir = ''
summaryMessage = ''
ctfPlotStepNum = 9
detect3DstepNum = 10
ctfCorrStepNum = 11
eraseStepNum = 12
trimStepNum = 20
nadStepNum = 21
redFiltStepNum = 21.5
cleanStepNum = 22
finalStepNum = 22

# Defaults
dfltDarkExcludeFraction = 0.33
dfltDarkExcludeRatio = 0.17
dfltMatchAtoBratio = 0.9
dfltCombineNumScales = 4
dfltCombineBoxSize = 32
dfltSampleThickness = 500
positionBinningTarget = 512
dfltPosThicknesses = (250, 400, 500, 600)
sizesForPosThicknesses = (0, 512, 1024, 2048)
dfltExtraWarpTargets = '0.4,0.45'
dfltFinalPatchSize = 'E'
solvematchMinFids = 8
solvematchMinEachSide = 3
solvematchMinSideRatio = 0.2
useFallbackRatio = 0.4
useVolMatch = -1
fb3dOptimalBinnedSize = 5.
fb3dMinBinnedSize = 4.
cryoPosExtraThick = 25
dfltPTmaxTiltAdjust = 20.
dfltPTmaxAdjustedAngle = 78.
expandFactor = 1.

logFile = None
finalRetVal = 0
renamingOnly = False
queueCommand = None
finishSetAndQuit = False
latestMessages = []
edfChanged = False
needGPUrelease = False
edfLines = []
makeSymLinks = None
originalStackExt = ''
standardTypeExts = standardTypeExtensions()
possibleStackExts = allowedRawStackExtensions()
direcFileErrorNames = ['batch defaults', 'scope template', 'system template',
                       'user template', 'batch directive']

# List of coms where the runner is handling messages, and the standard tags for messages
handlingMessages = ['autofidseed', 'transferfid', 'track', 'align', 'solvematch',
                    'autopatchfit', 'dualvolmatch']
standardTags = [('ERROR:', 0), ('WARNING:', 0)]

# This handles the problem with # being lost as a comment in standard input
PipForbidComments('CPUMachineList', 'cpus', False)

(opts, nonopts) = PipReadOrParseOptions(sys.argv, options, progname, 1, 1, 0)
os.environ['PIP_PRINT_ENTRIES'] = '0'
os.environ['SKIP_PHYSICAL_DPI'] = '1'

doPID = PipGetBoolean('PID', 0)
printPID(doPID)
bypassEtomo = PipGetBoolean('BypassEtomo', 0)

# Intercept stack extension renaming call and process completely
fromExtension = PipGetString('StackExtension', '')
if fromExtension:
   if fromExtension not in possibleStackExts:
      exitError(fromExtension + ' is not an allowed stack extension')
   renamingOnly = True
   (nameStyle, typeExtension) = getNamingStyle(None)
   if nameStyle < 0:
      exitError('The naming style must be entered for a stack renaming run')
   if PipNumberOfEntries('RootName') != 1:
      exitError('One root name must be entered for a stack renaming run')
   setName = PipGetString('RootName', '')
   dualAxis = PipGetInteger('AxisOfExtension', 0)
   if dualAxis < 0 or dualAxis > 2:
      exitError('Axis entry must be 0, 1, or 2')
   stack = setName + '.'
   if dualAxis:
      stack = setName + 'a.'

   stackExtension = ''
   extIfAlreadySet = ''
   if dualAxis == 2:
      checkRenameStack(setName + 'b.')
   checkRenameStack(stack)
   if dualAxis == 1:
      suppressAbort = True
      checkRenameStack(setName + 'b.')
   sys.exit(0)


# Get current directory
startingDir = os.getcwd()

# If doing parallel batch, get possible translation paths
parallelRoot = PipGetString('ParallelBatchRootName', '')
fromParallelPaths = []
toParallelPaths = []
numTranslations = 0
if parallelRoot:
   proChunkCheckFile = 'processchunks-' + str(myPID) + '.cmds'
   numTranslations = PipNumberOfEntries('TranslatePathsFrom')
   if numTranslations != PipNumberOfEntries('TranslatePathsTo'):
      exitError('There are not the same number of entries for TranslatePathsFrom ' +\
                'and TranslatePathsTo')
   for ind in range(numTranslations):
      fromParallelPaths.append(PipGetString('TranslatePathsFrom', ''))
      toParallelPaths.append(PipGetString('TranslatePathsTo', ''))

# Get the directive file names
numByArg = PipNumberOfEntries('DirectiveFile')
numSets = nonopts + numByArg
dirFiles = []
direcTranslateInds = []
if not numSets:
   exitError('You must enter at least one directive file')
for ind in range(numSets):
   if ind < numByArg:
      dfile = translateParallelPath(PipGetString('DirectiveFile', ''))
   else:
      dfile = translateParallelPath(PipGetNonOptionArg(ind))
   if not os.path.exists(dfile):
      exitError('Directive file ' + dfile + ' does not exist')
   dirFiles.append(dfile)
   direcTranslateInds.append(lastTranslateInd)

# Get current directory and delivery directory
currentDirs = []
curDirTranslateInds = []
numCurrent = PipNumberOfEntries('CurrentLocation')
for ind in range(numCurrent):
   currentDirs.append(translateParallelPath(PipGetString('CurrentLocation', '')))
   curDirTranslateInds.append(lastTranslateInd)
   
numDeliver = PipNumberOfEntries('DeliverToDirectory')
deliverDirs = []
deliverTranslateInds = []
for ind in range(numDeliver):
      deliverDirs.append(translateParallelPath(PipGetString('DeliverToDirectory', '')))
      deliverTranslateInds.append(lastTranslateInd)
      
if numDeliver == 1 and not os.path.isdir(deliverDirs[0]):
   if os.path.exists(deliverDirs[0]):
      exitError('The -deliver entry ' + deliverDirs[0] + ' must be a directory; there ' +\
                   'is a file with this name instead')
   exitError('The location in which to create dataset directories (' + deliverDirs[0] + \
                ') must already exist; it does not')
   
makeSubDir = PipGetBoolean('MakeSubDirectory', 0)
if numDeliver and makeSubDir:
   exitError('You cannot enter both -deliver and -make options')
doDelivery = numDeliver or makeSubDir
if numDeliver and numDeliver != 1 and numDeliver != numSets:
   exitError('If the -deliver option is used, it must be entered either once, ' +\
                'or once per data set')
if doDelivery and numCurrent != 1 and numCurrent != numSets:
   exitError('The -current option must be entered either once, or once per data set, ' +\
                'if -deliver or -make is entered')

# Get the root names if any
numRootOpts = PipNumberOfEntries('RootName')
rootByOption = []
for ind in range(numRootOpts):
   rootByOption.append(PipGetString('RootName', ''))

# Check validity of these entries and boost the directive file list
if numRootOpts and numSets > 1 and numSets != numRootOpts:
   exitError('If you enter root names, you must enter either one directive file or' + \
             ' one per root name')
if numRootOpts:
   for ind in range(numSets, numRootOpts):
      dirFiles.append(dirFiles[0])
      direcTranslateInds.append(direcTranslateInds[0])
   numSets = numRootOpts

if not doDelivery and (numRootOpts or numCurrent) and numCurrent != numSets:
   exitError('The -current option must be entered for each data set')

# Get other options

# Get the CPU list from environment or option
if os.getenv('MULTI_PROC_CPU_LIST') != None:
   cpuList = os.environ['MULTI_PROC_CPU_LIST']
   if cpuList == 'None':
      cpuList = ''
else:
   cpuList = PipGetString('CPUMachineList', '')

# Get cluster queue options and alternative of cores when running on a node
coresPerClusterJob = 0
(queueCommand, runningOnQueue, maxQueueJobs) = \
   getQueueOptions('QueueCommand', 'MULTI_PROC_QUEUE_COMMAND', 'MaxJobsOnQueue',
                   'MULTI_PROC_MAX_QUEUE_JOBS', '')
   
coresPerClusterJob = getClusterJobOption('CoresPerClusterJob', 'MULTI_PROC_JOB_CORES')
if coresPerClusterJob:
   if queueCommand:
      exitError('You cannot enter CoresPerClusterJob along with QueueCommand')
   if cpuList:
      exitError('You cannot enter CoresPerClusterJob along with a CPU machine list')
   cpuList = str(coresPerClusterJob)

# For CPU's, validate the number if it is just a number, then parse the list and add up
# the values after # if any, set flag for parallel processing = # of cores
localName = (platform.node().split('.'))[0]
parallelCPU = 0
mostCPUs = 0
topCPUmachine = '1'
topCpuLimit = 1
localCpuLimit = 1
firstCpuLimit = 1
if '#' in cpuList:
   warning('The # sign should not be used in the -cpus entry; use : instead', False)

if cpuList != '':
   localCpuLimit = 0
   firstCpuLimit = 0
   try:
      parallelCPU = int(cpuList)
      localCpuLimit = parallelCPU
      firstCpuLimit = parallelCPU
      topCpuLimit = parallelCPU
      if parallelCPU < 1 or parallelCPU > 128:
         exitError('A number of cores entered with -cpus must be between 1 and 128')
   except ValueError:
      for machine in cpuList.split(','):
         msplit = machine.replace('#', ':').split(':')
         if len(msplit) > 2:
            exitError('A machine name cannot be followed by two : or # signs')
         if len(msplit) < 2:
            numCPU = 1
         else:
            try:
               numCPU = int(msplit[1])
               if numCPU < 1:
                  exitError('The value after : or # is less than 1 in ' + machine)
            except ValueError:
               exitError('Failed to convert value after : or # to integer in ' + machine)
         parallelCPU += numCPU
         if msplit[0] == 'localhost' or (msplit[0].split('.'))[0] == localName:
            localCpuLimit += numCPU
         if numCPU > mostCPUs:
            mostCPUs = numCPU
            topCPUmachine = msplit[0]
            topCpuLimit = numCPU
         if not firstCpuLimit:
            firstCpuLimit = numCPU

# After all that, replace the #'s in the list that will get used, and accept an option
# to set or revise the limit on local cores, and if it is high enough, then revise the
# machine with most CPUs
cpuList = cpuList.replace('#', ':')

# Get thread limit from variable if running in parallel, or from option
if queueCommand:
   threadLimit = localLimit = firstCpuLimit = topCpuLimit = 1
elif coresPerClusterJob:
   threadLimit = localLimit = firstCpuLimit = topCpuLimit = coresPerClusterJob

elif os.getenv('MULTI_PROC_THREAD_LIMIT') != None:
   try:
      temp = os.environ['MULTI_PROC_THREAD_LIMIT']
      localLimit = int(temp)
   except ValueError:
      exitError('Converting MULTI_PROC_THREAD_LIMIT value of ' + temp + ' to an integer')
else:
   localLimit = PipGetInteger('LimitLocalThreads', 0)
if localLimit < 0 or localLimit > 128:
   exitError('The limit on number of local threads must be between 1 and 128')
if localLimit:
   localCpuLimit = localLimit
   if localLimit > topCpuLimit:
      topCPUmachine = 'localhost'
      topCpuLimit = localLimit

# For GPU's, set useGPU if entered, insist a number is 1, and set flag for parallel GPU
# if there is an actual machine list to number of GPUs, adding up entries 
# with : separators.  Here keep parallelGPU = 1 if
# it is running without splitting on a single remote machine
gpuList = PipGetString('GPUMachineList', '')
if os.getenv('MULTI_PROC_GPU_POOL'):
   gpuList = os.environ['MULTI_PROC_GPU_POOL']
   if gpuList == 'None':
      gpuList = ''
useGPU = gpuList != ''
parallelGPU = 0

# Determine if a GPU is being done with a separate queue
(gpuQueueCommand, nonsense, maxGpuQueueJobs) = \
   getQueueOptions('GPUQueueCommand', 'MULTI_PROC_GPU_QUEUE', 'MaxGPUJobsOnQueue', 
                   'MULTI_PROC_MAX_GPU_JOBS', 'GPU ')
if gpuQueueCommand:
   if useGPU:
      exitError('You cannot enter both a GPU queue command and a GPU machine list')
   useGPU = 1
   if maxGpuQueueJobs > 1:
      parallelGPU = maxGpuQueueJobs
   gpuList = '1'

# And find out if there are multiple GPUs for the current multicore node
gpusPerClusterJob = getClusterJobOption('GPUsPerClusterJob', 'MULTI_PROC_JOB_GPUS')
if gpusPerClusterJob:
   if queueCommand or gpuQueueCommand:
      exitError('You cannot enter GPUs per job with a queue command')
   if not coresPerClusterJob:
      exitError('You cannot enter GPUs per job without cores per job')

   useGPU = gpusPerClusterJob
   gpuList = '1'
   if useGPU > 1:
      gpuList = 'localhost'
      parallelGPU = useGPU
      for ind in range(useGPU):
         gpuList += ':' + str(ind + 1)

# Now process non-queue case for a gpu machine list
maxParallelGPUs = 0
if useGPU and not (gpuQueueCommand or gpusPerClusterJob):
   if queueCommand:
      exitError('You cannot enter a GPU machine list with a queue command')
   try:
      useGPU = int(gpuList)
      if useGPU != 1:
         exitError('The entry for a GPU list must be 1 to use just the local GPU')

   except ValueError:
      if coresPerClusterJob:
         exitError('With cores per node specified, an entry for the GPU machine ' + \
                   'list must be a single positive number')
      
      for machine in gpuList.split(','):
         msplit = machine.split(':')
         parallelGPU += max(1, len(msplit) - 1)

if parallelRoot and useGPU:
   maxParallelGPUs = PipGetInteger('MaxGPUsInParallelBatch', 4)
   if maxParallelGPUs < 1:
      exitError('Maximum GPUs for parallel batch runs must be positive')
   if gpuQueueCommand:
      parallelGPU = min(parallelGPU, maxParallelGPUs)
      maxGpuQueueJobs = min(maxGpuQueueJobs, maxParallelGPUs)
      maxParallelGPUs =0
   elif not (queueCommand or coresPerClusterJob):
      hostname = socket.gethostname().split('.')[0]

niceness = PipGetInteger('NiceValue', 15)
remoteStartDir = PipGetString('RemoteDirectory', '')
doOneAxis = PipGetInteger('ProcessOneAxis', 0)
startingStep = PipGetFloat('StartingStep', 0.)
endingStep = PipGetFloat('EndingStep', 10000000.)
firstStart = PipGetBoolean('StartForFirstSetOnly', 0)
exitOnError = PipGetBoolean('ExitOnError', 0)
etomoDebug = PipGetInteger('EtomoDebug', 0)
useFirstCPUforSingle = PipGetBoolean('SingleOnFirstCPU',0)
if not firstStart and startingStep > 100:
   warning('Entering a starting step above 100 has no effect without ' +\
              '-first option', False)
   while startingStep > 100:
      startingStep -= 100

# Get naming style and extension information, and set the output format if necessary
(nameStyleDflt, typeExtension) = defaultNamingStyle()
(nameStyle, typeExtension) = getNamingStyle(typeExtension)
if nameStyle < 0:
   nameStyle = nameStyleDflt
comExt = comExtensionFromOption(0)

setupBname = 'laterbsetup' + comExt

setOutputFormatIfNeeded(typeExtension)

# TESTING ONLY: set test variables for running etomo
if typeExtension and not bypassEtomo:
   os.environ['TEST_NAMING_STYLE'] = str(nameStyle)
if comExt == '.pcm' and not bypassEtomo:
   os.environ['TEST_USE_PCM_FOR_COM'] = '1'

skipTiltalign = 0
if not (startingStep <= 6.005 or (startingStep > 99.995 and startingStep <= 106.005)):
   skipTiltalign = PipGetBoolean('UseExistingAlignment', 1)

validation = PipGetInteger('ValidationType', 0)
remoteDataDir = ''
emailAddress = PipGetString('EmailAddress', '')
SMTPserver = PipGetString('SMTPserver', 'localhost')
if emailAddress:
   if pyVersion < 250:
      emailAddress = ''
      warning('No emails will be sent; the Python version must be 2.5 or ' + \
                 'higher to send emails', False)
   else:
      import smtplib
      from email.mime.text import MIMEText

# Set up check file
checkFile = PipGetString('CheckFile', '')
if checkFile:
   checkFile = imodAbsPath(checkFile)
else:
   checkFile = imodAbsPath(os.path.join(startingDir, progname + '.' +
                                    str(os.getpid()) + '.input'))
if validation <= 0:
   prnstr('To quit all processing, place a Q in the file: ' +
          translateParallelPath(checkFile))
   if parallelRoot:
      prnstr('Running on host ' + platform.node())
if os.path.exists(checkFile):
   cleanupFiles([checkFile])

baseComDict = {}
validComDict = {}
validRunDict = {}
validOtherDict = {}
fileType = 'directive'
defaultsFile = os.path.join(IMOD_DIR, 'com/batchDefaults.adoc')

if validation >= 0:
   validateFile = os.path.join(os.path.join(IMOD_DIR, 'com'), \
       'directives.csv')
   if not os.path.exists(validateFile):
      exitError('Cannot find file for validating directives, ' + validateFile)

   processValidationFile()
   if validation > 1:
      fileType = 'template'

# Start looping on the data sets
for dfileInd in range(numSets):

   # These need to be initialized for abortSet to work right
   dualAxis = False
   axisLet = None
   setName = '# ' + str(dfileInd + 1)

   closeLogFileWriteEdf()
   dfile = dirFiles[dfileInd]
   axisNum = 0
   dsetDirTranslateInd = -1
   
   # Check for an F in the check file before going on
   # No need to remove processchunks check file, it takes care of it itself
   checkForQuit()
   if finishSetAndQuit:
      prnstr('Exiting after finishing dataset as requested   [brt6]')
      sys.exit(0)

   os.chdir(startingDir)
   allDirectives = [{}, {}, {}, {}, {}]
   absDirectiveFile = imodAbsPath(dfile)
   (err, directLines, cptExtDirective, rewriteBatch, 
    dirForRewrite) = readDirectiveOrTemplate(dfile, batDictInd)
   if err:
      continue
   
   err = readDirectiveOrTemplate(defaultsFile, 0)
   if err[0]:
      continue

   dfilePrn = reverseTranslatePath(dfile, direcTranslateInds[dfileInd])
   
   prnstr(fmtstr('Beginning to process {} file # {} : {}   [brt7]', fileType,
                 dfileInd + 1, dfilePrn))
   startTime = time.time()

   # If validating a single template file, do not process it, just run the checks
   if validation > 1:
      if not checkAllDirectives('template'):
         prnstr('Directives all seem OK in that file')
      continue
   
   # Get essential setup items from directives
   if scanSetupDirectives():
      continue

   if validation >= 0:
      if checkAllDirectives():
         continue
      if validation > 0:
         prnstr('Directives all seem OK in that file')
         continue

   checkDefaultsInBatchFile()

   # If no fiducial size, add directive for 0 now
   if fidSizeNm == None:
      directLines.append(copyPrefix + 'gold = 0.')
      rewriteBatch = True;

   # Etomo will look for 'set '
   prnstr('Starting data set ' + setName + '   at ' + \
             datetime.datetime.now().strftime('%H:%M:%S') + '   [brt1]')
   prnstr('', flush = True)
      
   stackExtension = ''
   doBfirst = False
   fromExtension = lookupDirective(setupPrefix, 'currentStackExt', 0, STRING_VALUE)
   if not fromExtension and dualAxis and doOneAxis != 1:
      fromExtension = lookupDirective(setupPrefix, 'currentBStackExt', 0, STRING_VALUE)
      if fromExtension:
         doBfirst = True

   extIfAlreadySet = ''
   if fromExtension:
      extIfAlreadySet = 'st'
      if typeExtension:
         extIfAlreadySet = fromExtension

   # Deliver stacks to dataset directory, from the source current dir or the one for this
   # data set
   if doDelivery:
      dirInd = min(dfileInd, len(currentDirs) - 1)
      deliverFromDir = currentDirs[dirInd]
      delFromTranslateInd = curDirTranslateInds[dirInd]
      dsetDirTranslateInd = delFromTranslateInd
      if deliverDirs:
         dsetDirTranslateInd = deliverTranslateInds[min(dfileInd, len(deliverDirs) - 1)]
      if doBfirst and deliverStack('b'):
         continue
      if dualAxis and doOneAxis < 2 and deliverStack('a'):
         continue
      if not doBfirst and dualAxis and doOneAxis != 1 and deliverStack('b'):
         continue
      if not dualAxis and deliverStack(''):
         continue

   # If there is a remote directory, need to shift it to the dataset dir
   # but translate the path of the starting dir since the dataset dir was also
   if remoteStartDir or numTranslations:
      dirToTranslate = datasetDir
      if not remoteStartDir:
         remoteStartDir = translateParallelPath(startingDir)
         dirToTranslate = translateParallelPath(datasetDir)
      (remoteDataDir, errMess) = transferRemoteDirectory \
                                 (remoteStartDir, translateParallelPath(startingDir), 
                                  dirToTranslate)
      
      if not remoteDataDir:
         abortSet(errMess)
         continue

   try:
      os.chdir(datasetDir)
   except OSError:
      abortSet('Error changing to directory ' + datasetDir)
      continue

   # Open log file now
   try:
      logFile = open('batchruntomo.log', 'a')

      # Etomo needs the comma
      prnLog('\nBatchruntomo started on data set ' + setName + ', ' + \
                datetime.datetime.now().ctime() + '  [brt11]')
      prnLog('In: ' + reverseTranslatePath(imodAbsPath(datasetDir), dsetDirTranslateInd)
             + '   [brt13]')

   except IOError:
      warning('Failed to open or write to log file in data set directory', False)
      try:
         logFile.close()
         logFile = None
      except:
         pass
   
   # Check existence of file(s)
   stack = setName + '.'
   numAxes = 1
   if dualAxis:
      numAxes = 2
      stack = setName + 'a.'
   if doBfirst and checkRenameStack(setName + 'b.'):
      continue
   if checkRenameStack(stack):
      continue
   if not doBfirst and dualAxis and doOneAxis != 1 and checkRenameStack(setName + 'b.'):
      continue
   stack += stackExtension[1:]

   # Add or correct a stackext entry in the batch file for etomo setup
   if cptExtDirective[0] < 0 or cptExtDirective[1] != stackExtension[1:]:
      rewriteBatch = True
      extLine = cpTomoExtText + ' = ' + stackExtension[1:]
      if cptExtDirective[0] < 0:
         directLines.append(extLine)
      else:
         directLines[cptExtDirective[0]] = extLine

   # Rewrite the batch file with any corrections in the data directory
   if rewriteBatch:
      absDirectiveFile = os.path.join(dirForRewrite, localBatchCopy)
      if writeTextFileReportErr(absDirectiveFile, directLines):
         continue

   # Get pixel size if scanHeader set
   if scanHeader and not pixelSize:
      try:
         px = getmrcpixel(stack)
         pixelSize = px / 10.
      except ImodpyError:
         reportImodError('Error getting pixel size from ' + stack)
         continue

   if startingStep <= 0 and not (dualAxis and doOneAxis > 1):
      if bypassEtomo:

         # run copytomocoms directly for testing
         # The adoc will need rotation and brotation
         copyargs = [fmtstr('pixel {}', pixelSize),
                     'name ' + setName,
                     'stackext ' + stackExtension[1:]]
         if comExt == '.pcm':
            copyargs.append('pcm 1')
         if typeExtension:
            copyargs.append('style ' + str(nameStyle))            
         if not lookupDirective(copyPrefix, 'userawtlt', 0, BOOL_VALUE):
            copyargs.append('extract 1')
         if not lookupDirective(copyPrefix, 'buserawtlt', 0, BOOL_VALUE):
            copyargs.append('bextract 1')
         if not lookupDirective(copyPrefix, 'rotation', 0, FLOAT_VALUE):
            try:
               anglesEtc = getmrc(stack, angleLineValues = True)
               if anglesEtc[0] != None:
                  copyargs.append('rotation ' + str(anglesEtc[0]))
               if dualAxis and not lookupDirective(copyPrefix, 'brotation', 0,
                                                   FLOAT_VALUE):
                  anglesEtc = getmrc(setName + 'b' + stackExtension, 
                                     angleLineValues = True)
                  if anglesEtc[0] != None:
                     copyargs.append('brotation ' + str(anglesEtc[0]))
            except ImodpyError:
               reportImodError('Error trying to get rotation angle')

         for key in allDirectives[batDictInd]:
            if copyPrefix in key:
               arg = key[len(copyPrefix):]
               val = allDirectives[batDictInd][key][0]
               if val:
                  copyargs.append(arg + ' ' + val)
         if scopeTmplText in allDirectives[batDictInd]:
            copyargs.append('change ' + allDirectives[batDictInd][scopeTmplText][0])
         if sysTmplText in allDirectives[batDictInd]:
            copyargs.append('change ' + allDirectives[batDictInd][sysTmplText][0])
         if userTmplText in allDirectives[batDictInd]:
            copyargs.append('change ' + allDirectives[batDictInd][userTmplText][0])
         copyargs.append('change ' + absDirectiveFile)
         try:
            runcmd('copytomocoms -StandardInput', copyargs, 'stdout')
         except ImodpyError:
            reportImodError('Error running copytomocoms')
            continue
      else:

         # Run etomo for setup and try to report errors
         failed = False
         errlog = ''
         try:
            comline = 'etomo --fromBRT --directive "' + absDirectiveFile + '"' + \
                      ' --namingstyle ' + str(nameStyle)
            if etomoDebug:
               comline += ' --debug ' + str(etomoDebug)
            if parallelCPU > 1:
               comline += ' --cpus "' + cpuList + '"'
            if gpuList != '':
               comline += ' --gpus "' + gpuList + '"'

            etlines = runcmd(comline)
            for l in etlines:
               prnLog(l.rstrip())
               ind = l.find('with log in')
               if ind > 0:
                  errlog = l[ind + 11:].strip()
         except ImodpyError:
            failed = True
            for l in getErrStrings():
               ind = l.find('check:')
               if ind > 0:
                  errlog = l[ind + 7:].strip()

         if errlog:
            if '/' in errlog or '\\' in errlog:
               try:
                  shutil.copy(errlog, '.')
               except Exception:
                  warning('Failed to copy ' + errlog + ' to dataset directory')

            loglines = readTextFile(errlog, None, True)
            if isinstance(loglines, str):
               warning('Error ' + loglines)
            else:
               printTaggedMessages(loglines, [('INFO:', 1), ('LOG:', 1)])
               err = 0
               for line in loglines:
                  if 'Pixel spacing was set in FEI file' in line or \
                     'Pixel spacing was set in file' in line:
                     err = fixAllSuppliedModelHeaders()
                     break

               if err:
                  prnLog('Cannot proceed with this data set')
                  continue
               
         else:
            if failed:
               prnLog('Running etomo failed, no etomo error log available')
            else:
               prnLog('Cannot access an error log from running etomo')
         if failed:
            reportImodError('etomo setup failed, cannot proceed with this data set')
            continue

         # Report on excluded views
         if excludedViewsA and lookupDirective(runtimePrefix + 'Preprocessing',
                                               'removeExcludedViews', 0, BOOL_VALUE):
            mess = 'Removed excluded views ' + excludedViewsA
            if dualAxis:
               mess += ' for axis A'
            prnLog(mess + '  [brt10]')
            
         if excludedViewsB and dualAxis:
            axisNum = 1
            if lookupDirective(runtimePrefix + 'Preprocessing', 'removeExcludedViews', 0,
                               BOOL_VALUE):
               prnLog('Removed excluded views ' + excludedViewsB + ' for axis B  [brt10]')
            axisNum = 0

   reportReachedStep(0)
   if fidSizeNm:
      fidSizePix = fidSizeNm / pixelSize
   axisFailed = False

   if originalStackExt:
      edfcom = edfDelAndAdd('OrigImageStackExt', originalStackExt)
      if modifyEdfLines(edfcom):
         continue
   
   # Loop on the axes.  axisInd is 0 or 1, axisNum is 0 for single, 1/2 for a/b
   for axisInd in range(numAxes):
      if endingStep <= 0 or (dualAxis and ((axisInd and doOneAxis == 1) or 
                                           (not axisInd and doOneAxis > 1))):
         continue
      axisLet = ''
      axisEdfLet = 'a'
      axisUpperLet = 'A'
      if dualAxis:
         axisNum = axisInd + 1
         axisLet = 'a'
         addToMess = 'A'
         if axisInd:
            axisLet = 'b'
            axisEdfLet = 'b'
            axisUpperLet = 'B'
            addToMess = 'B' + logSuffixTag 
         prnLog('Starting axis ' + addToMess)
         if axisInd:
            prnLog('')

      if axisInd and doOneAxis > 1 and os.path.exists(setupBname):
         if runOneProcess(setupBname, message =
                          'Doing setup tasks now that second axis is present'):
            startingStep = 0
            axisFailed = True
            continue
            
      dataName = setName + axisLet
      axisCom = axisLet + comExt
      setRootAndExtension(dataName, typeExtension)
      while startingStep > 100:
         startingStep -= 100.
         continue
      if runOneAxis():

         # When an axis fails, mark failure to prevent going on to combine
         axisFailed = True
         if firstStart:
            startingStep = 0
            skipTiltalign = 0
         continue
      if firstStart:
         startingStep = 0
         skipTiltalign = 0
      if dualAxis:
         message = 'Completed axis ' + axisLet.upper() + ' of dataset ' + setName
      else:
         message = 'Completed dataset ' + setName
      if endingStep < finalStepNum:
         endPrint = int(round(endingStep))
         if math.fabs(endPrint - endingStep) > 0.005:
            endPrint = endingStep
         message += fmtstr(' through step {}', endPrint)
      if not dualAxis:
         (minutes, seconds, frac) = elapsedTimeComponents(startTime);
         message += fmtstr('  in {:02d}:{:02d}.{}   [brt4]', minutes, seconds, frac)
      prnLog(message)
      summaryMessage += message + '\n'

   # Try to combine dual axis dataset
   if dualAxis and doOneAxis != 1:
      axisNum = 0
      axisLet = 'a'
      axisInd = 2
      axisCom = comExt
      dataName = setName
      setRootAndExtension(dataName, typeExtension)
      if axisFailed:
         abortSet('One of the axes failed')
      elif not runCombine():
         (minutes, seconds, frac) = elapsedTimeComponents(startTime);
         prnLog(fmtstr('Completed dataset {}  in {:02d}:{:02d}.{}   [brt4]', setName,
                       minutes, seconds, frac))
      
   closeLogFileWriteEdf()
   prnLog('', flush = True);

sendEmail('Batchruntomo finished all data sets', summaryMessage)
mess = "Batch run finished; "
if finalRetVal:
   mess += fmtstr("failures occurred for {} datasets", finalRetVal);
else:
   mess += "no failures occurred"
prnstr(mess)
sys.exit(0)
