CL2D Class Reference

#include <mpi_classify_CL2D.h>

Collaboration diagram for CL2D:

Public Member Functions
	CL2D ()
	~CL2D ()
	CL2D (const CL2D &)=delete
	CL2D (const CL2D &&)=delete
CL2D &	operator= (const CL2D &)=delete
CL2D &	operator= (const CL2D &&)=delete
void	readImage (Image< double > &I, size_t objId, bool applyGeo) const
	Read Image. More...
void	initialize (MetaDataDb &_SF, std::vector< MultidimArray< double > > &_codes0)
	Initialize. More...
void	shareAssignments (bool shareAssignment, bool shareUpdates, bool shareNonCorr)
	Share assignments. More...
void	shareSplitAssignments (Matrix1D< int > &assignment, CL2DClass *node1, CL2DClass *node2) const
	Share split assignment. More...
void	write (const FileName &fnODir, const FileName &fnRoot, int level) const
	Write the nodes. More...
void	lookNode (MultidimArray< double > &I, int oldnode, int &newnode, CL2DAssignment &bestAssignment)
void	transferUpdates ()
void	run (const FileName &fnODir, const FileName &fnOut, int level)
int	cleanEmptyNodes ()
void	splitNode (CL2DClass *node, CL2DClass *&node1, CL2DClass *&node2, std::vector< size_t > &finalAssignment) const
void	splitFirstNode ()

Public Attributes
size_t	Nimgs =0
	Number of images. More...
MetaDataDb *	SF =nullptr
	Pointer to input metadata. More...
std::vector< CL2DClass * >	P
	List of nodes. More...

Detailed Description

Class for a CL2D

Definition at line 160 of file mpi_classify_CL2D.h.

Constructor & Destructor Documentation

CL2D() [1/3]

CL2D::CL2D ( )

inline

Empty constructor

Definition at line 172 of file mpi_classify_CL2D.h.

{}

~CL2D()

CL2D::~CL2D

(

)

Destructor

Definition at line 723 of file mpi_classify_CL2D.cpp.

{
    int qmax=P.size();
    for (int q=0; q<qmax; q++)
        delete P[q];
}

CL2D() [2/3]

CL2D::CL2D ( const CL2D &  )

delete

CL2D() [3/3]

CL2D::CL2D ( const CL2D &&  )

delete

Member Function Documentation

cleanEmptyNodes()

int CL2D::cleanEmptyNodes

(

)

Clean empty nodes. The number of nodes removed is returned.

Definition at line 1255 of file mpi_classify_CL2D.cpp.

{
    int retval = 0;
    std::vector<CL2DClass *>::iterator ptr = P.begin();
    while (ptr != P.end())
        if ((*ptr)->currentListImg.size() == 0)
        {
            ptr = P.erase(ptr);
            retval++;
        }
        else
            ptr++;
    return retval;
}

initialize()

void CL2D::initialize	(	MetaDataDb &	_SF,
		std::vector< MultidimArray< double > > &	_codes0
	)

Initialize.

Definition at line 748 of file mpi_classify_CL2D.cpp.

{
    if (prm->node->rank == 1)
        std::cout << "Initializing ...\n";

    SF = &_SF;
    Nimgs = SF->size();

    // Start with _Ncodes0 codevectors
    CL2DAssignment assignment;
    assignment.corr = 1;
    bool initialCodesGiven = _codes0.size() > 0;
    for (int q = 0; q < prm->Ncodes0; q++)
    {
        P.push_back(new CL2DClass());
        if (initialCodesGiven)
        {
            P[q]->Pupdate = _codes0[q];
            P[q]->nextListImg.push_back(assignment);
            P[q]->transferUpdate(false);
        }
    }

    // Estimate sigma and if no previous classes have been given,
    // assign randomly
    prm->sigma = 0;
    if (prm->node->rank == 1)
        init_progress_bar(Nimgs);
    Image<double> I;
    MultidimArray<double> Iaux, Ibest;
    bool oldUseCorrelation = prm->useCorrelation;
    prm->useCorrelation = true; // Since we cannot make the assignment before calculating sigma
    CL2DAssignment bestAssignment;
    size_t idx=0;
    SF->fillConstant(MDL_REF,"-1");
    for (size_t _ : prm->SF.ids())
    {
        if ((idx+1)%prm->node->size==prm->node->rank)
        {
            int q = -1;
            if (!initialCodesGiven)
                q = idx % (prm->Ncodes0);
            size_t objId = prm->objId[idx];
            readImage(I, objId, true);

            // Measure the variance of the signal outside a circular mask
            double avg, stddev;
            I().computeAvgStdev_within_binary_mask(prm->mask, avg, stddev);
            prm->sigma += stddev;

            // Put it randomly in one of the classes
            bestAssignment.objId = assignment.objId = objId;
            if (!initialCodesGiven)
            {
                bestAssignment.corr = 1;
                P[q]->updateProjection(I(), bestAssignment);
            }
            else
            {
                bestAssignment.corr = 0;
                q = -1;
                for (int qp = 0; qp < prm->Ncodes0; qp++)
                {
                    Iaux = I();
                    P[qp]->fit(Iaux, assignment);
                    if (assignment.corr > bestAssignment.corr)
                    {
                        bestAssignment = assignment;
                        Ibest = Iaux;
                        q = qp;
                    }
                }
                if (q != -1)
                    P[q]->updateProjection(Ibest, bestAssignment);
            }
            SF->setValue(MDL_REF, q + 1, objId);
            if (idx % 100 == 0 && prm->node->rank == 1)
                progress_bar(idx);
        }
        idx++;
    }
    prm->node->barrierWait();
    if (prm->node->rank == 1)
        progress_bar(Nimgs);
    prm->useCorrelation = oldUseCorrelation;

    // Put sigma in common
    MPI_Allreduce(MPI_IN_PLACE, &(prm->sigma), 1, MPI_DOUBLE, MPI_SUM,
                  MPI_COMM_WORLD);
    prm->sigma *= sqrt(2.0) / Nimgs;

    // Share all assignments
    shareAssignments(true, true, false);

    // Now compute the histogram of corr values
    if (!prm->classicalMultiref)
    {
        if (prm->node->rank == 1)
        {
            std::cout << "Computing histogram of correlation values\n";
            init_progress_bar(Nimgs);
        }

        CL2DAssignment inClass, outClass;
        size_t idx=0;
        for (size_t _ : prm->SF.ids())
        {
            if ((idx+1)%prm->node->size==prm->node->rank)
            {
                size_t objId = prm->objId[idx];
                readImage(I, objId, false);

                int q;
                SF->getValue(MDL_REF, q, objId);
                if (q == -1)
                    continue;
                q -= 1;

                outClass.objId = inClass.objId = objId;
                if (q != -1)
                {
                    Iaux = I();
                    P[q]->fit(Iaux, inClass);
                    P[q]->updateProjection(Iaux, inClass);
                    if (prm->Ncodes0 > 1)
                        for (int qp = 0; qp < prm->Ncodes0; qp++)
                        {
                            if (qp == q)
                                continue;
                            Iaux = I();
                            P[qp]->fit(Iaux, outClass);
                            P[qp]->updateNonProjection(outClass.corr);
                        }
                }
                if (prm->node->rank == 1 && idx % 100 == 0)
                    progress_bar(idx);
            }
            idx++;
        }
        prm->node->barrierWait();
        if (prm->node->rank == 1)
            progress_bar(Nimgs);

        shareAssignments(false, true, true);
    }
}

