#include <mpi_classify_CL2D_core_analysis.h>

Inheritance diagram for ProgClassifyCL2DCore:

Collaboration diagram for ProgClassifyCL2DCore:

Public Member Functions
	ProgClassifyCL2DCore (int argc, char **argv)
	Empty constructor. More...

void	readParams ()
	Read argument from command line. More...

void	show ()
	Show. More...

void	defineParams ()
	Usage. More...

void	produceSideInfo ()
	Produce side info. More...

void	produceClassInfo ()
	Produce class input format. More...

void	computeCores ()
	Remove outliers. More...

void	computeStableCores ()
	Compute cores. More...

void	gatherResults (int firstLevel, const String &suffix)
	Gather results. More...

void	run ()

Public Member Functions inherited from XmippProgram
const char *	getParam (const char *param, int arg=0)

const char *	getParam (const char param, const char subparam, int arg=0)

int	getIntParam (const char *param, int arg=0)

int	getIntParam (const char param, const char subparam, int arg=0)

double	getDoubleParam (const char *param, int arg=0)

double	getDoubleParam (const char param, const char subparam, int arg=0)

float	getFloatParam (const char *param, int arg=0)

float	getFloatParam (const char param, const char subparam, int arg=0)

void	getListParam (const char *param, StringVector &list)

int	getCountParam (const char *param)

bool	checkParam (const char *param)

bool	existsParam (const char *param)

void	addParamsLine (const String &line)

void	addParamsLine (const char *line)

ParamDef *	getParamDef (const char *param) const

virtual void	quit (int exit_code=0) const

virtual int	tryRun ()

void	initProgress (size_t total, size_t stepBin=60)

void	setProgress (size_t value=0)

void	endProgress ()

void	processDefaultComment (const char param, const char left)

void	setDefaultComment (const char param, const char comment)

virtual void	initComments ()

void	setProgramName (const char *name)

void	addUsageLine (const char *line, bool verbatim=false)

void	clearUsage ()

void	addExampleLine (const char *example, bool verbatim=true)

void	addSeeAlsoLine (const char *seeAlso)

void	addKeywords (const char *keywords)

const char *	name () const

virtual void	usage (int verb=0) const

virtual void	usage (const String &param, int verb=2)

int	version () const

virtual void	show () const

virtual void	read (int argc, const char **argv, bool reportErrors=true)

virtual void	read (int argc, char **argv, bool reportErrors=true)

void	read (const String &argumentsLine)

	XmippProgram ()

	XmippProgram (int argc, const char **argv)

virtual	~XmippProgram ()

Public Attributes
FileName	fnRoot

FileName	fnODir

int	NPCA

double	thPCAZscore

int	tolerance

ClassifyCL2DCoreAction	action

std::shared_ptr< MpiNode >	node

std::unique_ptr< MpiTaskDistributor >	taskDistributor

std::vector< CL2DBlock >	blocks

int	maxLevel

size_t	Ydim

size_t	Xdim

Public Attributes inherited from XmippProgram
bool	doRun

bool	runWithoutArgs

int	verbose
	Verbosity level. More...

int	debug

Additional Inherited Members
Protected Member Functions inherited from XmippProgram
void	defineCommons ()

Protected Attributes inherited from XmippProgram
int	errorCode

ProgramDef *	progDef
	Program definition and arguments parser. More...

std::map< String, CommentList >	defaultComments

int	argc
	Original command line arguments. More...

const char **	argv

Detailed Description

Core analysis parameters.

Definition at line 53 of file mpi_classify_CL2D_core_analysis.h.

Constructor & Destructor Documentation

◆ ProgClassifyCL2DCore()

ProgClassifyCL2DCore::ProgClassifyCL2DCore	(	int	argc,
		char **	argv
	)

Empty constructor.

Definition at line 38 of file mpi_classify_CL2D_core_analysis.cpp.

 {
     node=std::make_shared<MpiNode>(argc, argv);
     if (!node->isMaster())
         verbose=0;
     taskDistributor=nullptr;
     maxLevel=-1;
     tolerance=0;
     thPCAZscore=3;
 }

Member Function Documentation

◆ computeCores()

void ProgClassifyCL2DCore::computeCores ( )

Remove outliers.

