pnorm.hpp

dysii是一款非常出色的滤波函数库

#ifndef INDII_ML_AUX_PNORM_HPP
#define INDII_ML_AUX_PNORM_HPP

#include "Norm.hpp"

namespace indii {
  namespace ml {
    namespace aux {
/**
 * Vector p-norm, \f$\|\cdot\|_p\f$.
 *
 * @author Lawrence Murray <[e-mail]>
 * @version $Rev: 404 $
 * @date $Date: 2008-03-05 14:52:55 +0000 (Wed, 05 Mar 2008) $
 *
 * @param P Degree of the p-norm.
 */
template <unsigned int P>
class PNorm : public Norm {
public:
  virtual double operator()(const vector& x) const;

  virtual vector sample(const unsigned int N) const;
  
};

    }
  }
}

#include "Random.hpp"
#include "vector.hpp"

#include <set>

using namespace indii::ml::aux;

template<unsigned int P>
inline double indii::ml::aux::PNorm<P>::operator()(const vector& x) const {
  return norm<vector,P>(x);
}

template<unsigned int P>
indii::ml::aux::vector indii::ml::aux::PNorm<P>::sample(
    const unsigned int N) const {   
  unsigned int i;
  std::multiset<double> u;
  std::multiset<double>::iterator iter, end;
  vector x(N);
  
  /* sample independent uniform variates on [0,1] and sort */
  for (i = 0; i < N - 1; i++) {
    u.insert(Random::uniform(0.0, 1.0));
  }
  
  /* calculate increments */
  if (!u.empty()) {
    iter = u.begin();
    end = u.end();
    x(0) = *iter;
    iter++;
    for (i = 1; iter != end; iter++, i++) {
      x(i) = *iter - x(i - 1);
    }
    x(N - 1) = 1.0 - x(N - 2);
  } else {
    x(0) = 1.0;
  }
  assert (sum(x) == 1.0);
  
  /* normalise to unit vector in given norm */
  for (i = 0; i < N; i++) {
    x(i) = pow(x(i), 1.0 / P);
    if (round(Random::uniform()) == 0) {
      x(i) *= -1;
    }
  }

  return x;
  
  /* alternative angular coordinates approach (for 2-norm only)? See
   * Wikipedia article on "n-sphere" */
  //unsigned int i;
  //double y;
  //vector x(N), phi(N - 1);
  
  /* sample angular coordinates (radians), last has range of 2*PI,
   * others PI */
  //if (N >= 2) {
  //  for (i = 0; i < N - 2; i++) {
  //    phi(i) = Random::uniform(0.0, M_PI);
  //  }
  //  phi(N - 2) = Random::uniform(0.0, 2.0*M_PI);
  //}
  
  /* calculate cartesian coordinates */
  //y = 1.0;
  //for (i = 0; i < N - 1; i++) {
  //  x(i) = y * cos(phi(i));
  //  y *= sin(phi(i));
  //}
  //x(N - 1) = y;
}

#endif