src.main/CjMatrixImpl_1D.hh

/* -*- Mode: C; indent-tabs-mode: t; c-basic-offset: 4; tab-width: 4 -*- */
/*
 * mitie
 * Copyright (C) Jason Spencer 2009 <mitie@jasonspencer.org>
 * 
 * mitie 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 3 of the License, or
 * (at your option) any later version.
 * 
 * mitie 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, see <http://www.gnu.org/licenses/>.
 */

#ifndef _CJMATRIXIMPL_1D_H_
#define _CJMATRIXIMPL_1D_H_ 1

#include "CjMatrix.hh"
#include "CjMatrixImpl.hh"
#include <iostream>     // for stdout debugging only
#include <boost/tokenizer.hpp>
#include <sstream>      // std::ostringstream
// for use in file load ctor - TODO can remove once it's made a free function
#include <vector>
#include <string>
#include <fstream>
#include <iterator>     // for back_inserter
#include <algorithm>    // for std::swap
#include <cassert>
#include "JLOG.hh"
#include <boost/lexical_cast.hpp>

template <typename> class CjMatrix;     // fwd decl
template <typename> class CjMatrixImpl; // fwd decl

template <typename T> class CjMatrixImpl_1D: public CjMatrixImpl<T> {

protected:
T * m_data;
unsigned int m_width;
unsigned int m_height;
unsigned int m_size;
T m_defaultValue;
bool m_defaultValueSet;

public:


CjMatrixImpl_1D ( const unsigned int width, const unsigned int height, const T * const defaultValue = 0 );

explicit CjMatrixImpl_1D (const CjMatrixImpl_1D<T> & cjm);


CjMatrixImpl_1D (const CjMatrix<T> & cjm, const unsigned int width, const unsigned int height, const T * const defaultValue = 0); // throws: E_CjMatrix_MemoryAllocFail, E_CjMatrix_SizeError


template <typename V> CjMatrixImpl_1D (V m, const unsigned int w, const unsigned int h, const unsigned int width, const unsigned int height, const T * const defaultValue = 0); // throws: E_CjMatrix_MemoryAllocFail, E_CjMatrix_SizeError

~CjMatrixImpl_1D();


void swap(CjMatrixImpl_1D & other);


CjMatrixImpl<T> * clone() const;

void resizeImpl(const unsigned int width, const unsigned int height);
void resize_with_copyImpl( const unsigned int width, const unsigned int height, const CjMatrix<T> * const cjm = (CjMatrix<T> *) 0);
void setColumnValuesImpl( const unsigned int col, const T & val );
void setRowValuesImpl( const unsigned int row, const T & val );
void setAllValuesImpl(const T & val);
void setDiagonalImpl(const T & val);
void setDefaultValueImpl(const T & val);
T getDefaultValueImpl() const;
bool getDefaultValueSetImpl() const;
bool equalElements(const CjMatrix<T> &) const;
void transposeImpl();
bool isSymmetricImpl() const;
unsigned int getWidthImpl() const;
unsigned int getHeightImpl() const;
T & atImpl(unsigned int row, unsigned int col);
const T & atImpl(unsigned int row, unsigned int col) const;
T & getImpl (unsigned int row, unsigned int col);
const T & getImpl (unsigned int row, unsigned int col) const;

// additional methods:
CjMatrixImpl_1D<T> & operator=(const CjMatrixImpl_1D<T> & );

private:

void resize_with_copyImpl(const unsigned int width, const unsigned int height, const CjMatrixImpl_1D<T> * const cjmi );

};

template <typename T> inline CjMatrixImpl<T> * CjMatrixImpl_1D<T>::clone() const {
        CjMatrixImpl_1D<T> * retval = new CjMatrixImpl_1D<T>(*this);
        return retval;
}


