mp_dict.c

LastWave

/*
 * Change the dx of a dictionary : change dx for dict->channels as well as all sub-dictionaries.
 */
void SetDictDX(DICT dict,LWFLOAT dx)
{
  unsigned char channel;
  SIGNAL signal;
  unsigned short i;
  SUBDICT subDict;
  STFT stft;
  MAXIMADICT maximaDict;
  //WTRANS wtrans;

  if(dict == NULL)         Errorf("SetDictDX : NULL dict");
  if(dict->channels==NULL) Errorf("SetDictDX : NULL dict->channels");
  if (dx <= 0)             Errorf("SetDictDX : dx %g <= 0\n",dx);
  
  /* First, the dx of all channels */
  for(channel=0;channel<dict->nChannels;channel++) {
    signal = GetChannel(dict,channel);
    signal->dx = dx;
  }
  /* Then, that of the sub-dictionaries */
  for(i = 0; i < dict->size; i++) {
    subDict = dict->subDicts[i];
    /* Case of STFT subDicts */
    if(GetTypeValue(subDict->dataContainer)==stftType) {
      stft=(STFT)(subDict->dataContainer);
      stft->dx = dx;
      continue;
    } 
    /* Case of MAXIMADICT subDicts */
    if(GetTypeValue(subDict->dataContainer)==maximaDictType) {
      maximaDict=(MAXIMADICT)subDict->dataContainer;
      Errorf("SetMaximaDictDX not implemented");
      //SetMaximaDictDX(maximaDict,dx); TODO !
      continue;
    }
    Errorf("SetDictDX : (Weired)");
  }
}

/* Change the dx of a dictionary : change x0 for dict->channels as well as all sub-dictionaries. */
void SetDictX0(DICT dict,LWFLOAT x0)
{
  unsigned char channel;
  SIGNAL signal;
  unsigned short i;
  SUBDICT subDict;
  STFT stft;
  MAXIMADICT maximaDict;
  //WTRANS wtrans;
  
  /* Checking arguments */
  if(dict == NULL)         Errorf("SetDictX0 : NULL dict");
  if(dict->channels==NULL) Errorf("SetDictX0 : NULL dict->channels");

  /* First, the x0 of all channels */
  for(channel=0;channel<dict->nChannels;channel++) {
    signal = GetChannel(dict,channel);
    signal->x0 = x0;
  }

  /* Then, that of the sub-dictionaries */
  for(i = 0; i < dict->size; i++) {
    subDict = dict->subDicts[i];
    /* Case of STFT subDicts */
    if(GetTypeValue(subDict->dataContainer)==stftType) {
      stft=(STFT)(subDict->dataContainer);
      stft->x0 = x0;
      continue;
    } 
    /* Case of MAXIMADICT subDicts */
    if(GetTypeValue(subDict->dataContainer)==maximaDictType) {
      maximaDict=(MAXIMADICT)subDict->dataContainer;
      Errorf("SetMaximaDictX0 not implemented");
      //SetMaximaDictX0(maximaDict,x0); TODO!
      continue;
    }
    Errorf("SetDictX0 : (Weired)");
  }
}

/*************************************************/
/*
 *	The main functionalities of a dictionary
 */				    
/*************************************************/
/*
 * Update all the 'main' sub-dictionaries, and the necessary 'aux' sub-dictionaries,
 * using the 'removedMolecule' if necessary, and resets [updateTimeIdMin,updateTimeIdMax].
 * The state of the dictionary after a call to this function is as follows :
 * 0/ removedMolecule = NULL.
 * 1/if there is no sub-dictionary of local maxima : all sub-dictionaries are up to date
 * and [updaTimeIdMin,updateTimeIdMax] = [dict->channels[0]->size-1,0] 
 * (when a molecule is removed from the dictionary, updateTimeIdMin/Max decreases/increases)
 * 2/if there is a sub-dictionary of local maxima  : this sub-dictionary is up to date,
 * and the state of [updaTimeIdMin,updateTimeIdMax] depends on whether the update of the 
 * sub-dictionary of local maxima involved an 'initialization' or not.
 */