Definition at line 145 of file mpi_classify_CL2D_core_analysis.cpp.

 {
     if (verbose && node->rank==0)
         std::cerr << "Computing cores ...\n";
     ProgAnalyzeCluster analyzeCluster;
     analyzeCluster.verbose=0;
     analyzeCluster.NPCA=NPCA;
     analyzeCluster.Niter=10;
     analyzeCluster.distThreshold=thPCAZscore;
     analyzeCluster.dontMask=false;
 
     MetaDataVec MD;
     size_t first, last;
     size_t Nblocks=blocks.size();
     if (verbose && node->rank==0)
         init_progress_bar(Nblocks);
     while (taskDistributor->getTasks(first, last))
         for (size_t idx=first; idx<=last; ++idx)
         {
             // Remove outliers in the PCA projection
             analyzeCluster.SFin.clear();
             analyzeCluster.fnSel=blocks[idx].block+"@"+blocks[idx].fnLevel;
             analyzeCluster.fnOut=blocks[idx].fnLevel.insertBeforeExtension((String)"_core_"+blocks[idx].block);
             analyzeCluster.run();
 
             // Remove outliers from file
             MD.read(analyzeCluster.fnOut);
             MD.removeDisabled();
             MD.write(analyzeCluster.fnOut,MD_APPEND);
 
             if (verbose && node->rank==0)
                 progress_bar(idx);
         }
     taskDistributor->wait();
     if (verbose && node->rank==0)
         progress_bar(Nblocks);
 
     // Gather all results
     gatherResults(0,"core");
 }

◆ computeStableCores()

void ProgClassifyCL2DCore::computeStableCores ( )

Compute cores.

Definition at line 187 of file mpi_classify_CL2D_core_analysis.cpp.

 {
     if (verbose && node->rank==0)
         std::cerr << "Computing stable cores ...\n";
     MetaDataDb thisClass, anotherClass, commonImages, thisClassCore;
     size_t first, last;
     Matrix2D<unsigned char> coocurrence;
     Matrix1D<unsigned char> maximalCoocurrence;
     int Nblocks=blocks.size();
     taskDistributor->reset();
     std::vector<size_t> commonIdx;
     std::map<String,size_t> thisClassOrder;
     String fnImg;
     while (taskDistributor->getTasks(first, last))
         for (size_t idx=first; idx<=last; ++idx)
         {
             // Read block
             CL2DBlock &thisBlock=blocks[idx];
             if (thisBlock.level<=tolerance)
                 continue;
             if (!existsBlockInMetaDataFile(thisBlock.fnLevelCore, thisBlock.block))
                 continue;
             thisClass.read(thisBlock.block+"@"+thisBlock.fnLevelCore);
             thisClassCore.clear();
 
             // Add MDL_ORDER
             if (thisClass.size()>0)
             {
                 size_t order=0;
                 thisClassOrder.clear();
                 for (size_t objId : thisClass.ids())
                 {
                     thisClass.getValue(MDL_IMAGE,fnImg,objId);
                     thisClassOrder[fnImg]=order++;
                 }
 
                 // Calculate coocurrence within all blocks whose level is inferior to this
                 size_t NthisClass=thisClass.size();
                 if (NthisClass>0)
                 {
                     try {
                        coocurrence.initZeros(NthisClass,NthisClass);
                     } catch (XmippError &e)
                     {
                        std::cerr << e.what() << std::endl;
                        std::cerr << "There is a memory allocation error. Most likely there are too many images in this class ("
                                  << NthisClass << " images). Consider increasing the number of initial and final classes\n";
                        REPORT_ERROR(ERR_MEM_NOTENOUGH,"While computing stable class");
                     }
                     for (int n=0; n<Nblocks; n++)
                     {
                         CL2DBlock &anotherBlock=blocks[n];
                         if (anotherBlock.level>=thisBlock.level)
                             break;
                         if (!existsBlockInMetaDataFile(anotherBlock.fnLevelCore, anotherBlock.block))
                             continue;
                         anotherClass.read(anotherBlock.block+"@"+anotherBlock.fnLevelCore);
                         anotherClass.intersection(thisClass,MDL_IMAGE);
                         commonImages.join1(anotherClass, thisClass, MDL_IMAGE,LEFT);
                         commonIdx.resize(commonImages.size());
                         size_t idx=0;
                         for (size_t objId : commonImages.ids())
                         {
                             commonImages.getValue(MDL_IMAGE,fnImg,objId);
                             commonIdx[idx++]=thisClassOrder[fnImg];
                         }
                         size_t Ncommon=commonIdx.size();
                         for (size_t i=0; i<Ncommon; i++)
                         {
                             size_t idx_i=commonIdx[i];
                             for (size_t j=i+1; j<Ncommon; j++)
                             {
                                 size_t idx_j=commonIdx[j];
                                 MAT_ELEM(coocurrence,idx_i,idx_j)+=1;
                             }
                         }
                     }
                 }
 
                 // Take only those elements whose coocurrence is maximal
                 maximalCoocurrence.initZeros(NthisClass);
                 int aimedCoocurrence=thisBlock.level-tolerance;
                 FOR_ALL_ELEMENTS_IN_MATRIX2D(coocurrence)
                 if (MAT_ELEM(coocurrence,i,j)==aimedCoocurrence)
                     VEC_ELEM(maximalCoocurrence,i)=VEC_ELEM(maximalCoocurrence,j)=1;
 
                 // Now compute core
                 for (size_t objId : thisClass.ids())
                 {
                     thisClass.getValue(MDL_IMAGE,fnImg,objId);
                     size_t idx=thisClassOrder[fnImg];
                     if (VEC_ELEM(maximalCoocurrence,idx))
                         thisClassCore.addRow(*thisClass.getRow(objId));
                 }
             }
             thisClassCore.write(thisBlock.fnLevel.insertBeforeExtension((String)"_stable_core_"+thisBlock.block),MD_APPEND);
         }
     taskDistributor->wait();
 
     // Gather all results
     gatherResults(tolerance+1,"stable_core");
 }