template <typename T> inline void CjMatrixImpl_1D<T>::swap(CjMatrixImpl_1D<T> & that) {
        std::swap(that.m_data, m_data);
        std::swap(that.m_width, m_width);
        std::swap(that.m_height, m_height);
        std::swap(that.m_size, m_size);
        std::swap(that.m_defaultValue, m_defaultValue);
        std::swap(that.m_defaultValueSet, m_defaultValueSet);
}


template <typename T> inline CjMatrixImpl_1D<T>::CjMatrixImpl_1D (const unsigned int width, const unsigned int height, const T * const defaultValue):m_data((T *)0),m_width(0),m_height(0),m_size(0),m_defaultValueSet(false) { // creates a matrix filled with the default value
JLOG(4, std::cout << "CjMatrixImpl_1D<T>::CjMatrixImpl_1D (const unsigned int (" << width << "), const unsigned int (" << height << "), T defaultValue ) called." << std::endl; )
        if(defaultValue) setDefaultValueImpl(*defaultValue);
        resizeImpl(width,height);
}


template <typename T> inline CjMatrixImpl_1D<T>::CjMatrixImpl_1D (const CjMatrixImpl_1D<T> & cjmi):m_data((T *)0),m_width(0),m_height(0),m_size(0),m_defaultValueSet(false) {
JLOG(4, std::cout << "CjMatrixImpl_1D<T>::CjMatrixImpl_1D (const CjMatrixImpl_1D &) called." << std::endl; )
        operator=(cjmi);        // could throw E_CjMatrix_MemoryAllocFail
}


template <typename T> inline CjMatrixImpl_1D<T> & CjMatrixImpl_1D<T>::operator=(const CjMatrixImpl_1D<T> & cjmi) { // could throw E_CjMatrix_MemoryAllocFail
JLOG(4, std::cout << "CjMatrixImpl_1D<T>::operator=(const CjMatrixImpl_1D<T> &) called." << std::endl;)
        if (this != &cjmi) {
                this->resize_with_copyImpl(cjmi.getWidthImpl(), cjmi.getHeightImpl(), &cjmi );  // could throw E_CjMatrix_MemoryAllocFail
        }
        return *this;
}

template <typename T> inline bool CjMatrixImpl_1D<T>::getDefaultValueSetImpl() const { return m_defaultValueSet; }

// create a new CjMatrixImpl_1D of 'height' and 'width' and copies the intersection with cjm into the upper left corner. If defaultValue is specified the non-intersection elements will be a copy of this dereferenced value
// throws: E_CjMatrix_MemoryAllocFail, E_CjMatrix_SizeError
template <typename T> CjMatrixImpl_1D<T>::CjMatrixImpl_1D (const CjMatrix<T> & cjm, const unsigned int width, const unsigned int height, const T * const defaultValue):m_data((T *)0),m_width(0),m_height(0),m_size(0),m_defaultValueSet(false) {
JLOG(4, std::cout << "CjMatrixImpl_1D<T>::CjMatrixImpl_1D (const CjMatrixImpl_1D & cjm, unsigned int width, unsigned int height, T * defaultValue = 0) called." << std::endl; )
        if(defaultValue != 0) setDefaultValueImpl(*defaultValue);
        resize_with_copyImpl(width, height, &cjm );  // could throw E_CjMatrix_MemoryAllocFail
}

