vectors.cpp

Gaussian Mixture Algorithm

  int n = npairs;
  Pair *p = new Pair[n+1];
  for (int i=0; i < n; i++)
    p[i] = pairs[i];
  p[n].i = -1;
  Rep *nr = new Rep;
  delete [] nr->pairs;
  nr->pairs = p;
  nr->size = size;
  nr->npairs = nr->mpairs = n;
  return nr;
}


inline double
SVector::Rep::qset(int i, double v)
{
  assert(i >= size);
  if (npairs >= mpairs)
    {
      resize(max(16, mpairs+mpairs));
      assert(npairs < mpairs);
    }
  Pair *p = &pairs[npairs++];
  size = i + 1;
  p->i = i;
  p->v = v;
  return v;
}


SVector::SVector()
{
  trim();
}


SVector::SVector(const FVector &v)
{
  int m = v.size();
  const VFloat *f = v;
  Rep *r = rep();
  r->resize(m);
  for (int i=0; i<m; i++)
    if (f[i] != 0)
      r->qset(i,f[i]);
  trim();
}


namespace {
  
  SVector::Pair *
  search(SVector::Pair *pairs, int npairs, int i)
  {
    int lo = 0;
    int hi = npairs - 1;
    while (lo <= hi)
      {
        int d = (lo + hi + 1) / 2;
        if (i == pairs[d].i)
          return &pairs[d];
        else if (i < pairs[d].i)
          hi = d-1;
        else
          lo = d+1;
      }
    return 0;
  }

}


double
SVector::get(int i) const
{
  const Rep *r = rep();
  if (i >= r->size)
    return 0;
  // binary search
  assert(i >= 0);
  SVector::Pair *pair = search(r->pairs, r->npairs, i);
  if (pair)
    return pair->v;
  return 0;
}


double 
SVector::set(int i, double v)
{
  w.detach();
  Rep *r = rep();
  if (v)
    {
      if (i >= r->size)
        return r->qset(i, v);
      assert(i >= 0);
      SVector::Pair *p = search(r->pairs, r->npairs, i);
      if (p)
        return p->v = v;
      if (r->npairs >= r->mpairs)
        {
          r->resize(max(16, r->mpairs + min(r->mpairs, 4096)));
          assert(r->npairs < r->mpairs);
        }
      SVector::Pair *s = r->pairs;
      p = s + r->npairs;
      r->npairs += 1;
      for (; p > s && p[-1].i > i;  p--)
        p[0] = p[-1];
      p[0].i = i;
      p[0].v = v;
    }
  else
    {
      SVector::Pair *s = r->pairs;
      SVector::Pair *p = search(s, r->npairs, i);
      if (p)
        {
          r->npairs -= 1;
          for (; p->i >= 0; p++)
            p[0] = p[1];
        }
    }
  return v;
}


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

void 
SVector::trim()
{
  w.detach();
  Rep *r = rep();
  r->resize(r->npairs);
}


SVector
SVector::slice(int fi, int ti) const
{
  assert(ti >= 0);
  assert(ti >= fi);
  SVector y;
  for(Pair *p = rep()->pairs; p->i >= 0 && p->i <= ti; p++)
    if (p->i >= fi)
      y.set(p->i, p->v);
  return y;
}


void 
SVector::add(const SVector &v2)
{
  operator=( ::combine(*this, 1, v2, 1) );
}


void 
SVector::add(const SVector &v2, double c2)
{
  operator=( ::combine(*this, 1, v2, c2) );
}


void 
SVector::combine(double c1, const SVector &v2, double c2)
{
  operator=( ::combine(*this, c1, v2, c2) );
}


void 
SVector::scale(double c1)
{
  if (c1)
    {
      w.detach();
      Rep *r = rep();
      Pair *pairs = r->pairs;
      int npairs = r->npairs;
      for (int i=0; i<npairs; i++)
        pairs[i].v *= c1;
    }
  else
    {
      clear();
    }
}


std::ostream& 
operator<<(std::ostream &f, const SVector &v)
{
  const SVector::Rep *r = v.rep();
  const SVector::Pair *pairs = r->pairs;
  int npairs = r->npairs;
  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<npairs; i++)
    {
      f << " " << pairs[i].i << ":";
      VFloat v = pairs[i].v;
      short iv = (int)v;
      if (v == (VFloat)iv)
        f << iv;
      else
        f << v;
    }
  f << std::endl;
  f.precision(oldprec);
  f.flags(oldflags);
  return f;
}


std::istream& 
operator>>(std::istream &f, SVector &v)
{
  v.clear();
  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;
      v.set(i,x);
    }
  v.trim();
  return f;
}