◆ defineParams()

void ProgClassifyCL2DCore::defineParams ( )

virtual

Usage.

Reimplemented from XmippProgram.

Definition at line 82 of file mpi_classify_CL2D_core_analysis.cpp.

 {
     addUsageLine("Compute the core of a CL2D clustering");
     addParamsLine("    --root <rootname>          : Rootname of the CL2D");
     addParamsLine("    --dir <dir>                : Output directory of the CL2D");
     addParamsLine("    --computeCore <thPCAZscore=3> <NPCA=2>: The class cores are computed by thresholding the Zscore of");
     addParamsLine("                               : the class images and their projections onto a nD PCA space (by default, n=2)");
     addParamsLine("or  --computeStableCore <tolerance=1> : The stable core is formed by all the images in the class that have been");
     addParamsLine("                               : in the same class (except in tolerance levels) in the whole hierarchy.");
     addParamsLine("                               : If tolerance=0, then the stable core is formed by the set of images that");
     addParamsLine("                               : have been always together. For this reason, the stable core can be computed only");
     addParamsLine("                               : for the level such that level>tolerance");
     addExampleLine("mpirun -np 4 `which xmipp_mpi_classify_CL2D_core_analysis` -i 2D/CL2D/run_001/results --computeCore");
     addExampleLine("mpirun -np 4 `which xmipp_mpi_classify_CL2D_core_analysis` -i 2D/CL2D/run_001/results --computeStableCore");
 }

◆ gatherResults()

void ProgClassifyCL2DCore::gatherResults	(	int	firstLevel,
		const String &	suffix
	)

Gather results.

Definition at line 290 of file mpi_classify_CL2D_core_analysis.cpp.

 {
     node->barrierWait();
     if (node->rank==0)
     {
         FileName fnBlock, fnClass, fnSummary, fnSummaryOriginal;
         Image<double> classAverage;
         // Compute class averages
         MetaDataVec classes, MD, MDoriginal;
         int Nblocks=blocks.size();
         for (int level=firstLevel; level<=maxLevel; level++)
         {
             classes.clear();
             fnSummary=formatString("%s/level_%02d/%s_classes_%s",fnODir.c_str(),level,fnRoot.c_str(),suffix.c_str());
             for (int idx=0; idx<Nblocks; idx++)
             {
                 if (blocks[idx].level!=level)
                     continue;
                 fnBlock=fnSummary+"_"+blocks[idx].block+".xmd";
                 if (fileExists(fnBlock))
                 {
                     MD.read(fnBlock);
                     MDoriginal.read(formatString("%s@%s/level_%02d/%s_classes.xmd",blocks[idx].block.c_str(),fnODir.c_str(),level,fnRoot.c_str()));
                     int classNo=textToInteger(blocks[idx].block.substr(6,6));
                     size_t classSize=MD.size();
                     fnClass.compose(classNo,fnSummary,"stk");
                     if (classSize>0)
                         getAverageApplyGeo(MD, classAverage());
                     else
                         classAverage().initZeros(Ydim,Xdim);
                     classAverage.write(fnClass);
 
                     size_t id=classes.addObject();
                     classes.setValue(MDL_REF,classNo,id);
                     classes.setValue(MDL_IMAGE,fnClass,id);
                     classes.setValue(MDL_CLASS_COUNT,classSize,id);
                     classes.setValue(MDL_MODELFRAC,((double)classSize)/MDoriginal.size(),id);
                 }
             }
             classes.write((String)"classes@"+fnSummary+".xmd");
         }
 
         // Write the rest of blocks
         for (int idx=0; idx<Nblocks; idx++)
         {
             fnSummary=formatString("%s/level_%02d/%s_classes_%s",fnODir.c_str(),blocks[idx].level,fnRoot.c_str(),suffix.c_str());
             fnBlock=fnSummary+"_"+blocks[idx].block+".xmd";
             if (fileExists(fnBlock))
             {
                 MD.read(fnBlock);
                 unlink(fnBlock.c_str());
                 MD.write(blocks[idx].block+"@"+fnSummary+".xmd",MD_APPEND);
             }
         }
     }
     node->barrierWait();
 }