lookNode()

void CL2D::lookNode	(	MultidimArray< double > &	I,
		int	oldnode,
		int &	newnode,
		CL2DAssignment &	bestAssignment
	)

Look for a node suitable for this image. The image is rotationally and translationally aligned with the best node.

Definition at line 944 of file mpi_classify_CL2D.cpp.

{
#ifdef DEBUG
    std::cout << "Looking for node. Oldnode=" << oldnode << std::endl;
#endif
    int Q = P.size();
    int bestq = -1;
    MultidimArray<double> bestImg, Iaux;
    Matrix1D<double> corrList;
    corrList.resizeNoCopy(Q);
    CL2DAssignment assignment;
    bestAssignment.likelihood = bestAssignment.corr = 0;
    size_t objId = bestAssignment.objId;
    for (int q = 0; q < Q; q++)
    {
        // Check if q is neighbour of the oldnode
        bool proceed = false;
        if (oldnode >= 0)
        {
            if (oldnode<Q)
            {
                int imax = P[oldnode]->neighboursIdx.size();
                for (int i = 0; i < imax; i++)
                {
                    if (P[oldnode]->neighboursIdx[i] == q)
                    {
                        proceed = true;
                        break;
                    }
                }
                if (!proceed)
                {
                    double threshold = 3.0 * P[oldnode]->currentListImg.size();
                    threshold = XMIPP_MAX(threshold,1000);
                    threshold = (double) (XMIPP_MIN(threshold,Nimgs)) / Nimgs;
                    proceed = (rnd_unif(0, 1) < threshold);
                }
            }
            else
            {
                // In some strange conditions we may loose the reference to its neighbours
                proceed = true;
            }
        }
        else
            proceed = true;

        if (proceed) {
            // Try this image
            Iaux = I;
            P[q]->fit(Iaux, assignment);
            VEC_ELEM(corrList,q) = assignment.corr;
#ifdef DEBUG
    std::cout << "   Proceeding with node " << q << " corr=" << assignment.corr << std::endl;
#endif
            if ((!prm->classicalMultiref && assignment.likelihood > bestAssignment.likelihood) ||
                (prm->classicalMultiref && assignment.corr > bestAssignment.corr) ||
                 prm->classifyAllImages && bestAssignment.corr==0) {
                bestq = q;
                bestImg = Iaux;
                bestAssignment = assignment;
            }
        }
    }

    I = bestImg;
    newnode = bestq;
    bestAssignment.objId = objId;

    // Assign it to the new node and remove it from the rest
    // of nodes if it was among the best
#ifdef DEBUG
    std::cout << "   New node=" << newnode << std::endl;
#endif
    if (newnode != -1)
    {
        P[newnode]->updateProjection(I, bestAssignment);
        if (!prm->classicalMultiref)
            for (int q = 0; q < Q; q++)
                if (q != newnode && corrList(q) > 0)
                    P[q]->updateNonProjection(corrList(q));
    }
}

operator=() [1/2]

CL2D& CL2D::operator= ( const CL2D &  )

delete

operator=() [2/2]

CL2D& CL2D::operator= ( const CL2D &&  )

delete

readImage()

void CL2D::readImage	(	Image< double > &	I,
		size_t	objId,
		bool	applyGeo
	)		const

Read Image.

Definition at line 731 of file mpi_classify_CL2D.cpp.

{
    if (applyGeo)
        I.readApplyGeo(*SF, objId);
    else
    {
        FileName fnImg;
        SF->getValue(MDL_IMAGE, fnImg, objId);
        I.read(fnImg);
    }
    I().setXmippOrigin();
    if (prm->normalizeImages)
        I().statisticsAdjust(0., 1.);
}

run()

