@@ -38,8 +38,7 @@ checkpointOutFile="gaps_checkpoint.out", checkpointInterval=1000,

 checkpointInFile=NULL, transposeData=FALSE, ...)

 {

     # parse parameters from ...

-    # TODO remove args from ... as they're processed, check if null

-    fullParams <- list("gaps"=params,

+    allParams <- list("gaps"=params,

         "nCores"=nCores,

         "messages"=messages,

         "outputFrequency"=outputFrequency,

@@ -47,10 +46,9 @@ checkpointInFile=NULL, transposeData=FALSE, ...)

         "checkpointInterval"=checkpointInterval,

         "checkpointInFile"=checkpointInFile,

         "transposeData"=transposeData,

-        "whichMatrixFixed"="N" # internal parameter

+        "whichMatrixFixed"="" # internal parameter

     )

-    fullParams <- parseOldParams(fullParams, list(...))

-    fullParams <- parseDirectParams(fullParams, list(...))

+    allParams <- parseExtraParams(allParams, list(...))

     # check file extension

     if (class(data) == "character" & !(file_ext(data) %in% c("tsv", "csv", "mtx")))

@@ -97,6 +95,10 @@ checkpointInFile=NULL, transposeData=FALSE, ...)

         else

             message("Running CoGAPS from a checkpoint")

     }

+    else

+    {

+        allParams$checkpointInFile <- ""

+    }

     # run cogaps

     gapsReturnList <- dispatchFunc(data, allParams, uncertainty)

