gtr.c

QccPack implementation in C

/*
 * 
 * QccPack: Quantization, compression, and coding libraries
 * Copyright (C) 1997-2005  James E. Fowler
 * 
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Library General Public
 * License as published by the Free Software Foundation; either
 * version 2 of the License, or (at your option) any later version.
 * 
 * This library 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
 * Library General Public License for more details.
 * 
 * You should have received a copy of the GNU Library General Public
 * License along with this library; if not, write to the
 * Free Software Foundation, Inc., 675 Mass Ave, Cambridge,
 * MA 02139, USA.
 * 
 */


#include "libQccPack.h"


static void QccAVQgtrCalcProb(QccVQCodebook *codebook, const QccVector cnt)
{
  int codeword;
  double sum_cnt = 0;
  
  for (codeword = 0; codeword < codebook->num_codewords; codeword++)
    sum_cnt += cnt[codeword];
  for (codeword = 0; codeword < codebook->num_codewords; codeword++)
    codebook->codeword_probs[codeword] = cnt[codeword]/sum_cnt;
}


static int QccAVQgtrSendSideInfo(QccAVQSideInfo *sideinfo, 
                                 QccAVQSideInfoSymbol *updatevector_symbol,
                                 int update_flag)
{
  QccAVQSideInfoSymbol flag_symbol;
  
  flag_symbol.type = QCCAVQSIDEINFO_FLAG;
  flag_symbol.flag = update_flag;
  flag_symbol.vector = NULL;
  flag_symbol.vector_indices = NULL;

  if (QccAVQSideInfoSymbolWrite(sideinfo, &flag_symbol))
    {
      QccErrorAddMessage("(QccAVQgtrSendSideInfo): Error calling QccAVQSideInfoSymbolWrite()");
      return(1);
    }
  
  if (update_flag == QCCAVQSIDEINFO_FLAG_UPDATE)
    {
      updatevector_symbol->type = QCCAVQSIDEINFO_VECTOR;
      if (QccAVQSideInfoSymbolWrite(sideinfo, updatevector_symbol))
        {
          QccErrorAddMessage("(QccAVQgtrSendSideInfo): Error calling QccAVQSideInfoSymbolWrite()");
          return(1);
        }
    }
  
  return(0);
}


static void QccAVQgtrCalcUpdateCosts(double *update_distortion,
                                     double *update_rate,
                                     const
                                     QccAVQSideInfoSymbol *new_codeword_symbol,
                                     QccVector current_vector,
                                     QccAVQSideInfo *sideinfo)
{
  int offset;
  int component;
  double current_prob;

  offset = sideinfo->codebook_coder.num_levels/2;

  *update_distortion = 
    QccVectorSquareDistance(current_vector,
                            new_codeword_symbol->vector,
                            sideinfo->vector_dimension);

  *update_rate = 0;

  switch (sideinfo->codebook_coder.type)
    {
    case QCCSQSCALARQUANTIZER_GENERAL:
      for (component = 0; component < sideinfo->vector_dimension; component++)
        {
          current_prob = sideinfo->codebook_coder.probs
               [new_codeword_symbol->vector_indices[component]];
          if (current_prob != 0.0)
            *update_rate += -QccMathLog2(current_prob);
        }
      break;

    case QCCSQSCALARQUANTIZER_UNIFORM:
    case QCCSQSCALARQUANTIZER_DEADZONE:
      for (component = 0; component < sideinfo->vector_dimension; component++)
        {
          current_prob = sideinfo->codebook_coder.probs
            [new_codeword_symbol->vector_indices[component] + offset];
          if (current_prob != 0.0)
            *update_rate += -QccMathLog2(current_prob);
        }
      break;
    }
}