◆ produceClassInfo()

void ProgClassifyCL2DCore::produceClassInfo ( )

Produce class input format.

Definition at line 348 of file mpi_classify_CL2D_core_analysis.cpp.

 {
     node->barrierWait();
     if (node->rank==0)
     {
         FileName fnSummaryOriginal;
         MetaDataVec MDoriginal;
         int Nblocks=blocks.size();
 
         // Read and write reference to 2D classes blocks
         for (int idx=0; idx<Nblocks; idx++)
         {
              int classNo=textToInteger(blocks[idx].block.substr(6,6));
              fnSummaryOriginal=formatString("%s@%s/level_%02d/%s_classes.xmd",blocks[idx].block.c_str(),fnODir.c_str(),
                                             blocks[idx].level,fnRoot.c_str());
              MDoriginal.read(fnSummaryOriginal);
              MDoriginal.fillConstant(MDL_REF, integerToString(classNo));
              MDoriginal.write(fnSummaryOriginal,MD_APPEND);
         }
     }
     node->barrierWait();
 }

◆ produceSideInfo()

void ProgClassifyCL2DCore::produceSideInfo ( )

Produce side info.

Definition at line 99 of file mpi_classify_CL2D_core_analysis.cpp.

 {
     // Get maximum CL2D level
     maxLevel=0;
     FileName fnLevel;
     do
         fnLevel=formatString("%s/level_%02d/%s_classes.xmd",fnODir.c_str(),maxLevel++,fnRoot.c_str());
     while (fnLevel.exists());
     maxLevel-=2;
     if (maxLevel==-1)
         REPORT_ERROR(ERR_ARG_MISSING,"Cannot find any CL2D analysis in the directory given");
 
     // Read all the blocks available in all MetaData
     StringVector blocksAux;
     CL2DBlock block;
     for (int level=0; level<=maxLevel; level++)
     {
         fnLevel=formatString("%s/level_%02d/%s_classes.xmd",fnODir.c_str(),level,fnRoot.c_str());
         getBlocksInMetaDataFile(fnLevel,blocksAux);
         block.level=level;
         block.fnLevel=fnLevel;
         block.fnLevelCore=fnLevel.insertBeforeExtension("_core");
         for (size_t i=0; i<blocksAux.size(); i++)
         {
             if (blocksAux[i].find("class")!=std::string::npos && 
                 blocksAux[i].find("images")!=std::string::npos)
             {
                 block.block=blocksAux[i];
                 blocks.push_back(block);
             }
         }
     }
 
     // Create a task file distributor for all blocks
     size_t Nblocks=blocks.size();
     taskDistributor=std::make_unique<MpiTaskDistributor>(Nblocks,1,node);
 
     // Get image dimensions
     if (Nblocks>0)
     {
         size_t Zdim, Ndim;
         getImageSizeFromFilename(blocks[0].block+"@"+blocks[0].fnLevel,Xdim,Ydim,Zdim,Ndim);
     }
 }

◆ readParams()

void ProgClassifyCL2DCore::readParams ( )

virtual

Read argument from command line.

Reimplemented from XmippProgram.

Definition at line 50 of file mpi_classify_CL2D_core_analysis.cpp.

 {
     fnRoot = getParam("--root");
     fnODir = getParam("--dir");
     if (checkParam("--computeCore"))
     {
         thPCAZscore = getDoubleParam("--computeCore",0);
         NPCA = getIntParam("--computeCore",1);
         action=COMPUTE_CORE;
     }
     else if (checkParam("--computeStableCore"))
     {
         tolerance = getIntParam("--computeStableCore",0);
         action=COMPUTE_STABLE_CORE;
     }
 }

◆ run()

void ProgClassifyCL2DCore::run ( )

virtual

Run.

Reimplemented from XmippProgram.

Definition at line 372 of file mpi_classify_CL2D_core_analysis.cpp.

 {
     show();
     produceSideInfo();
     produceClassInfo();
     if (action==COMPUTE_CORE)
         computeCores();
     else
         computeStableCores();
 }

◆ show()

void ProgClassifyCL2DCore::show ( )

Show.

Definition at line 68 of file mpi_classify_CL2D_core_analysis.cpp.

 {
     if (!verbose)
         return;
     std::cout << "CL2D rootname:        " << fnRoot << std::endl
     << "CL2D output dir:      " << fnODir << std::endl;
     if (action==COMPUTE_STABLE_CORE)
         std::cout << "Tolerance:                " << tolerance << std::endl;
     else
         std::cout << "Threshold PCA Zscore:     " << thPCAZscore << std::endl
                   << "Number of PCA dimensions: " << NPCA << std::endl;
 }