SpecImage.hpp

Go to the documentation of this file.

/*
        NAG Copyright 2016.
*/
#ifndef SPECIMAGE_HPP
#define SPECIMAGE_HPP
#include <string>
#include <memory>
#include <vector>


namespace nag
{

        typedef enum source
        {
                SIMS, 
                MALDI, 
                none 
        } source_t;

        typedef enum pca_mat
        {
                Covariance, 
                Correlation 
        } pca_mat_t;
        
        typedef enum kmeans_metric
        {
                Euclidean, 
                Cosine 
        } kmeans_metric_t;

        typedef enum matrix_ordering
        {
                PixelsByChannels, 
                ChannelsByPixels 
        } matrix_ordering_t;

        typedef enum wavelet_name
        {
                MexicanHat, 
                Haar 
        } wavelet_t;

        typedef enum wavelet_extension
        {
                Periodic, 
                Zero, 
                Symmetric 
        } wavelet_extension_t;

        template<class FPType> struct csr_mat
        {
                size_t nnz; 
                size_t m; 
                size_t n; 
                std::vector<size_t> row_ptr;
                std::vector<uint32_t> col_indx; 
                std::vector<FPType> values; 
        };
        template<class FPType> using csr_mat_t = struct csr_mat<FPType>;

        struct cfree
        {
                void operator()(void *ptr) { free(ptr); }
        };

        template<class FPType> struct pca_results
        {
                int k; 
                size_t m; 
                size_t n; 
                FPType pca_tot; 
                std::unique_ptr<FPType[], cfree> loadings; 
                std::unique_ptr<FPType[], cfree> scores; 
                std::vector<FPType> eigenvalues; 

                // We won't get a default copy constructor since we can't copy unique_ptrs, so define a default move constructor
                // - we seem to get away without this on Linux but not Windows
                pca_results(pca_results&& a) : loadings(std::move(a.loadings)), scores(std::move(a.scores)), eigenvalues(std::move(a.eigenvalues)) { k = a.k; m = a.m; n = a.n; pca_tot = a.pca_tot; };
                // Since we've defined a constructor the implicit default constructor is deleted, so explicitly request one
                pca_results() = default;
                pca_results& operator=(pca_results&& a) { loadings = std::move(a.loadings); scores = std::move(a.scores); eigenvalues = std::move(a.eigenvalues); k = a.k; m = a.m; n = a.n; pca_tot = a.pca_tot; return *this; };
                // note we don't have a copy assignment operator since we're working with unique_ptrs.
        };
        template<class FPType> using pca_results_t = struct pca_results<FPType>;

        template<class FPType> struct tsne_results
        {
                int k; 
                size_t m; 
                size_t d; 
                FPType kl_divergence; 
                FPType theta; 
                FPType perplexity; 
                std::unique_ptr<FPType[], cfree> embedding; 
                int total_iter; 
                
                // We won't get a default copy constructor since we can't copy unique_ptrs, so define a default move constructor
                // - we seem to get away without this on Linux but not Windows
                tsne_results(tsne_results&& a) : embedding(std::move(a.embedding) ) { k = a.k; m = a.m; d = a.d; kl_divergence = a.kl_divergence; theta = a.theta; perplexity = a.perplexity; total_iter = a.total_iter;};
                // Since we've defined a constructor the implicit default constructor is deleted, so explicitly request one
                tsne_results() = default;
                tsne_results& operator=(tsne_results&& a) { embedding = std::move(a.embedding); k = a.k; m = a.m; kl_divergence = a.kl_divergence; theta = a.theta; perplexity = a.perplexity; total_iter = a.total_iter; return *this; };
                // note we don't have a copy assignment operator since we're working with unique_ptrs.
        };
        template<class FPType> using tsne_results_t = struct tsne_results<FPType>;

        template<class FPType> struct kmeans_results
        {
                int k; 
                size_t m; 
                size_t n;
                kmeans_metric_t metric; 
                std::unique_ptr<FPType[], cfree> cmeans; 
                std::vector<int> inc; 
                std::vector<int> nic; 
                std::vector<FPType> wcs; 

                kmeans_results(kmeans_results&& a) : cmeans(std::move(a.cmeans)), inc(std::move(a.inc)), nic(std::move(a.nic)), wcs(std::move(a.wcs)) { k = a.k; m = a.m; n = a.n; metric = a.metric; };
                kmeans_results() = default;
                kmeans_results& operator=(kmeans_results&& a) { cmeans = std::move(a.cmeans); inc = std::move(a.inc); nic = std::move(a.nic); wcs = std::move(a.wcs); k = a.k; m = a.m; n = a.n; metric = a.metric; return *this; };
        };
        template<class FPType> using kmeans_results_t = struct kmeans_results<FPType>;