void UpdateDict(DICT dict) {
  unsigned short i;
  SUBDICT subDict;

  /* If already up to date, do nothing */
  if(dict->updateTimeIdMin > dict->updateTimeIdMax) return;
  
  /* When there is a sub-dictionary of local maxima */
  if(subDict=GetMaximaDictSubDict(dict)) {
    /* This takes care of the proper management of [updateTimeIdMin,updateTimeIdMax] */
    UpdateSubDict(subDict);
  }
  /* When there is no sub-dictionary of local maxima */
  else {
    for(i = 0; i < dict->size; i++) {
      subDict = dict->subDicts[i];
      /* Display a clock for waiting */
      switch(i%4) {
//      case 0 : Printf("\\"); Flush(); break;
//      case 1 : Printf("|"); Flush(); break;
//      case 2 : Printf("/"); Flush(); break;
//      case 3 : Printf("-"); Flush(); break;
      }
      UpdateSubDict(subDict);
//      Printf("\b");
    }
    dict->updateTimeIdMin = dict->channels[0]->size-1;
    dict->updateTimeIdMax = 0;
  }
  /* Deletes the removedMolecule if exists */
  if(dict->removedMolecule) dict->removedMolecule = DeleteMolecule(dict->removedMolecule);
}

/*
 * Returns a MOLECULE corresponding to the maximum of a DICT over a search range.
 * If the specified search range is empty or the  maximum value is zero, we return NULL.
 * The returned MOLECULE is not optimized (interpolate,chirped,...).
 * The maximum is looked for only in 'main' sub-dictionaries, which should be up to date
 * (else an error is generated).
 * An error is also generated if [updateTimeIdMin,updateTimeIdMax] is not empty in the case
 * when there is no sub-dictionary of local maxima.
 */
MOLECULE GetMaxDict(DICT dict,LISTV searchRange) 
{
  unsigned short i;
  SUBDICT subDict,subDictMax=NULL;
  static LISTV tmpParamMax = NULL;
  LISTV paramMax = NULL;
  char  maxFound = NO;
  LWFLOAT tmpMaxValue = 0.0,maxValue = 0.0;
  MOLECULE molecule = NULL;
  MAXIMADICT maximaDict;

  /* Allocate (once only) the list of parameters that will be used     */
  /* to locate quickly the &mol after the first exhaustive search */
  if(tmpParamMax == NULL)     tmpParamMax = NewListv();
  else                        ClearListv(tmpParamMax);

  /* First, locate the subdict that contains the maximum. */
  if((subDictMax=GetMaximaDictSubDict(dict))!=NULL) {
    /* 
     * If there is a sub-dictionary of local maxima, this is the only sub-dictionary where to look.
     * and there is a maximum as soon as it is not empty 
     */
    maximaDict = (MAXIMADICT)(subDictMax->dataContainer);
    if(maximaDict->nMaxima>0) {
      maxFound = YES;
    }
  } else {
    /*
     * Check that[updateTimeIdMin,updateTimeIdMax] is empty
     * (an individual check of each sub-dictionary being up to date will be done during the search for the max)
     */
    if(dict->updateTimeIdMin <= dict->updateTimeIdMax) 
      Errorf("GetMaxDict : dict should be updated first in [%d %d]",dict->updateTimeIdMin,dict->updateTimeIdMax);
    for(i = 0; i < dict->size; i++) {
      subDict = dict->subDicts[i];
      /* Skip auxiliary subDicts */
      if(subDict->flagMain==NO) continue;
      /* If the search range is empty for the subDict, we skip to the next subDict */
      if(GetMaxSubDict(subDict,searchRange,&tmpMaxValue,tmpParamMax)==NO) continue;
      /* The search range was non empty, hence we can compare the local maximum to the best we have so far */
      if(maxFound == NO || tmpMaxValue > maxValue) {
	maxValue = tmpMaxValue;
	paramMax = TNewListv();
	CopyListv(tmpParamMax,paramMax);
	subDictMax = subDict;
      }
      maxFound = YES;
    }
  } 

  /* If no maximum has been found (i.e. the search range was empty or the max was zero) we return NULL */
  if(maxFound == NO) return(NULL);

  /* Now we have to fill the maximum variable */
  molecule = TNewMolecule();
  GetMaxSubDict(subDictMax,paramMax,NULL,molecule);
  return(molecule);
}