@@ -34,7 +34,7 @@ setMethod("initialize", "CogapsParams",

         .Object@maxGibbsMassP <- 100

         .Object@seed <- getMilliseconds(as.POSIXlt(Sys.time()))

         .Object@singleCell <- FALSE

-        .Object@distributed <- "A"

+        .Object@distributed <- ""

         .Object@nSets <- 3

         .Object@cut <- .Object@nPatterns

         .Object@minNS <- ceiling(.Object@nSets / 2)

@@ -23,7 +23,7 @@ distributedCogaps <- function(data, allParams, uncertainty)

     # set fixed matrix

     initialA <- initialP <- matrix(0)

-    if (allParams$modelParams@distributed == "genome-wide")

+    if (allParams$gaps@distributed == "genome-wide")

     {

         allParams$whichMatrixFixed <- "P"

         initialP <- consensusMatrix

@@ -122,42 +122,28 @@ patternMatch <- function(allPatterns, allParams)

     cc <- corcut(allPatterns, allParams)

     ### split by maxNS

-    indx <- which(unlist(lapply(cc$PatsByClust, function(x) nrow(x) > maxNS)))

-    i <- 1

+    indx <- which(sapply(cc$PatsByClust, function(x) ncol(x) > allParams$gaps@maxNS))

     while (length(indx) > 0)

     { 

-        icc <- corcut(cc$AByClust[[indx[1]]], minNS, 2, cluster.method)

-        if (length(icc[[2]]) == 0)

-        {

-            indx <- indx[-1]

-            next

-        }

-        else

-        {

-            cc$AByClust[[indx[1]]] <- icc[[2]][[1]]

-            cc$RtoMeanPattern[[indx[1]]] <- icc[[1]][[1]]

-            if (length(icc[[2]]) > 1)

-            {

-                cc$AByClust<-append(cc$AByClust,icc[[2]][2])

-                cc$RtoMeanPattern<-append(cc$RtoMeanPattern,icc[[1]][2])

-            } 

-            indx <- which(unlist(lapply(cc$AByClust,function(x) nrow(x) > maxNS)))

-        }

+        allParams$gaps@cut <- 2

+        icc <- corcut(cc$PatsByClust[[indx[1]]], allParams)

+

+        cc$PatsByClust[[indx[1]]] <- icc$PatsByClust[[1]]

+        cc$RtoMeanPattern[[indx[1]]] <- icc$RtoMeanPattern[[1]]

+        cc$PatsByClust <- append(cc$PatsByClust, icc$PatsByClust[2])

+        cc$RtoMeanPattern <- append(cc$RtoMeanPattern, icc$RtoMeanPattern[2])

+

+        indx <- which(sapply(cc$PatsByClust, function(x) ncol(x) > allParams$gaps@maxNS))

     }

-    #weighted.mean(AByClustDrop[[1]],RtoMPDrop[[1]])

-    AByCDSWavg <- t(sapply(1:length(cc$AByClust), function(z)

-        apply(cc$AByClust[[z]], 1, function(x)

-            weighted.mean(x, (cc$RtoMeanPattern[[z]])^3))))

-    rownames(AByCDSWavg) <- lapply(1:length(cc$AByClust), function(x)

-        paste("Pattern",x))

-    #scale As

-    Amax <- apply(AByCDSWavg, 1, max)

-    AByCDSWavgScaled <- t(sapply(1:dim(AByCDSWavg)[1], function(x)

-        AByCDSWavg[x,] / Amax[x]))

-    rownames(AByCDSWavgScaled) <- rownames(AByCDSWavg)

-

-    return(AByCDSWavgScaled)

+    # created matrix of mean patterns - weighted by R closeness of fit

+    PatsByCDSWavg <- t(sapply(1:length(cc$PatsByClust), function(z)

+        apply(cc$PatsByClust[[z]], 1, function(x) weighted.mean(x, (cc$RtoMeanPattern[[z]])^3))))

+    rownames(PatsByCDSWavg) <- lapply(1:length(cc$PatsByClust), function(x) paste("Pattern", x))

+    return(PatsByCDSWavg)

+

+    # scale

+    return(apply(PatsByCDSWavg, 1, function(row) row / max(row)))

 }

 corcut <- function(allPatterns, allParams)

@@ -5,6 +5,16 @@

 #include <string>

+static std::vector<unsigned> convertRVec(const Rcpp::NumericVector &rvec)

+{

+    std::vector<unsigned> vec;

+    for (unsigned i = 0; i < rvec.size(); ++i)

+    {

+        vec.push_back(rvec[i]);

+    }

+    return vec;

+}

+

 static RowMatrix convertRMatrix(const Rcpp::NumericMatrix &rmat)

 {

     RowMatrix mat(rmat.nrow(), rmat.ncol());

@@ -37,55 +47,86 @@ static bool isNull(const RowMatrix &mat)

     return mat.nRow() == 1 && mat.nCol() == 1;

 }

+static bool isNull(const Rcpp::NumericMatrix &mat)

+{

+    return mat.nrow() == 1 && mat.ncol();

+}

+

+static bool isNull(const Rcpp::NumericVector &vec)

+{

+    return vec.size() == 1;

+}

+

 static bool isNull(const std::string &path)

 {

     return path.empty();

 }

 template <class T>

-static Rcpp::List cogapsRun(const T &data, Rcpp::S4 params, const T &unc,

-const RowMatrix &fixedMatrix, const std::string &checkpointInFile)

+static Rcpp::List cogapsRun(const T &data, const Rcpp::List &allParams,

+const T &uncertainty, const Rcpp::NumericVector &indices,

+const Rcpp::NumericMatrix &initialA, const Rcpp::NumericMatrix &initialP)

 {

     GapsDispatcher dispatcher;

     // check if we're initializing with a checkpoint or not

-    if (!isNull(checkpointInFile))

+    if (!isNull(allParams["checkpointInFile"]))

     {

-        dispatcher.initialize(data, checkpointInFile);

+        dispatcher.initialize(data, allParams["transposeData"],

+            allParams["checkpointInFile"]);

     }

     else

     {

-        dispatcher.initialize(data, params.slot("nPatterns"), params.slot("seed"));

+        // check to see if we're subsetting the data

+        if (!isNull(allParams["gaps"].slot("distributed")))

+        {

+            dispatcher.initialize(data, allParams["transposeData"],

+                allParams["gaps"].slot("distributed") == "genome-wide",

+                convertRVec(indices), allParams["gaps"].slot("nPatterns"),

+                allParams["gaps"].slot("seed"));

+        }

+        else

+        {

+            dispatcher.initialize(data, allParams["transposeData"],

+                allParams["gaps"].slot("nPatterns"), allParams["gaps"].slot("seed"));

+        }

         // set optional parameters

-        dispatcher.setMaxIterations(params.slot("nIterations"));

-        dispatcher.setOutputFrequency(params.slot("outputFrequency"));

-

-        dispatcher.setSparsity(params.slot("alphaA"), params.slot("alphaP"),

-            params.slot("singleCell"));

-        

-        dispatcher.setMaxGibbsMass(params.slot("maxGibbsMassA"),

-            params.slot("maxGibbsMassP"));

-

-        dispatcher.printMessages(params.slot("messages"));

+        dispatcher.setMaxIterations(allParams["gaps"].slot("nIterations"));

+        dispatcher.setSparsity(allParams["gaps"].slot("alphaA"),

+            allParams["gaps"].slot("alphaP"), allParams["gaps"].slot("singleCell"));

+        dispatcher.setMaxGibbsMass(allParams["gaps"].slot("maxGibbsMassA"),

+            allParams["gaps"].slot("maxGibbsMassP"));

+

+        // set initial values for A and P matrix

+        if (!isNull(initialA))

+        {

+            dispatcher.setAMatrix(convertRMatrix(initialA));

+        }

+        if (!isNull(initialP))

+        {

+            dispatcher.setPMatrix(convertRMatrix(initialP));

+        }

         // check if running with a fixed matrix

-        if (!isNull(fixedMatrix))

+        if (!isNull(allParams["whichMatrixFixed"]))

         {

-            dispatcher.setFixedMatrix(params.slot("whichMatrixFixed"), fixedMatrix);

+            dispatcher.setFixedMatrix(allParams["whichMatrixFixed"]);

         }

     }

     // set the uncertainty matrix

-    if (!isNull(unc))

+    if (!isNull(uncertainty))

     {

-        dispatcher.setUncertainty(unc);

+        dispatcher.setUncertainty(uncertainty);

     }

     // set parameters that aren't saved in the checkpoint

     dispatcher.setNumCoresPerSet(params.slot("nCores"));

-    dispatcher.setCheckpointInterval(params.slot("checkpointInterval"));

+    dispatcher.printMessages(params.slot("messages"));

+    dispatcher.setOutputFrequency(params.slot("outputFrequency"));

     dispatcher.setCheckpointOutFile(params.slot("checkpointOutFile"));

+    dispatcher.setCheckpointInterval(params.slot("checkpointInterval"));

     // run the dispatcher and return the GapsResult in an R list

     GapsResult result(dispatcher.run());

@@ -105,24 +146,23 @@ const RowMatrix &fixedMatrix, const std::string &checkpointInFile)

 Rcpp::List cogaps_cpp_from_file(const std::string &data,

 const Rcpp::List &allParams,

 const std::string &uncertainty,

-const Rcpp::NumericVector &indices=Rcpp::NumericVector(),

-const Rcpp::NumericMatrix &initialA=Rcpp::NumericMatrix(),

-const Rcpp::NumericMatrix &initialP=Rcpp::NumericMatrix())

+const Rcpp::NumericVector &indices=Rcpp::NumericVector(1),

+const Rcpp::NumericMatrix &initialA=Rcpp::NumericMatrix(1,1),

+const Rcpp::NumericMatrix &initialP=Rcpp::NumericMatrix(1,1))

 {

-    return cogapsRun(data, params, unc, convertRMatrix(fixedMatrix),

-        checkpointInFile);

+    return cogapsRun(data, allParams, uncertainty, indices, initialA, initialP);

 }

 // [[Rcpp::export]]

 Rcpp::List cogaps_cpp(const Rcpp::NumericMatrix &data,

 const Rcpp::List &allParams,

 const Rcpp::NumericMatrix &uncertainty,

-const Rcpp::NumericVector &indices=Rcpp::NumericVector(),

-const Rcpp::NumericMatrix &initialA=Rcpp::NumericMatrix(),

-const Rcpp::NumericMatrix &initialP=Rcpp::NumericMatrix())