// constructs a width*height array and copies the intersection with w*h elements from m
// throws: E_CjMatrix_MemoryAllocFail, E_CjMatrix_SizeError
// type V must be dereferencable by [][]
template <typename T> template <typename V> CjMatrixImpl_1D<T>::CjMatrixImpl_1D (V m, const unsigned int w, const unsigned int h, const unsigned int width, const unsigned int height, const T * const defaultValue):m_data((T *)0),m_width(0),m_height(0),m_size(0),m_defaultValueSet(false) {
        //cannot be implemented in terms of operator= or copy constructor as they call this function
        if((0==width)||(0==height)) throw E_CjMatrix_SizeError("Neither width nor height can be zero in CjMatrixImpl_1D (const T * const * const m, const unsigned int w, const unsigned int h, const unsigned int width, const unsigned int height, const T * const defaultValue = 0)");
        // find width and height of overlap between input matrix and new matrix size - this is the size we will copy
        unsigned int min_h = std::min(h,height);
        unsigned int min_w = std::min(w,width);

        CjMatrixImpl_1D <T> tmp(width,height);
        // copy min_h rows and min_w columns
        JLOG(4, std::cout << " min_w is " << min_w << " and min_h is " << min_h << std::endl; )
        for(unsigned int row=0; row < min_h; ++row)
                std::copy( &m[row][0], &m[row][min_w], tmp.m_data+(row*width) );

        if(defaultValue != 0) tmp.setDefaultValueImpl(*defaultValue);
        // pad remainder with defaultValue
        if(m_defaultValueSet) {
                for(unsigned int col=min_w; col < width; ++col) tmp.setColumnValuesImpl(col, m_defaultValue);
                for(unsigned int row=min_h; row < height; ++row) tmp.setRowValuesImpl(row, m_defaultValue);
                }
        swap(tmp);
}

template <typename T> void CjMatrixImpl_1D<T>::resize_with_copyImpl(const unsigned int width, const unsigned int height, const CjMatrix<T> * const cjm ) {
        //std::cout << "Current size is " << m_width << "x" << m_height << " - resizing to " << width << "x" << height << std::endl;
        //cannot be implemented in terms of operator= or copy constructor as they call this function
        if((0==width)||(0==height)) throw E_CjMatrix_SizeError("Neither width nor height can be zero in CjMatrixImpl_1D<T>::resize_with_copyImpl");
        // find dimension of input matrix (either this or cjm)
        unsigned int old_w = (cjm==0)?m_width:cjm->getWidth();
        unsigned int old_h = (cjm==0)?m_height:cjm->getHeight();
        // find width and height of overlap between input matrix and new matrix size - this is the size we will copy
        unsigned int min_h = std::min(old_h,height);
        unsigned int min_w = std::min(old_w,width);

        CjMatrixImpl_1D <T> tmp(width,height);
        // copy min_h rows and min_w columns
        JLOG(4, std::cout << " min_w is " << min_w << " and min_h is " << min_h << std::endl; )
        for(unsigned int row=0; row < min_h; ++row)
                if(cjm) std::copy( cjm->row(row).cbegin(), cjm->row(row).cbegin()+min_w, tmp.m_data+(row*width) );
                else std::copy( m_data+(row*m_width), m_data+(row*m_width)+min_w, tmp.m_data+(row*width) );
        // pad remainder with defaultValue
        if(m_defaultValueSet) {
                for(unsigned int col=min_w; col < width; ++col) tmp.setColumnValuesImpl(col, m_defaultValue);
                for(unsigned int row=min_h; row < height; ++row) tmp.setRowValuesImpl(row, m_defaultValue);
                }
        swap(tmp);
}