static void QccAVQgtrUpdateCodebookCoderProbs(QccAVQSideInfo *sideinfo,
                                              const QccAVQSideInfoSymbol 
                                              *new_codeword_symbol,
                                              double omega)
{
  int cnt_total = 0;
  int component;
  int offset;
  int codeword;

  offset = sideinfo->codebook_coder.num_levels/2;

  for (codeword = 0; codeword < sideinfo->codebook_coder.num_levels; 
       codeword++)
    {
      sideinfo->codebook_coder.probs[codeword] =
        (int)(sideinfo->codebook_coder.probs[codeword] * omega + 0.5);
      cnt_total +=
        sideinfo->codebook_coder.probs[codeword];
    }

  switch (sideinfo->codebook_coder.type)
    {
    case QCCSQSCALARQUANTIZER_GENERAL:
      for (component = 0; component < sideinfo->vector_dimension; component++)
        sideinfo->codebook_coder.probs
          [new_codeword_symbol->vector_indices[component]] += 1;
      break;

    case QCCSQSCALARQUANTIZER_UNIFORM:
    case QCCSQSCALARQUANTIZER_DEADZONE:
      for (component = 0; component < sideinfo->vector_dimension; component++)
        sideinfo->codebook_coder.probs
          [new_codeword_symbol->vector_indices[component] + offset] += 1;
      break;
    }

  cnt_total += sideinfo->vector_dimension;

  for (codeword = 0; codeword < sideinfo->codebook_coder.num_levels; 
       codeword++)
    sideinfo->codebook_coder.probs[codeword] /= cnt_total;

}