bool 
SVector::save(std::ostream &f) const
{
  const Rep *r = rep();
  const Pair *pairs = r->pairs;
  int npairs = r->npairs;
  f.write((const char*)&npairs, sizeof(int));
  f.write((const char*)pairs, sizeof(Pair)*npairs);
  return f.good();
}


bool 
SVector::load(std::istream &f)
{
  clear();
  int npairs = 0;
  f.read((char*)&npairs, sizeof(int));
  if (npairs < 0)
    f.setstate(std::ios::badbit);
  if (!f.good())
    return false;
  rep()->resize(npairs);
  for (int i=0; i<npairs; i++)
    {
      Pair pair;
      f.read((char*)&pair, sizeof(Pair));
      if (f.good())
        set(pair.i, pair.v);
    }
  trim();
  return f.good();
}


double 
dot(const FVector &v1, const FVector &v2)
{
  int m = min(v1.size(), v2.size());
  const VFloat *f1 = v1;
  const VFloat *f2 = v2;
  double sum = 0.0;
  while (--m >= 0)
    sum += (double)(*f1++) * (double)(*f2++);
  return sum;
}


double 
dot(const FVector &v1, const SVector &v2)
{
  int m = v1.size();
  const VFloat *f = v1;
  const SVector::Pair *p = v2;
  double sum = 0;
  if (p)
    for( ; p->i >= 0 && p->i < m; p++)
      sum += (double)p->v * (double)f[p->i];
  return sum;
}


double 
dot(const SVector &v1, const FVector &v2)
{
  int m = v2.size();
  const VFloat *f = v2;
  const SVector::Pair *p = v1;
  double sum = 0;
  if (p)
    for( ; p->i >= 0 && p->i < m; p++)
      sum += (double)p->v * (double)f[p->i];
  return sum;
}


double 
dot(const SVector &v1, const SVector &v2)
{
  const SVector::Pair *p1 = v1;
  const SVector::Pair *p2 = v2;
  double sum = 0;
  if (p1 && p2)
    while (p1->i >= 0 && p2->i >= 0)
      {
        if (p1->i < p2->i)
          p1++;
        else if (p1->i > p2->i)
          p2++;
        else
          sum += (double)(p1++)->v * (double)(p2++)->v;
      }
  return sum;
}


SVector 
combine(const SVector &v1, double a1, const SVector &v2, double a2)
{
  const SVector::Pair *p1 = v1;
  const SVector::Pair *p2 = v2;
  SVector ans;
  SVector::Rep *r = ans.rep();
  r->resize(v1.npairs() + v2.npairs());
  SVector::Pair *p = r->pairs;
  while (p1->i >= 0 && p2->i >= 0)
     {
       if (p1->i < p2->i)
         {
           double v = (double)p1->v * a1;
           if (v)
             {
               p->i = p1->i;
               p->v = v;
               p++;
             }
           p1++;
         }
       else if (p1->i > p2->i)
         {
           double v = (double)p2->v * a2;
           if (v)
             {
               p->i = p2->i;
               p->v = v;
               p++;
             }
           p2++;
         }
       else
         {
           double v = (double)p1->v * a1 + (double)p2->v * a2;
           if (v)
             {
               p->i = p1->i;
               p->v = v;
               p++;
             }
           p1++;
           p2++;
         }
     }
  while  (p1->i >= 0)
    {
      double v = (double)p1->v * a1;
      if (v)
        {
          p->i = p1->i;
          p->v = v;
          p++;
        }
      p1++;
    }
  while  (p2->i >= 0)
    {
      double v = (double)p2->v * a2;
      if (v)
        {
          p->i = p2->i;
          p->v = v;
          p++;
        }
      p2++;
    }
  r->npairs = p - r->pairs;
  r->size = (r->npairs > 0) ? p[-1].i + 1 : 0;
  ans.trim();
  return ans;
}


FVector 
combine(const FVector &v1, double a1, const SVector &v2, double a2)
{
  FVector r = v1;
  r.combine(a1, v2, a2);
  return r;
}


FVector 
combine(const SVector &v1, double a1, const FVector &v2, double a2)
{
  FVector r = v2;
  r.combine(a2, v1, a1);
  return r;
}


FVector 
combine(const FVector &v1, double a1, const FVector &v2, double a2)
{
  FVector r = v1;
  r.combine(a1, v2, a2);
  return r;
}



/* -------------------------------------------------------------
   Local Variables:
   c++-font-lock-extra-types: ( "\\sw+_t" "[A-Z]\\sw*[a-z]\\sw*" "std::\\sw+")
   End:
   ------------------------------------------------------------- */