        template<class FPType> struct nmf_results
        {
                int k; 
                size_t m; 
                size_t n; 
                std::unique_ptr<FPType[], cfree> W; 
                std::unique_ptr<FPType[], cfree> H; 
                FPType tol; 

                nmf_results(nmf_results&& a) : W(std::move(a.W)), H(std::move(a.H)) { k = a.k; m = a.m; n = a.n; tol = a.tol; };
                nmf_results() = default;
                nmf_results& operator=(nmf_results&& a) { W = std::move(a.W); H = std::move(a.H); k = a.k; m = a.m; n = a.n; tol = a.tol; return *this; };
        };
        template<class FPType> using nmf_results_t = struct nmf_results<FPType>;

        template<class FPType> struct cwt_results
        {
                size_t n; 
                int nscales; 
                size_t npixels; 
                std::vector<uint32_t> scales; 
                std::vector<size_t> pixels; 
                std::unique_ptr<FPType[], cfree> coefficients; 

                cwt_results(cwt_results&& a) : scales(std::move(a.scales)), pixels(std::move(a.pixels)), coefficients(std::move(a.coefficients)) { n = a.n; nscales = a.nscales; npixels = a.npixels; }
                cwt_results() = default;
                cwt_results& operator=(cwt_results&& a) { scales = std::move(a.scales); pixels = std::move(a.pixels); coefficients = std::move(a.coefficients); n = a.n; nscales = a.nscales; npixels = a.npixels; return *this; }
        };
        template<class FPType> using cwt_results_t = struct cwt_results<FPType>;

        class Error
        {
                public:
                        std::string message; 
                        Error(std::string msg) : message(msg) {};
        };

        
        template<class FPType> class SpecImage
        {

                protected:
                        void *_imgh; 
                        source_t _source_type; 
        
                public:
                        SpecImage() : _imgh(NULL) {};

                        SpecImage(const csr_mat_t<FPType>& matrix, const FPType *mzvals = NULL);

                        SpecImage(size_t m, size_t n, const FPType *matrix, const FPType *mzvals = NULL);

                        virtual ~SpecImage();
        

                        /*
                                Preprocessing functions
                        */
                        void peakSelect();


                        /*
                                Analysis functions
                        */
                        pca_results_t<FPType> pca(uint32_t k, pca_mat_t t = Covariance);

                        kmeans_results_t<FPType> kmeans(uint32_t k, FPType* initial_cmeans = NULL, kmeans_metric_t metric = Cosine);

                        nmf_results_t<FPType> nmf(uint32_t k, FPType *initial_W = NULL, FPType *initial_H = NULL, 
                                                        FPType tol = 1.e-2, uint32_t max_iter = 200);

                        cwt_results_t<FPType> cwt(const std::vector<uint32_t>& scales, const std::vector<size_t>& pixels, 
                                                wavelet_t wav = Haar, wavelet_extension_t ext = Periodic);

                        tsne_results_t<FPType> tsne(uint32_t k, pca_results_t<FPType> *data = NULL, FPType *embedding = NULL,
                                                    FPType perplexity = 50, FPType theta = 0.3, uint32_t max_iter = 1000,
                                                    uint32_t verbosity = 50, FPType min_grad_norm = 1.0e-5,
                                                    uint32_t n_unimproved_iter = 30, FPType learning_rate = 100,
                                                    FPType early_exaggeration = 4.0, uint32_t exaggerated_iter = 50);
        

                        /*
                                Accessor functions. No set counterparts; setting 
                                should be done by constructing a new object.

                                Limits can be set only since they are provided by the user and stored internally
                                as matrix columns or rows, depending on the matrix ordering.
                        */
                        /*
                                Get matrix functions.
        
                                Perform conversion if matrix is currently in different format, or requires transpose.
                                Sparse matrix is row major (so e.g. pixelsByChannels stores channels contiguously).
                        */
                        csr_mat_t<FPType> getSparseMatrix(matrix_ordering_t order = PixelsByChannels);

                        std::unique_ptr<FPType[]> getDenseMatrix(matrix_ordering_t order = PixelsByChannels);


                        std::vector<FPType> getSinglePeak(FPType minmz, FPType maxmz);

        
                        size_t getM();

                        size_t getN();

                        size_t getNNZ();
        
                        int getXDim();

                        int getYDim();

                        int getZDim();

                        // rely on compiler elision / move semantics here
                        std::vector<FPType> getMZVals();
        
                        void setLimits(const std::vector<FPType>& limits);

                        source_t getSource() { return _source_type; };
        
                        //TODO?
                        // Functions to add spectra incrementally? Or remove spectra?
                        // Matrix operators: +, *, '(transpose) ?
        };

}
#endif