vectors.cpp

Gaussian Mixture Algorithm

// -*- C++ -*-
// Little library of vectors and sparse vectors
// Copyright (C) 2007- Leon Bottou

// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// 
// This program 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, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111, USA



// $Id: vectors.cpp,v 1.1 2008/10/27 23:26:31 cp71 Exp $


#include "vectors.h"
#include <ios>
#include <iomanip>
#include <iostream>
#include <cassert>
#include <cctype>


namespace
{
  template <typename T> inline T min(T a, T b) {
    return (a < b) ? a  : b;
  }
  template <typename T> inline T max(T a, T b) {
    return (a < b) ? b  : a;
  }
}



void 
FVector::Rep::resize(int n)
{
  if (size != n)
    {
      if (n > 0)
        {
          VFloat *newdata = new VFloat[n];
          int i = 0;
          int s = min(size,n);
          for (; i<s; i++)
            newdata[i] = data[i];
          for (; i<n; i++)
            newdata[i] = 0;
          delete [] data;
          data = newdata;
          size = n;
        }
      else
        {
          delete [] data;
          data = 0;
          size = 0;
        }
    }
}


FVector::Rep *
FVector::Rep::copy()
{
  Rep *newrep = new Rep;
  newrep->resize(size);
  for (int i=0; i<size; i++)
    newrep->data[i] = data[i];
  return newrep;
}


FVector::FVector()
{
}


FVector::FVector(int n)
{
  Rep *r = rep();
  r->resize(n);
}


FVector::FVector(const SVector &v)
{
  Rep *r = rep();
  r->resize(v.size());
  int npairs = v.npairs();
  const SVector::Pair *pairs = v;
  for (int i=0; i<npairs; i++, pairs++)
    r->data[pairs->i] = pairs->v;
}


double
FVector::get(int i) const
{
  const Rep *r = rep();
  if (i >=0 && i<r->size)
    return r->data[i];
  assert(i>=0);
  return 0;
}


double 
FVector::set(int i, double v)
{
  w.detach();
  Rep *r = rep();
  if (i >= r->size)
    r->resize(i+1);
  assert(i>=0);
  r->data[i] = v;
  return v;
}


void 
FVector::clear()
{
  w.detach();
  rep()->resize(0);
}


void 
FVector::resize(int n)
{
  w.detach();
  assert(n >= 0);
  rep()->resize(n);
}


void 
FVector::touch(int i)
{
  if (i >= rep()->size)
    resize(i+1);
}


FVector
FVector::slice(int fi, int ti) const
{
  assert(ti >= 0);
  assert(ti >= fi);
  FVector y;
  int s = size();
  if (s > 0)
    {
      fi = min(fi, s-1);
      ti = max(ti, s-1);
      int n = ti - fi + 1;
      y.resize(n);
      VFloat *yp = y.rep()->data;
      VFloat *xp = rep()->data + fi;
      for (int i=0; i<n; i++)
        yp[i] = xp[i];
    }
  return y;
}


void 
FVector::add(double c1)
{
  w.detach();
  Rep *r = rep();
  VFloat *d = r->data;
  for (int i=0; i<r->size; i++)
    d[i] += c1;
}


void 
FVector::add(const FVector &v2)
{
  w.detach();
  Rep *r = rep();
  int m = max(r->size, v2.size());
  if (m > r->size)
    r->resize(m);
  VFloat *d = r->data;
  const VFloat *s = (const VFloat*) v2;
  for (int i=0; i<m; i++)
    d[i] += s[i];
}


void 
FVector::add(const SVector &v2)
{
  w.detach();
  Rep *r = rep();
  int m = max(r->size, v2.size());
  if (m > r->size)
    r->resize(m);
  VFloat *d = r->data;
  int npairs = v2.npairs();
  const SVector::Pair *pairs = v2;
  for (int i=0; i<npairs; i++, pairs++)
    d[pairs->i] += pairs->v;
}


void 
FVector::add(const FVector &v2, double c2)
{
  w.detach();
  Rep *r = rep();
  int m = max(r->size, v2.size());
  if (m > r->size)
    r->resize(m);
  VFloat *d = r->data;
  const VFloat *s = (const VFloat*) v2;
  for (int i=0; i<m; i++)
    d[i] += s[i] * c2;
}


void 
FVector::add(const SVector &v2, double c2)
{
  w.detach();
  Rep *r = rep();
  int m = max(r->size, v2.size());
  if (m > r->size)
    r->resize(m);
  VFloat *d = r->data;
  int npairs = v2.npairs();
  const SVector::Pair *pairs = v2;
  for (int i=0; i<npairs; i++, pairs++)
    d[pairs->i] += pairs->v * c2;
}


