[Author Prev][Author Next][Thread Prev][Thread Next][Author Index][Thread Index]

[freehaven-cvs] Not yet useful for anybody else, but heres some code...



Update of /home/freehaven/cvsroot/doc/e2e-traffic/src
In directory moria.mit.edu:/tmp/cvs-serv15356

Added Files:
	.cvsignore Makefile comb.h rng.cpp rng.h sim.cpp sim.h 
	simmain.cpp test.cpp vec.cpp vec.h 
Log Message:
Not yet useful for anybody else, but heres some code to run traffic analysis sims.

--- NEW FILE: .cvsignore ---
sim
test
gmon.out
boost*


--- NEW FILE: Makefile ---
(This appears to be a binary file; contents omitted.)

--- NEW FILE: comb.h ---

#ifndef _COMB_H
#define _COMB_H

#include <assert.h>

inline long fact(int n) {
  assert(n>=0);
  int r = 1;
  for (int i = 2; i <= n; ++i) {
    r *= i;
  }
  return r;  
}

inline long comb(int n, int k) {
  /* COMB(n,k) = n!/(k!(n-k)!) */
  assert(n>=0 && k>=0 && k<=n);
  int t = n-k;
  if (k>t) { int tmp = t; t = k; k = tmp; }
  int r = 1;
  for (int i = t+1; i <= n; ++i) {
    r *= i;
  }
  for (int i = 2; i <= k; ++i) {
    r /= i;
  }
  return r;
}

inline long bins(int n, int k) {
  return comb(n+k, k);
}

#endif

--- NEW FILE: rng.cpp ---

#include <stdlib.h>
#include "rng.h"

#include "boost/random.hpp"

//
// If we're prototyping, use a lagged finonacci generator: it's fastest.
//
#ifdef FAST_RANDOM
#define RNGCLASS boost::lagged_fibonacci19937
#else
#define RNGCLASS boost:mt19937
#endif

class RNGState {
public:
  RNGCLASS _rng;
  boost::uniform_01<RNGCLASS> uniform;
  boost::normal_distribution<RNGCLASS> norm;

  RNGState() : _rng(), uniform(_rng), norm(_rng) {}
};

static RNGState state;

double
rng()
{
  return state.uniform();
}

double
normal_rng()
{
  return state.norm();
}

CumulativeDist::CumulativeDist(const std::vector<double> &dist, 
                               bool isCumulative)
  : prob(dist.size()), cumDist(dist.size())
{
  int sz = dist.size();
  if (isCumulative) {
    cumDist = dist;
    
    for (int i = 0; i < sz-1; ++i) {
      prob[i] = cumDist[i+1] - cumDist[i];
    }
    prob[sz-1] = 1.0 - cumDist[sz-1];
  } else {
    prob = dist;

    double tot = 0.0;
    for (int i = 0; i < sz; ++i) {
      tot += dist[i];
      cumDist[i] = tot;
    }
    for (int i = 0; i < sz; ++i) {
      cumDist[i] /= tot;
    }
  }
}

int
CumulativeDist::get() const
{
  std::vector<double>::const_iterator it = 
    std::lower_bound(cumDist.begin(), cumDist.end(), rng());
  return it - cumDist.begin();
}


--- NEW FILE: rng.h ---

#ifndef _RNG_H
#define _RNG_H

#include <vector>

template <class C>
class Dist
{
 public:
  virtual C get() const = 0;
  virtual double getP(const C &v) const = 0;
  virtual ~Dist() {}
};

template <class C>
class ConstDist : public Dist<C>
{
 private:
  C val;
 public:
  ConstDist(const C &vv) : val(vv) {}
  C get() const { return val; }
  double getP(const C &v) const { return v == val ? 1.0 : 0.0; }
  ~ConstDist() {}
};

class CumulativeDist : public Dist<int>
{
 private:
  std::vector<double> prob;
  std::vector<double> cumDist;
 public:
  CumulativeDist(const std::vector<double> &dist, bool isCumulative=0);
  CumulativeDist(const CumulativeDist &d) : cumDist(d.cumDist) {}
  CumulativeDist & operator=(const CumulativeDist &d)
    { cumDist = d.cumDist; return *this; }
  int get() const;
  double getP(const int &v) const
    { return (0 <= v && v <= (int)prob.size()) ? prob[v] : 0.0; }
  ~CumulativeDist() {}
};

double rng();
double normal_rng();