void CL2D::run	(	const FileName &	fnODir,
		const FileName &	fnOut,
		int	level
	)

Quantize with the current number of codevectors

Definition at line 1039 of file mpi_classify_CL2D.cpp.

{
    size_t Q = P.size();

    if (prm->node->rank == 1)
        std::cout << "Quantizing with " << Q << " codes...\n";
#ifdef DEBUG
    Image<double> save;
    for (int q=0; q<Q; q++)
    {
        save()=P[q]->P;
        save.write(formatString("PPPclass%04d.xmp",q));
    }
    std::cout << "Classes written. Press any key" << std::endl;
    char c; std::cin >> c;
#endif

    std::vector<int> oldAssignment, newAssignment;

    int iter = 1;
    bool goOn = true;
    MetaDataVec MDChanges;
    Image<double> I;
    int progressStep = XMIPP_MAX(1,Nimgs/60);
    CL2DAssignment assignment;
    FileName fnResultsDir=formatString("%s/level_%02d",fnODir.c_str(),level);
    fnResultsDir.makePath(0755);
    while (goOn)
    {
        if (prm->node->rank == 1)
        {
            std::cout << "Iteration " << iter << " ...\n";
            init_progress_bar(Nimgs);
        }

        size_t K = std::min((size_t)prm->Nneighbours+1,Q);
        if (K == 0)
            K = Q;
        for (size_t q = 0; q < Q; q++)
            P[q]->lookForNeighbours(P, K);

        int node;
        double corrSum = 0;
        SF->getColumnValues(MDL_REF, oldAssignment);
        int *ptrOld = &(oldAssignment[0]);
        for (size_t n = 0; n < Nimgs; ++n, ++ptrOld)
            *ptrOld -= 1;
        SF->fillConstant(MDL_REF, "-1");
        size_t idx=0;
        for (size_t _ : prm->SF.ids())
        {
            if ((idx+1)%prm->node->size==prm->node->rank)
            {
                size_t objId = prm->objId[idx];
                readImage(I, objId, false);
                LOG(((String)"Processing image: "+I.name()).c_str());

                assignment.objId = objId;
                lookNode(I(), oldAssignment[idx], node, assignment);
                SF->setValue(MDL_REF, node + 1, objId);
                corrSum += assignment.corr;
                if (prm->node->rank == 1 && idx % progressStep == 0)
                    progress_bar(idx);
            }
            idx++;
        }
        prm->node->barrierWait();

        // Gather all pieces computed by nodes
        MPI_Allreduce(MPI_IN_PLACE, &corrSum, 1, MPI_DOUBLE, MPI_SUM,
                      MPI_COMM_WORLD);
        shareAssignments(true, true, true);

        // Some report
        size_t idMdChanges=0;
        int finish=0;
        if (prm->node->rank == 1)
        {
            progress_bar(Nimgs);
            double avgSimilarity = corrSum / Nimgs;
            if (avgSimilarity==0)
            {
                std::cerr << "The average correlation is 0.0, make sure that the maximum allowed shift is enough\n";
                   
                //MPI_Abort(MPI_COMM_WORLD,ERR_UNCLASSIFIED);
                finish=1;
            }
            std::cout << "\nAverage correlation with input vectors=" << avgSimilarity << std::endl;
            idMdChanges = MDChanges.addObject();
            MDChanges.setValue(MDL_ITER, iter, idMdChanges);
            MDChanges.setValue(MDL_CL2D_SIMILARITY, avgSimilarity, idMdChanges);
        }
        MPI_Bcast(&finish,1, MPI_INT, 1, MPI_COMM_WORLD);
        if (finish)
        {
            MPI_Finalize();
            exit(ERR_UNCLASSIFIED);
        }


        // Count changes
        SF->getColumnValues(MDL_REF, newAssignment);
        int *ptrNew = &(newAssignment[0]);
        for (size_t n = 0; n < Nimgs; ++n, ++ptrNew)
            *ptrNew -= 1;
        size_t Nchanges = 0;
        if (iter > 1)
        {
            int *ptrNew = &(newAssignment[0]);
            ptrOld = &(oldAssignment[0]);
            for (size_t n = 0; n < Nimgs; ++n, ++ptrNew, ++ptrOld)
                if (*ptrNew != *ptrOld)
                    ++Nchanges;
        }
        if (prm->node->rank == 1)
        {
            FileName fnClasses=formatString("%s/%s_classes.xmd",fnResultsDir.c_str(),fnOut.c_str());
            if (iter==1)
            {
                MetaDataVec dummy;
                dummy.write(formatString("classes@%s",fnClasses.c_str()));
            }

            std::cout << "Number of assignment changes=" << Nchanges
            << std::endl;
            MDChanges.setValue(MDL_CL2D_CHANGES, (int)Nchanges, idMdChanges);
            MDChanges.write(formatString("info@%s",fnClasses.c_str()),MD_APPEND);
        }

        // Check if there are empty nodes
        if (Q > 1)
        {
            bool smallNodes;
            do
            {
                smallNodes = false;
                int largestNode = -1, sizeLargestNode = -1, smallNode = -1;
                size_t sizeSmallestNode = Nimgs + 1;
                for (size_t q = 0; q < Q; q++)
                {
                    if (P[q]->currentListImg.size() < sizeSmallestNode)
                    {
                        smallNode = q;
                        sizeSmallestNode = P[q]->currentListImg.size();
                    }
                    if ((int) (P[q]->currentListImg.size()) > sizeLargestNode)
                    {
                        sizeLargestNode = P[q]->currentListImg.size();
                        largestNode = q;
                    }
                }
                if (sizeLargestNode==0)
                    REPORT_ERROR(ERR_UNCLASSIFIED,"All classes are of size 0: normally this happens when images are too noisy");
                if (largestNode == -1 || smallNode == -1)
                    break;
                if (sizeSmallestNode < prm->PminSize * Nimgs / Q * 0.01 && sizeSmallestNode<0.25*sizeLargestNode)
                {
                    if (prm->node->rank == 1 && prm->verbose)
                        std::cout << "Splitting node " << largestNode << " (size=" << sizeLargestNode << ") "
                        << " by overwriting " << smallNode << " (size=" << sizeSmallestNode << ")" << std::endl;
                    smallNodes = true;

                    // Clear the old assignment of the images in the small node
                    std::vector<CL2DAssignment> &currentListImgSmall =
                        P[smallNode]->currentListImg;
                    size_t iimax = currentListImgSmall.size();
                    for (size_t ii = 0; ii < iimax; ii++)
                        SF->setValue(MDL_REF, -1, currentListImgSmall[ii].objId);
                    delete P[smallNode];

                    // Clear the old assignment of the images in the large node
                    std::vector<CL2DAssignment> &currentListImgLargest =
                        P[largestNode]->currentListImg;
                    iimax = currentListImgLargest.size();
                    for (size_t ii = 0; ii < iimax; ii++)
                        SF->setValue(MDL_REF, -1,
                                     currentListImgLargest[ii].objId);

                    // Now split the largest node
                    auto *node1 = new CL2DClass();
                    auto *node2 = new CL2DClass();
                    std::vector<size_t> splitAssignment;
                    splitNode(P[largestNode], node1, node2, splitAssignment);
                    delete P[largestNode];

                    // Keep the results of the split
                    P[largestNode] = node1;
                    P[smallNode] = node2;
                    iimax = splitAssignment.size();
                    for (size_t ii = 0; ii < iimax; ii += 2)
                    {
                        if (splitAssignment[ii + 1] == 1)
                            SF->setValue(MDL_REF, largestNode + 1,
                                         splitAssignment[ii]);
                        else
                            SF->setValue(MDL_REF, smallNode + 1,
                                         splitAssignment[ii]);
                    }
                }
            }
            while (smallNodes);
        }

        if (prm->node->rank == 1)
            write(fnODir,fnOut,level);

        if ((iter > 1 && Nchanges < 0.005 * Nimgs && Q > 1) || iter >= prm->Niter)
            goOn = false;
        iter++;
    }

    std::sort(P.begin(), P.end(), SDescendingClusterSort());
}

