git.gag.com Git - debian/gnuradio/blob - gr-trellis/src/lib/trellis_siso_combined_f.cc

   1 /* -*- c++ -*- */
   2 /*
   3  * Copyright 2004 Free Software Foundation, Inc.
   4  *
   5  * This file is part of GNU Radio
   6  *
   7  * GNU Radio is free software; you can redistribute it and/or modify
   8  * it under the terms of the GNU General Public License as published by
   9  * the Free Software Foundation; either version 3, or (at your option)
  10  * any later version.
  11  *
  12  * GNU Radio is distributed in the hope that it will be useful,
  13  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15  * GNU General Public License for more details.
  16  *
  17  * You should have received a copy of the GNU General Public License
  18  * along with GNU Radio; see the file COPYING.  If not, write to
  19  * the Free Software Foundation, Inc., 51 Franklin Street,
  20  * Boston, MA 02110-1301, USA.
  21  */
  22
  23 #ifdef HAVE_CONFIG_H
  24 #include "config.h"
  25 #endif
  26
  27 #include <trellis_siso_combined_f.h>
  28 #include <gr_io_signature.h>
  29 #include <stdexcept>
  30 #include <assert.h>
  31 #include <iostream>
  32
  33 static const float INF = 1.0e9;
  34
  35 trellis_siso_combined_f_sptr
  36 trellis_make_siso_combined_f (
  37     const fsm &FSM,
  38     int K,
  39     int S0,
  40     int SK,
  41     bool POSTI,
  42     bool POSTO,
  43     trellis_siso_type_t SISO_TYPE,
  44     int D,
  45     const std::vector<float> &TABLE,
  46     trellis_metric_type_t TYPE)
  47 {
  48   return trellis_siso_combined_f_sptr (new trellis_siso_combined_f (FSM,K,S0,SK,POSTI,POSTO,SISO_TYPE,D,TABLE,TYPE));
  49 }
  50
  51 trellis_siso_combined_f::trellis_siso_combined_f (
  52     const fsm &FSM,
  53     int K,
  54     int S0,
  55     int SK,
  56     bool POSTI,
  57     bool POSTO,
  58     trellis_siso_type_t SISO_TYPE,
  59     int D,
  60     const std::vector<float> &TABLE,
  61     trellis_metric_type_t TYPE)
  62   : gr_block ("siso_combined_f",
  63                           gr_make_io_signature (1, -1, sizeof (float)),
  64                           gr_make_io_signature (1, -1, sizeof (float))),
  65   d_FSM (FSM),
  66   d_K (K),
  67   d_S0 (S0),
  68   d_SK (SK),
  69   d_POSTI (POSTI),
  70   d_POSTO (POSTO),
  71   d_SISO_TYPE (SISO_TYPE),
  72   d_D (D),
  73   d_TABLE (TABLE),
  74   d_TYPE (TYPE)//,
  75   //d_alpha(FSM.S()*(K+1)),
  76   //d_beta(FSM.S()*(K+1))
  77 {
  78     int multiple;
  79     if (d_POSTI && d_POSTO)
  80         multiple = d_FSM.I()+d_FSM.O();
  81     else if(d_POSTI)
  82         multiple = d_FSM.I();
  83     else if(d_POSTO)
  84         multiple = d_FSM.O();
  85     else
  86         throw std::runtime_error ("Not both POSTI and POSTO can be false.");
  87     //printf("constructor: Multiple = %d\n",multiple);
  88     set_output_multiple (d_K*multiple);
  89     //what is the meaning of relative rate for a block with 2 inputs?
  90     //set_relative_rate ( multiple / ((double) d_FSM.I()) );
  91     // it turns out that the above gives problems in the scheduler, so
  92     // let's try (assumption O>I)
  93     //set_relative_rate ( multiple / ((double) d_FSM.O()) );
  94     // I am tempted to automate like this
  95     if(d_FSM.I() <= d_D)
  96       set_relative_rate ( multiple / ((double) d_D) );
  97     else
  98       set_relative_rate ( multiple / ((double) d_FSM.I()) );
  99 }
 100
 101
 102 void
 103 trellis_siso_combined_f::forecast (int noutput_items, gr_vector_int &ninput_items_required)
 104 {
 105   int multiple;
 106   if (d_POSTI && d_POSTO)
 107       multiple = d_FSM.I()+d_FSM.O();
 108   else if(d_POSTI)
 109       multiple = d_FSM.I();
 110   else if(d_POSTO)
 111       multiple = d_FSM.O();
 112   else
 113       throw std::runtime_error ("Not both POSTI and POSTO can be false.");
 114   //printf("forecast: Multiple = %d\n",multiple);
 115   assert (noutput_items % (d_K*multiple) == 0);
 116   int input_required1 =  d_FSM.I() * (noutput_items/multiple) ;
 117   int input_required2 =  d_D * (noutput_items/multiple) ;
 118   //printf("forecast: Output requirements:  %d\n",noutput_items);
 119   //printf("forecast: Input requirements:  %d   %d\n",input_required1,input_required2);
 120   unsigned ninputs = ninput_items_required.size();
 121   assert(ninputs % 2 == 0);
 122   for (unsigned int i = 0; i < ninputs/2; i++) {
 123     ninput_items_required[2*i] = input_required1;
 124     ninput_items_required[2*i+1] = input_required2;
 125   }
 126 }
 127
 128 inline float min(float a, float b)
 129 {
 130   return a <= b ? a : b;
 131 }
 132
 133 inline float min_star(float a, float b)
 134 {
 135   return (a <= b ? a : b)-log(1+exp(a <= b ? a-b : b-a));
 136 }
 137
 138 void siso_algorithm_combined(int I, int S, int O,
 139              const std::vector<int> &NS,
 140              const std::vector<int> &OS,
 141              const std::vector< std::vector<int> > &PS,
 142              const std::vector< std::vector<int> > &PI,
 143              int K,
 144              int S0,int SK,
 145              bool POSTI, bool POSTO,
 146              float (*p2mymin)(float,float),
 147              int D,
 148              const std::vector<float> &TABLE,
 149              trellis_metric_type_t TYPE,
 150              const float *priori, const float *observations, float *post//,
 151              //std::vector<float> &alpha,
 152              //std::vector<float> &beta
 153              )
 154 {
 155   float norm,mm,minm;
 156   std::vector<float> alpha(S*(K+1));
 157   std::vector<float> beta(S*(K+1));
 158   float *prioro = new float[O*K];
 159
 160
 161   if(S0<0) { // initial state not specified
 162       for(int i=0;i<S;i++) alpha[0*S+i]=0;
 163   }
 164   else {
 165       for(int i=0;i<S;i++) alpha[0*S+i]=INF;
 166       alpha[0*S+S0]=0.0;
 167   }
 168
 169   for(int k=0;k<K;k++) { // forward recursion
 170       calc_metric(O, D, TABLE, &(observations[k*D]), &(prioro[k*O]),TYPE); // calc metrics
 171       norm=INF;
 172       for(int j=0;j<S;j++) {
 173           minm=INF;
 174           for(unsigned int i=0;i<PS[j].size();i++) {
 175               //int i0 = j*I+i;
 176               mm=alpha[k*S+PS[j][i]]+priori[k*I+PI[j][i]]+prioro[k*O+OS[PS[j][i]*I+PI[j][i]]];
 177               minm=(*p2mymin)(minm,mm);
 178           }
 179           alpha[(k+1)*S+j]=minm;
 180           if(minm<norm) norm=minm;
 181       }
 182       for(int j=0;j<S;j++)
 183           alpha[(k+1)*S+j]-=norm; // normalize total metrics so they do not explode
 184   }
 185
 186   if(SK<0) { // final state not specified
 187       for(int i=0;i<S;i++) beta[K*S+i]=0;
 188   }
 189   else {
 190       for(int i=0;i<S;i++) beta[K*S+i]=INF;
 191       beta[K*S+SK]=0.0;
 192   }
 193
 194   for(int k=K-1;k>=0;k--) { // backward recursion
 195       norm=INF;
 196       for(int j=0;j<S;j++) {
 197           minm=INF;
 198           for(int i=0;i<I;i++) {
 199               int i0 = j*I+i;
 200               mm=beta[(k+1)*S+NS[i0]]+priori[k*I+i]+prioro[k*O+OS[i0]];
 201               minm=(*p2mymin)(minm,mm);
 202           }
 203           beta[k*S+j]=minm;
 204           if(minm<norm) norm=minm;
 205       }
 206       for(int j=0;j<S;j++)
 207           beta[k*S+j]-=norm; // normalize total metrics so they do not explode
 208   }
 209
 210
 211   if (POSTI && POSTO)
 212   {
 213     for(int k=0;k<K;k++) { // input combining
 214         norm=INF;
 215         for(int i=0;i<I;i++) {
 216             minm=INF;
 217             for(int j=0;j<S;j++) {
 218                 mm=alpha[k*S+j]+prioro[k*O+OS[j*I+i]]+beta[(k+1)*S+NS[j*I+i]];
 219                 minm=(*p2mymin)(minm,mm);
 220             }
 221             post[k*(I+O)+i]=minm;
 222             if(minm<norm) norm=minm;
 223         }
 224         for(int i=0;i<I;i++)
 225             post[k*(I+O)+i]-=norm; // normalize metrics
 226     }
 227
 228
 229     for(int k=0;k<K;k++) { // output combining
 230         norm=INF;
 231         for(int n=0;n<O;n++) {
 232             minm=INF;
 233             for(int j=0;j<S;j++) {
 234                 for(int i=0;i<I;i++) {
 235                     mm= (n==OS[j*I+i] ? alpha[k*S+j]+priori[k*I+i]+beta[(k+1)*S+NS[j*I+i]] : INF);
 236                     minm=(*p2mymin)(minm,mm);
 237                 }
 238             }
 239             post[k*(I+O)+I+n]=minm;
 240             if(minm<norm) norm=minm;
 241         }
 242         for(int n=0;n<O;n++)
 243             post[k*(I+O)+I+n]-=norm; // normalize metrics
 244     }
 245   }
 246   else if(POSTI)
 247   {
 248     for(int k=0;k<K;k++) { // input combining
 249         norm=INF;
 250         for(int i=0;i<I;i++) {
 251             minm=INF;
 252             for(int j=0;j<S;j++) {
 253                 mm=alpha[k*S+j]+prioro[k*O+OS[j*I+i]]+beta[(k+1)*S+NS[j*I+i]];
 254                 minm=(*p2mymin)(minm,mm);
 255             }
 256             post[k*I+i]=minm;
 257             if(minm<norm) norm=minm;
 258         }
 259         for(int i=0;i<I;i++)
 260             post[k*I+i]-=norm; // normalize metrics
 261     }
 262   }
 263   else if(POSTO)
 264   {
 265     for(int k=0;k<K;k++) { // output combining
 266         norm=INF;
 267         for(int n=0;n<O;n++) {
 268             minm=INF;
 269             for(int j=0;j<S;j++) {
 270                 for(int i=0;i<I;i++) {
 271                     mm= (n==OS[j*I+i] ? alpha[k*S+j]+priori[k*I+i]+beta[(k+1)*S+NS[j*I+i]] : INF);
 272                     minm=(*p2mymin)(minm,mm);
 273                 }
 274             }
 275             post[k*O+n]=minm;
 276             if(minm<norm) norm=minm;
 277         }
 278         for(int n=0;n<O;n++)
 279             post[k*O+n]-=norm; // normalize metrics
 280     }
 281   }
 282   else
 283     throw std::runtime_error ("Not both POSTI and POSTO can be false.");
 284
 285   delete [] prioro;
 286
 287 }
 288
 289
 290
 291
 292
 293
 294 int
 295 trellis_siso_combined_f::general_work (int noutput_items,
 296                         gr_vector_int &ninput_items,
 297                         gr_vector_const_void_star &input_items,
 298                         gr_vector_void_star &output_items)
 299 {
 300   assert (input_items.size() == 2*output_items.size());
 301   int nstreams = output_items.size();
 302   //printf("general_work:Streams:  %d\n",nstreams);
 303   int multiple;
 304   if (d_POSTI && d_POSTO)
 305       multiple = d_FSM.I()+d_FSM.O();
 306   else if(d_POSTI)
 307       multiple = d_FSM.I();
 308   else if(d_POSTO)
 309       multiple = d_FSM.O();
 310   else
 311       throw std::runtime_error ("Not both POSTI and POSTO can be false.");
 312
 313   assert (noutput_items % (d_K*multiple) == 0);
 314   int nblocks = noutput_items / (d_K*multiple);
 315   //printf("general_work:Blocks:  %d\n",nblocks);
 316   //for(int i=0;i<ninput_items.size();i++)
 317       //printf("general_work:Input items available:  %d\n",ninput_items[i]);
 318
 319   float (*p2min)(float, float) = NULL;
 320   if(d_SISO_TYPE == TRELLIS_MIN_SUM)
 321     p2min = &min;
 322   else if(d_SISO_TYPE == TRELLIS_SUM_PRODUCT)
 323     p2min = &min_star;
 324
 325
 326   for (int m=0;m<nstreams;m++) {
 327     const float *in1 = (const float *) input_items[2*m];
 328     const float *in2 = (const float *) input_items[2*m+1];
 329     float *out = (float *) output_items[m];
 330     for (int n=0;n<nblocks;n++) {
 331       siso_algorithm_combined(d_FSM.I(),d_FSM.S(),d_FSM.O(),
 332         d_FSM.NS(),d_FSM.OS(),d_FSM.PS(),d_FSM.PI(),
 333         d_K,d_S0,d_SK,
 334         d_POSTI,d_POSTO,
 335         p2min,
 336         d_D,d_TABLE,d_TYPE,
 337         &(in1[n*d_K*d_FSM.I()]),&(in2[n*d_K*d_D]),
 338         &(out[n*d_K*multiple])//,
 339         //d_alpha,d_beta
 340         );
 341     }
 342   }
 343
 344   for (unsigned int i = 0; i < input_items.size()/2; i++) {
 345     consume(2*i,d_FSM.I() * noutput_items / multiple );
 346     consume(2*i+1,d_D * noutput_items / multiple );
 347   }
 348
 349   return noutput_items;
 350 }