+const Rcpp::NumericVector &indices=Rcpp::NumericVector(1),

+const Rcpp::NumericMatrix &initialA=Rcpp::NumericMatrix(1,1),

+const Rcpp::NumericMatrix &initialP=Rcpp::NumericMatrix(1,1))

 {

-    return cogapsRun(convertRMatrix(data), params, convertRMatrix(unc),

-        convertRMatrix(fixedMatrix), checkpointInFile);

+    return cogapsRun(convertRMatrix(data), allParams,

+        convertRMatrix(uncertainty), indices, initialA, initialP);

 }

 // [[Rcpp::export]]

@@ -54,7 +54,11 @@ private:

     GapsDispatcher& operator=(const GapsDispatcher &p); // don't allow copies

     template <class DataType>

-    void loadData(const DataType &data);

+    void loadData(const DataType &data, bool transposeData);

+

+    template <class DataType>

+    void loadData(const DataType &data, bool transposeData, bool partitionRows,

+        const std:vector<unsigned> &indices);

 public:

@@ -62,10 +66,13 @@ public:

     ~GapsDispatcher();

     template <class DataType>

-    void initialize(const DataType &data, unsigned nPatterns, uint32_t seed=0);

+    void initialize(const DataType &data, bool transposeData,

+        bool partitionRows, const std::vector<unsigned> &indices,

+        unsigned nPatterns, uint32_t seed);

     template <class DataType>

-    void initialize(const DataType &data, const std::string &cptFile);

+    void initialize(const DataType &data, bool transposeData,

+        const std::string &cptFile);

     template <class DataType>

     void setUncertainty(const DataType &unc);

@@ -75,7 +82,10 @@ public:

     void setOutputFrequency(unsigned n);

     void setSparsity(float alphaA, float alphaP, bool singleCell);

     void setMaxGibbsMass(float maxA, float maxP);

