Addressed jigsaw slowness issues, both in iterations and error propagation,... (250bcdef) · Commits · aflab / astrogeology / ISIS3

isis/src/base/objs/SparseBlockMatrix/SparseBlockMatrix.cpp

+31 −1

Original line number	Diff line number	Diff line
		@@ -23,6 +23,14 @@ using namespace boost::numeric::ublas;

		namespace Isis {

		/**
		* Default constructor.
		*/
		SparseBlockColumnMatrix::SparseBlockColumnMatrix() {
		m_startColumn = 0;
		}


		/**
		* Destructor. See description of wipe method below.
		*/
		@@ -71,6 +79,8 @@ namespace Isis {
		// insert matrix into map
		this->insert(it.key(),m);
		}

		m_startColumn = src.startColumn();
		}


		@@ -123,6 +133,26 @@ namespace Isis {
		}


		/**
		* Sets starting column for block in full matrix.
		*
		* @param nStartColumn value for starting column in full matrix for this block columns
		*/
		void SparseBlockColumnMatrix::setStartColumn(int nStartColumn) {
		m_startColumn = nStartColumn;
		}


		/**
		* Sets starting column for block in full matrix.
		*
		* @return int returns the starting column in the full matrix
		*/
		int SparseBlockColumnMatrix::startColumn() const {
		return m_startColumn;
		}


		/**
		* Returns total number of matrix elements in map (NOTE: NOT the number of matrix blocks). The sum
		* of all the elements in all of the matrix blocks.

isis/src/base/objs/SparseBlockMatrix/SparseBlockMatrix.h

+12 −1

Original line number	Diff line number	Diff line
		@@ -62,12 +62,15 @@ namespace Isis {
		* ISIS coding standards.
		* @history 2016-08-10 Jeannie Backer - Replaced boost matrix with Isis::LinearAlgebra::Matrix.
		* References #4163.
		* @history 2017-05-09 Ken Edmundson - Added m_startColumn member and mutator/accessor methods
		* to SparseBlockColumnMatrix. Done to eliminate lengthy computation of
		* leading colums and rows. References #4664.
		*/
		class SparseBlockColumnMatrix :
		public QMap< int, LinearAlgebra::Matrix * > {

		public:
		SparseBlockColumnMatrix(){} // default constructor
		SparseBlockColumnMatrix(); // default constructor
		~SparseBlockColumnMatrix(); // destructor

		// copy constructor
		@@ -80,11 +83,19 @@ namespace Isis {

		void zeroBlocks();
		bool insertMatrixBlock(int nColumnBlock, int nRows, int nCols);

		void setStartColumn(int nStartColumn);
		int startColumn() const;
		int numberOfElements();
		int numberOfRows();
		int numberOfColumns();
		void print(std::ostream& outstream);
		void printClean(std::ostream& outstream);

		protected:
		int m_startColumn; /**< starting column for this Block Column in full matrix
		e.g. for Block Column 4, if the preceding Block Columns each have 6
		columns, then the starting column for Block Column 4 is 24 */
		};

		// operators to read/write SparseBlockColumnMatrix to/from binary disk file

isis/src/base/objs/SparseBlockMatrix/SparseBlockMatrix.truth

+22 −0

Original line number	Diff line number	Diff line
		@@ -149,6 +149,28 @@ Printing SparseBlockColumnMatrix...

		Empty SparseBlockColumnMatrix...

		----- Insert boost block matrix of zeros in each column, set values, and startColumn
		# block columns: 3
		# matrix blocks: 3
		# diagonal matrix blocks: 1
		# off-diagonal matrix blocks: 2
		block 0 startColumn: 0
		block 1 startColumn: 3
		block 2 startColumn: 6
		Printing SparseBlockMatrix...
		Printing SparseBlockColumnMatrix...
		0
		[3,3]((0,1,2),(3,4,5),(6,7,8))

		Printing SparseBlockColumnMatrix...
		2
		[3,3]((0,1,2),(3,4,5),(6,7,8))

		Printing SparseBlockColumnMatrix...
		3
		[3,3]((0,1,2),(3,4,5),(6,7,8))


		----- copy method
		# matrix blocks: 1
		# diagonal matrix blocks: 1

isis/src/base/objs/SparseBlockMatrix/unitTest.cpp

+41 −0

Original line number	Diff line number	Diff line
		@@ -346,6 +346,47 @@ int main(int argc, char *argv[]) {
		e.print();
		}

		try {
		cerr << endl << "----- Insert boost block matrix of zeros in each column, set values, and startColumn" << endl;
		SparseBlockMatrix sbm;
		sbm.setNumberOfColumns(3);
		sbm.insertMatrixBlock(0, 0, 3, 3);
		sbm.insertMatrixBlock(1, 2, 3, 3);
		sbm.insertMatrixBlock(2, 3, 3, 3);

		for ( int i = 0; i < 3; i++ ) {
		for ( int j = 0; j < 3; j++ ) {
		((sbm[0])[0])(i,j) = 3*i+j;
		}
		}
		for ( int i = 0; i < 3; i++ ) {
		for ( int j = 0; j < 3; j++ ) {
		((sbm[1])[2])(i,j) = 3*i+j;
		}
		}
		for ( int i = 0; i < 3; i++ ) {
		for ( int j = 0; j < 3; j++ ) {
		((sbm[2])[3])(i,j) = 3*i+j;
		}
		}

		sbm.at(0)->setStartColumn(0);
		sbm.at(1)->setStartColumn(3);
		sbm.at(2)->setStartColumn(6);

		cerr << " # block columns: " << sbm.size() << endl;
		cerr << " # matrix blocks: " << sbm.numberOfBlocks() << endl;
		cerr << " # diagonal matrix blocks: " << sbm.numberOfDiagonalBlocks() << endl;
		cerr << "# off-diagonal matrix blocks: " << sbm.numberOfOffDiagonalBlocks() << endl;
		cerr << " block 0 startColumn: " << sbm.at(0)->startColumn() << endl;
		cerr << " block 1 startColumn: " << sbm.at(1)->startColumn() << endl;
		cerr << " block 2 startColumn: " << sbm.at(2)->startColumn() << endl;
		sbm.print(std::cerr);
		}
		catch(IException &e) {
		e.print();
		}

		try {
		cerr << endl << "----- copy method" << endl;
		SparseBlockMatrix sbm;

isis/src/control/objs/BundleAdjust/BundleAdjust.cpp

+90 −66

Original line number	Diff line number	Diff line
		@@ -327,6 +327,9 @@ namespace Isis {
		* @todo answer comments with questions, TODO, ???, and !!!
		*/
		void BundleAdjust::init(Progress *progress) {

		m_previousNumberImagePartials = 0;

		// initialize
		//
		// JWB
		@@ -490,6 +493,8 @@ namespace Isis {
		// initializations for cholmod
		initializeCHOLMODLibraryVariables();

		// initialize normal equations matrix
		initializeNormalEquationsMatrix();
		}


		@@ -546,12 +551,11 @@ namespace Isis {

		/**
		* Initializations for CHOLMOD sparse matrix package.
		* Calls cholmod_start, sets m_cholmodCommon options, and resizes m_sparseNormals.
		* Calls cholmod_start, sets m_cholmodCommon options.
		*
		* @return @b bool If the CHOLMOD library variables were successfully initialized.
		*/
		bool BundleAdjust::initializeCHOLMODLibraryVariables() {

		if ( m_rank <= 0 ) {
		return false;
		}
		@@ -567,12 +571,48 @@ namespace Isis {
		m_cholmodCommon.nmethods = 1;
		m_cholmodCommon.method[0].ordering = CHOLMOD_AMD;

		// set size of sparse block normal equations matrix
		return true;
		}


		/**
		* Initialize Normal Equations matrix (m_sparseNormals).
		*
		* Ken NOTE: Currently we are explicitly setting the start column for each block in the normal
		* equations matrix below. I think it should be possible (and relatively easy) to make
		* the m_sparseNormals smart enough to set the start column of a column block
		* automatically when it is added to the matrix.
		*
		* @return @b bool.
		*/
		bool BundleAdjust::initializeNormalEquationsMatrix() {

		int nBlockColumns = m_bundleObservations.size();

		if (m_bundleSettings->solveTargetBody())
		nBlockColumns += 1;

		m_sparseNormals.setNumberOfColumns(nBlockColumns);

		m_sparseNormals.at(0)->setStartColumn(0);

		int nParameters = 0;
		if (m_bundleSettings->solveTargetBody()) {
		m_sparseNormals.setNumberOfColumns(m_bundleObservations.size()+1);
		nParameters += m_bundleSettings->numberTargetBodyParameters();
		m_sparseNormals.at(1)->setStartColumn(nParameters);

		int observation = 0;
		for (int i = 2; i < nBlockColumns; i++) {
		nParameters += m_bundleObservations.at(observation)->numberParameters();
		m_sparseNormals.at(i)->setStartColumn(nParameters);
		observation++;
		}
		}
		else {
		m_sparseNormals.setNumberOfColumns(m_bundleObservations.size());
		for (int i = 0; i < nBlockColumns; i++) {
		m_sparseNormals.at(i)->setStartColumn(nParameters);
		nParameters += m_bundleObservations.at(i)->numberParameters();
		}
		}

		return true;
		@@ -1031,8 +1071,10 @@ namespace Isis {
		// update number of observations
		int numObs = m_bundleResults.numberObservations();
		m_bundleResults.setNumberObservations(numObs + 2);

		formMeasureNormals(N22, N12, n1, n2, coeffTarget, coeffImage, coeffPoint3D, coeffRHS,
		measure->observationIndex());

		} // end loop over this points measures

		formPointNormals(N22, N12, n2, m_RHS, point);
		@@ -1044,7 +1086,6 @@ namespace Isis {
		} // end loop over 3D points

		// finally, form the reduced normal equations

		formWeightedNormals(n1, m_RHS);

		// update number of unknown parameters
		@@ -1113,9 +1154,9 @@ namespace Isis {
		N11TargetImage.clear();
		N11TargetImage = prod(trans(coeffTarget),coeffImage);

		m_sparseNormals.insertMatrixBlock(observationIndex+1, 0,
		m_sparseNormals.insertMatrixBlock(blockIndex, 0,
		numTargetPartials, coeffImage.size2());
		((m_sparseNormals[observationIndex+1])[0]) += N11TargetImage;
		((m_sparseNormals[blockIndex])[0]) += N11TargetImage;

		// form N12 target portion
		static matrix<double> N12Target(numTargetPartials, 3);
		@@ -1153,14 +1194,7 @@ namespace Isis {

		N11 = prod(trans(coeffImage), coeffImage);

		int t = 0;
		//testing
		for (int a = 0; a < observationIndex; a++) {
		BundleObservationQsp observation = m_bundleObservations.at(a);
		t += observation->numberParameters();
		}
		// account for target parameters
		t += numTargetPartials;
		int t = m_sparseNormals.at(blockIndex)->startColumn();

		// insert submatrix at column, row
		m_sparseNormals.insertMatrixBlock(blockIndex, blockIndex,
		@@ -1261,7 +1295,6 @@ namespace Isis {
		bundleControlPoint->setAdjustedSurfacePoint(SurfacePoint);

		// form Q (this is N22{-1} * N12{T})
		// Q = prod(N22, trans(N12));
		productATransB(N22, N12, Q);

		// form product of N22(inverse) and n2; store in NIC
		@@ -1271,7 +1304,6 @@ namespace Isis {
		productAB(N12, Q);

		// accumulate -nj
		// nj -= prod(trans(Q),n2);
		accumProductAlphaAB(-1.0, Q, n2, nj);

		return true;
		@@ -1372,7 +1404,7 @@ namespace Isis {

		int columnIndex = Qit.key();

		subrangeStart = m_sparseNormals.getLeadingColumnsForBlock(columnIndex);
		subrangeStart = m_sparseNormals.at(columnIndex)->startColumn();
		subrangeEnd = subrangeStart + Qit.value()->size2();

		v2 += alpha * prod(*(Qit.value()),subrange(v1,subrangeStart,subrangeEnd));
		@@ -1473,7 +1505,7 @@ namespace Isis {
		return;
		}

		int numTargetParameters = m_bundleSettings->numberTargetBodyParameters();
		int numParams;

		QMapIterator<int, LinearAlgebra::Matrix*> Qit(Q);

		@@ -1485,22 +1517,7 @@ namespace Isis {

		LinearAlgebra::Vector blockProduct = prod(trans(*Qblock),n2);

		int numParams = 0;
		for (int observationIndex = 0; observationIndex < columnIndex; observationIndex++) {
		if (numTargetParameters > 0 && observationIndex == 0) {
		numParams += numTargetParameters;
		}
		else {
		if (numTargetParameters > 0 ) {
		BundleObservationQsp observation = m_bundleObservations.at(observationIndex-1);
		numParams += observation->numberParameters();
		}
		else {
		BundleObservationQsp observation = m_bundleObservations.at(observationIndex);
		numParams += observation->numberParameters();
		}
		}
		}
		numParams = m_sparseNormals.at(columnIndex)->startColumn();

		for (unsigned i = 0; i < blockProduct.size(); i++) {
		nj(numParams+i) += alpha*blockProduct(i);
		@@ -1622,7 +1639,7 @@ namespace Isis {
		return false;
		}

		int numLeadingColumns = m_sparseNormals.getLeadingColumnsForBlock(columnIndex);
		int numLeadingColumns = normalsColumn->startColumn();

		QMapIterator< int, LinearAlgebra::Matrix * > it(*normalsColumn);

		@@ -1631,7 +1648,9 @@ namespace Isis {

		int rowIndex = it.key();

		int numLeadingRows = m_sparseNormals.getLeadingRowsForBlock(rowIndex);
		// note: as the normal equations matrix is symmetric, the # of leading rows for a block is
		// equal to the # of leading columns for a block column at the "rowIndex" position
		int numLeadingRows = m_sparseNormals.at(rowIndex)->startColumn();

		LinearAlgebra::Matrix *normalsBlock = it.value();
		if ( !normalsBlock ) {
		@@ -1816,7 +1835,14 @@ namespace Isis {
		BundleObservationQsp observation = measure.parentBundleObservation();

		int numImagePartials = observation->numberParameters();

		// we're saving the number of image partials in m_previousNumberImagePartials
		// to compare to the previous computePartials call to avoid unnecessary resizing of the
		// coeffImage matrix
		if (numImagePartials != m_previousNumberImagePartials) {
		coeffImage.resize(2,numImagePartials);
		m_previousNumberImagePartials = numImagePartials;
		}

		// clear partial derivative matrices and vectors
		if (m_bundleSettings->solveTargetBody()) {
		@@ -2093,8 +2119,6 @@ namespace Isis {
		t += numParameters;
		}

		// TODO: CHECK - do we need point index in case of rejected points????

		// TODO: Below code should move into BundleControlPoint->updateParameterCorrections
		// except, what about the productAlphaAV method?

		@@ -2116,7 +2140,6 @@ namespace Isis {
		boost::numeric::ublas::bounded_vector< double, 3 > &corrections = point->corrections();

		// subtract product of Q and nj from NIC
		// NIC -= prod(Q, m_imageSolution);
		productAlphaAV(-1.0, NIC, Q, m_imageSolution);

		// get point parameter corrections
		@@ -2599,20 +2622,19 @@ namespace Isis {
		int columnIndex = 0;
		int numColumns = 0;
		int numBlockColumns = m_sparseNormals.size();

		for (i = 0; i < numBlockColumns; i++) {

		// columns in this column block
		SparseBlockColumnMatrix *normalsColumn = m_sparseNormals.at(i);
		if (i == 0) {
		numColumns = normalsColumn->numberOfColumns();
		int numRows = m_sparseNormals.at(i)->numberOfRows();
		int numRows = normalsColumn->numberOfRows();
		inverseMatrix.insertMatrixBlock(i, numRows, numColumns);
		inverseMatrix.zeroBlocks();
		}
		else {
		if (normalsColumn->numberOfColumns() == numColumns) {
		int numRows = m_sparseNormals.at(i)->numberOfRows();
		int numRows = normalsColumn->numberOfRows();
		inverseMatrix.insertMatrixBlock(i, numRows, numColumns);
		inverseMatrix.zeroBlocks();
		}
		@@ -2710,7 +2732,9 @@ namespace Isis {
		}

		// get corresponding Q matrix
		SparseBlockRowMatrix Q = point->cholmodQMatrix();
		// NOTE: we are getting a reference to the Q matrix stored
		// in the BundleControlPoint for speed (without the & it is dirt slow)
		SparseBlockRowMatrix &Q = point->cholmodQMatrix();

		T.clear();