shareAssignments()

void CL2D::shareAssignments	(	bool	shareAssignment,
		bool	shareUpdates,
		bool	shareNonCorr
	)

Share assignments.

Definition at line 536 of file mpi_classify_CL2D.cpp.

{
    if (shareAssignment)
    {
        // Put assignment in common
        std::vector<int> nodeRef;
        if (SF->containsLabel(MDL_REF))
            SF->getColumnValues(MDL_REF, nodeRef);
        else
            nodeRef.resize(SF->size(), -1);

        MPI_Allreduce(MPI_IN_PLACE, &(nodeRef[0]), nodeRef.size(), MPI_INT,
                      MPI_MAX, MPI_COMM_WORLD);
        SF->setColumnValues(MDL_REF, nodeRef);
#ifdef DEBUG_WITH_LOG
        FileName fnImg;
        for (size_t objId : SF->id())
        {
            int classRef;
            SF->getValue(MDL_IMAGE,fnImg,objId);
            SF->getValue(MDL_REF,classRef,objId);
            LOG(formatString("Image %s: class: %d",fnImg.c_str(),classRef).c_str());
        }
#endif
    }

    // Share code updates
    if (shareUpdates)
    {
        std::vector<CL2DAssignment> auxList;
        std::vector<double> auxList2;
        int Q = P.size();
        for (int q = 0; q < Q; q++)
        {
            MPI_Allreduce(MPI_IN_PLACE, MULTIDIM_ARRAY(P[q]->Pupdate),
                          MULTIDIM_SIZE(P[q]->Pupdate), MPI_DOUBLE, MPI_SUM,
                          MPI_COMM_WORLD);

            // Share nextClassCorr and nextNonClassCorr
            std::vector<double> receivedNonClassCorr;
            std::vector<CL2DAssignment> receivedNextListImage;
            int listSize;
            for (size_t rank = 0; rank < prm->node->size; rank++)
            {
                if (rank == prm->node->rank)
                {
                    listSize = P[q]->nextListImg.size();
                    MPI_Bcast(&listSize, 1, MPI_INT, rank, MPI_COMM_WORLD);
                    MPI_Bcast(&(P[q]->nextListImg[0]),
                              P[q]->nextListImg.size() * sizeof(CL2DAssignment),
                              MPI_CHAR, rank, MPI_COMM_WORLD);
                    if (shareNonCorr)
                    {
                        listSize = P[q]->nextNonClassCorr.size();
                        MPI_Bcast(&listSize, 1, MPI_INT, rank, MPI_COMM_WORLD);
                        MPI_Bcast(&(P[q]->nextNonClassCorr[0]),
                                  P[q]->nextNonClassCorr.size(), MPI_DOUBLE,
                                  rank, MPI_COMM_WORLD);
                    }
                }
                else
                {
                    MPI_Bcast(&listSize, 1, MPI_INT, rank, MPI_COMM_WORLD);
                    auxList.resize(listSize);
                    MPI_Bcast(&(auxList[0]), listSize * sizeof(CL2DAssignment),
                              MPI_CHAR, rank, MPI_COMM_WORLD);
                    receivedNextListImage.reserve(
                        receivedNextListImage.size() + listSize);
                    for (int j = 0; j < listSize; j++)
                        receivedNextListImage.push_back(auxList[j]);
                    if (shareNonCorr)
                    {
                        MPI_Bcast(&listSize, 1, MPI_INT, rank, MPI_COMM_WORLD);
                        auxList2.resize(listSize);
                        MPI_Bcast(&(auxList2[0]), listSize, MPI_DOUBLE, rank,
                                  MPI_COMM_WORLD);
                        receivedNonClassCorr.reserve(
                            receivedNonClassCorr.size() + listSize);
                        for (int j = 0; j < listSize; j++)
                            receivedNonClassCorr.push_back(auxList2[j]);
                    }
                }
            }
            // This is important to ensure that all nodes have all images in the same order
            std::sort(receivedNextListImage.begin(),receivedNextListImage.end(),CL2DAssignmentComparator);

            // Copy the received elements
            listSize = receivedNextListImage.size();
            P[q]->nextListImg.reserve(P[q]->nextListImg.size() + listSize);
            for (int j = 0; j < listSize; j++)
                P[q]->nextListImg.push_back(receivedNextListImage[j]);
            listSize = receivedNonClassCorr.size();
            P[q]->nextNonClassCorr.reserve(
                P[q]->nextNonClassCorr.size() + listSize);
            for (int j = 0; j < listSize; j++)
                P[q]->nextNonClassCorr.push_back(receivedNonClassCorr[j]);
        }

        transferUpdates();
    }
}