template <typename T> void CjMatrixImpl_1D<T>::resize_with_copyImpl(const unsigned int width, const unsigned int height, const CjMatrixImpl_1D<T> * const cjmi ) {
        //std::cout << "Current size is " << m_width << "x" << m_height << " - resizing to " << width << "x" << height << std::endl;
        //cannot be implemented in terms of operator= or copy constructor as they call this function
        if(cjmi == this) return;
        if((0==width)||(0==height)) throw E_CjMatrix_SizeError("Neither width nor height can be zero in CjMatrixImpl_1D<T>::resize_with_copyImpl");
        // find dimension of input matrix (either this or cjm)
        unsigned int old_w = (cjmi==0)?m_width:cjmi->getWidthImpl();
        unsigned int old_h = (cjmi==0)?m_height:cjmi->getHeightImpl();
        // find width and height of overlap between input matrix and new matrix size - this is the size we will copy
        unsigned int min_h = std::min(old_h,height);
        unsigned int min_w = std::min(old_w,width);

        CjMatrixImpl_1D <T> tmp(width,height);
        // copy min_h rows and min_w columns
        JLOG(4, std::cout << " min_w is " << min_w << " and min_h is " << min_h << std::endl; )
        for(unsigned int row=0; row < min_h; ++row)
                if(cjmi) std::copy( cjmi->m_data+(row*width), cjmi->m_data+(row*width)+min_w, tmp.m_data+(row*width) );
                else std::copy( m_data+(row*m_width), m_data+(row*m_width)+min_w, tmp.m_data+(row*width) );
        // pad remainder with defaultValue
        if(m_defaultValueSet) {
                for(unsigned int col=min_w; col < width; ++col) tmp.setColumnValuesImpl(col, m_defaultValue);
                for(unsigned int row=min_h; row < height; ++row) tmp.setRowValuesImpl(row, m_defaultValue);
                }
        swap(tmp);
}

template <typename T> void CjMatrixImpl_1D<T>::setColumnValuesImpl(const unsigned int col, const T & val) {
        for(unsigned int i=0;i<m_height;++i) m_data[i*m_width+col] = val;
}

template <typename T> void CjMatrixImpl_1D<T>::setRowValuesImpl(const unsigned int row, const T & val) {
        std::fill_n( m_data+(row*m_width), m_width, val );
}

template <typename T> void CjMatrixImpl_1D<T>::setAllValuesImpl(const T & val) {
        std::fill( m_data, m_data+m_size, val );
}

template <typename T> void CjMatrixImpl_1D<T>::setDiagonalImpl(const T & val) {
        for( unsigned int i=0 ; i < std::min(m_width, m_height) ; ++i ) m_data[i*m_width+i] = val;
}

template <typename T> void CjMatrixImpl_1D<T>::resizeImpl(const unsigned int width, const unsigned int height) {
// this function cannot be implemented in terms of any of the constructors because many of them call this function or another function which calls this function
unsigned int new_size = width*height;
if((0==width)||(0==height)) throw E_CjMatrix_SizeError("Neither width nor height can be zero in CjMatrixImpl_1D<T>::resizeImpl");

if((m_size==new_size)&&(m_data != (T *) 0)) {
        m_width = width;
        m_height = height;
} else {
        try {
                T * new_arr = new T[new_size];
                std::swap(new_arr,m_data);
                m_size = new_size;
                m_width = width;
                m_height = height;
                delete [] new_arr;
                new_arr = (T *)0;
                
        } catch (std::bad_alloc &) { /* std::cerr << "** ALLOC ERROR WHEN RESIZING 1D MATRIX **" << std::endl; */throw E_CjMatrix_MemoryAllocFail(); }
        }
if (m_defaultValueSet) setAllValuesImpl(m_defaultValue);
}

template <typename T> void CjMatrixImpl_1D<T>::setDefaultValueImpl(const T & val) { m_defaultValue = val; m_defaultValueSet = true; }
template <typename T> T CjMatrixImpl_1D<T>::getDefaultValueImpl() const { return m_defaultValue; }

template <typename T> unsigned int CjMatrixImpl_1D<T>::getWidthImpl() const { return m_width; };
template <typename T> unsigned int CjMatrixImpl_1D<T>::getHeightImpl() const { return m_height; };

template <typename T> inline T & CjMatrixImpl_1D<T>::atImpl(const unsigned int row, const unsigned int col) { // throws E_CjMatrix_BoundsError
        if((row>=m_height)||(col>=m_width)) {
                throw E_CjMatrix_BoundsError( "CjMatrixImpl_1D<T>::atImpl(const unsigned int row, const unsigned int col): Bounds error - (" + boost::lexical_cast<std::string>(row) + "," + boost::lexical_cast<std::string>(col) + ") was requested but height = " + boost::lexical_cast<std::string>(m_height) + " width = " + boost::lexical_cast<std::string>(m_width) );
        }
        return getImpl(row,col);
}