-    void setFixedMatrix(char which, const RowMatrix &mat);

+

+    void setAMatrix(const RowMatrix &mat);

+    void setPMatrix(const RowMatrix &mat);

+    void setFixedMatrix(char which);

     void setNumCoresPerSet(unsigned n);

     void setCheckpointInterval(unsigned n);

@@ -85,14 +95,15 @@ public:

 };

 template <class DataType>

-void GapsDispatcher::initialize(const DataType &data, unsigned nPatterns,

+void GapsDispatcher::initialize(const DataType &data, bool transposeData,

+bool partitionRows, const std::vector<unsigned> &indices, unsigned nPatterns,

 uint32_t seed)

 {

     mSeed = seed;

     mNumPatterns = nPatterns;

     gaps::random::setSeed(mSeed);

-    loadData(data);

+    loadData(data, transposeData, partitionRows, indices);

     mInitialized = true;

 }

@@ -105,7 +116,7 @@ void GapsDispatcher::initialize(const DataType &data, const std::string &cptFile

     ar >> mSeed >> mNumPatterns >> mMaxIterations >> mPrintMessages >>

         mCheckpointsCreated >> mPhase;

-    loadData(data);

+    loadData(data, transposeData);

     ar >> *mRunners[0];

     mInitialized = true;

 }

@@ -118,11 +129,22 @@ void GapsDispatcher::setUncertainty(const DataType &unc)

 }

 template <class DataType>

-void GapsDispatcher::loadData(const DataType &data)

+void GapsDispatcher::loadData(const DataType &data, bool transposeData)

+{

+    gaps_printf("Loading Data...");

+    gaps_flush();

+    mRunners.push_back(new GapsRunner(data, transposeData, mNumPatterns));

+    gaps_printf("Done!\n");

+}

+

+template <class DataType>

+void GapsDispatcher::loadData(const DataType &data, bool transposeData,

+bool partitionRows, const std:vector<unsigned> &indices)

 {

     gaps_printf("Loading Data...");

     gaps_flush();

-    mRunners.push_back(new GapsRunner(data, mNumPatterns));

+    mRunners.push_back(new GapsRunner(data, transposeData, partitionRows,

+        indices, mNumPatterns));

     gaps_printf("Done!\n");

 }

@@ -37,7 +37,11 @@ private:

 public:

     template <class DataType>

-    GapsRunner(const DataType &data, unsigned nPatterns);

+    GapsRunner(const DataType &data, bool transposeData, unsigned nPatterns);

+

+    template <class DataType>

+    GapsRunner(const DataType &data, bool transposeData, bool partitionRows,

+        const std::vector<unsigned> &indices, unsigned nPatterns);

     template <class DataType>

     void setUncertainty(const DataType &unc);

@@ -70,9 +74,10 @@ public:

 };

 template <class DataType>

-GapsRunner::GapsRunner(const DataType &data, unsigned nPatterns)

+GapsRunner::GapsRunner(const DataType &data, bool transposeData, unsigned nPatterns)

     :

-mASampler(data, nPatterns), mPSampler(data, nPatterns),

+mASampler(data, transposeData, nPatterns),

+mPSampler(data, transposeData, nPatterns),

 mStatistics(mASampler.dataRows(), mPSampler.dataCols(), nPatterns),

 mSamplePhase(false), mNumUpdatesA(0), mNumUpdatesP(0), mFixedMatrix('N')

 {

@@ -81,10 +86,31 @@ mSamplePhase(false), mNumUpdatesA(0), mNumUpdatesP(0), mFixedMatrix('N')

 }

 template <class DataType>

-void GapsRunner::setUncertainty(const DataType &unc)

+GapsRunner::GapsRunner(const DataType &data, bool transposeData,

+bool partitionRows, const std::vector<unsigned> &indices, unsigned nPatterns)

+    :

+mASampler(data, transposeData, partitionRows, indices, nPatterns),

+mPSampler(data, transposeData, partitionRows, indices, nPatterns),

+mStatistics(mASampler.dataRows(), mPSampler.dataCols(), nPatterns),

+mSamplePhase(false), mNumUpdatesA(0), mNumUpdatesP(0), mFixedMatrix('N')

+{

+    mASampler.sync(mPSampler);

+    mPSampler.sync(mASampler);

+}

+

+template <class DataType>

+void GapsRunner::setUncertainty(const DataType &unc, bool transposeData)

+{

+    mASampler.setUncertainty(unc, transposeData);

+    mPSampler.setUncertainty(unc, transposeData);

+}

+

+template <class DataType>