inline int rng(int m) { return static_cast<int>(rng()*m); }

template<class C>
const C &rng_pick(const std::vector<C> &v) { return v[rng(v.size())]; }
template<class C>
C &rng_pick(std::vector<C> &v) { return v[rng(v.size())]; }

template<class C> void 
rng_shuffle(std::vector<C> &v, int n=-1) {
  int sz = v.size();
  if (n == -1) { n = v.size()-1; }
  for (int i = n; i; ++i) {
    int swap = i+rng(sz-i);
    if (i != swap) {
      C tmp = v[swap];
      v[swap] = v[i];
      v[i] = tmp;
    }
  }
}

#endif // _RNG_H

--- NEW FILE: sim.cpp ---

#include "sim.h"
#include "rng.h"

Alice::~Alice() {}
Background::~Background() {}

// ======================================================================

UniformAlice::UniformAlice(const std::vector<int> &r, 
                           Dist<int> *mD, Dist<int> *dD)
  : recipients(r), msgDist(mD), dummyDist(dD)
{
}

void 
UniformAlice::getTraffic(vec<int> &v_out, int &nOut)
{
  nOut = msgDist->get();
  for (int i = 0; i < nOut; ++i) {
    ++ v_out[rng_pick(recipients)];
  }
  if (dummyDist)
    nOut += dummyDist->get();
}

UniformBackground::UniformBackground(int nR, int nPR)
  : nRecipients(nR), nPerRound(nPR)
{
}

void
UniformBackground::getNTraffic(vec<int> &v_out, int nMessages)
{
  assert(nMessages >= 0);
  while (nMessages--) {
    ++v_out[rng(nRecipients)];
  }
}

void
UniformBackground::getTraffic(vec<int> &v_out)
{
  assert(nPerRound >= 0);
  getNTraffic(v_out, nPerRound);
}

BatchMix::BatchMix(int b)
  : Mixnet(), batchSize(b)
{
}

void
BatchMix::reset()
{
}

void
BatchMix::addRound(const vec<int> &alice,
                   int nAlice, const vec<int> &background,
                   FwdAttacker *a)
{
  assert(alice.size() == background.size());
  int nOther = background.total(0);
  vec<int> total(alice);
  total += background;
  a->addRound(nAlice, nOther, total);
}

SDAttacker::SDAttacker(vec<double> &bg)
  : nRounds(0), nAlice(0), nOther(0), background(bg), observed(bg.size(), 0)
{
}

void
SDAttacker::reset()
{
  nRounds = 0;
  nAlice = 0;
  nOther = 0;
  observed.reset(0);
}

void
SDAttacker::addRound(int nA, int nO, const vec<int> &rcvd)
{
  ++nRounds;
  if (nA) {
    nAlice += nA;
    nOther += nO;
    observed += rcvd;
  }
}

void
SDAttacker::guessAlice(vec<double> &r)
{
  r = vec<double>(observed);
  //cerr << "1." << r << endl;
  //cerr << "1b." << background << endl;
  //  cerr << "1c. " << nAlice << ", " << nOther << endl;
  r -= background*nOther;
  //cerr << "2." << r << endl;
  r /= nAlice;
  //cerr << "3." << r << endl;
}

SDTrial::SDTrial(int nR, int nAR, int b, int g)
{
  aMsgs = new ConstDist<int>(1);
  aDummies = new ConstDist<int>(0);

  truth = std::vector<int>(nAR, 0);
  for (int i=0; i<nAR; ++i) { truth[i] = i; }
  
  alice = new UniformAlice(truth, aMsgs, aDummies);
  background = new UniformBackground(nR);
  vec<double> uniform(nR, 1.0/b);
  attacker = new SDAttacker(uniform);
  mixnet = new BatchMix(b);
  nRecips = nR;
  nBatch = b;
  granularity = g;
}

SDTrial::~SDTrial()
{
  delete mixnet;
  delete attacker;
  delete background;
  delete alice;
  delete aMsgs;
  delete aDummies;
}