template <typename T> inline const T & CjMatrixImpl_1D<T>::atImpl(const unsigned int row, const unsigned int col) const { // throws E_CjMatrix_BoundsError
        if((row>=m_height)||(col>=m_width)) {
                throw E_CjMatrix_BoundsError( "CjMatrixImpl_1D<T>::atImpl(const unsigned int row, const unsigned int col) const: Bounds error - (" + boost::lexical_cast<std::string>(row) + "," + boost::lexical_cast<std::string>(col) + ") was requested but height = " + boost::lexical_cast<std::string>(m_height) + " width = " + boost::lexical_cast<std::string>(m_width) );
        }
        return getImpl(row,col);
}

template <typename T> inline T & CjMatrixImpl_1D<T>::getImpl (const unsigned int row, const unsigned int col) {
        return *(m_data+col+(row*m_width));
}

template <typename T> inline const T & CjMatrixImpl_1D<T>::getImpl (const unsigned int row, const unsigned int col) const {
        return *(m_data+col+(row*m_width));
}

template <typename T> bool CjMatrixImpl_1D<T>::equalElements(const CjMatrix<T> & cjm) const {
        if ( ( cjm.getWidth() != m_width ) || ( cjm.getHeight() != m_height ) ) return false;
        bool different = false;
        JLOG(4, std::cout << "CjMatrixImpl_1D<T>::equalElements: comparing " << m_width << "x" << m_height << " against " << cjm.getWidth() << "x" << cjm.getHeight() << std::endl; )
        for(unsigned int row=0;(row<m_height) && !different;++row)
                for(unsigned int col=0;col<m_width;++col)
                        if(atImpl(row,col)!=cjm(row,col)) {
                                JLOG(4, std::cout << "CjMatrixImpl_1D<T>::equalElements: Element mismatch at " << row << "," << col << ": " << atImpl(row,col) << " vs " << cjm(row,col) << std::endl; )
                                different = true;
                                break;
                        }
        JLOG(4, std::cout << "CjMatrixImpl_1D<T>::equalElements: end." << std::endl; )
        return !different;
}


// the most elegant way to tranpose is just to set a flag which swaps row and column in operator().. but we want that to be absolutely as fast as possible - so we'll actually perform the tranpose.
template <typename T> void CjMatrixImpl_1D<T>::transposeImpl(){ // could throw E_CjMatrix_MemoryAllocFail
// we could make a temp matrix and use column iterators to copy to rows, but this will be faster and closer to the data structures

if((m_height==1)||(m_width==1)) { std::swap(m_width,m_height); return; }

if(m_height==m_width) {         // we're square - no need to realloc memory - just move values around..
        for(unsigned int i=0;i<m_height;++i)
                for(unsigned int j=i+1;j<m_width;++j) {
                        std::swap( *(m_data+j+(m_width*i)), *(m_data+i+(m_width*j)) );
                        }
        return;
} else {
        // ok, so we're not square or single row/column
        // There might be a fast way of transposing inplace and without the use of scratch memory but let's keep it simple.. make a new matrix and copy the elements across, then delete the old one.
        try {   
                        T * new_arr = new T[m_size]; // could throw std::bad_alloc
                        for(unsigned int i=0;i<m_height;++i)
                                for(unsigned int j=0;j<m_width;++j) {
                                        std::swap( m_data[j+(m_width*i)], new_arr[i+(m_height*j)] );
                                        }
                        std::swap(m_width,m_height);
                        std::swap(new_arr,m_data);
                        delete [] new_arr;
                } catch (std::bad_alloc &) { throw E_CjMatrix_MemoryAllocFail(); }
}

}