+void GapsRunner::setUncertainty(const DataType &unc, bool transposeData,

+bool partitionRows, const std::vector<unsigned> &indices)

 {

-    mASampler.setUncertainty(unc);

-    mPSampler.setUncertainty(unc);

+    mASampler.setUncertainty(unc, transposeData, partitionRows, indices);

+    mPSampler.setUncertainty(unc, transposeData, partitionRows, indices);

 }

 #endif // __COGAPS_GAPS_RUNNER_H__

\ No newline at end of file

@@ -174,27 +174,75 @@ public:

 /******************** IMPLEMENTATION OF TEMPLATED FUNCTIONS *******************/

 template <class DataType>

-AmplitudeGibbsSampler::AmplitudeGibbsSampler(const DataType &data, unsigned nPatterns)

+AmplitudeGibbsSampler::AmplitudeGibbsSampler(const DataType &data,

+bool transposeData, unsigned nPatterns)

     :

-GibbsSampler(data, true, nPatterns)

+GibbsSampler(data, transposeData, true, nPatterns)

 {

     mQueue.setDimensionSize(mBinSize, mNumCols);

 }

 template <class DataType>

-PatternGibbsSampler::PatternGibbsSampler(const DataType &data, unsigned nPatterns)

+AmplitudeGibbsSampler::AmplitudeGibbsSampler(const DataType &data,

+bool transposeData, bool partitionRows, const std::vector<unsigned> &indices,

+unsigned nPatterns)

     :

-GibbsSampler(data, false, nPatterns)

+GibbsSampler(data, transposeData, partitionRows, indices, true, nPatterns)

+{

+    mQueue.setDimensionSize(mBinSize, mNumCols);

+}

+

+template <class DataType>

+PatternGibbsSampler::PatternGibbsSampler(const DataType &data,

+bool transposeData, unsigned nPatterns)

+    :

+GibbsSampler(data, transposeData, false, nPatterns)

+{

+    mQueue.setDimensionSize(mBinSize, mNumRows);

+}

+

+template <class DataType>

+PatternGibbsSampler::PatternGibbsSampler(const DataType &data,

+bool transposeData, bool partitionRows, const std::vector<unsigned> &indices,

+unsigned nPatterns)

+    :

+GibbsSampler(data, transposeData, partitionRows, indices, false, nPatterns)

 {

     mQueue.setDimensionSize(mBinSize, mNumRows);

 }

+// normal constructor

+template <class T, class MatA, class MatB>

+template <class DataType>

+GibbsSampler<T, MatA, MatB>::GibbsSampler(const DataType &data,

+bool transposeData, bool amp, unsigned nPatterns)

+    :

+mDMatrix(data, transposeData), mSMatrix(mDMatrix.pmax(0.1f, 0.1f)), 

+mAPMatrix(mDMatrix.nRow(), mDMatrix.nCol()),

+mMatrix(amp ? mDMatrix.nRow() : nPatterns, amp ? nPatterns : mDMatrix.nCol()),

+mOtherMatrix(NULL), mQueue(mMatrix.nRow() * mMatrix.nCol()), mLambda(0.f),

+mMaxGibbsMass(100.f), mAnnealingTemp(1.f), mNumRows(mMatrix.nRow()),

+mNumCols(mMatrix.nCol()), mAvgQueue(0.f), mNumQueues(0.f)

+{

+    // calculate atomic domain size

+    mBinSize = std::numeric_limits<uint64_t>::max()

+        / static_cast<uint64_t>(mNumRows * mNumCols);

+    mQueue.setDomainSize(mBinSize * mNumRows * mNumCols);

+    mDomain.setDomainSize(mBinSize * mNumRows * mNumCols);

+

+    // default sparsity parameters

+    setSparsity(0.01, false);

+}

+

+// constructor that allows for subsetting the data

 template <class T, class MatA, class MatB>

 template <class DataType>

 GibbsSampler<T, MatA, MatB>::GibbsSampler(const DataType &data,

+bool transposeData, bool partitionRows, const std::vector<unsigned> &indices,

 bool amp, unsigned nPatterns)

     :

-mDMatrix(data), mSMatrix(mDMatrix.pmax(0.1f, 0.1f)), 

+mDMatrix(data, transposeData, partitionRows, indices),

+mSMatrix(mDMatrix.pmax(0.1f, 0.1f)),

 mAPMatrix(mDMatrix.nRow(), mDMatrix.nCol()),

 mMatrix(amp ? mDMatrix.nRow() : nPatterns, amp ? nPatterns : mDMatrix.nCol()),

 mOtherMatrix(NULL), mQueue(mMatrix.nRow() * mMatrix.nCol()), mLambda(0.f),

@@ -3,86 +3,224 @@

 #include "../file_parser/MatrixElement.h"

 #include "../GapsAssert.h"

-template<class MatA, class MatB>