void 
FVector::add(const FVector &v2, double c2, const FVector &q2)
{
  w.detach();
  Rep *r = rep();
  int m = r->size;
  m = max(m, v2.size());
  m = min(m, q2.size());
  if (m > r->size)
    r->resize(m);
  VFloat *d = r->data;
  const VFloat *s = (const VFloat*) v2;
  const VFloat *q = (const VFloat*) q2;
  for (int i=0; i<m; i++)
    d[i] += s[i] * q[i] * c2;
}


void 
FVector::add(const SVector &v2, double c2, const FVector &q2)
{
  w.detach();
  Rep *r = rep();
  int m = r->size;
  m = max(m, v2.size());
  m = min(m, q2.size());
  if (m > r->size)
    r->resize(m);
  VFloat *d = r->data;
  const VFloat *q = (const VFloat*) q2;
  int npairs = v2.npairs();
  const SVector::Pair *pairs = v2;
  for (int i=0; i<npairs && pairs->i < m; i++, pairs++)
    d[pairs->i] += pairs->v * q[pairs->i] * c2;
}


void 
FVector::scale(double c1)
{
  w.detach();
  Rep *r = rep();
  VFloat *d = r->data;
  for (int i=0; i<r->size; i++)
    d[i] *= c1;
}


void 
FVector::combine(double c1, const FVector &v2, double c2)
{
  w.detach();
  Rep *r = rep();
  int m = max(r->size, v2.size());
  if (m > r->size)
    r->resize(m);
  VFloat *d = r->data;
  const VFloat *s = (const VFloat*) v2;
  for (int i=0; i<m; i++)
    d[i] = d[i] * c1 + s[i] * c2;
}


void 
FVector::combine(double c1, const SVector &v2, double c2)
{
  w.detach();
  Rep *r = rep();
  int m = max(r->size, v2.size());
  if (m > r->size)
    r->resize(m);
  VFloat *d = r->data;
  int npairs = v2.npairs();
  const SVector::Pair *pairs = v2;
  int j = 0;
  for (int i=0; i<npairs; i++, pairs++)
    {
      if (c1 != 1)
        for (; j < pairs->i; j++)
          d[j] = d[j] * c1;
      j = pairs->i;
      d[j] = d[j] * c1 + pairs->v * c2;
      j++;
    }
  if (c1 != 1)
    for (; j < pairs->i; j++)
      d[j] = d[j] * c1;
}


std::ostream& 
operator<<(std::ostream &f, const FVector &v)
{
  std::ios::fmtflags oldflags = f.flags();
  std::streamsize oldprec = f.precision();
  f << std::scientific << std::setprecision(sizeof(VFloat)==4 ? 7 : 16);
  for (int i=0; i<v.size(); i++)
    if (v[i] || (i+1 == v.size()))
      {
        f << " " << i << ":";
        VFloat x = v[i];
        short ix = (int)x;
        if (x == (VFloat)ix)
          f << ix;
        else
          f << x;
      }
  f << std::endl;
  f.precision(oldprec);
  f.flags(oldflags);
  return f;
}


std::istream& 
operator>>(std::istream &f, FVector &v)
{
  int sz = 0;
  int msz = 1024;
  v.clear();
  v.resize(msz);
  for(;;)
    {
      int c = f.get();
      if (!f.good() || (c=='\n' || c=='\r'))
        break;
      if (::isspace(c))
        continue;
      int i;
      f.unget();
      f >> std::skipws >> i >> std::ws;
      if (f.get() != ':')
        {
          f.unget();
          f.setstate(std::ios::badbit);
          break;
        }
      double x;
      f >> std::skipws >> x;
      if (!f.good())
        break;
      if (i >= sz)
        sz = i + 1;
      if (i >= msz)
        {
          while(i >= msz)
            msz += msz;
          v.resize(msz);
        }
      v.set(i,x);
    }
  v.resize(sz);
  return f;
}

bool 
FVector::save(std::ostream &f) const
{
  int i = size();
  const VFloat *d = rep()->data;
  f.write((const char*)&i, sizeof(int));
  f.write((const char*)d, sizeof(VFloat)*i);
  return f.good();
}


bool
FVector::load(std::istream &f)
{
  int i = 0;
  clear();
  f.read((char*)&i, sizeof(int));
  if (i<0)
    f.setstate(std::ios::badbit);
  if (!f.good())
    return false;
  resize(i);
  VFloat *d = rep()->data;
  f.read((char*)d, sizeof(VFloat)*i);
  return f.good();
}



// ----------------------------------------


void
SVector::Rep::resize(int n)
{
  if (n != mpairs)
    {
      Pair *p = new Pair[n+1];
      int m = min(n, npairs);
      int i = 0;
      for (; i < m; i++)
        p[i] = pairs[i];
      for (; i <= n; i++)
        p[i].i = -1;
      delete [] pairs;
      pairs = p;
      npairs = m;
      mpairs = n;
      size = (m>0) ? p[m-1].i + 1 : 0;
    }
}


SVector::Rep *
SVector::Rep::copy()
{