int
SDTrial::attempt()
{
  mixnet->reset();
  attacker->reset();
  
  vec<int> aTraffic(nRecips, 0);
  vec<int> bTraffic(nRecips, 0);
  vec<double> recips(nRecips, 0.0);

  int n = 0;
  while (1) {
    aTraffic.reset(0);
    bTraffic.reset(0);
    int nAlice;
    alice->getTraffic(aTraffic, nAlice);
    background->getNTraffic(bTraffic, nBatch - aTraffic.total(0));
    //cout << n << " "<<bTraffic << endl;
    mixnet->addRound(aTraffic, nAlice, bTraffic, attacker);
    if (!(n % granularity)) {
      attacker->guessAlice(recips);
      // cout << n << " " << recips << endl;
      //cout << n << " ------ " << endl;
      //pvec(cout, truth);
      //pvec(cout, recips.topN(truth.size()));
      //cout << n << " ------ " << endl;
      if (truth == recips.topN(truth.size()))
	return n;
    }
    ++n;
  }
}

// ======================================================================
// Unknown background

void
DistBackground::getNTraffic(vec<int> &v_out, int nMessages) 
{
  while (nMessages-- > 0) {
    ++ v_out[ cdist.get() ];
  }
}

void 
DistBackground::getTraffic(vec<int> &v_out)
{
  getNTraffic(v_out, nMessages->get());
}

UnkBGBatchAttacker::UnkBGBatchAttacker(int nR)
  : vObservations(nR, 0), uObservations(nR, 0)
{
  reset();
}

void
UnkBGBatchAttacker::reset()
{
  nAlice = nOther = nBg = 0;
  vObservations.reset(0);
  uObservations.reset(0);
}

void
UnkBGBatchAttacker::addRound(int nA, int nO, const vec<int> &nReceived)
{
  if (nA) {
    nAlice += nA;
    nOther += nO;
    vObservations += nReceived;
  } else {
    nBg += nO;
    uObservations += nReceived;
  }
}

void 
UnkBGBatchAttacker::guessAlice(vec<double> &r)
{
  vec<double> u(uObservations);
  u *= (nOther / nBg);
  r = vec<double>(vObservations);
  r -= u;
  r /= nAlice;
}

// ======================================================================

DelayMix::DelayMix(int nR, int mD, 
                   Dist<int> *d)
  : maxDelay(mD), poolIdx(0), pools(mD), delayDist(d)
{
  for (int i = 0; i < mD; ++i) {
    pools[i] = new vec<int>(nR, 0);
  }
}

void
DelayMix::reset()
{
  poolIdx = 0;
  for (int i = 0; i < maxDelay; ++i) {
    pools[i]->reset(0);
  }
}

DelayMix::~DelayMix()
{
  for (int i = 0; i < maxDelay; ++i) {
    delete pools[i];
  }
}

void
DelayMix::addRound(const vec<int> &alice,
                   int nAlice,
                   const vec<int> &background,
                   FwdAttacker *a) 
{
  int nOther = background.total(0);
  vec<int> total(alice);
  total += background;
  int sz = alice.size();
  for (int i = 0; i < sz; ++i) {
    for (int j = 0; j < total[i]; ++j) {
      ++ pools[(poolIdx+getDelay())%maxDelay];
    }
  }
  a->addRound(nAlice, nOther, *pools[poolIdx]);
  pools[poolIdx]->reset(0);
  ++poolIdx;
  if (poolIdx >= maxDelay)
    poolIdx = 0;
}

DelayMixAttacker::DelayMixAttacker(int nR, int mD,  Dist<int> *dDist)
  : nRecips(nR), maxDelay(mD), nAliceIdx(0), nAliceHist(mD, 0), 
    nOtherHist(mD, 0),
    background(nR, 0.0),
    observed(nR, 0.0), 
    nObservedOther(0.0), nObservedAlice(0.0),
    delayDist(dDist)
{
}

double
DelayMixAttacker::expectedAliceMsgs()
{
  double tot = 0.0;
  for (int i = 0; i < maxDelay; ++i) {
    tot += delayDist->getP(i)*aHist(i);
  }
  return tot;
}

double
DelayMixAttacker::expectedOtherMsgs()
{
  double tot = 0.0;
  for (int i = 0; i < maxDelay; ++i) {
    tot += delayDist->getP(i)*oHist(i);
  }
  return tot;
}

void
DelayMixAttacker::reset()
{
  nAliceIdx = 0;
  nAliceHist = std::vector<int>(nRecips, 0);
  nOtherHist = std::vector<int>(nRecips, 0);  
  background.reset(0.0);
  observed.reset(0.0);
  nObservedOther = nObservedAlice = 0.0;
}