shareSplitAssignments()

void CL2D::shareSplitAssignments	(	Matrix1D< int > &	assignment,
		CL2DClass *	node1,
		CL2DClass *	node2
	)		const

Share split assignment.

Definition at line 639 of file mpi_classify_CL2D.cpp.

{
    // Share assignment
    int imax = VEC_XSIZE(assignment);
    MPI_Allreduce(MPI_IN_PLACE, MATRIX1D_ARRAY(assignment), imax, MPI_INT,
                  MPI_MAX, MPI_COMM_WORLD);

    // Share code updates
    std::vector<CL2DAssignment> auxList;
    std::vector<double> auxList2;
    CL2DClass *node = NULL;
    for (int q = 0; q < 2; q++)
    {
        if (q == 0)
            node = node1;
        else
            node = node2;
        MPI_Allreduce(MPI_IN_PLACE, MULTIDIM_ARRAY(node->Pupdate),
                      MULTIDIM_SIZE(node->Pupdate), MPI_DOUBLE, MPI_SUM,
                      MPI_COMM_WORLD);

        // Share nextClassCorr and nextNonClassCorr
        std::vector<double> receivedNonClassCorr;
        std::vector<CL2DAssignment> receivedNextListImage;
        int listSize;
        for (size_t rank = 0; rank < prm->node->size; rank++)
        {
            if (rank == prm->node->rank)
            {
                // Transmit node 1 next list of images
                listSize = node->nextListImg.size();
                MPI_Bcast(&listSize, 1, MPI_INT, rank, MPI_COMM_WORLD);
                MPI_Bcast(&(node->nextListImg[0]),
                          node->nextListImg.size() * sizeof(CL2DAssignment),
                          MPI_CHAR, rank, MPI_COMM_WORLD);
                // Transmit node 1 non class corr
                listSize = node->nextNonClassCorr.size();
                MPI_Bcast(&listSize, 1, MPI_INT, rank, MPI_COMM_WORLD);
                MPI_Bcast(&(node->nextNonClassCorr[0]),
                          node->nextNonClassCorr.size(), MPI_DOUBLE, rank,
                          MPI_COMM_WORLD);
            }
            else
            {
                // Receive node 1 next list of images
                MPI_Bcast(&listSize, 1, MPI_INT, rank, MPI_COMM_WORLD);
                auxList.resize(listSize);
                MPI_Bcast(&(auxList[0]), listSize * sizeof(CL2DAssignment),
                          MPI_CHAR, rank, MPI_COMM_WORLD);
                receivedNextListImage.reserve(
                    receivedNextListImage.size() + listSize);
                for (int j = 0; j < listSize; j++)
                    receivedNextListImage.push_back(auxList[j]);
                // Receive node 1 non class corr
                MPI_Bcast(&listSize, 1, MPI_INT, rank, MPI_COMM_WORLD);
                auxList2.resize(listSize);
                MPI_Bcast(&(auxList2[0]), listSize, MPI_DOUBLE, rank,
                          MPI_COMM_WORLD);
                receivedNonClassCorr.reserve(
                    receivedNonClassCorr.size() + listSize);
                for (int j = 0; j < listSize; j++)
                    receivedNonClassCorr.push_back(auxList2[j]);
            }
        }
        // This is important to ensure that all nodes have all images in the same order
        std::sort(receivedNextListImage.begin(),receivedNextListImage.end(),CL2DAssignmentComparator);

        // Copy the received elements
        listSize = receivedNextListImage.size();
        node->nextListImg.reserve(node->nextListImg.size() + listSize);
        for (int j = 0; j < listSize; j++)
            node->nextListImg.push_back(receivedNextListImage[j]);
        listSize = receivedNonClassCorr.size();
        node->nextNonClassCorr.reserve(node->nextNonClassCorr.size() + listSize);
        for (int j = 0; j < listSize; j++)
            node->nextNonClassCorr.push_back(receivedNonClassCorr[j]);
    }

    node1->transferUpdate();
    node2->transferUpdate();
}

splitFirstNode()

void CL2D::splitFirstNode

(

)

Split the widest node

