Max-Lloyd-Algorithmus

[Rüdiger K. (sleipnir)]

Hiho! Ich wollte den Max-Lloyd-Algorithmus zur Optimalquantisierung in 
C++ implementieren.

Allerdings wollte ich nicht über die Integration der 
Wahrscheinlichkeitsdichteverteilung gehen sondern habe den Algorithmus 
wie folgt abgewandelt:
 Wiederhole:
   * ordne jedem Trainigssample den Rekonstruktionswert zu, der ihm am 
nächsten liegt.
   * berechne die Rekonstruktionswerte als Mittelwert der ihnen 
zugeordneten Samples.

Die Ergebnisse sehen für laplaceverteilte Daten schon ganz gut aus, 
allerdings nicht mehr für einen Sinus als Eingangssignal. Dort müßten 
sich die Rekonstruktionswerte eigentlich im Bereich der großen 
Amplituden konzentrieren, was sie aber nicht wirklich tun.

template<class T>  bool MaxLloydQuantizer<T>::train(const itpp::Vec<T> &trainingdata,double epsilon,unsigned int max_iter)
{
#ifdef DEBUG
  if(!initialized) cerr << "MaxLloydQuantizer<T>::train() called before initializing the quantizer!" << endl;
#endif
  if(trainingdata.length()==0) return false;
  T minv,maxv;
  minv=maxv=trainingdata(0);
  for(const T* ptr=trainingdata._data();ptr < trainingdata._data()+trainingdata.length();ptr++)
  {
    register T val = *ptr;
    minv = ((val < minv) ? val : minv);
    maxv = ((val > maxv) ? val : maxv);
  }

  // berechne nun die Intervallgrenzen eines gleichverteilten Quantisierers
  T step=(maxv-minv)/N;
  T w=minv;
  for(unsigned int k=0;k<=N;k++,w+=step) d(k)=w;
  w=minv+step/2;
  for(unsigned int k=0;k<N;k++,w+=step) r(k)=w;

  itpp::Vec<T> r_neu(N);
  itpp::Vec<unsigned int> r_anz(N);

  unsigned int iteration=0;
  for(;iteration < max_iter;iteration++)
  {
    // ordne jedem Trainingswert einen Rekonstruktionswert zu
    for(int k=0;k<trainingdata.length();k++)
    {
      T val=trainingdata(k);
      T min_dist=std::abs(r(0)-val);
      unsigned min_index=0;
      for(unsigned int l=1;l<N;l++)
      {
        T dist=std::abs(T(r(l)-val));
        if(dist < min_dist)
        {
          min_dist=dist;
          min_index=l;
        }
      }
      r_neu(min_index)+=val;
      r_anz(min_index)++;
    }
    double change=0.0;
    for(unsigned int k=0;k<N;k++)
    {
      double r_neuw = (r_anz(k) != 0 ? r_neu(k)/r_anz(k) : r(k));
      double tmp=r_neuw-r(k);
      change+=tmp*tmp;
      r(k)=r_neuw;
    }

    if(change/N <= epsilon) break;
    r_neu.clear();
    r_anz.clear();
  }
  // berechne Intervallgrenzen neu

  for(unsigned int i=1;i<N;i++)
  {
    T di_neu=(r(i-1)+r(i))/2;

    d(i)=di_neu;
  }
  trained=true;
  return (iteration!=max_iter);
}