#define BG_THRESHOLD 0.05
void 
DelayMixAttacker::addRound(int nAlice, int nOther, 
                           const vec<int> &r)
{
  nAliceHist[nAliceIdx] = nAlice;
  nOtherHist[nAliceIdx] = nOther;
  
  double exAlice = expectedAliceMsgs();
  double exOther = expectedAliceMsgs();
  if (exAlice < BG_THRESHOLD) {
    for (int i = 0; i < nRecips; ++i)
      background[i] += r[i]*exOther;
  } else {
    for (int i = 0; i < nRecips; ++i)
      observed[i] += r[i]*exAlice;
    nObservedOther += exOther;
    nObservedAlice += exAlice;
  }

  ++nAliceIdx;
  if (nAliceIdx >= maxDelay)
    nAliceIdx = 0;
}

void
DelayMixAttacker::guessAlice(vec<double> &res)
{
  vec<double> u(background);
  u *= (nObservedOther / u.total(0.0));
  res = observed - u;
  res -= u;
  res /= nObservedAlice;
}

--- NEW FILE: sim.h ---

#ifndef _SIM_H
#define _SIM_H

#include "vec.h"
#include "rng.h"

class Alice {
public:
  Alice() {}
  virtual void getTraffic(vec<int> &v_out, int &nOut) = 0;
  virtual ~Alice();
};

class Background {
public:
  Background() {}
  virtual void getNTraffic(vec<int> &v_out, int nMessages) = 0;
  virtual void getTraffic(vec<int> &v_out) = 0;
  virtual ~Background();
};

class FwdAttacker {
 public:
  FwdAttacker() {}

  virtual void reset() = 0;  
  virtual void addRound(int nSentAlice, int nSentOther, 
                        const vec<int> &nReceived) = 0;
  virtual void guessAlice(vec<double> &recipients) = 0;

  virtual ~FwdAttacker() {}
};

class Mixnet {
public:
  Mixnet() {}
  virtual void reset() = 0;
  virtual void addRound(const vec<int> &alice,
                        int nAlice,
                        const vec<int> &background,
                        FwdAttacker *r) = 0;
  virtual ~Mixnet() {}
};

class Trial {
public:
  Trial() {}
  virtual int attempt() = 0;
  virtual ~Trial() {}
};

// ======================================================================
// Original statistical disclosure attack

class UniformAlice : public Alice {
private:
  std::vector<int> recipients;
  Dist<int> *msgDist;
  Dist<int> *dummyDist;
public:
  UniformAlice(const std::vector<int> &r,
               Dist<int> *msgDist, Dist<int> *dummyDist);
  void getTraffic(vec<int> &v_out, int &nOut);
  ~UniformAlice() {}
};

class UniformBackground : public Background {
private:
  int nRecipients;
  int nPerRound;
public:
  UniformBackground(int nR, int nPR=-1);
  void getNTraffic(vec<int> &v_out, int nMessages);
  void getTraffic(vec<int> &v_out);
  ~UniformBackground() {}
};

class BatchMix : public Mixnet {
private:
  int batchSize;
public:
  BatchMix(int b);
  void reset();
  void addRound(const vec<int> &alice, 
                int nAlice, 
                const vec<int> &background,
                FwdAttacker *f);
  ~BatchMix() {}
};

class SDAttacker : public FwdAttacker {
private:
  int nRounds;
  int nAlice;
  int nOther;
  vec<double> background;
  vec<int> observed;
public:
  SDAttacker(vec<double> &background);
  void reset();
  void addRound(int nSentAlice, int nSentOther,
                const vec<int> &nReceived);
  void guessAlice(vec<double> &recipients);
  ~SDAttacker() {}
};


// Simple statistical disclosure attack.
class SDTrial : public Trial {
protected:
  Dist<int> *aMsgs;
  Dist<int> *aDummies;
  Alice *alice;
  Background *background;
  BatchMix *mixnet;
  FwdAttacker *attacker;
  std::vector<int> truth;
  int nRecips;
  int nBatch;
  int granularity;
public:
  SDTrial(int nRecipients, int nAliceRecipients, int batchSize, 
	  int granularity=5);
  int attempt();
  ~SDTrial();
};

// ======================================================================
// Attack with unknown background

class DistBackground : public Background {
 private:
  std::vector<double> dist;
  CumulativeDist cdist;
  Dist<int> *nMessages;
 public:
  DistBackground(const std::vector<double> &dist, Dist<int> *nMsgs)
    : dist(dist), cdist(dist), nMessages(nMsgs) {}
  void getNTraffic(vec<int> &v_out, int nMessages);
  void getTraffic(vec<int> &v_out);
  ~DistBackground() {}
};

