def _getIterClasses(self, it, clean=False):
        """ Return the file with the classes for this iteration.
        If the file doesn't exists, it will be created. 
        """
        from convert import readSetOfClasses3D
        
        dataClasses = self.protocol._getFileName('classes_scipion', iter=it)
        
        if clean:
            cleanPath(dataClasses)
        
        if not exists(dataClasses):
            fileparList = []
            volumeList = []
            for ref in self._refsList:
                filepar = self.protocol._getFileName('output_par_class', iter=it, ref=ref)
                volFn = self.protocol._getFileName('iter_vol_class', iter=it, ref=ref)
                fileparList.append(filepar)
                volumeList.append(volFn)
            
            clsSet = em.SetOfClasses3D(filename=dataClasses)
            clsSet.setImages(self.inputParticles.get())
            readSetOfClasses3D(clsSet, fileparList, volumeList)
            clsSet.write()
            clsSet.close()

        return dataClasses

    def createOutputStep(self):
        fnImgs = self._getExtraPath('images.stk')
        if os.path.exists(fnImgs):
            cleanPath(fnImgs)

        outputSet = self._createSetOfParticles()
        imgSet = self.inputSet.get()
        imgFn = self._getExtraPath("anglesCont.xmd")
        self.newAssignmentPerformed = os.path.exists(self._getExtraPath("angles.xmd"))
        self.samplingRate = self.inputSet.get().getSamplingRate()
        if isinstance(imgSet, SetOfClasses2D):
            outputSet = self._createSetOfClasses2D(imgSet)
            outputSet.copyInfo(imgSet.getImages())
        elif isinstance(imgSet, SetOfAverages):
            outputSet = self._createSetOfAverages()
            outputSet.copyInfo(imgSet)
        else:
            outputSet = self._createSetOfParticles()
            outputSet.copyInfo(imgSet)
            if not self.newAssignmentPerformed:
                outputSet.setAlignmentProj()
        outputSet.copyItems(imgSet,
                            updateItemCallback=self._processRow,
                            itemDataIterator=md.iterRows(imgFn, sortByLabel=md.MDL_ITEM_ID))
        self._defineOutputs(outputParticles=outputSet)
        self._defineSourceRelation(self.inputSet, outputSet)

 def getBestVolumesStep(self):
     volumes = []
     inliers = []
     
     for n in range(self.nRansac.get()):
         fnAngles = self._getTmpPath("angles_ransac%05d"%n+".xmd")
         md=xmipp.MetaData("[email protected]"+fnAngles)
         numInliers=md.getValue(xmipp.MDL_WEIGHT,md.firstObject())
         volumes.append(fnAngles)
         inliers.append(numInliers)
     
     index = sorted(range(inliers.__len__()), key=lambda k: inliers[k])
     fnBestAngles = ''
     threshold=self.getCCThreshold()
  
     i = self.nRansac.get()-1
     indx = 0
     while i >= 0 and indx < self.numVolumes:
         fnBestAngles = volumes[index[i]]
         fnBestAnglesOut = self._getPath("proposedVolume%05d"%indx+".xmd")
         copyFile(fnBestAngles, fnBestAnglesOut)
         self._log.info("Best volume %d = %s" % (indx, fnBestAngles))
         if not self.useAll:
             self.runJob("xmipp_metadata_utilities","-i %s -o %s --query select \"maxCC>%f \" --mode append" %(fnBestAnglesOut,fnBestAnglesOut,threshold))
             if not isMdEmpty(fnBestAnglesOut):
                 indx += 1
         else:
             indx += 1
         i -= 1
         
     # Remove unnecessary files
     for n in range(self.nRansac.get()):
         fnAngles = self._getTmpPath("angles_ransac%05d"%n+".xmd")
         cleanPath(fnAngles)