// template <typename T> CjMatrixImpl_1D<T>::~CjMatrixImpl_1D() { if(m_data!=(T *)0) { delete [] m_data; m_data = (T *)0; } }
template <typename T> CjMatrixImpl_1D<T>::~CjMatrixImpl_1D() { delete [] m_data; m_data = (T *)0; }

template <typename T> bool CjMatrixImpl_1D<T>::isSymmetricImpl() const {
        if(m_width != m_height) return false;
        for(unsigned int i = 1; i < m_width; ++i)
                for(unsigned int j = i+1; j < m_height; ++j)
                        if(getImpl(i,j) != getImpl(j,i)) return false;
        return true;
}

// free functions

template <typename T> CjMatrixImpl_1D <T> * makeCjMatrixImpl_1D(const char * const inputFileName, const char * const delim = ", \t", unsigned int max_width = 0, unsigned int max_height = 0, T * default_value = 0 ) {
std::cout << "makeCjMatrixImpl_1D(const char *, const char *) called." << std::endl;
std::string inputText;  // lines are read from file, one at a time into this string
std::string buf;
typedef std::vector<T> vectorT_t;
typedef std::vector< vectorT_t > vec_VectorT_t;
vec_VectorT_t rows;     // holds pointers to the vectors containing converted input data

T val;

int lineNo = 0;
int ncols = -1;
unsigned int j;

        std::ifstream inFile(inputFileName);

        if( !inFile ) {
                throw E_CjMatrix_FileError(std::string("makeCjMatrixImpl_1D(const char *, const char *): Couldn't open file for reading: ").append(inputFileName) );
                }
        std::vector <T> v;
        while( inFile.good() ) {
                std::getline(inFile, inputText);
                lineNo++;
                std::vector<std::string> tokens;
                boost::char_separator<char> sep (delim);
                boost::tokenizer< boost::char_separator<char> > tok(inputText, sep);
                std::copy( tok.begin(), tok.end(), std::back_inserter(tokens));

                if(tokens.size()==0) {
//                      std::cout << "Skipping blank line number " << lineNo << std::endl;
                        continue;       // empty line - skip it
                        }
                if(tokens.front().compare("#") == 0) {
//                      std::cout << "Skipping comment line " << lineNo << std::endl;
                        continue;       // comment line - skip it
                        }
                j = 1;
                v.clear();
                for(std::vector<std::string>::iterator it = tokens.begin(); it != tokens.end(); ++it, ++j) {
                        std::istringstream ss(*it);
                        ss >> std::noskipws >> val;
                        if ( (!ss) && (!ss.eof()) ) {
//                              std::cout << "Could not parse entry " << j << " on line " << lineNo << std::endl;
                                std::ostringstream oss;
                                oss << "makeCjMatrixImpl_1D(const char *, const char *): Invalid input in entry " << j << " on line " << lineNo << " in file " << inputFileName;
                                throw E_CjMatrix_FileError(oss.str());
                                }
                        v.push_back(val);      // val is cloned into the back of words
//                      std::cout << "Just pushed " << val << std::endl;
                        }

//              std::cout << " just read " << v.size() << " values from line '" << inputText << "'" << std::endl;
                if( ncols == -1 ) ncols = v.size();     // set the number of nodes
                else if( (int)v.size() != ncols) {
                        std::ostringstream oss;
                        oss << "makeCjMatrixImpl_1D(const char *, const char *): Incorrect number of entries in line " << lineNo << " - was expecting " << ncols << " but got " << v.size();
                        throw E_CjMatrix_FileError(oss.str());
                        }
                rows.push_back(v);
        }
        std::cout << "Just read " << ncols << " by " << rows.size() << " from file " << inputFileName << std::endl;

        return new CjMatrixImpl_1D<T>(rows, ncols, rows.size(), (max_width!=0)?max_width:(unsigned int)ncols, (max_height!=0)?max_height:rows.size(), default_value);
}

#endif // _CJMATRIXIMPL_1D_H_