class UnkBGBatchAttacker : public FwdAttacker {
 private:
  int nAlice;
  int nOther;
  vec<int> vObservations;
  int nBg;
  vec<int> uObservations;
 public:
  UnkBGBatchAttacker(int nRecips);
  void reset();
  void addRound(int nSentAlice, int nSentOther, const vec<int> &nReceived);
  void guessAlice(vec<double> &nRecipients);
  ~UnkBGBatchAttacker() {}
};

// ======================================================================
// Timed mixes and mix-nets.

class DelayMix : public Mixnet
{
private:
  int maxDelay;
  int poolIdx;
  std::vector< vec<int>* > pools;
  Dist<int> *delayDist;

protected:
  virtual int getDelay() { return delayDist->get(); }

public:
  DelayMix(int nRecips, int maxDelay, Dist<int> *delayDist);
  void reset();
  void addRound(const vec<int> &alice,
                int nAlice,
                const vec<int> &background,
                FwdAttacker *a);
  ~DelayMix();
};

class DelayMixAttacker : public FwdAttacker
{
private:
  int nRecips;
  int maxDelay;

  int nAliceIdx;
  std::vector<int> nAliceHist;
  std::vector<int> nOtherHist;

  vec<double> background;
  vec<double> observed;
  double nObservedOther;
  double nObservedAlice;

  Dist<int> *delayDist;

  int aHist(int rds) { return rds>maxDelay ? 0 : 
      nAliceHist[(maxDelay+nAliceIdx-rds)%maxDelay]; }
  int oHist(int rds) { return rds>maxDelay ? 0 : 
      nOtherHist[(maxDelay+nAliceIdx-rds)%maxDelay]; }
  
  double expectedAliceMsgs();
  double expectedOtherMsgs();

public:
  DelayMixAttacker(int nRecips, int maxDelay, Dist<int> *delayDist);
  void reset();
  void addRound(int nSentAlice, int nSentOther, const vec<int> &nReceived);
  void guessAlice(vec<double> &nRecipients);
  ~DelayMixAttacker() {}
};

#endif /* _SIM_H */

--- NEW FILE: simmain.cpp ---

#include "sim.h"
#include "vec.h"
#include <vector>
#include <algorithm>
#include <iostream>

using std::cout;
using std::ostream;
using std::endl;

template<class C>
void pvec(ostream &o, const std::vector<C> &v)
{
  o << "{";
  for (typename std::vector<C>::const_iterator i = v.begin(); i != v.end(); ++i) {
    o << *i << " ";
  }
  o << "}" << endl;
}

class XTrial : public Trial {
public:
  int attempt() 
  {
    SDTrial x(100, 10, 1000);
    return x.attempt();
  }
};
  
struct percentiles {
  int min, p25, p50, p75, p90, p95, max;
};


void
getPercentile(Trial &test, int nTrials, percentiles &p)
{
  std::vector<int> results(nTrials);

  for (int i = 0; i < nTrials; ++i) {
    results[i] = test.attempt();
    cout << results[i] << endl;
  }
  std::sort(results.begin(), results.end());

  p.min = results[0];
  p.p25 = results[(int)(nTrials *.25)];
  p.p50 = results[(int)(nTrials *.50)];
  p.p75 = results[(int)(nTrials *.75)];
  p.p90 = results[(int)(nTrials *.90)];
  p.p95 = results[(int)(nTrials *.95)];
  p.max = results[nTrials-1];
}


int 
main(int c, char **v)
{
  SDTrial trial(10000, 10, 3000);
  // XTrial trial;
  percentiles p;
  getPercentile(trial, 5, p);
  cout << "======" 
       << p.min << " " << p.p50 << " "
       << p.p90 << " " << p.max << "=====" << endl;

  return 0;
}


--- NEW FILE: test.cpp ---

#include <iostream>
#include "vec.h"
#include "comb.h"

using namespace std;

#define test(x) if (!(x)) { cerr << __LINE__ << ": !" << (x) << endl; return; }
#define test_eq(x,y) if ((x)!=(y)) { \
     cerr << __LINE__ << ": " << (x) << "!=" << (y) << endl; return; }