int QccAVQgtrEncode(const QccDataset *dataset, 
                    QccVQCodebook *codebook, 
                    QccChannel *channel,
                    QccAVQSideInfo *sideinfo,
                    double lambda,
                    double omega,
                    QccVQDistortionMeasure distortion_measure)
{
  int return_value;
  int block_size;
  int vector_dimension;
  int current_vector;
  double winning_distortion;
  double update_distortion;
  double update_rate;
  int winner;
  QccDataset tmp_dataset;
  QccVector cnt = NULL;
  int codeword;
  double newcnt;
  double delta_J;
  QccAVQSideInfoSymbol *new_codeword_symbol = NULL;

  if ((dataset == NULL) || (codebook == NULL) || (channel == NULL) ||
      (sideinfo == NULL))
    return(0);

  if (dataset->vectors == NULL)
    return(0);
  
  vector_dimension = dataset->vector_dimension;
  if ((vector_dimension != codebook->codeword_dimension) ||
      (vector_dimension != sideinfo->vector_dimension))
    {
      QccErrorAddMessage("(QccAVQgtrEncode): Codebook %s, dataset %s, and sideinfo %s have different vector dimensions",
                         codebook->filename, dataset->filename, sideinfo->filename);
      goto QccAVQError;
    }
  
  if ((new_codeword_symbol = 
       QccAVQSideInfoSymbolAlloc(vector_dimension)) == NULL)
    {
      QccErrorAddMessage("(QccAVQgtrEncode): Error calling QccAVQSideInfoSymbolAlloc()");
      goto QccAVQError;
    }
  new_codeword_symbol->type = QCCAVQSIDEINFO_VECTOR;

  /*
    precision = 
    (int)(log((double)sideinfo->codebook_coder.num_levels)/log(2.0) + 0.5);
  */

  block_size = QccDatasetGetBlockSize(dataset);
  
  QccDatasetInitialize(&tmp_dataset);
  tmp_dataset.num_vectors = 1;
  tmp_dataset.vector_dimension = dataset->vector_dimension;
  tmp_dataset.access_block_size = QCCDATASET_ACCESSWHOLEFILE;
  tmp_dataset.num_blocks_accessed = 0;
  if (QccDatasetAlloc(&tmp_dataset))
    {
      QccErrorAddMessage("(QccAVQgtrEncode): Error calling QccDatasetAlloc()");
      goto QccAVQError;
    }
  
  if ((cnt = 
       QccVectorAlloc(sideinfo->max_codebook_size)) == NULL)
    {
      QccErrorAddMessage("(QccAVQgtrEncode): Error allocating memory");
      goto QccAVQError;
    }
  
  for (current_vector = 0; current_vector < block_size; current_vector++)
    {
      QccVectorCopy(tmp_dataset.vectors[0],
                    dataset->vectors[current_vector],
                    vector_dimension);
      
      if (QccVQEntropyConstrainedVQ(&tmp_dataset, codebook, lambda, 
                                    (QccVector)&winning_distortion, &winner,
                                    distortion_measure))
        {
          QccErrorAddMessage("(QccAVQgtrEncode): Error calling QccVQEntropyConstrainedVQ()");
          goto QccAVQError;
        }

      QccVectorCopy(new_codeword_symbol->vector,
                    dataset->vectors[current_vector],
                    vector_dimension);
      QccAVQSideInfoCodebookCoder(sideinfo, new_codeword_symbol);

      /*
        update_distortion = 
        QccVectorSquareDistance(dataset->vectors[current_vector],
        new_codeword_symbol->vector,
        vector_dimension);
      */

      QccAVQgtrCalcUpdateCosts(&update_distortion,
                               &update_rate,
                               new_codeword_symbol,
                               dataset->vectors[current_vector],
                               sideinfo);

      for (codeword = 0; codeword < codebook->num_codewords; codeword++)
        cnt[codeword] = 
          (int)(codebook->codeword_probs[codeword] * omega + 0.5);

      /*
        delta_J = (update_distortion - winning_distortion) +
        lambda*vector_dimension*precision;
      */
      delta_J = (update_distortion - winning_distortion) +
        lambda * update_rate;
      
      if (delta_J < 0.0)
        {
          /*  Update codebook  */
          if (codebook->num_codewords < sideinfo->max_codebook_size)
            {
              if (QccVQCodebookAddCodeword(codebook,
                                           (QccVQCodeword)
                                           new_codeword_symbol->vector))
                {
                  QccErrorAddMessage("(QccAVQgtrEncode): Error calling QccVQCodebookAddCodeword()");
                  goto QccAVQError;
                }
              winner = codebook->num_codewords - 1;
              newcnt = 1;
            }
          else
            {
              QccVectorCopy((QccVector)codebook->codewords
                            [codebook->num_codewords - 1], 
                            new_codeword_symbol->vector,
                            vector_dimension);
              newcnt = cnt[winner] + 1;
            }
          
          if (QccVQCodebookMoveCodewordToFront(codebook, 
                                               codebook->num_codewords - 1))
            {
              QccErrorAddMessage("(QccAVQgtrEncode): Error calling QccVQCodebookMoveCodewordToFront()");
              goto QccAVQError;
            }
          
          if (winner != (codebook->num_codewords - 1))
            {
              cnt[winner] = newcnt/2;
              cnt[codebook->num_codewords - 1] = cnt[winner];
            }
          else
            cnt[winner] = newcnt;
          
          QccAVQgtrCalcProb(codebook, cnt);
          QccVectorMoveComponentToFront(codebook->codeword_probs,
                                        codebook->num_codewords,
                                        codebook->num_codewords - 1);
          
          channel->channel_symbols[current_vector] = QCCCHANNEL_NULLSYMBOL;
          
          if (QccAVQgtrSendSideInfo(sideinfo, new_codeword_symbol, 
                                    QCCAVQSIDEINFO_FLAG_UPDATE))
            {
              QccErrorAddMessage("(QccAVQgtrEncode): Error calling QccAVQgtrSendSideInfo()");
              goto QccAVQError;
            }
          QccAVQgtrUpdateCodebookCoderProbs(sideinfo,
                                            new_codeword_symbol,
                                            omega);
        }
      else
        {
          /*  No codebook update  */
          QccVQCodebookMoveCodewordToFront(codebook, winner);
          cnt[winner]++;
          QccAVQgtrCalcProb(codebook, cnt);
          if (QccVectorMoveComponentToFront(codebook->codeword_probs,