---

mmatrix.cc

/***************************************************************************
 *   This file is part of the MagiC++ library.                             *
 *                                                                         *
 *   Copyright (C) 1998-2001 Marko Grönroos <magi@iki.fi>                  *
 *                                                                         *
 ***************************************************************************
 *                                                                         *
 *  This library is free software; you can redistribute it and/or          *
 *  modify it under the terms of the GNU Library General Public            *
 *  License as published by the Free Software Foundation; either           *
 *  version 2 of the License, or (at your option) any later version.       *
 *                                                                         *
 *  This library 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      *
 *  Library General Public License for more details.                       *
 *                                                                         *
 *  You should have received a copy of the GNU Library General Public      *
 *  License along with this library; see the file COPYING.LIB.  If         *
 *  not, write to the Free Software Foundation, Inc., 59 Temple Place      *
 *  - Suite 330, Boston, MA 02111-1307, USA.                               *
 *                                                                         *
 ***************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include "magic/mmath.h"
#include "magic/mmatrix.h"
#include "magic/mstream.h"
#include "magic/mlist.h"

void Matrix::load (const String& filename) {
    FILE* in = fopen (filename, "r");
    if (!in)
        throw file_not_found (format("Matrix file '%s' not found", (CONSTR)filename));
    load (in);
    fclose (in);
}

void Matrix::load (FILE* in) {
    List<Vector*> rows;     // Will contain the rows in the input file
    int rowCount=0;         // Number of rows in the input file

    // Read the datafile
    String buffer;
    String item;        // One numeric or text item
    item.reserve (20);

    // Read one line at a time
    int rowCnt = 0;
    while (fgetS (in, buffer) != -1) {
        buffer.chop(); // Remove trailing newline
        Vector* prow = new Vector (cols);
        int itemIndex=0;    // Index of the above item within the row
        for (int i=0; i<buffer.length()+1; i++) {
            if (i==buffer.length() || buffer[i]==' ' || buffer[i]=='\t') {
                if (itemIndex<cols)
                    (*prow)[itemIndex++] = (item=="x")? UNDEFINED_FLOAT : item.toDouble();
                item = "";
            } else
                item += buffer[i];
        }
        rows.add (prow);
        rowCnt++;
    }

    // Adjust our size
    make (rowCnt,cols);

    // Put the pattern list into the set
    int r=0;
    for (ListIter<Vector*> l (rows); !l.exhausted(); l.next(), r++) {
        for (int c=0; c<cols; c++)
            get(r,c) = (*l.get())[c];
    }
}

void Matrix::save (FILE* out) const {
    for (int r=0; r<rows; r++) {
        for (int c=0; c<cols; c++)
            fprintf (out, "%g ", get(r,c));
        fprintf (out, "\n");
    }
}

const Matrix& Matrix::transpose () {
    for (int i=0; i<rows; i++)
        for (int j=i; j<cols; j++)
            swap (get(i,j), get(j,i));
    return *this;
}
 
const Matrix& Matrix::multiplyToSum (double s) {
    double cs = sum ();
    ASSERTWITH (fabs(cs)>1E-10, "Tried to normalize a zero matrix");

    operator*= (1.0/cs);
    return *this;
}

double Matrix::sum () const {
    int size=rows*cols;
    double res=0.0;
    for (int i=0; i<size; i++)
        res += data[i];
    return res;
}

const Matrix& Matrix::operator= (double x) {
    int size=rows*cols;
    for (int i=0; i<size; i++)
        data[i]= x;
    return *this;
}

const Matrix& Matrix::operator= (const Matrix& orig) {
    make (orig.rows, orig.cols);
    int size=rows*cols;
    for (int i=0; i<size; i++)
        data[i]= orig.data[i];
    return *this;
}

const Matrix& Matrix::operator+= (const Matrix& other) {
    int size=rows*cols;
    for (int i=0; i<size; i++)
        data[i] += other.data[i];
    return *this;
}

const Matrix& Matrix::operator+= (double x) {
    int size=rows*cols;
    for (int i=0; i<size; i++)
        data[i] += x;
    return *this;
}

const Matrix& Matrix::operator*= (const Matrix& other) {
    int size=rows*cols;
    for (int i=0; i<size; i++)
        data[i] *= other.data[i];
    return *this;
}

const Matrix& Matrix::operator*= (double x) {
    int size=rows*cols;
    for (int i=0; i<size; i++)
        data[i] *= x;
    return *this;
}

void Matrix::operator>> (OStream& out) const {
    out.printf ("{");
    for (int i=0; i<rows; i++) {
        if (i>0)
            out.printf (" ");
        out.printf ("{");
        for (int j=0; j<cols; j++) {
            out.printf ("%2f", get (i,j));
            if (j<cols-1)
                out.printf (",");
        }
        out.print ("}");
        if (i<rows-1)
            out.printf ("\n");
    }
    out.printf ("}\n");
}

void Matrix::splitVertical (Matrix& a, Matrix& b, int column) const {
    a.make (rows, column);
    b.make (rows, cols-column);
    for (int r=0; r<rows; r++) {
        for (int c=0; c<column; c++)
            a.get(r,c) = get(r,c);
        for (int c=0; c<cols-column; c++)
            b.get(r,c) = get(r,c+column);
    }
}

void Matrix::splitHorizontal (Matrix& a, Matrix& b, int row) const {
    a.make (row, cols);
    b.make (rows-row, cols);
    for (int c=0; c<cols; c++) {
        for (int r=0; r<row; r++)
            a.get(r,c) = get(r,c);
        for (int r=0; r<rows-row; r++)
            b.get(r,c) = get(r+row,c);
    }
}

void Matrix::joinVertical (const Matrix& a, const Matrix& b) {
    ASSERT (a.rows == b.rows);
    make (a.rows, a.cols + b.cols);
    for (int r=0; r<rows; r++) {
        for (int c=0; c<a.cols; c++)
            get(r,c) = a.get(r,c);
        for (int c=0; c<b.cols; c++)
            get(r,c+a.cols) = a.get(r,c);
    }
}

void Matrix::joinHorizontal (const Matrix& a, const Matrix& b) {
    ASSERT (a.cols == b.cols);
    make (a.rows+b.rows, a.cols);
    for (int c=0; c<cols; c++) {
        for (int r=0; r<a.rows; r++)
            get(r,c) = a.get(r,c);
        for (int r=0; r<b.rows; r++)
            get(r+a.rows,c) = a.get(r,c);
    }
}

Matrix Matrix::sub (int row0, int row1, int col0, int col1) const {
    ASSERT (row0>=0 && row0<rows);
    ASSERT ((row1>=-1 && row1<rows) || row1==-1);
    ASSERT (col0>=0 && col0<cols);
    ASSERT ((col1>=-1 && col1<cols) || col1==-1);

    if (row1==-1)
        row1 = rows-1;
    if (col1==-1)
        col1 = cols-1;

    Matrix result;
    result.make (row1-row0+1, col1-col0+1);
    for (int r=0; r<result.rows; r++)
        for (int c=0; c<result.cols; c++)
            result.get(r,c) = (*this).get(r+row0, c+col0);
    return result;
}

void Matrix::insertColumn (int cols) {
    
}


/*
SubMatrix::SubMatrix (Matrix& m, int sr, int sc, int er, int ec) : matrix (m) {
    ASSERT (sr>0 && sc>0 && er>0 && ec>0);
    ASSERT (sr<m.rows);
    ASSERT (sc<m.cols);
    ASSERT (er<m.rows);
    ASSERT (ec<m.cols);
    ASSERT (sr<=er);
    ASSERT (sc<=ec);
    srow=sr, scol=sc, erow=er, ecol=ec;
    rows = erow-srow+1;
    cols = ecol-scol+1;

}

double& SubMatrix::get (int row, int col) {
    return matrix.PackTable<double>::get (row-srow, col-scol);
}
*/

---