-static void copyMatrix(MatA &dest, const MatB &source)

+void RowMatrix::allocate()

 {

-    GAPS_ASSERT(dest.nRow() == source.nRow());

-    GAPS_ASSERT(dest.nCol() == source.nCol());

-    for (unsigned i = 0; i < source.nRow(); ++i)

+    for (unsigned i = 0; i < mNumRows; ++i)

+    {

+        mData.push_back(Vector(mNumCols));

+    }

+}

+

+float& RowMatrix::operator()(unsigned r, unsigned c)

+{

+    return mData[r][c];

+}

+

+float RowMatrix::operator()(unsigned r, unsigned c) const

+{

+    return mData[r][c];

+}

+

+Vector& getRow(unsigned row) {return mData[row];}

+const Vector& getRow(unsigned row) const {return mData[row];}

+

+float* rowPtr(unsigned row) {return mData[row].ptr();}

+const float* rowPtr(unsigned row) const {return mData[row].ptr();}

+

+void RowMatrix::allocate()

+{

+    for (unsigned i = 0; i < mNumRows; ++i)

+    {

+        mRows.push_back(Vector(mNumCols));

+    }

+}

+

+float& operator()(unsigned r, unsigned c) {return mData[c][r];}

+float operator()(unsigned r, unsigned c) const {return mData[c][r];}

+

+Vector& getCol(unsigned col) {return mData[col];}

+const Vector& getCol(unsigned col) const {return mData[col];}

+

+float* colPtr(unsigned col) {return mData[col].ptr();}

+const float* colPtr(unsigned col) const {return mData[col].ptr();}

+

+

+template <class MatA, class MatB>

+static void copyMatrix(MatA &dest, const MatB &source, bool transpose=false)

+{

+    GAPS_ASSERT(dest.nRow() == (transpose ? source.nCol() : source.nRow()));

+    GAPS_ASSERT(dest.nCol() == (transpose ? source.nRow() : source.nCol()))

+    for (unsigned i = 0; i < dest.nRow(); ++i)

     {

-        for (unsigned j = 0; j < source.nCol(); ++j)

+        for (unsigned j = 0; j < dest.nCol(); ++j)

         {

-            dest(i,j) = source(i,j);

+            dest(i,j) = transpose ? source(j,i) : source(i,j);

         }

     }

 }

-// if partitionRows is false, partition columns instead

-// rows of matrix should be partition dimension, i.e. need to transpose

-// is partitionRows is false

-/*

 template <class Matrix>

-static void fill(Matrix &mat, FileParser &p, bool partitionRows, std::vector<unsigned> whichIndices)

+static void createFromFile(Matrix &mat, const std::string &path, bool transpose)

 {

-    unsigned rowSelect = partitionRows ? 0 : 1;

-    unsigned colSelect = partitionRows ? 1 : 0;

+    FileParser p(path);

+    mNumRows = transpose ? p.nCol() : p.nRow();

+    mNumCols = transpose ? p.nRow() : p.nCol();

+    allocate();

     while (p.hasNext())

     {

         MatrixElement e(p.getNext());

-        std::vector<unsigned>::iterator newRowIndex = std::find(whichIndices.begin(), whichIndices.end(), e[rowSelect]);

-        if (newRowIndex != whichIndices.end())

+        unsigned row = transpose ? e.col() : e.row();

+        unsigned col = transpose ? e.row() : e.col();

+        this->operator()(row, col) = e.value();

+    }

+}

+

+// apply transpose first, then partition either rows or columns

+static void createFromMatrix(const RowMatrix &mat, bool transpose, bool partitionRows,

+const std::vector<unsigned> &indices)

+{

+    // create matrix

+    mNumRows = partitionRows ? indices.size() : (transpose ? mat.nCol() : mat.nRow());

+    mNumCols = partitionRows ? (transpose ? mat.nRow() : mat.nCol()) : indices.size();

+    allocate();

+

+    // fill matrix

+    for (unsigned i = 0; i < mat.nRow(); ++i)

+    {

+        for (unsigned j = 0; j < mat.nCol(); ++j)

         {

-            unsigned row = std::distance(whichIndices.begin(), newRowIndex);

-            mat.operator()(row, e[colSelect]) = e.value();

+            unsigned dataIndex = transpose ? (partitionRows ? j : i) : (partitionRows ? i : j);

+            std::vector<unsigned>::iterator pos = std::find(indices.begin(), indices.end(), dataIndex);

+            if (pos != indices.end())

+            {

+                unsigned index = std::distance(indices.begin(), pos);

+                unsigned row = partitionRows ? index : (transpose ? j : i);

+                unsigned col = partitionRows ? (transpose ? i : j) : index;

+                this->operator()(row,col) = mat(i,j)

+            }

         }

     }

 }

-*/

 /****************************** ROW MATRIX *****************************/