Definition at line 1633 of file mpi_classify_CL2D.cpp.

{
    std::sort(P.begin(), P.end(), SDescendingClusterSort());
    int Q = P.size();
    P.push_back(new CL2DClass());
    P.push_back(new CL2DClass());
    std::vector<size_t> splitAssignment;
    splitNode(P[0], P[Q], P[Q + 1], splitAssignment);
    delete P[0];
    P[0] = NULL;
    P.erase(P.begin());
}

splitNode()

void CL2D::splitNode	(	CL2DClass *	node,
		CL2DClass *&	node1,
		CL2DClass *&	node2,
		std::vector< size_t > &	finalAssignment
	)		const

Split node

Definition at line 1272 of file mpi_classify_CL2D.cpp.

{
    const CL2DClass* firstNode=node;
    std::vector<CL2DClass *> toDelete;
    Matrix1D<int> newAssignment, oldAssignment, firstSplitAssignment;
    Image<double> I;
    MultidimArray<double> Iaux1, Iaux2, corrList;
    MultidimArray<int> idx;
    CL2DAssignment assignment, assignment1, assignment2;
    CL2DClass *firstSplitNode1 = NULL;
    CL2DClass *firstSplitNode2 = NULL;
    size_t minAllowedSize = 0;

    bool oldclassicalMultiref = prm->classicalMultiref;
    prm->classicalMultiref = false;
    bool finish;
    bool success = true;
    int splitTrials=0;
    do
    {
        minAllowedSize = (size_t)(prm->PminSize/2 * 0.01 * node->currentListImg.size());

        // Sort the currentListImg to make sure that all nodes have the same order
        std::sort(node->currentListImg.begin(),node->currentListImg.end(),CL2DAssignmentComparator);
#ifdef DEBUG
            std::cout << "Splitting node " << node << "(" << node->currentListImg.size() << ") into " << node1 << " and " << node2 << std::endl;
            Image<double> save;
            save()=node->P;
            save.write("PPPnodeToSplit.xmp");
#endif
        finish = true;
        node2->neighboursIdx = node1->neighboursIdx = node->neighboursIdx;
        node2->P = node1->P = node->P;

        size_t imax = node->currentListImg.size();
        if (imax < minAllowedSize)
        {
            toDelete.push_back(node1);
            toDelete.push_back(node2);
#ifdef DEBUG
            std::cout << "Pushing to delete " << node1 << std::endl;
            std::cout << "Pushing to delete " << node2 << std::endl;
#endif
            success = false;
            break;
        }

        // Compute the corr histogram
        if (prm->node->rank == 1 && prm->verbose >= 2)
            std::cerr << "Calculating corr distribution at split ..."
            << std::endl;
        corrList.initZeros(imax);
        for (size_t i = 0; i < imax; i++)
        {
            if ((i + 1) % (prm->node->size) == prm->node->rank)
            {
                readImage(I, node->currentListImg[i].objId, false);
                node->fit(I(), assignment);
                A1D_ELEM(corrList,i) = assignment.corr;
            }
            if (prm->node->rank == 1 && i % 25 == 0 && prm->verbose >= 2)
                progress_bar(i);
        }
        if (prm->node->rank == 1 && prm->verbose >= 2)
            progress_bar(imax);
        MPI_Allreduce(MPI_IN_PLACE, MULTIDIM_ARRAY(corrList), imax, MPI_DOUBLE,
                      MPI_MAX, MPI_COMM_WORLD);
        newAssignment.initZeros(imax);

        // Compute threshold
        corrList.indexSort(idx);
        double corrThreshold = corrList(idx(XSIZE(idx)/2)-1);
#ifdef DEBUG
        corrList.printStats();
        std::cout << " corrThreshold= " << corrThreshold << std::endl;
#endif
        if (corrThreshold == 0)
        {
            if (firstSplitNode1 != NULL)
            {
                toDelete.push_back(node1);
                toDelete.push_back(node2);
#ifdef DEBUG
            std::cout << "Pushing to delete " << node1 << std::endl;
            std::cout << "Pushing to delete " << node2 << std::endl;
#endif
                success = false;
                break;
            }
            else
            {
                // Split at random
                for (size_t i = 0; i < imax; i++)
                {
                    assignment.objId = node->currentListImg[i].objId;
                    readImage(I, assignment.objId, false);
                    node->fit(I(), assignment);
                    if ((i + 1) % 2 == 0)
                    {
                        node1->updateProjection(I(), assignment,true);
                        VEC_ELEM(newAssignment,i) = 1;
                        node2->updateNonProjection(assignment.corr, true);
                    }
                    else
                    {
                        node2->updateProjection(I(), assignment,true);
                        VEC_ELEM(newAssignment,i) = 2;
                        node1->updateNonProjection(assignment.corr, true);
                    }
                }
                node1->transferUpdate();
                node2->transferUpdate();
#ifdef DEBUG
            std::cout << "Splitting at random " << node1->currentListImg.size() << " " << node2->currentListImg.size() << std::endl;
            save()=node1->P;
            save.write("PPPnode1.xmp");
            save()=node2->P;
            save.write("PPPnode2.xmp");
            std::cout << "Press any key" << std::endl;
            char c; std::cin >> c;
#endif
                success = true;
                break;
            }
        }
        else
        {
            // Split according to corr
            if (prm->node->rank == 1 && prm->verbose >= 2)
                std::cerr << "Splitting by corr threshold ..." << std::endl;
            for (size_t i = 0; i < imax; i++)
            {
                if ((i + 1) % (prm->node->size) == prm->node->rank)
                {
                    assignment.objId = node->currentListImg[i].objId;
                    readImage(I, assignment.objId, false);
                    node->fit(I(), assignment);
                    if (assignment.corr < corrThreshold)
                    {
                        node1->updateProjection(I(), assignment);
                        VEC_ELEM(newAssignment,i) = 1;
                        node2->updateNonProjection(assignment.corr);
                    }
                    else
                    {
                        node2->updateProjection(I(), assignment);
                        VEC_ELEM(newAssignment,i) = 2;
                        node1->updateNonProjection(assignment.corr);
                    }
                }
                if (prm->node->rank == 1 && i % 25 == 0 && prm->verbose >= 2)
                    progress_bar(i);
            }
            if (prm->node->rank == 1 && prm->verbose >= 2)
                progress_bar(imax);
            shareSplitAssignments(newAssignment, node1, node2);
        }

#ifdef DEBUG
    std::cout << "After first split Node1: " << node1->currentListImg.size() << " Node2 " << node2->currentListImg.size() << std::endl;
    save()=node1->P;
    save.write("PPPnode1.xmp");
    save()=node2->P;
    save.write("PPPnode2.xmp");
    std::cout << "Press any key" << std::endl;
    char c; std::cin >> c;
#endif

    // Backup the first split in case it fails
        if (firstSplitNode1 == NULL)
        {
#ifdef DEBUG
    std::cout << "Creating backup\n";
#endif
            firstSplitAssignment = newAssignment;
            firstSplitNode1 = new CL2DClass(*node1);
            firstSplitNode2 = new CL2DClass(*node2);
        }

        // Split iterations
        for (int it = 0; it < prm->Niter; it++)
        {
            if (prm->node->rank == 1 && prm->verbose >= 2)
            {
                init_progress_bar(imax);
            }

            oldAssignment = newAssignment;
            newAssignment.initZeros();
            for (size_t i = 0; i < imax; i++)
            {
                if ((i + 1) % (prm->node->size) == prm->node->rank)
                {
                    // Read image
                    assignment1.objId = assignment2.objId = assignment.objId
                                                            = node->currentListImg[i].objId;
                    readImage(I, assignment.objId, false);

                    Iaux1 = I();
                    node1->fit(Iaux1, assignment1);
                    Iaux2 = I();
                    node2->fit(Iaux2, assignment2);

                    //std::cout << "Image " << i << " Likelihood: " << assignment1.likelihood << " " << assignment2.likelihood << std::endl;
                    //std::cout << "Image " << i << " Corr: " << assignment1.corr << " " << assignment2.corr << std::endl;
                    if (assignment1.likelihood > assignment2.likelihood && !prm->classicalSplit)
                    {
                        node1->updateProjection(Iaux1, assignment1);
                        VEC_ELEM(newAssignment,i) = 1;
                        node2->updateNonProjection(assignment2.corr);
                        //std::cout << "Assigned to 1" << std::endl;
                    }
                    else if (assignment2.likelihood > assignment1.likelihood && !prm->classicalSplit)
                    {
                        node2->updateProjection(Iaux2, assignment2);
                        VEC_ELEM(newAssignment,i) = 2;
                        node1->updateNonProjection(assignment1.corr);
                        //std::cout << "Assigned to 2" << std::endl;
                    }
                    else if (prm->classifyAllImages || prm->classicalSplit)
                    {
                        if (assignment1.corr > assignment2.corr)
                        {
                            node1->updateProjection(Iaux1, assignment1);
                            VEC_ELEM(newAssignment,i) = 1;
                            node2->updateNonProjection(assignment2.corr);
                            //std::cout << "Assigned to 1" << std::endl;
                        }
                        else if (assignment2.corr > assignment1.corr)
                        {
                            node2->updateProjection(Iaux2, assignment2);
                            VEC_ELEM(newAssignment,i) = 2;
                            node1->updateNonProjection(assignment1.corr);
                            //std::cout << "Assigned to 2" << std::endl;
                        }
                    }
                }
                if (prm->node->rank == 1 && i % 25 == 0 && prm->verbose >= 2)
                    progress_bar(i);
            }
            if (prm->node->rank == 1 && prm->verbose >= 2)
                progress_bar(imax);
            shareSplitAssignments(newAssignment, node1, node2);

#ifdef DEBUG
    std::cout << "After refinement iteration: " << it << " Node1: " << node1->currentListImg.size() << " Node2 " << node2->currentListImg.size() << std::endl;
    save()=node1->P;
    save.write("PPPnode1.xmp");
    save()=node2->P;
    save.write("PPPnode2.xmp");
    std::cout << "Press any key" << std::endl;
    std::cin >> c;
#endif

            int Nchanges = 0;
            FOR_ALL_ELEMENTS_IN_MATRIX1D(newAssignment)
            if (newAssignment(i) != oldAssignment(i))
                Nchanges++;
            if (prm->node->rank == 1 && prm->verbose >= 2)
                std::cout << "Number of assignment split changes=" << Nchanges
                << std::endl;

            // Check if one of the nodes is too small
            if (node1->currentListImg.size() < minAllowedSize
                || node2->currentListImg.size() < minAllowedSize
                || Nchanges < 0.005 * imax)
                break;
        }

        if (node1->currentListImg.size() < minAllowedSize && node2->currentListImg.size() < minAllowedSize)
        {
            if (prm->node->rank == 1 && prm->verbose >= 2)
                std::cout << "Removing both nodes, they are too small. Current size: "
                << node1->currentListImg.size() << " Minimum allowed: "
                << minAllowedSize << "...\n";
            if (node1 != node)
                delete node1;
            if (node2 != node)
                delete node2;
            success = false;
            finish = true;
        }
        else if (node1->currentListImg.size() < minAllowedSize)
        {
            if (prm->node->rank == 1 && prm->verbose >= 2)
                std::cout << "Removing node1, it's too small "
                << node1->currentListImg.size() << " "
                << minAllowedSize << "...\n";
            if (node1 != node)
                delete node1;
            node1 = new CL2DClass();
            toDelete.push_back(node2);
#ifdef DEBUG
            std::cout << "Pushing to delete " << node2 << std::endl;
#endif
            node = node2;
            node2 = new CL2DClass();
            finish = false;
        }
        else if (node2->currentListImg.size() < minAllowedSize)
        {
            if (prm->node->rank == 1 && prm->verbose >= 2)
                std::cout << "Removing node2, it's too small "
                << node2->currentListImg.size() << " "
                << minAllowedSize << "...\n";
            if (node2 != node)
                delete node2;
            node2 = new CL2DClass();
            toDelete.push_back(node1);
#ifdef DEBUG
            std::cout << "Pushing to delete " << node1 << std::endl;
#endif
            node = node1;
            node1 = new CL2DClass();
            finish = false;
        }
        splitTrials++;
        if (splitTrials>=prm->NSplitTrials)
        {
            success = false;
            finish = true;  
        }
    }
    while (!finish);
    for (size_t i = 0; i < toDelete.size(); i++)
        if (toDelete[i] != firstNode)
        {
#ifdef DEBUG
            std::cout << "deleting from list " << toDelete[i] << std::endl;
#endif
            delete toDelete[i];
        }

    if (success)
    {
        for (size_t i = 0; i < node1->currentListImg.size(); i++)
        {
            splitAssignment.push_back(node1->currentListImg[i].objId);
            splitAssignment.push_back(1);
        }
        for (size_t i = 0; i < node2->currentListImg.size(); i++)
        {
            splitAssignment.push_back(node2->currentListImg[i].objId);
            splitAssignment.push_back(2);
        }
        delete firstSplitNode1;
        delete firstSplitNode2;
    }
    else
    {
        node1 = firstSplitNode1;
        node2 = firstSplitNode2;
        splitAssignment.reserve(VEC_XSIZE(firstSplitAssignment));
        FOR_ALL_ELEMENTS_IN_MATRIX1D(firstSplitAssignment)
        splitAssignment.push_back(VEC_ELEM(firstSplitAssignment,i));
    }
    prm->classicalMultiref = oldclassicalMultiref;
}

