From 5ebcb1f52a169a2c2b190025144ca9608b786491 Mon Sep 17 00:00:00 2001
From: drowe67 <drowe67@01035d8c-6547-0410-b346-abe4f91aad63>
Date: Tue, 7 Feb 2012 01:21:37 +0000
Subject: [PATCH] compiling OK, making a few mods

git-svn-id: https://svn.code.sf.net/p/freetel/code@312 01035d8c-6547-0410-b346-abe4f91aad63
---
 codec2-dev/unittest/vq_train_jvm.c | 448 +++++++++++++++++++++++++++++
 1 file changed, 448 insertions(+)
 create mode 100755 codec2-dev/unittest/vq_train_jvm.c

diff --git a/codec2-dev/unittest/vq_train_jvm.c b/codec2-dev/unittest/vq_train_jvm.c
new file mode 100755
index 00000000..a3f5b51b
--- /dev/null
+++ b/codec2-dev/unittest/vq_train_jvm.c
@@ -0,0 +1,448 @@
+/*---------------------------------------------------------------------------*\
+                                                                          
+  FILE........: vq_train_jvm.c                                                  
+  AUTHOR......: Jean-Marc Valin                                            
+  DATE CREATED: 21 Jan 2012
+                                                               
+  Multi-stage Vector Quantoser training program developed by Jean-Marc at 
+  linux.conf.au 2012.  Minor mods by David Rowe
+                                                              
+\*---------------------------------------------------------------------------*/
+
+/*
+  Copyright (C) 2012 Jean-Marc Valin 
+
+  All rights reserved.
+
+  This program is free software; you can redistribute it and/or modify
+  it under the terms of the GNU Lesser General Public License version 2, as
+  published by the Free Software Foundation.  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 Lesser General Public License
+  along with this program; if not, see <http://www.gnu.org/licenses/>.
+*/
+
+
+#include <valgrind/memcheck.h>
+
+#include <assert.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#define MIN(a,b) ((a)<(b)?(a):(b))
+#define COEF 0.70f
+#define MAX_ENTRIES 16384
+
+void compute_weights(const float *x, float *w, int ndim)
+{
+  int i;
+  w[0] = MIN(x[0], x[1]-x[0]);
+  for (i=1;i<ndim-1;i++)
+    w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]);
+  w[ndim-1] = MIN(x[ndim-1]-x[ndim-2], M_PI-x[ndim-1]);
+  
+  for (i=0;i<ndim;i++)
+    w[i] = 1./(.01+w[i]);
+  w[0]*=3;
+  w[1]*=2;
+}
+
+int find_nearest(const float *codebook, int nb_entries, float *x, int ndim)
+{
+  int i, j;
+  float min_dist = 1e15;
+  int nearest = 0;
+  
+  for (i=0;i<nb_entries;i++)
+  {
+    float dist=0;
+    for (j=0;j<ndim;j++)
+      dist += (x[j]-codebook[i*ndim+j])*(x[j]-codebook[i*ndim+j]);
+    if (dist<min_dist)
+    {
+      min_dist = dist;
+      nearest = i;
+    }
+  }
+  return nearest;
+}
+
+int find_nearest_weighted(const float *codebook, int nb_entries, float *x, const float *w, int ndim)
+{
+  int i, j;
+  float min_dist = 1e15;
+  int nearest = 0;
+  
+  for (i=0;i<nb_entries;i++)
+  {
+    float dist=0;
+    for (j=0;j<ndim;j++)
+      dist += w[j]*(x[j]-codebook[i*ndim+j])*(x[j]-codebook[i*ndim+j]);
+    if (dist<min_dist)
+    {
+      min_dist = dist;
+      nearest = i;
+    }
+  }
+  return nearest;
+}
+
+int quantize_lsp(const float *x, const float *codebook1, const float *codebook2, 
+		 const float *codebook3, int nb_entries, float *xq, int ndim)
+{
+  int i, n1, n2, n3;
+  float err[ndim], err2[ndim], err3[ndim];
+  float w[ndim], w2[ndim], w3[ndim];
+  
+  w[0] = MIN(x[0], x[1]-x[0]);
+  for (i=1;i<ndim-1;i++)
+    w[i] = MIN(x[i]-x[i-1], x[i+1]-x[i]);
+  w[ndim-1] = MIN(x[ndim-1]-x[ndim-2], M_PI-x[ndim-1]);
+  
+  /*
+  for (i=0;i<ndim;i++)
+    w[i] = 1./(.003+w[i]);
+  w[0]*=3;
+  w[1]*=2;*/
+  compute_weights(x, w, ndim);
+  
+  for (i=0;i<ndim;i++)
+    err[i] = x[i]-COEF*xq[i];
+  n1 = find_nearest(codebook1, nb_entries, err, ndim);
+  
+  for (i=0;i<ndim;i++)
+  {
+    xq[i] = COEF*xq[i] + codebook1[ndim*n1+i];
+    err[i] -= codebook1[ndim*n1+i];
+  }
+  for (i=0;i<ndim/2;i++)
+  {
+    err2[i] = err[2*i];  
+    err3[i] = err[2*i+1];
+    w2[i] = w[2*i];  
+    w3[i] = w[2*i+1];
+  }
+  n2 = find_nearest_weighted(codebook2, nb_entries, err2, w2, ndim/2);
+  n3 = find_nearest_weighted(codebook3, nb_entries, err3, w3, ndim/2);
+  
+  for (i=0;i<ndim/2;i++)
+  {
+    xq[2*i] += codebook2[ndim*n2/2+i];
+    xq[2*i+1] += codebook3[ndim*n3/2+i];
+  }
+  return 0;
+}
+
+void split(float *codebook, int nb_entries, int ndim)
+{
+  int i,j;
+  for (i=0;i<nb_entries;i++)
+  {
+    for (j=0;j<ndim;j++)
+    {
+      float delta = .01*(rand()/(float)RAND_MAX-.5);
+      codebook[i*ndim+j] += delta;
+      codebook[(i+nb_entries)*ndim+j] = codebook[i*ndim+j] - delta;
+    }
+  }
+}
+
+void update(float *data, int nb_vectors, float *codebook, int nb_entries, int ndim)
+{
+  int i,j;
+  int count[nb_entries];
+  int nearest[nb_vectors];
+  
+  for (i=0;i<nb_entries;i++)
+    count[i] = 0;
+  
+  for (i=0;i<nb_vectors;i++)
+  {
+    nearest[i] = find_nearest(codebook, nb_entries, data+i*ndim, ndim);
+  }
+  for (i=0;i<nb_entries*ndim;i++)
+    codebook[i] = 0;
+  
+  for (i=0;i<nb_vectors;i++)
+  {
+    int n = nearest[i];
+    count[n]++;
+    for (j=0;j<ndim;j++)
+      codebook[n*ndim+j] += data[i*ndim+j];
+  }
+
+  float w2=0;
+  for (i=0;i<nb_entries;i++)
+  { 
+    for (j=0;j<ndim;j++)
+      codebook[i*ndim+j] *= (1./count[i]);
+    w2 += (count[i]/(float)nb_vectors)*(count[i]/(float)nb_vectors);
+  }
+  fprintf(stderr, "%f / %d\n", 1./w2, nb_entries);
+}
+
+void update_weighted(float *data, float *weight, int nb_vectors, float *codebook, int nb_entries, int ndim)
+{
+  int i,j;
+  float count[MAX_ENTRIES][ndim];
+  int nearest[nb_vectors];
+  
+  for (i=0;i<nb_entries;i++)
+    for (j=0;j<ndim;j++)
+      count[i][j] = 0;
+  
+  for (i=0;i<nb_vectors;i++)
+  {
+    nearest[i] = find_nearest_weighted(codebook, nb_entries, data+i*ndim, weight+i*ndim, ndim);
+  }
+  for (i=0;i<nb_entries*ndim;i++)
+    codebook[i] = 0;
+  
+  for (i=0;i<nb_vectors;i++)
+  {
+    int n = nearest[i];
+    for (j=0;j<ndim;j++)
+    {
+      float w = sqrt(weight[i*ndim+j]);
+      count[n][j]+=w;
+      codebook[n*ndim+j] += w*data[i*ndim+j];
+    }
+  }
+
+  //float w2=0;
+  for (i=0;i<nb_entries;i++)
+  { 
+    for (j=0;j<ndim;j++)
+      codebook[i*ndim+j] *= (1./count[i][j]);
+    //w2 += (count[i]/(float)nb_vectors)*(count[i]/(float)nb_vectors);
+  }
+  //fprintf(stderr, "%f / %d\n", 1./w2, nb_entries);
+}
+
+void vq_train(float *data, int nb_vectors, float *codebook, int nb_entries, int ndim)
+{
+  int i, j, e;
+  e = 1;
+  for (j=0;j<ndim;j++)
+    codebook[j] = 0;
+  for (i=0;i<nb_vectors;i++)
+    for (j=0;j<ndim;j++)
+      codebook[j] += data[i*ndim+j];
+  for (j=0;j<ndim;j++)
+    codebook[j] *= (1./nb_vectors);
+  
+  
+  while (e< nb_entries)
+  {
+    split(codebook, e, ndim);
+    fprintf(stderr, "%d\n", e);
+    e<<=1;
+    for (j=0;j<ndim;j++)
+      update(data, nb_vectors, codebook, e, ndim);
+  }
+}
+
+void vq_train_weighted(float *data, float *weight, int nb_vectors, float *codebook, int nb_entries, int ndim)
+{
+  int i, j, e;
+  e = 1;
+  for (j=0;j<ndim;j++)
+    codebook[j] = 0;
+  for (i=0;i<nb_vectors;i++)
+    for (j=0;j<ndim;j++)
+      codebook[j] += data[i*ndim+j];
+  for (j=0;j<ndim;j++)
+    codebook[j] *= (1./nb_vectors);
+   
+  while (e<nb_entries)
+  {
+    split(codebook, e, ndim);
+    fprintf(stderr, "%d\n", e);
+    e<<=1;
+    for (j=0;j<ndim;j++)
+      update_weighted(data, weight, nb_vectors, codebook, e, ndim);
+  }
+}
+
+
+int main(int argc, char **argv)
+{
+  int i,j;
+  FILE *ftrain;
+  int nb_vectors, nb_entries, ndim;
+  float *data, *pred, *codebook, *codebook2, *codebook3;
+  float *weight, *weight2, *weight3;
+  float *delta, *delta2;
+  float tmp, err;
+  int ret;
+
+  printf("Jean-Marc Valin's Split VQ training program....\n");
+
+  if (argc != 5) {
+      printf("usage: %s TrainTextFile K(dimension) M(codebook size) VQFile\n", argv[0]);
+      exit(0);      
+  }
+  
+  ndim = atoi(argv[2]);
+  nb_vectors = atoi(argv[3]);
+  nb_entries = atoi(argv[3]);
+
+  /* determine size of training file */
+
+  ftrain = fopen(argv[1],"rt");  assert(ftrain != NULL);
+  nb_vectors = 0;
+  while (1) {
+    if (feof(ftrain))
+      break;
+    for (j=0;j<ndim;j++)
+    {
+	ret = fscanf(ftrain, "%f ", &tmp);
+    }
+    printf("\r%d lines",nb_vectors++);
+  }
+
+  rewind(ftrain);
+
+  printf("ndim %d nb_vectors %d nb_entries %d\n", ndim, nb_vectors, nb_entries);
+
+  data = malloc(nb_vectors*ndim*sizeof(*data));
+  weight = malloc(nb_vectors*ndim*sizeof(*weight));
+  weight2 = malloc(nb_vectors*ndim*sizeof(*weight2));
+  weight3 = malloc(nb_vectors*ndim*sizeof(*weight3));
+  pred = malloc(nb_vectors*ndim*sizeof(*pred));
+  codebook = malloc(nb_entries*ndim*sizeof(*codebook));
+  codebook2 = malloc(nb_entries*ndim*sizeof(*codebook2));
+  codebook3 = malloc(nb_entries*ndim*sizeof(*codebook3));
+  
+  for (i=0;i<nb_vectors;i++)
+  {
+    if (feof(ftrain))
+      break;
+    for (j=0;j<ndim;j++)
+    {
+	ret = fscanf(ftrain, "%f ", &data[i*ndim+j]);
+    }
+  }
+  nb_vectors = i;
+  VALGRIND_CHECK_MEM_IS_DEFINED(data, nb_entries*ndim);
+
+  /* determine weights for each training vector */
+
+  for (i=0;i<nb_vectors;i++)
+  {
+    compute_weights(data+i*ndim, weight+i*ndim, ndim);
+    for (j=0;j<ndim/2;j++)
+    {
+      weight2[i*ndim/2+j] = weight[i*ndim+2*j];
+      weight3[i*ndim/2+j] = weight[i*ndim+2*j+1];
+    }
+  }
+
+  /* initial predictor state */
+
+  for (i=0;i<ndim;i++)
+    pred[i] = data[i] - M_PI*(i+1)/(ndim+1);
+
+  /* generate predicted data for training */
+
+  for (i=1;i<nb_vectors;i++)
+  {
+    for (j=0;j<ndim;j++)
+      pred[i*ndim+j] = data[i*ndim+j] - COEF*data[(i-1)*ndim+j];
+  }
+
+  VALGRIND_CHECK_MEM_IS_DEFINED(pred, nb_entries*ndim);
+
+  /* train first stage */
+
+  vq_train(pred, nb_vectors, codebook, nb_entries, ndim);
+  
+  delta = malloc(nb_vectors*ndim*sizeof(*data));
+  err = 0;
+  for (i=0;i<nb_vectors;i++)
+  {
+    int nearest = find_nearest(codebook, nb_entries, &pred[i*ndim], ndim);
+    for (j=0;j<ndim;j++)
+    {
+      //delta[i*ndim+j] = data[i*ndim+j] - codebook[nearest*ndim+j];
+      //printf("%f ", delta[i*ndim+j]);
+      //err += (delta[i*ndim+j])*(delta[i*ndim+j]);
+      delta[i*ndim/2+j/2+(j&1)*nb_vectors*ndim/2] = pred[i*ndim+j] - codebook[nearest*ndim+j];
+      //printf("%f ", delta[i*ndim/2+j/2+(j&1)*nb_vectors*ndim/2]);
+      err += (delta[i*ndim/2+j/2+(j&1)*nb_vectors*ndim/2])*(delta[i*ndim/2+j/2+(j&1)*nb_vectors*ndim/2]);
+    }
+    //printf("\n");
+  }
+  fprintf(stderr, "LSP RMS error: %f\n", sqrt(err/nb_vectors/ndim));
+
+#if 1
+  vq_train(delta, nb_vectors, codebook2, nb_entries, ndim/2);
+  vq_train(delta+ndim*nb_vectors/2, nb_vectors, codebook3, nb_entries, ndim/2);
+#else
+  vq_train_weighted(delta, weight2, nb_vectors, codebook2, nb_entries, ndim/2);
+  vq_train_weighted(delta+ndim*nb_vectors/2, weight3, nb_vectors, codebook3, nb_entries, ndim/2);
+#endif
+  err = 0;
+  
+  delta2 = delta + nb_vectors*ndim/2;
+
+  for (i=0;i<nb_vectors;i++)
+  {
+    int n1, n2;
+    n1 = find_nearest(codebook2, nb_entries, &delta[i*ndim/2], ndim/2);
+    for (j=0;j<ndim/2;j++)
+    {
+      delta[i*ndim/2+j] = delta[i*ndim/2+j] - codebook2[n1*ndim/2+j];
+      err += (delta[i*ndim/2+j])*(delta[i*ndim/2+j]);
+    }
+
+    n2 = find_nearest(codebook3, nb_entries, &delta2[i*ndim/2], ndim/2);
+    for (j=0;j<ndim/2;j++)
+    {
+      delta[i*ndim/2+j] = delta[i*ndim/2+j] - codebook2[n2*ndim/2+j];
+      err += (delta2[i*ndim/2+j])*(delta2[i*ndim/2+j]);
+    }
+  }
+  fprintf(stderr, "LSP RMS error stage 2: %f\n", sqrt(err/nb_vectors/ndim));
+  
+  float xq[ndim];
+  for (i=0;i<ndim;i++)
+    xq[i] = M_PI*(i+1)/(ndim+1);
+  
+  for (i=0;i<nb_vectors;i++)
+  {
+    quantize_lsp(data+i*ndim, codebook, codebook2, 
+		 codebook3, nb_entries, xq, ndim);
+    /*for (j=0;j<ndim;j++)
+      printf("%f ", xq[j]);
+    printf("\n");*/
+  }
+  
+  for (i=0;i<nb_entries;i++)
+  {
+    for (j=0;j<ndim;j++)
+      printf("%f ", codebook[i*ndim+j]);
+    printf("\n");
+  }
+  printf("\n");
+  for (i=0;i<nb_entries;i++)
+  {
+    for (j=0;j<ndim/2;j++)
+      printf("%f ", codebook2[i*ndim/2+j]);
+    printf("\n");
+  }
+  printf("\n");
+  for (i=0;i<nb_entries;i++)
+  {
+    for (j=0;j<ndim/2;j++)
+      printf("%f ", codebook2[i*ndim/2+j]);
+    printf("\n");
+  }
+
+  return 0;
+}
-- 
2.25.1