-RowMatrix::RowMatrix(unsigned nrow, unsigned ncol)

-: mNumRows(nrow), mNumCols(ncol)

+void RowMatrix::allocate()

 {

     for (unsigned i = 0; i < mNumRows; ++i)

     {

         mRows.push_back(Vector(mNumCols));

     }

+}    

+

+RowMatrix::RowMatrix(unsigned nrow, unsigned ncol)

+: mNumRows(nrow), mNumCols(ncol)

+{

+    allocate();

 }

 RowMatrix::RowMatrix(const ColMatrix &mat)

 : mNumRows(mat.nRow()), mNumCols(mat.nCol())

 {

-    for (unsigned i = 0; i < mNumRows; ++i)

+    allocate();

+    copyMatrix(*this, mat);

+}

+

+RowMatrix::RowMatrix(const RowMatrix &mat, bool transpose)

+    :

+mNumRows(transpose ? mat.nCol() : mat.nRow()),

+mNumCols(transpose ? mat.nRow() : mat.nCol())

+{

+    allocate();

+    copyMatrix(*this, mat, true);

+}

+

+RowMatrix::RowMatrix(const std::string &path, bool transpose)

+{

+    createFromFile(path, transpose);

+}

+

+// apply transpose first, then partition either rows or columns

+RowMatrix::RowMatrix(const RowMatrix &mat, bool transpose, bool partitionRows,

+const std::vector<unsigned> &indices)

+{

+    // create matrix

+    mNumRows = partitionRows ? indices.size() : (transpose ? mat.nCol() : mat.nRow());

+    mNumCols = partitionRows ? (transpose ? mat.nRow() : mat.nCol()) : indices.size();

+    allocate();

+

+    // fill matrix

+    for (unsigned i = 0; i < mat.nRow(); ++i)

     {

-        mRows.push_back(Vector(mNumCols));

-        for (unsigned j = 0; j < mNumCols; ++j)

+        for (unsigned j = 0; j < mat.nCol(); ++j)

         {

-            this->operator()(i,j) = mat(i,j);

+            unsigned dataIndex = transpose ? (partitionRows ? j : i) : (partitionRows ? i : j);

+            std::vector<unsigned>::iterator pos = std::find(indices.begin(), indices.end(), dataIndex);

+            if (pos != indices.end())

+            {

+                unsigned index = std::distance(indices.begin(), pos);

+                unsigned row = partitionRows ? index : (transpose ? j : i);

+                unsigned col = partitionRows ? (transpose ? i : j) : index;

+                this->operator()(row,col) = mat(i,j)

+            }

         }

     }

 }

-RowMatrix::RowMatrix(const std::string &path)

+// apply transpose first, then partition either rows or columns

+RowMatrix::RowMatrix(const std::string &path, bool transpose, bool partitionRows,

+const std::vector<unsigned> &indices)

 {

+    // create matrix

     FileParser p(path);

-    mNumRows = p.nRow();

-    mNumCols = p.nCol();

-

-    // allocate matrix

-    for (unsigned i = 0; i < mNumRows; ++i)

+    mNumRows = partitionRows ? indices.size() : (transpose ? p.nCol() : p.nRow());

+    mNumCols = partitionRows ? (transpose ? p.nRow() : p.nCol()) : indices.size();

+    allocate();

+

+    // fill

+    unsigned indexLoc = transpose ? (partitionRows ? 1 : 0) : (partitionRows ? 0 : 1);

+    unsigned rowLoc = transpose ? 1 : 0;

+    unsigned colLoc = transpose ? 0 : 1;

+    while (p.hasNext())

     {

-        mRows.push_back(Vector(mNumCols));

+        MatrixElement e(p.getNext());

+

+        unsigned dataIndex = e[indexLoc];

+        std::vector<unsigned>::iterator pos = std::find(indices.begin(), indices.end(), dataIndex);

+        if (pos != indices.end())

+        {

+            unsigned index = std::distance(indices.begin(), pos);

+            unsigned row = partitionRows ? index : e[rowLoc];

+            unsigned col = partitionRows ? e[colLoc] : index;

+            this->operator()(row, col) = e.value();

+        }

     }

+}

+

+

+// if partitionRows is false, partition columns instead

+// rows of matrix should be partition dimension, i.e. need to transpose

+// is partitionRows is false