def exportData(emxDir, inputSet, ctfSet=None, xmlFile='data.emx', binaryFile=None):
    """ Export micrographs, coordinates or particles to  EMX format. """
    cleanPath(emxDir)
    makePath(emxDir) 
    emxData = emxlib.EmxData()
    micSet=None
    
    if binaryFile is None:
        binaryFile = xmlFile.replace('.emx', '.mrc')

    if isinstance(inputSet, SetOfMicrographs):
        _micrographsToEmx(emxData, inputSet, emxDir, ctfSet)
        
    elif isinstance(inputSet, SetOfCoordinates):
        micSet = inputSet.getMicrographs()
        _micrographsToEmx(emxData, micSet, emxDir, ctfSet)
        _particlesToEmx(emxData, inputSet, None, micSet)
        
    elif isinstance(inputSet, SetOfParticles):
        print ("SetOfParticles-----------------------------------------")
        if inputSet.hasCoordinates():
            micSet = inputSet.getCoordinates().getMicrographs()
            _micrographsToEmx(emxData, micSet, emxDir, writeData=False)
        fnMrcs = join(emxDir, binaryFile)
        _particlesToEmx(emxData, inputSet, fnMrcs, micSet)
        
    fnXml = join(emxDir, xmlFile)
    emxData.write(fnXml)

    def __createTemporaryCtfs(self, obj, setOfMics):
        """ Create a temporary .sqlite file to visualize CTF while the
             protocol has not finished yet.
            """
        cleanPath(obj._getPath("ctfs_temporary.sqlite"))
        ctfSet = self.protocol._createSetOfCTF("_temporary")

        for mic in setOfMics:
            micFn = mic.getFileName()
            micDir = obj._getExtraPath(removeBaseExt(mic.getFileName()))
            samplingRate = mic.getSamplingRate() * self.protocol.ctfDownFactor.get()
            mic.setSamplingRate(samplingRate)
            out = self.protocol._getCtfOutPath(micDir)
            psdFile = self.protocol._getPsdPath(micDir)

            if exists(out) and exists(psdFile):
                ctfModel = em.CTFModel()
                readCtfModel(ctfModel, out,
                             ctf4=self.protocol.useCtffind4.get())
                ctfModel.setPsdFile(psdFile)
                ctfModel.setMicrograph(mic)
                ctfSet.append(ctfModel)

        if not ctfSet.isEmpty():
            ctfSet.write()
            ctfSet.close()

        return ctfSet

    def _calculateAuxiliaryFile(self):
        """create new ctf_set with ctf that satisfies the
        constraint and persist it
        """
        try:
            self.setOfCTFsConst
        except AttributeError:
            pass
        else:
            #TODO close the mapper, if not the object cannot be reused (Although it should be able)
            self.setOfCTFsConst.close()

        #metadata file with protocol output
        #get temporary fileName for metadata file
        self.targetFile = self.protocol._getTmpPath(self.tmpMetadataFile)
        resolutionThreshold = self.resolutionThreshold.get()
        print "TODO: this should be closer to the mapper. Here it does not make any sense. ROB"
        #TODO check if this is necessary
        cleanPath(self.targetFile)
        
        #metadata with selected CTFs
        self.setOfCTFsConst  = data.SetOfCTF(filename=self.targetFile)
        #object read metadata file
        ctfs  = data.SetOfCTF(filename=self.pairsFile)
        #condition to be satisfized for CTFs
        for ctf in ctfs:
            if ctf.resolution < resolutionThreshold:
                self.setOfCTFsConst.append(ctf)
        #new file with selected CTFs
        self.setOfCTFsConst.write()
        #check if empty
        if self.setOfCTFsConst.getSize() < 1:
            print "WARNING: Empty set of CTFs."

    def processMovieStep(self, movieId, movieFn, *args):
        movieFolder = self._getMovieFolder(movieId)
        movieName = basename(movieFn)
        #export SCIPION_DEBUG=1 # passwd=a
        #startDebugger()

        if self._filterMovie(movieId, movieFn):
            makePath(movieFolder)
            createLink(movieFn, join(movieFolder, movieName))
            toDelete = [movieName]
    
            if movieName.endswith('bz2'):
                movieMrc = movieName.replace('.bz2', '') # we assume that if compressed the name ends with .mrc.bz2
                toDelete.append(movieMrc)
                if not exists(movieMrc):
                    self.runJob('bzip2', '-d -f %s' % movieName, cwd=movieFolder)
            else:
                movieMrc = movieName
            
            self.info("Processing movie: %s" % movieMrc)
            
            if movieMrc.endswith('.em'):
                movieMrc = movieMrc + ":ems"

            self._processMovie(movieId, movieMrc, movieFolder, *args)
            
            if self.cleanMovieData:
                cleanPath(movieFolder)
            else:
                self.info('Clean movie data DISABLED. Movie folder will remain in disk!!!')

    def reformatPdbOutputStep(self, numberOfModes):
        self._enterWorkingDir()
        
        makePath('modes')
        Natoms = self._countAtoms("atoms.pdb")
        fhIn = open('diagrtb.eigenfacs')
        fhAni = open('vec_ani.txt','w')
        
        for n in range(numberOfModes):
            # Skip two lines
            fhIn.readline()
            fhIn.readline()
            fhOut=open('modes/vec.%d'%(n+1),'w')
            for i in range(Natoms):
                line=fhIn.readline()
                fhOut.write(line)
                fhAni.write(line.rstrip().lstrip()+" ")
            fhOut.close()
            if n!=(numberOfModes-1):
                fhAni.write("\n")
        fhIn.close()
        fhAni.close()
        self.runJob("nma_prepare_for_animate.py","",env=getNMAEnviron())
        cleanPath("vec_ani.txt")
        moveFile('vec_ani.pkl', 'extra/vec_ani.pkl')

        self._leaveWorkingDir()

 def _createCluster(self):
     """ Create the cluster with the selected particles
     from the cluster. This method will be called when
     the button 'Create Cluster' is pressed.
     """
     # Write the particles
     prot = self.protocol
     project = prot.getProject()
     inputSet = prot.getInputParticles()
     fnSqlite = prot._getTmpPath('cluster_particles.sqlite')
     cleanPath(fnSqlite)
     partSet = SetOfParticles(filename=fnSqlite)
     partSet.copyInfo(inputSet)
     for point in self.getData():
         if point.getState() == Point.SELECTED:
             particle = inputSet[point.getId()]
             partSet.append(particle)
     partSet.write()
     partSet.close()
             
     from protocol_batch_cluster import BatchProtNMACluster
     newProt = project.newProtocol(BatchProtNMACluster)
     clusterName = self.clusterWindow.getClusterName()
     if clusterName:
         newProt.setObjLabel(clusterName)
     newProt.inputNmaDimred.set(prot)
     newProt.sqliteFile.set(fnSqlite)
     
     project.launchProtocol(newProt)

 def cleanPrime(self):
     self._enterDir(self._getExtraPath())
     cleanPath("cmdline.txt")
     cleanPattern("*.txt")
     cleanPattern("startvol_state*.spi")
     # Get last iteration
     for i in range(1, self.getLastIteration()):
         cleanPattern("recvol_state*_iter%d.spi"%i)
     self._leaveDir()

    def visualizeClasses(self, e=None):
        prot = self.protocol
        classDir = prot.getClassDir()
        classAvg = 'classavg'
        classVar = 'classvar'
        classDoc = 'docclass'
        
        params = {'[class_dir]': classDir,
                  '[desired-classes]': self.numberOfClasses.get(),
                  '[particles]': prot._params['particles'] + '@******',
                  '[class_doc]': join(classDir, classDoc + '***'), 
                  '[class_avg]': join(classDir, classAvg + '***'),
                  '[class_var]': join(classDir, classVar + '***'),        
                  }
        
        prot.runTemplate('mda/classavg.msa', prot.getExt(), params)

        particles = prot.inputParticles.get()
        particles.load()
        sampling = particles.getSamplingRate()
        
        setFn = prot._getTmpPath('classes2D.sqlite')
        cleanPath(setFn)
        classes2D = SetOfClasses2D(filename=setFn)
        classes2D.setImages(particles)

        # We need to first create a map between the particles index and
        # the assigned class number
        classDict = {}
        for classId in range(1, self.numberOfClasses.get()+1):
            docClass = prot._getPath(classDir, classDoc + '%03d.stk' % classId)
            doc = SpiderDocFile(docClass)
            for values in doc.iterValues():
                imgIndex = int(values[0])
                classDict[imgIndex] = classId
            doc.close()

        updateItem = lambda p, i: p.setClassId(classDict[i])

        def updateClass(cls):
            rep = cls.getRepresentative()
            rep.setSamplingRate(particles.getSamplingRate())
            avgFn = prot._getPath(classDir,
                                  classAvg + '%03d.stk' % cls.getObjId())
            rep.setLocation(1, avgFn)

        particlesRange = range(1, particles.getSize()+1)
        classes2D.classifyItems(updateItemCallback=updateItem,
                                updateClassCallback=updateClass,
                                itemDataIterator=iter(particlesRange))

        classes2D.write()
        classes2D.close()

        return [ClassesView(self.getProject(), prot.strId(),
                            classes2D.getFileName(), particles.strId())]

    def visualizeClasses(self, e=None):
        prot = self.protocol
        classDir = prot.getClassDir()
        classAvg = 'classavg'
        classVar = 'classvar'
        classDoc = 'docclass'
        ext = prot.getExt()
        
        params = {'[class_dir]': classDir,
                  '[desired-classes]': self.numberOfClasses.get(),
                  '[particles]': prot._params['particles'] + '@******',
                  '[class_doc]': join(classDir, classDoc + '***'), 
                  '[class_avg]': join(classDir, classAvg + '***'),
                  '[class_var]': join(classDir, classVar + '***'),        
                  }
        
        prot.runTemplate('mda/classavg.msa', prot.getExt(), params)

        particles = prot.inputParticles.get()
        particles.load()
        sampling = particles.getSamplingRate()
        
        setFn = self._getPath('classes2D.sqlite')
        cleanPath(setFn)
        classes2D = SetOfClasses2D(filename=setFn)
        classes2D.setImages(particles)
            
        for classId in range(1, self.numberOfClasses.get()+1):
            class2D = Class2D()
            class2D.setObjId(classId)
            
            avgImg = Particle()
            avgImg.setSamplingRate(sampling)
            avgFn = prot._getPath(classDir, classAvg + '%03d.stk' % classId)
            avgImg.setLocation(1, avgFn)
            #avgImg.setLocation(classId, 'classavg.stk')
            
            class2D.setRepresentative(avgImg)
            classes2D.append(class2D)
            
            docClass = prot._getPath(classDir, classDoc + '%03d.stk' % classId)
            doc = SpiderDocFile(docClass)
            
            for values in doc.iterValues():
                imgId = int(values[0])
                img = particles[imgId]
                class2D.append(img)
                
            classes2D.update(class2D)            
        classes2D.write()
        classes2D.close()

        return [ClassesView(self.getProject(),
                            prot.strId(), classes2D.getFileName(), 
                            prot.inputParticles.get().strId())]
                              

 def evaluateResiduals(self):
     # Evaluate each image
     fnAutoCorrelations = self._getExtraPath("autocorrelations.xmd")
     stkAutoCorrelations = self._getExtraPath("autocorrelations.stk")
     stkResiduals = self._getExtraPath("residuals.stk")
     anglesOutFn=self._getExtraPath("anglesCont.xmd")
     self.runJob("xmipp_image_residuals", " -i %s -o %s --save_metadata_stack %s" % (stkResiduals, stkAutoCorrelations, fnAutoCorrelations), numberOfMpi=1)
     self.runJob("xmipp_metadata_utilities", '-i %s --operate rename_column "image imageResidual"' % fnAutoCorrelations, numberOfMpi=1)
     self.runJob("xmipp_metadata_utilities", '-i %s --set join %s imageResidual' % (anglesOutFn, fnAutoCorrelations), numberOfMpi=1)
     cleanPath(fnAutoCorrelations)

