ml2d.h

Go to the documentation of this file.

/***************************************************************************
 *
 * Authors:    Sjors Scheres           scheres@cnb.csic.es (2004)
 *
 * Unidad de  Bioinformatica of Centro Nacional de Biotecnologia , CSIC
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
 * 02111-1307  USA
 *
 *  All comments concerning this program package may be sent to the
 *  e-mail address 'xmipp@cnb.csic.es'
 ***************************************************************************/

#ifndef _ML2D_H
#define _ML2D_H

#include <sys/time.h>
#include <pthread.h>
#include <random>
#include <vector>
#include <core/xmipp_fftw.h>
#include <core/xmipp_fft.h>
#include <core/args.h>
#include <core/xmipp_funcs.h>
#include <core/metadata_extension.h>
#include <core/xmipp_image.h>
#include <core/geometry.h>
#include <data/filters.h>
#include <data/mask.h>
#include <data/ctf.h>
#include <core/xmipp_threads.h>
#include <core/xmipp_program.h>
#include <core/xmipp_log.h>


#define SIGNIFICANT_WEIGHT_LOW 1e-8
#define SMALLANGLE 1.75
#define DATALINELENGTH 12

#define FN_EXTRA(file) formatString("%sextra/%s", fn_root.c_str(), file)
#define FN_IMAGES_MD(base)   (base) + "images.xmd"
#define FN_CLASSES_MD(base)  formatString("classes@%sclasses.xmd", (base).c_str())
#define FN_LOGS_MD(base)  formatString("info@%sclasses.xmd", (base).c_str())
#define FN_CLASS_IMAGES_MD(base, ref)  formatString("class%06d_images@%sclasses.xmd", (ref), (base).c_str())
#define FN_CLASSES_STK(base) (base) + "classes.stk"

#define FOR_ALL_GLOBAL_IMAGES() \
    for (size_t imgno = 0; imgno < nr_images_global; ++imgno)

#define FOR_ALL_LOCAL_IMAGES() \
    for (size_t imgno = myFirstImg; imgno <= myLastImg; ++imgno)

#define IMG_LOCAL_INDEX (imgno - myFirstImg)
#define IMG_REAL_INDEX(imgno) (img_id[(imgno)])
#define IMG_BLOCK(imgno) (imgno) % blocks

//For MPI
#define IS_MASTER (rank == 0)
//threads tasks
typedef enum { TH_EXIT, TH_ESI_REFNO, TH_ESI_UPDATE_REFNO, TH_RR_REFNO, TH_RRR_REFNO, TH_PFS_REFNO } ThreadTask;
//output types constants
typedef enum { OUT_BLOCK, OUT_ITER, OUT_FINAL, OUT_REFS, OUT_IMGS } OutputType;


class ModelML2D
{
public:
    int n_ref;
    std::vector < Image<double> > Iref, WsumMref;
    double sigma_noise, wsum_sigma_noise;
    double sigma_offset, wsum_sigma_offset;
    std::vector<double> sumwsc, sumw_mirror;
    std::vector<double> alpha_k, mirror_fraction, scale;
    double sumw_allrefs, sumw_allrefs2;
    double avePmax, sumfracweight;
    double LL;

    int dim;
    bool do_student, do_student_sigma_trick, do_norm;

    ModelML2D();
    ModelML2D(int n_ref);

    void initData();
    void setNRef(int n_ref);
    void combineModel(const ModelML2D &model, int sign);
    void addModel(const ModelML2D &model);
    void substractModel(const ModelML2D &model);

    void update();
    void print(int tabs = 0) const;
}
;//close class ModelML2D




class ML2DBaseProgram: public XmippProgram
{
public:
    FileName fn_img, fn_ref, fn_root, fn_frac, fn_sig, fn_doc, fn_oext, fn_scratch, fn_control;
    String cline;
    //double sigma_noise2;
    bool do_mirror;
    bool fix_fractions;
    bool fix_sigma_offset;
    bool fix_sigma_noise;
    int istart;
    int iter;
    int Niter;
    size_t dim, dim2, hdim;
    double ddim2;
    size_t nr_psi;
    size_t nr_flip;
    double psi_step;
    double psi_max;
    size_t nr_nomirror_flips;
    size_t nr_images_global;
    size_t nr_images_local;
    size_t myFirstImg, myLastImg;
    double eps;
    MetaDataDb MDimg, MDref, MDlog;
    std::vector<Matrix2D<double> > F;
    std::vector < Image<double> > Iold;
    MultidimArray<double> P_phi, Mr2;
    bool fast_mode;
    double C_fast;
    std::vector<Matrix1D<double> > Vtrans;
    bool zero_offsets;
    bool limit_rot;
    double search_rot;
    bool save_mem1, save_mem2, save_mem3;
    std::vector<double> imgs_oldphi, imgs_oldtheta;
    bool do_ML3D;
    bool do_generate_refs;
    std::vector<std::vector<double> > imgs_offsets;
    double trymindiff_factor;
    double df, df2, dfsigma2;

    bool do_norm;
    std::vector<double> imgs_scale, imgs_bgmean, imgs_trymindiff;
    std::vector<int> imgs_optrefno;
    double average_scale;
    int factor_nref;
    int refs_per_class;
    std::vector<double> conv;
    //Just for debbuging
    size_t current_image;

    ModelML2D model, *current_model;
    //Some incremental stuff
    //Number of blocks for IEM
    size_t blocks;
    //Current processing block
    size_t current_block;

    //Vector of image IDs in the MetaData object (change order for randomize)
    std::vector<size_t> img_id;
    MultidimArray<double> docfiledata;
    bool do_restart;

    bool referenceExclusive, allowFastOption, allowThreads, allowIEM, allowRestart;
    String defaultRoot;
    int defaultNiter, defaultStartingIter;

    int threads;

    String outRefsMd;

    MultidimArray<double>  spectral_signal;

    void initSamplingStuff();

public:
    /* Constructor for some initializations */
    ML2DBaseProgram();
    virtual void produceSideInfo() = 0;
    virtual void produceSideInfo2() = 0;

    /* This function will perform an ML iteration
     * usually composed by 'expectation' and 'maximization'
     */
    virtual void iteration() = 0;

    virtual void expectation() = 0;

    virtual void endIteration();

    virtual void maximization() = 0;

    virtual void setNumberOfLocalImages();

    virtual void randomizeImagesOrder();

    virtual void createThreads()
    {}
    ;//by default do nothing

    virtual void destroyThreads()
    {}
    ;//by default do nothing

    virtual void writeOutputFiles(const ModelML2D &model, OutputType outputType = OUT_FINAL) = 0;


    virtual bool checkConvergence();

    virtual void addPartialDocfileData(const MultidimArray<double> &data, size_t first, size_t last) = 0;
    virtual void run();

    virtual void defineBasicParams(XmippProgram * prog);
    virtual void defineAdditionalParams(XmippProgram * prog, const char * sectionLine);
    virtual void defineHiddenParams(XmippProgram *prog);

    int seed;
}
;//end of class ML2DBaseProgram

FileName getIterExtraPath(const FileName &fn_root, int iter, bool makePath=true);


#endif