+/*

+template <class Matrix>

+static void fill(Matrix &mat, FileParser &p, bool partitionRows, std::vector<unsigned> whichIndices)

+{

+    unsigned rowSelect = partitionRows ? 0 : 1;

+    unsigned colSelect = partitionRows ? 1 : 0;

-    // populate matrix

     while (p.hasNext())

     {

         MatrixElement e(p.getNext());

-        this->operator()(e.row(), e.col()) = e.value();

+        std::vector<unsigned>::iterator newRowIndex = std::find(whichIndices.begin(), whichIndices.end(), e[rowSelect]);

+        if (newRowIndex != whichIndices.end())

+        {

+            unsigned row = 

+            mat.operator()(row, e[colSelect]) = e.value();

+        }

     }

 }

+*/

+

 /*

 RowMatrix::RowMatrix(FileParser &p, bool partitionRows, std::vector<unsigned> whichIndices)

@@ -8,6 +8,174 @@

 #include <algorithm>

 #include <vector>

+/****************************** MATRIX INTERFACE ******************************/

+

+class RowMatrix;

+class ColMatrix;

+typedef Matrix RowMatrix; // when we don't care about row/col major order

+

+template <class T>

+class GenericMatrix

+{

+protected:

+

+    std::vector<Vector> mData;

+    unsigned mNumRows, mNumCols;

+

+    T* impl();

+

+public:

+

+    GenericMatrix(unsigned nrow, unsigned ncol);

+

+    GenericMatrix(const Matrix &mat, bool transpose=false);

+    GenericMatrix(const std::string &path, bool transpose=false);

+

+    GenericMatrix(const Matrix &mat, bool transpose, bool partitionRows,

+        const std::vector<unsigned> &indices);

+    GenericMatrix(const std::string &path, bool transpose, bool partitionRows,

+        const std::vector<unsigned> &indices);

+

+    unsigned nRow() const;

+    unsigned nCol() const;

+

+    T operator*(float val) const;

+    T operator/(float val) const;

+    T pmax(float scale, float max) const;

+

+    friend Archive& operator<< <T>(Archive &ar, GenericMatrix &mat);

+    friend Archive& operator>> <T>(Archive &ar, GenericMatrix &mat);

+};

+

+class RowMatrix : public Matrix<RowMatrix>

+{

+private:

+

+    friend class GenericMatrix;

+    void allocate();

+

+public:

+

+    RowMatrix(unsigned nrow, unsigned ncol);

+

+    RowMatrix(const Matrix &mat, bool transpose=false);

+    RowMatrix(const std::string &path, bool transpose=false);

+

+    RowMatrix(const Matrix &mat, bool transpose, bool partitionRows,

+        const std::vector<unsigned> &indices);

+    RowMatrix(const std::string &path, bool transpose, bool partitionRows,

+        const std::vector<unsigned> &indices);

+

+    // simplifies assigning to/from standard Matrix (RowMatrix)

+    RowMatrix& operator=(const ColMatrix &mat);

+

+    // for single element access - do not loop over elements with this

+    float& operator()(unsigned r, unsigned c) {return mData[r][c];}

+    float operator()(unsigned r, unsigned c) const {return mData[r][c];}

+

+    // for convenience when doing non-performance critical math

+    Vector& getRow(unsigned row) {return mData[row];}

+    const Vector& getRow(unsigned row) const {return mData[row];}

+

+    // raw pointer allows for high performance looping over rows with SIMD

+    float* rowPtr(unsigned row) {return mData[row].ptr();}

+    const float* rowPtr(unsigned row) const {return mData[row].ptr();}

+};

+

+class ColMatrix : public Matrix<ColMatrix>

+{

+private:

+

+    friend class GenericMatrix;

+    void allocate();

+

+public:

+

+    ColMatrix(unsigned nrow, unsigned ncol);

+

+    ColMatrix(const Matrix &mat, bool transpose=false);

+    ColMatrix(const std::string &path, bool transpose=false);

+

+    ColMatrix(const Matrix &mat, bool transpose, bool partitionRows,

+        const std::vector<unsigned> &indices);

+    ColMatrix(const std::string &path, bool transpose, bool partitionRows,

+        const std::vector<unsigned> &indices);

+

+    // simplifies assigning to/from standard Matrix (RowMatrix)

+    ColMatrix& operator=(const RowMatrix &mat);

+

+    // for single element access - do not loop over elements with this

+    float& operator()(unsigned r, unsigned c) {return mData[c][r];}

+    float operator()(unsigned r, unsigned c) const {return mData[c][r];}

+

+    // for convenience when doing non-performance critical math

+    Vector& getCol(unsigned col) {return mData[col];}

+    const Vector& getCol(unsigned col) const {return mData[col];}

+

+    // raw pointer allows for high performance looping over rows with SIMD

+    float* colPtr(unsigned col) {return mData[col].ptr();}

+    const float* colPtr(unsigned col) const {return mData[col].ptr();}

+};