#define test_neq(x,y) if ((x)==(y)) { \
     cerr << __LINE__ << ": " << (x) << "==" << (y) << endl; return; }

vector<int> vec<double>::topN(int n) const;

void 
test_vec() {
    vec<int> vi1(5);
    vec<int> vi2(5, 7);
    vec<double> vd(6, 2.0);
    vec<int> vi3(vd);
    vec<double> vd2(vi2);

    test_eq(vi3[0], 2);

    vi1[0] = 3;
    test_eq(3, vi1[0]);
    vi1[1] = 4;
    vi1[2] = 5;
    vi1[3] = 6;
    vi1[4] = 7;
    vi3 = vi1+vi2;
    test_eq(vi3[0], 10);
    test_eq(vi3[4], 14);
    vi1 += vi2;
    test_eq(vi1[0], 10);
    test_eq(vi1[4], 14);

    vi1 *= 2;
    test_eq(vi1[0], 20);
    test_eq(vi1[4], 28);

    vd /= 3;
    test_eq(vd[0], 2.0/3);
    test_eq(vd[5], 2.0/3);

    test_eq(vi1.minidx(), 0);
    test_eq(vi1.maxidx(), 4);


    vec<double> vdx(10, 0.0);
    vdx[7] = 3;
    vdx[3] = 10;
    vdx[9] = 60;
    vdx[1] = 1;
    std::vector<int> top3 = vdx.topN(3);
    test_eq(top3[0], 3);
    test_eq(top3[1], 7);
    test_eq(top3[2], 9);
    std::vector<int> top4 = vdx.topN(4);
    test_eq(top4[0], 8);
    test_eq(top4[1], 3);
    test_eq(top4[2], 7);
    test_eq(top4[3], 9);

    cerr << "Vectors OK" << endl;
}

void 
test_comb() {
    test_eq(1, fact(0));
    test_eq(1, fact(1));
    test_eq(2, fact(2));
    test_eq(120, fact(5));

    test_eq(1, comb(10, 0));
    test_eq(1, comb(1, 1));
    test_eq(5, comb(5, 1));
    test_eq(10, comb(5, 2));
    test_eq(10, comb(5, 3));
    test_eq(5, comb(5, 4));
    test_eq(1, comb(5, 5));

    test_eq(1, comb(6, 0));
    test_eq(6, comb(6, 1));
    test_eq(15, comb(6, 2));
    test_eq(20, comb(6, 3));
    test_eq(15, comb(6, 4));
    test_eq(6, comb(6, 5));
    test_eq(1, comb(6, 6));

    cerr << "Combinatorics OK" << endl;
}

int 
main(int c, char **v) {
    test_vec();
    test_comb();

    return 0;
}  


--- NEW FILE: vec.cpp ---

#include "vec.h"
#include <functional>

template<class C> 
vec<C>::vec(int n) : rep(n) 
{}

template class vec<int>;
template class vec<double>;

--- NEW FILE: vec.h ---

#ifndef _VEC_H 
#define _VEC_H

#include <assert.h>
#include <vector>
#include <algorithm>
#include <functional>
#include <iostream>

template <class C>
class vec {
 private:
    vec();

    class Item {
     public:
      int index;
      const vec<C> *v;
      Item() {}
      Item(const vec<C> &vec, int i) : index(i), v(&vec) {}
      const C& getValue() const { return (*v)[index]; }
      bool operator<(const vec<C>::Item &o) const { 
        return getValue() < o.getValue(); 
      }
      bool operator>(const vec<C>::Item &o) const { 
        return getValue() > o.getValue(); 
      }
      bool operator==(const vec<C>::Item &o) const { 
        return getValue() == o.getValue(); 
      }
      ~Item() {}
    };

 protected:
    std::vector<C> rep;
    
 public:
    explicit vec(int len);
    vec(int len, C val) : rep(len, val) {}
    vec(const vec<C> &other) : rep(other.rep) {}
    template <class D> vec(const vec<D> &other) : rep(other.size()) {
	int sz = rep.size();
	for (int i = 0; i < sz; ++i) {
	    rep[i] = static_cast<C>(other[i]);
	}  
    }

    void reset(const C& v) {
      int sz = size();
      for(int i = 0; i < sz; ++i) {
        rep[i] = v;
      }
    }

    int size() const {
	return rep.size();
    }
 
    vec<C> &operator=(const vec<C>& other) {
	rep = other.rep;
	return *this;
    }
    