// TODO : clean this up
void AddShiftedSignals(SIGNAL input,SIGNAL output,long shift)
{
  long iOut,iOutMin,iOutMax;
  
  /* Some Tests */
  if (output == input)   Errorf("AddShiftedSignals : you must have output != input\n");
  if (output->size == 0) Errorf("AddShiftedSignals : output should not be empty\n");
  
  iOutMin = MAX(shift,0);
  iOutMax = MIN(output->size,shift+input->size);
  
  for (iOut = iOutMin; iOut < iOutMax; iOut++)
    output->Y[iOut] += input->Y[iOut-shift];
}


static void SubstractShiftedSignals(SIGNAL input,SIGNAL output,long shift)
{
  long iOut,iOutMin,iOutMax;

  /* Some Tests */
  if (output == input)   Errorf("SubstractShiftedSignals : you must have output != input\n");
  if (output->size == 0) Errorf("SubstractShiftedSignals : output should not be empty\n");

  iOutMin = MAX(shift,0);
  iOutMax = MIN(output->size,shift+input->size);
  
  for (iOut = iOutMin; iOut < iOutMax; iOut++)
    output->Y[iOut] -= input->Y[iOut-shift];
}


/*
 * Builds an atom and removes it from a given channel of a dictionary.
 * The [updateTimeIdMin,updateTimeIdMax] of the dict is enlarged if necessary.
 * All sub-dictionaries are marked as out of date.
 */
static void RemoveAtomFromDictChannel(DICT dict,unsigned char channel,ATOM atom) 
{
  static ATOM copy = NULL;
  static SIGNAL atomSignal = NULL;
  // These are signed integers for they may be negative
  // TODO : change to unsigned ?
  int atomTimeIdMin,atomTimeIdMax;

  LWFLOAT shift;

  // TODO : What borderType do we really want ? Shall we store it in dict ?
  char borderType = BorderPad0;
  LWFLOAT real1,imag1,real2,imag2;
  /* Checking arguments */
  CheckDictNotEmpty(dict);
  CheckAtom(atom);

  /* Allocating (just once) */
  if(atomSignal == NULL) atomSignal = NewSignal();
  
  // DEBUG : COMPRENDRE POURQUOI CA ARRIVE
  copy = CopyAtom(atom,copy);
  SCAtomInnerProduct(GetChannel(dict,channel),copy,borderType,&(copy->coeffR),&(copy->coeffI));
  CastAtomReal(copy);
  if(fabsf((atom->coeffR-copy->coeffR)/atom->coeffR)>1e-5 || fabsf((atom->coeffI-copy->coeffI)/atom->coeffI)>1e-5) {
    //    Printf("FreqId %g Phase %g %g\n",atom->freqId,atom->cosPhase,atom->sinPhase);
    //    Printf("Given (%g,%g) (%g,%g,%g) [%g,%g]\n",atom->coeffR,atom->coeffI,atom->coeff2,atom->cosPhase,atom->sinPhase,atom->realGG,atom->imagGG);
    //    Printf("Computed (%g,%g) (%g,%g,%g) [%g,%g]\n",copy->coeffR,copy->coeffI,copy->coeff2,copy->cosPhase,copy->sinPhase,copy->realGG,copy->imagGG);
  }
  /* Build a real valued atom  */
  BuildAtom(atom,atomSignal,NULL,borderType,YES);
  /* Substract it to the dict->signal */
  // TODO : check the shift is integer! and clean this up anyway
  shift = (atomSignal->x0-atom->x0)/atom->dx;
  if (shift < 0) shift = (int) (shift-.5);
  else shift = (int) (shift+.5);
  SubstractShiftedSignals(atomSignal,GetChannel(dict,channel),(long) shift);
    
  /* Compute the limits where the signal has changed */
  ComputeWindowSupport(atom->windowSize,atom->windowShape,atom->timeId,&atomTimeIdMin,&atomTimeIdMax);
  dict->updateTimeIdMin = MIN(dict->updateTimeIdMin,atomTimeIdMin);
  dict->updateTimeIdMax = MAX(dict->updateTimeIdMax,atomTimeIdMax);

  // TODO : treat other border effects (pass it as an argument to the function ?)
  // TODO : clean this up !!!
  switch(borderType) {
  case BorderPad0 :
    if(atomTimeIdMin < 0 || dict->updateTimeIdMax >= dict->signalSize)
      Errorf("RemoveAtomFromDictChannels : signal updated out of its limits!");
    break;
  default :
    Errorf("RemoveAtomFromDictChannels : border type '%s' not treated yet",BorderType2Name(borderType));
  }
}