transferUpdates()

void CL2D::transferUpdates

(

)

Transfer all updates

Definition at line 1030 of file mpi_classify_CL2D.cpp.

{
    int Q = P.size();
    for (int q = 0; q < Q; q++)
        P[q]->transferUpdate();
}

write()

void CL2D::write	(	const FileName &	fnODir,
		const FileName &	fnRoot,
		int	level
	)		const

Write the nodes.

Definition at line 898 of file mpi_classify_CL2D.cpp.

{
    int Q = P.size();
    MetaDataVec SFout;
    Image<double> I;
    FileName fnResultsDir=formatString("%s/level_%02d",fnODir.c_str(),level);
    FileName fnOut = formatString("%s/%s_classes.stk",fnResultsDir.c_str(),fnRoot.c_str()), fnClass;
    fnOut.deleteFile();
    for (int q = 0; q < Q; q++)
    {
        fnClass.compose(q + 1, fnOut);
        I() = P[q]->P;
        I.write(fnClass, q, true, WRITE_APPEND);
        size_t id = SFout.addObject();
        SFout.setValue(MDL_REF, q + 1, id);
        SFout.setValue(MDL_IMAGE, fnClass, id);
        SFout.setValue(MDL_CLASS_COUNT,P[q]->currentListImg.size(), id);
    }
    FileName fnSFout = formatString("%s/%s_classes.xmd",fnResultsDir.c_str(),fnRoot.c_str());
    SFout.write("classes@"+fnSFout, MD_APPEND);

    // Make the selfiles of each class
    FileName fnImg;

    for (int q = 0; q < Q; q++)
    {
        MetaDataVec SFq;
        std::vector<CL2DAssignment> &currentListImg = P[q]->currentListImg;
        int imax = currentListImg.size();
        for (int i = 0; i < imax; i++)
        {
            const CL2DAssignment &assignment = currentListImg[i];
            MDRowSql row = SF->getRowSql(assignment.objId);
            row.setValue(MDL_FLIP, assignment.flip);
            row.setValue(MDL_SHIFT_X, assignment.shiftx);
            row.setValue(MDL_SHIFT_Y, assignment.shifty);
            row.setValue(MDL_ANGLE_PSI, assignment.psi);
            SFq.addRow(row);
        }
        MetaDataVec SFq_sorted;
        SFq_sorted.sort(SFq, MDL_IMAGE);
        SFq_sorted.write(formatString("class%06d_images@%s",q+1,fnSFout.c_str()),MD_APPEND);
    }
}

Member Data Documentation

Nimgs

size_t CL2D::Nimgs =0

Number of images.

Definition at line 163 of file mpi_classify_CL2D.h.

P

std::vector<CL2DClass *> CL2D::P

List of nodes.

Definition at line 169 of file mpi_classify_CL2D.h.

SF

MetaDataDb* CL2D::SF =nullptr

Pointer to input metadata.

Definition at line 166 of file mpi_classify_CL2D.h.

---

The documentation for this class was generated from the following files:

• xmipp/libraries/parallel/mpi_classify_CL2D.h
• xmipp/libraries/parallel/mpi_classify_CL2D.cpp