def writeSqliteIterData(imgStar, imgSqlite, **kwargs):
    """ Given a Relion images star file (from some iteration)
    create the corresponding SetOfParticles (sqlite file)
    for this iteration. This file can be visualized sorted
    by the LogLikelihood.
    """
    cleanPath(imgSqlite)
    imgSet = em.SetOfParticles(filename=imgSqlite)
    readSetOfParticles(imgStar, imgSet, **kwargs)
    imgSet.write()

    def __createSet(self, SetClass, template, suffix):
        """ Create a set and set the filename using the suffix. 
        If the file exists, it will be delete. """
        setFn = self._getPath(template % suffix)
        # Close the connection to the database if
        # it is open before deleting the file
        cleanPath(setFn)

        SqliteDb.closeConnection(setFn)
        setObj = SetClass(filename=setFn)
        return setObj

 def evaluateStep(self, outImgsFn):
     # Evaluate each image
     fnAutoCorrelations = self._getExtraPath("autocorrelations.xmd")
     stkAutoCorrelations = self._getExtraPath("autocorrelations.stk")
     stkDiff = self._getExtraPath("diff.stk")
     args1 = " -i %s -o %s --save_metadata_stack %s"
     args2 = " -i %s --set merge %s"
     outClasses = '[email protected]' + outImgsFn
     self.runJob("xmipp_image_residuals", args1 % (stkDiff, stkAutoCorrelations, fnAutoCorrelations), numberOfMpi=1)
     self.runJob("xmipp_metadata_utilities", '-i %s --operate rename_column "image image1"' % fnAutoCorrelations, numberOfMpi=1)
     self.runJob("xmipp_metadata_utilities", args2 % (outClasses, fnAutoCorrelations), numberOfMpi=1)
     cleanPath(fnAutoCorrelations)

 def produceResiduals(self, fnVol, fnAngles, Ts):
     if fnVol.endswith(".mrc"):
         fnVol+=":mrc"
     anglesOutFn=self._getExtraPath("anglesCont.stk")
     residualsOutFn=self._getExtraPath("residuals.stk")
     projectionsOutFn=self._getExtraPath("projections.stk")
     xdim=self.inputVolume.get().getDim()[0]
     self.runJob("xmipp_angular_continuous_assign2", "-i %s -o %s --ref %s --optimizeAngles --optimizeGray --optimizeShift --max_shift %d --oresiduals %s --oprojections %s --sampling %f" %\
                 (fnAngles,anglesOutFn,fnVol,floor(xdim*0.05),residualsOutFn,projectionsOutFn,Ts))
     fnNewParticles=self._getExtraPath("images.stk")
     if os.path.exists(fnNewParticles):
         cleanPath(fnNewParticles)

 def computeModesStep(self, fnPseudoatoms, numberOfModes, cutoffStr):
     (baseDir,fnBase)=os.path.split(fnPseudoatoms)
     fnBase=fnBase.replace(".pdb","")
     fnDistanceHist=os.path.join(baseDir,'extra',fnBase+'_distance.hist')
     rc = self._getRc(fnDistanceHist)
     self._enterWorkingDir()
     self.runJob('nma_record_info.py', "%d %s.pdb %d" % (numberOfModes, fnBase, rc),env=getNMAEnviron())
     self.runJob("nma_pdbmat.pl","pdbmat.dat",env=getNMAEnviron())
     self.runJob("nma_diag_arpack","",env=getNMAEnviron())
     if not exists("fort.11"):
         self._printWarnings(redStr("Modes cannot be computed. Check the number of modes you asked to compute and/or consider increasing cut-off distance. The maximum number of modes allowed by the method for pseudoatomic normal mode analysis is 3 times the number of pseudoatoms but the protocol allows only up to 200 modes as 20-100 modes are usually enough.  If the number of modes is below the minimum between 200 and 3 times the number of pseudoatoms, consider increasing cut-off distance."))
     cleanPath("diag_arpack.in", "pdbmat.dat")
     self._leaveWorkingDir()

    def convertInputStep(self, particlesId):
        """ Create the input file in STAR format as expected by Relion.
        If the input particles comes from Relion, just link the file.
        Params:
            particlesId: use this parameters just to force redo of convert if 
                the input particles are changed.
        """
        imgSet = self._getInputParticles()
        imgStar = self._getFileName('input_star')

        self.info("Converting set from '%s' into '%s'" %
                           (imgSet.getFileName(), imgStar))
        
        # Pass stack file as None to avoid write the images files
        writeSetOfParticles(imgSet, imgStar, self._getExtraPath())
        
        if self.doCtfManualGroups:
            self._splitInCTFGroups(imgStar)
        
        if not self.IS_CLASSIFY:
            if self.realignMovieFrames:
                movieParticleSet = self.inputMovieParticles.get()
                
                auxMovieParticles = self._createSetOfMovieParticles(suffix='tmp')
                auxMovieParticles.copyInfo(movieParticleSet)
                # Discard the movie particles that are not present in the refinement set
                for movieParticle in movieParticleSet:
                    particle = imgSet[movieParticle.getParticleId()]
                    if particle is not None:
                        auxMovieParticles.append(movieParticle)
                writeSetOfParticles(auxMovieParticles,
                                    self._getFileName('movie_particles'), None, originalSet=imgSet,
                                    postprocessImageRow=self._postprocessImageRow)
                mdMovies = md.MetaData(self._getFileName('movie_particles'))
                mdParts = md.MetaData(self._getFileName('input_star'))

                if getVersion() == "1.4":
                    mdParts.renameColumn(md.RLN_IMAGE_NAME, md.RLN_PARTICLE_ORI_NAME)
                else:
                    mdParts.renameColumn(md.RLN_IMAGE_NAME, md.RLN_PARTICLE_NAME)
                mdParts.removeLabel(md.RLN_MICROGRAPH_NAME)
                
                detectorPxSize = movieParticleSet.getAcquisition().getMagnification() * movieParticleSet.getSamplingRate() / 10000
                mdAux = md.MetaData()
                mdMovies.fillConstant(md.RLN_CTF_DETECTOR_PIXEL_SIZE, detectorPxSize)
                
                mdAux.join2(mdMovies, mdParts, md.RLN_PARTICLE_ID, md.RLN_IMAGE_ID, md.INNER_JOIN)
                
                mdAux.write(self._getFileName('movie_particles'), md.MD_OVERWRITE)
                cleanPath(auxMovieParticles.getFileName())

    def createVolumesSqlite(self, files, path, samplingRate):
        from em import SetOfVolumes, Volume
        cleanPath(path)
        volSet = SetOfVolumes(filename=path)
        volSet.setSamplingRate(samplingRate)

        for volFn in files:
            vol = Volume()
            vol.setFileName(volFn)
            volSet.append(vol)
        volSet.write()
        volSet.close()
        
        return volSet