/*
 * Builds a molecule and removes it from the signal of a dictionary.
 * The signalEnergy is updated. An error is generated if it does not decrease.
 * The [updateTimeIdMin,updateTimeIdMax] of the dict is enlarged if necessary.
 * All sub-dictionaries are marked as out of date.
 * The molecule is memorized in the dict->removedMolecule field for possible
 * use at the next UpdateDict step.
 * If a molecule has already been memorized, an error is generated.
 */
void RemoveMoleculeFromDict(DICT dict,MOLECULE molecule) 
{
  unsigned char channel;
  unsigned short k;
  SIGNAL signal;
  ATOM atom;
  unsigned long i;
  LWFLOAT prevSigEnergy;

  SUBDICT subDict;
  // TODO : remove that when possible
  STFT stft;

  /* No other molecule should have been memorized */
  if(dict->removedMolecule) Errorf("RemoveMoleculeFromDict : a molecule has already been removed");

  
  for(channel = 0; channel < molecule->nChannels; channel++) {
    for(k = 0; k < molecule->dim; k++) {
      atom = GetMoleculeAtom(molecule,channel,k);
      RemoveAtomFromDictChannel(dict,channel,atom);
    }
  }

  /* Sub-dictionaries of the dict are now out of date */
  for(i = 0; i < dict->size; i++) {
    subDict = dict->subDicts[i];
    subDict->flagUpToDate = NO;
    // TODO : remove that when possible
    if(GetTypeValue(subDict->dataContainer)==stftType) {
      stft=(STFT)(subDict->dataContainer);
      stft->flagUpToDate=NO;
    }
  }

  /* Update the signalEnergy */
  // TODO fast update using [updateTimeIdMin updateTimeIdMax] ???
  prevSigEnergy = dict->signalEnergy;
  dict->signalEnergy = 0.0;
  for(channel=0; channel<dict->nChannels;channel++) {
    signal = GetChannel(dict,channel);
    for(i = 0; i < dict->signalSize; i++)
      dict->signalEnergy += signal->Y[i]*signal->Y[i];
  }
  /* The energy should strictly decrease */
  // TODO : remove when we are confident enough this will never happen
  if(isnan(dict->signalEnergy) || prevSigEnergy<dict->signalEnergy) {
    PrintInfoValue(molecule);
    Errorf("RemoveMoleculeFromDict : (WEIRED) increases the signalEnergy from %g to %g",prevSigEnergy,dict->signalEnergy);
  } 

  /* Memorize the molecule */
  dict->removedMolecule = molecule;
  AddRefValue(molecule);  

}

/*
 * If dict->channels is already set or the new signal(s) have inconsistent size/dx/x0, an error is generated
 * The [updateTimeIdMin,updateTimeIdMax] is set to [0,dict->signalSize-1]
 * The signalEnergy is initialized.