    vec<C> operator+(const vec<C>& other) const {
	return (vec<C>(*this) += other);
    }
    vec<C> operator-(const vec<C>& other) const {
	return (vec<C>(*this) -= other);
    }
    vec<C> operator*(const vec<C>& other) const {
	return (vec<C>(*this) *= other);
    }
    vec<C> operator/(const vec<C>& other) const {
	return (vec<C>(*this) /= other);
    }

    vec<C> operator+(const C& other) const {
	return (vec<C>(*this) += other);
    }
    vec<C> operator-(const C& other) const {
	return (vec<C>(*this) -= other);
    }
    vec<C> operator*(const C& other) const {
	return (vec<C>(*this) *= other);
    }
    vec<C> operator/(const C& other) const {
	return (vec<C>(*this) /= other);
    }

    const C &operator[](int i) const { return rep[i]; }
    C &operator[](int i) { return rep[i]; }

    template<class D>
    vec<C> &operator+=(const vec<D>& other) {
	int sz = rep.size();
	assert(other.size() == sz);
	for (int i = 0; i < sz; ++i) {
	    rep[i] += other.rep[i];
	}
	return *this;
    }
    template<class D>
    vec<C> &operator-=(const vec<D>& other) {
	int sz = rep.size();
	assert(other.size() == sz);
	for (int i = 0; i < sz; ++i) {
	    rep[i] -= other.rep[i];
	}
	return *this;
    }
    template<class D>
    vec<C> &operator*=(const vec<D>& other) {
	int sz = rep.size();
	assert(other.size() == sz);
	for (int i = 0; i < sz; ++i) {
	    rep[i] *= other.rep[i];
	}
	return *this;
    }
    template<class D>
    vec<C> &operator/=(const vec<D> other) {
	int sz = rep.size();
	assert(other.size() == sz);
	for (int i = 0; i < sz; ++i) {
	    rep[i] /= other.rep[i];
	}
	return *this;
    }
    vec<C> operator+=(const C& other) {
	int sz = rep.size();
	for (int i = 0; i < sz; ++i) { rep[i]+=other; }
	return *this;
    }
    vec<C> operator-=(const C& other) {
	int sz = rep.size();
	for (int i = 0; i < sz; ++i) { rep[i]-=other; }
	return *this;
    }
    vec<C> operator*=(const C& other) {
	int sz = rep.size();
	for (int i = 0; i < sz; ++i) { rep[i]*=other; }
	return *this;
    }
    vec<C> operator/=(const C& other) {
	int sz = rep.size();
	for (int i = 0; i < sz; ++i) { rep[i]/=other; }
	return *this;
    }

    C total(C start) const {
      int sz = rep.size();
      for (int i = 0; i < sz; ++i) {
        start += rep[i];
      }
      return start;
    }

    int maxidx() const {
	int sz = rep.size();
	C max = rep[0];
	int maxidx = 0;
	for (int i = 1; i < sz; ++i) {
	    if (rep[i] > max) {
		max = rep[i];
		maxidx = i;
	    }
	}
	return maxidx;
    }

    int minidx() const {
	int sz = rep.size();
	C min = rep[0];
	int minidx = 0;
	for (int i = 1; i < sz; ++i) {
	    if (rep[i] < min) {
		min = rep[i];
		minidx = i;
	    }
	}
	return minidx;
    }
    
    // Returns indices with higest values, in order of indices.
    std::vector<int> topN(int n) const
    {
      int sz = size();
      std::vector<vec<C>::Item> items( sz );
      for (int i = 0; i<sz; ++i) {
        items[i] = vec<C>::Item(*this, i);
      }
      
      std::partial_sort(items.begin(), items.begin()+n, items.end(),
                        std::greater<vec<C>::Item>() );
      
      std::vector<int> r(n);
      for (int i = 0; i < n; ++i) {
        r[i] = items[i].index;
      }
      std::sort(r.begin(), r.end());
      
      return r;
    }

    ~vec() {}
};

template<class C>
inline std::ostream &
operator<<(std::ostream &o, const vec<C> &v) {
  int sz = v.size();
  o << "[";
  for (int i = 0; i < sz; ++i) {
    o << v[i] << " ";
  }
  o << "]";
  return o;
}

#endif /* _VEC_H */

***********************************************************************
To unsubscribe, send an e-mail to majordomo@seul.org with
unsubscribe freehaven-cvs       in the body. http://freehaven.net/