/***************************************************************************
                          Melody.cpp  -  description
                             -------------------
    begin                : Don Mai 8 2003
    copyright            : (C) 2003 by klaus frieler
    email                : kgf@omniversum.de
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
#include <cmath>
#include <time.h>
#include <cstdlib>
#include <list>
#include "../libmelody/ZN.h"
#include "Melody.h"
#include "Types.h"
#include "EditDistance.cpp"
#include "Fourier.h"
#include "gauss.h"
#include "Scales.h"
#include "MetricTree.h"
#include "VectorManip.cpp"
#include "Ngrams.cpp"
#include "BarTime.h"
#include "Misc.h"
#include "freq.cpp"

static const char* genToneNames[] = { "C", "C#", "D", "Eb", "E", "F", "F#", "G", "G#", "A", "Bb", "B"};

#define sgn(x) (abs((x))/(x))

void dumpMelody(const char* text, PitchSequence &mel)
{
   uint pC = 0;
   if(!text) return;
   cout<<text<<endl;
   for(uint i = 0; i<mel.size(); ++i){
      pC = mel[i];
      if(pC>= C0) pC = (pC-C0)%12;
      else pC = (12+(pC-C0))%12;
      cout<<genToneNames[pC]<<" ";      
   }
   cout<<endl;

}

uint Melody::nGramLength = 4;

void Melody::getPitchRange(MIDIPitch &low, MIDIPitch &high) const
{
   low = 127; high = 0;
   for (uint i = 0; i<m_pitches.size(); ++i){
      if(m_pitches[i]<low) low = m_pitches[i];
      if (m_pitches[i]>high) high = m_pitches[i];
   }
}

PitchSequence Melody::getPitches(uint idx)
{
   PitchSequence mel;
   for (uint i = 0; i<m_pitches.size(); ++i){
      if(i>=idx)
         mel.push_back(m_pitches[i]);
   }
   return mel;   
}

data Melody::getPitchesAsDouble(uint idx) const
{
   data mel;
   for (uint i = 0; i<m_pitches.size(); ++i){
      if(i>=idx)
         mel.push_back((double)m_pitches[i]);
   }
   return mel;
}

uint Melody::getLength() const
{
   uint sum = 0;
   for(uint i = 0; i< m_durations.size(); ++i){
      sum += m_durations[i];
   }
   return sum;
}

double Melody::getMeanPitch() const
{
   return getMeanPitch(m_pitches);   
}

double Melody::getMeanPitch(const PitchSequence& pitches) const
{
   double sum = 0;
   uint len = pitches.size();
   for(uint i = 0; i< len; ++i){
      sum += pitches[i];
   }
   if(len) sum = sum / (double)len;

   return sum;
}
durations Melody::getDurations() const
{
   return m_durations;
}
data Melody::getDurationsAsDoubles() const
{
   data ret;
   for(uint i = 0; i<m_durations.size(); ++i){
      ret.push_back(m_durations[i]);
   }
   return ret;
}

vector<uint> Melody::getDurationModus(uint &modFreq) const
{
   return getModus<uint>(m_durations, modFreq);
}

vector<MIDIPitch> Melody::getPitchModus(uint &modFreq) const
{
   MIDIPitch low, high;
   getPitchRange(low, high);
   return getModus<MIDIPitch>(m_pitches, modFreq);
}

Melody Melody::normalize() const
{
   Melody ret;
   uint i = 0;
   
   //out()<<m_pitches.size()<<" "<<m_durations.size()<<" "<<m_timeLine.size()<<endl;
   if(m_pitches.size() != m_durations.size() ||  m_pitches.size() != m_timeLine.size()){
      out()<<m_pitches.size()<<" "<<m_durations.size()<<" "<<m_timeLine.size()<<endl;
      out()<<"Invalid Melody class object ("<<m_sName<<")"<<endl;
      return ret; 
   }
   for(i = 0; i< m_pitches.size(); ++i){
      ret.m_timeLine.push_back(m_timeLine[i]);
      for(uint j = 0; j<m_durations[i]; ++j){
         ret.m_pitches.push_back(m_pitches[i]);
         ret.m_durations.push_back((uint)1);                 
      }
   }
   ret.normalizeTimeLine();
   BarPosition barPos = ret.m_timeLine[0];
   ret.m_timeLine.clear();
   ret.m_timeLine.push_back(barPos);
   //m_signature.dump();
   uint timeBase = m_timeLine.getTimebase();
   if(ret.m_pitches.size()<2) return ret;
   for(i = 0; i< ret.m_pitches.size()-1; ++i){
      //barPos.dump();
      //barPos.addSemiQuaver(m_signature);
      barPos.addTicks(timeBase);
      ret.m_timeLine.push_back(barPos);
      //barPos.dump();
   }
   return ret;
}

Melody& Melody::transpose(const char semitone)
{
   for(uint i = 0; i< m_pitches.size(); ++i){
      m_pitches[i] += semitone;
   }
   return *this;
}

data Melody::contourize(uint mode) const
{
   data mel;
   indizes x;
   uint i = 0, j = 0;
   PitchSequence y;
   if(mode & STEINBECK){
      getContourPoints(x, y, eCont_steinbeck);
   }
   else if(mode & NATURAL){
      getContourPoints(x, y, eCont_all);
   }
   else {
      getContourPoints(x, y, eCont_classic);
   }
   for(i = 0, j =0; i<m_pitches.size()-1; ++i){
      if(i == x[j]) j++;
      double y1 = (double)y[j-1], y2 = (double)y[j];
      double x1 = (double)x[j-1], x2 = (double)x[j];
      double val = (y2-y1)/(x2-x1)*(i - x1) +y1;
      //out()<<"y1 = "<< y1<<" y2 = "<<y2<<" x1 = "<<x1<<" x2 = "<<x2<<endl;
      //out()<<"Pitch ["<<i<<"]: "<<(uint)m_pitches[i]<<" -> "<<val<<endl;
      mel.push_back(val);
   }
   mel.push_back(m_pitches[m_pitches.size()-1]);
   //out()<<"Pitch: "<<(uint)m_pitches[m_pitches.size()-1]<<" -> "<<(uint)m_pitches[m_pitches.size()-1]<<endl;
   if(mode & INTERVALS){
      data mel1;
      for(i = 0; i<mel.size()-1; ++i){
         mel1.push_back(mel[i+1]-mel[i]);
      }
      mel = mel1;
      //dumpData("Interval contour:", mel1);
   }
   return mel;
}


Melody Melody::intervals() const
{
   Melody ret;
   if(m_pitches.size() != m_durations.size()){
      out()<<"Invalid Melody class object."<<endl;
      return ret;
   }
   for(uint i = 0; i< m_pitches.size()-1; ++i){
      ret.m_pitches.push_back(m_pitches[i+1] - m_pitches[i]);
      ret.m_durations.push_back(m_durations[i]);
   }
   return ret;
}

uint Melody::getLeftOffset()
{
   uint len = getLength();
   uint maxLen = len/10 > 8 ? 8: len/10;
   uint sum = 0;
   
   for( uint i = 0; i<getSize(); ++i){
      sum += m_durations[i];
      if(sum>maxLen) return i;
   }

   return getSize()-1;
}

list<char> Melody::getTranspositionHints(const Melody& melody)
{
   list<char> trans;
   list<char>::iterator iter;
   uint modFreq = 0;
   vector<MIDIPitch> modus1 = getPitchModus(modFreq);
   vector<MIDIPitch> modus2 = melody.getPitchModus(modFreq);
   MIDIPitch mean1  = (MIDIPitch)floor(0.5+getMeanPitch());
   MIDIPitch mean2  = (MIDIPitch)floor(0.5+melody.getMeanPitch());
   MIDIPitch low1 = 0, high1 = 0;
   MIDIPitch low2 = 0, high2 = 0;

   
   //printMelody();
   //melody.printMelody();
   getPitchRange(low1, high1);
   melody.getPitchRange(low2, high2);
   if(low1 == high1 || low2 == high2){
      if(low1 == 0 || low2 == 0){
         trans.push_back(0);
         return trans;
      }
   }
   for(uint i = 0; i<modus1.size(); ++i){
      for(uint j = 0; j<modus2.size(); ++j){
         trans.push_back(modus1[i]-modus2[j]);
      }
   }
   trans.push_back(mean1-mean2);

   trans.push_back(low1-low2);
   trans.push_back(high1-high2);
   trans.sort();
   trans.unique();
   /*
   for(iter = trans.begin(); iter != trans.end(); ++iter){
      out()<<(int)*iter<<" ";
   }
   out()<<endl;
   */
   return trans;
}

char Melody::getBestTransposition(const Melody& melody)
{
   char trans = 0;
   double dummy = 0;
   list<char> transHints;
   double corr = 0, maxCorr = 0;
   double mean1 = getMeanPitch(), mean2  = melody.getMeanPitch();
   transHints = getTranspositionHints(melody);
   //return trans;//tmp
   for(char i = *(transHints.begin()); i<=*(transHints.rbegin()); ++i){
      if(fabs(mean2+i-mean1)>2.5) continue;
      Melody tmp(melody);
      //cout<<"Testing : "<<(int)i<<endl;
      tmp.transpose(i);
      //tmp.printMelodyNames(cout);
      corr = 0.5*(similarity(tmp, SimType(EDIT_DISTANCE,_VOID_), dummy)
                  + similarity(tmp, SimType(EDIT_DISTANCE|WEIGHTED), dummy));
      if(corr>maxCorr) {
         trans = i;
         maxCorr = corr;
      }

   }
   /*
   data dummyvec;
   uint ton1 = getRawTonalityKrumhansl(-1, dummyvec);
   uint ton2 = melody.getRawTonalityKrumhansl(-1, dummyvec);
   iData key1 = getKeys(0);
   iData key2 = melody.getKeys(0);
   cout<<"Keys 1:";
   for(uint i = 0; i<key1.size(); ++i){
      cout<<getKeyName(key1[i]).c_str()<<" ";
   }
   cout<<endl;
   cout<<"Keys 2:";
   for(uint i = 0; i<key2.size(); ++i){
      cout<<getKeyName(key2[i]).c_str()<<" ";
   }
   cout<<endl;
   cout<<"Melody1 = "<<getKeyName(ton1).c_str()<<" Melody2 = "<<getKeyName(ton2).c_str()<<" proposed transposition = "<<int(trans)<<endl;
   */
   return trans;
}

iData Melody::getKeys(uint mode) const
{
   iData harm;
   uint maxIdx = 0;
   double max = 0;
   uint nBar = getBarCount();
   data tonCache1, tonCache2, tonCurr;
   getRawTonalityKrumhansl(1, tonCache1);
   getMax<double>(tonCache1, maxIdx);
   if(nBar < 2){
      getMax<double>(tonCache1, maxIdx);
      harm.push_back(maxIdx);
      return harm;    
   }
   getRawTonalityKrumhansl(2, tonCache2);
   getMax<double>(tonCache2, maxIdx);

   data tmp = add<double>(mult<double>(tonCache1, 2.0), tonCache2);
   getMax<double>(tmp, maxIdx);
   harm.push_back(maxIdx);
   for(uint i = 2; i<nBar; ++i){      
      getRawTonalityKrumhansl(i+1, tonCurr);
      tmp = add<double>(tonCurr, add<double>(tonCache1, mult<double>(tonCache2, 2.0)));
      getMax<double>(tmp, maxIdx);
      harm.push_back(maxIdx);
      tonCache1 = tonCache2;
      tonCache2 = tonCurr;
   }
   if(nBar == 2){
      tmp = add<double>(mult<double>(tonCurr, 2.0), tonCache2);
      getMax<double>(tmp, maxIdx);
      harm.push_back(maxIdx);
   }
   if(nBar>2){
      tmp = add<double>(mult<double>(tonCurr, 2.0), tonCache2);
      getMax<double>(tmp, maxIdx);
      harm.push_back(maxIdx);
   }
   cout<<"Harmonies: ";
   for(uint i = 0; i<harm.size(); ++i){
      cout<<getKeyName(harm[i]).c_str()<<" ";
   }
   cout<<endl;
   uint dummy;
   iData ret = getModus<uint>(harm, dummy);
   return ret;
}
double  Melody::similarity(const Melody& melody, const CombiType& combiType, double &tmp)const
{
   CombiType::const_iterator   ct_iter;
   SimList::const_iterator     sl_iter;

   double sum = 0.0;
   
   for(ct_iter = combiType.begin(); ct_iter != combiType.end(); ++ct_iter){
      double prod = 1;
      for(sl_iter = ct_iter->begin(); sl_iter != ct_iter->end(); ++sl_iter){
         double sim = similarity(melody, *sl_iter, tmp);
         if(sim == -1) return -1.0; //invalid sim type!!
         prod *= sim;
      }
      sum += prod*(ct_iter->m_dWeight);
   }
   if(sum>1.0) sum = 1.0;
   if(sum<0.0) sum = 0.0;
   
   return sum;  
}

double  
Melody::similarity(const Melody& melody,  const SimType& simType, double &tmp)const
{
   EditDistance<data> editDistance;
   Melody mel1(*this), mel2(melody);
   uint lof1 = 0, lof2 = 0;
   data   dMel1, dMel2;
   double corr = -1.0;
   uint len1 = 0 , len2 = 0, maxLen = 0, diff = 0;
   const uint mode = simType.m_iMode;
   const uint submode = simType.m_iSubmode;
   const uint accentType = simType.m_iAccentFlags;
   tmp = 0;
   if(mel1.empty() || mel2.empty()){
      //out()<<"Found invalid melody!"<<endl;
      return -1.0;
   }
   if(mode & ACCENTS){
      return getAccentSimilarity(melody, mode, accentType);
   }

   if(mode & CONSTSIM)
      return 1.0;
   if(mode & WEIGHTED){
      mel1 = mel1.normalize();
      mel2 = mel2.normalize();
      if(mel1.empty() || mel2.empty()) return -1.0;
   }
   if(mode & HARM_CORR){
      return mel1.harmonicCorrelation(mel2, mode, submode, tmp);
   }
   lof1 = mel1.getLeftOffset();
   lof2 = mel2.getLeftOffset();
   
   if(mode & CONTOUR){
      dMel1 = mel1.contourize(mode);
      dMel2 = mel2.contourize(mode);
      //dumpData("Contour: ", dMel1);
      //dumpData("Contour: ", dMel2);
   }
   else
      if(mode & RANKS){
         dMel1 = mel1.getRanks();
         dMel2 = mel2.getRanks();
         //dumpData("Ranks: ", dMel1, out());
         //dumpData("Ranks: ", dMel2, out());
      }

   else
      if(mode & INTERVALS && !(mode & FOURIER)){
         dMel1 = mel1.getIntervals(false, mode);
         dMel2 = mel2.getIntervals(false, mode);
   }
   else
      if(mode & RHYTHM){
         if(mode & FUZZY){
            return fuzzyRhythmCorr(mel2, tmp, mode, submode);
         }
         dMel1 = mel1.m_onsets;
         dMel2 = mel2.m_onsets;
      }
   else{
         dMel1 = mel1.getPitchesAsDouble();
         dMel2 = mel2.getPitchesAsDouble();

   }
      
   len1  = dMel1.size();
   len2  = dMel2.size();
   maxLen = len1<len2? len2 : len1;
   diff = len1<len2? len2-len1 : len1-len2;

   if(mode & EDIT_DISTANCE){
      tmp = (double) editDistance(dMel1, dMel1.size(), dMel2, dMel2.size());      
      /*
      SpikeCost< pair<double, uint> > spikeCost(0.5, 0.1);
      EditDistance2< pair<double, uint> > ed(&spikeCost);
      vector <pair<double, uint> > onsetPitch1, onsetPitch2;
      onsetPitch1 = mel1.getAsOnsetPitchPairs();
      onsetPitch2 = mel2.getAsOnsetPitchPairs();
      double tmp = (double)   ed(onsetPitch1, onsetPitch2);
      double maxDist = ed.getMaxDist();
      //out()<<"ed = "<<tmp<<" maxDist = "<<maxDist<<endl;
      if(maxDist) corr = tmp/maxDist;
      else corr = 1;
      */
      corr = tmp /maxLen;
      return 1-corr;
   }

   if(len1<len2)
      corr = getMaxCorrelation(dMel1, dMel2, lof1, tmp, mode, submode);
   else
      corr = getMaxCorrelation(dMel2, dMel1, lof2, tmp, mode, submode);   

   return corr;
}
double  Melody::phraseSimilarity(const Melody& melody, const CombiType& combiType)const
{
   MelodyList mel1 = getPhrases();
   MelodyList mel2 = melody.getPhrases();
   SimilarityCost simCost(combiType);
   uint l = combiType.size();
   EditDistance2<PMelody>   ed(&simCost);
   if(mel1.empty() || mel2.empty()) return -1;
   double ret = ed(mel1, mel2);
      //delete mel2[i];
   uint max = mel1.size()> mel2.size() ? mel1.size() : mel2.size();
   ret /= (double) max;
   return 1-ret;
}

string Melody::getKeyName(const uint maxIdx) const
{
   string ret;
   if(maxIdx<24) {
      ret = genToneNames[maxIdx % 12];
      if(maxIdx>11) ret +="m";
   }
   return ret;
}
void Melody::printRatio(ostream& out, const uint count, const uint denom) const
{
   uint gcd = ggT(count, denom);
   cout<<count/gcd<<"/"<<denom/gcd;
}
void Melody::printInfo(ostream& out) const
{
   MIDIPitch low, high;
   getPitchRange(low, high);
   uint modFreq = 0;
   PitchSequence tmp = getPitchModus(modFreq);
   out<<endl;
   out<<"Melody "<<m_sName.c_str()<<endl;
   out<<"-----------------------------------"<<endl;
   out<<"No. of pitches        = "<<getSize()<<endl;
   out<<"Total Duration (1/16) = "<<getLength()<<endl;
   out<<"Signature:            = "<<m_signature.m_nCount<<"/"<<m_signature.m_nDenom<<endl;
   out<<"Tatum:                = ";
   printRatio(out, m_timeLine.getTimebase(), 4*m_signature.m_nTicks);
   out<<endl;
   out<<"Number of Bars        = "<<getBarCount()<<endl;
   out<<"Lowest pitch          = "<<(uint)low<<endl;
   out<<"Highest pitch         = "<<(uint)high<<endl;
   out<<"Mean pitch            = "<<getMeanPitch()<<endl;
   out<<"Most frequent pitches = ";
   for(uint i = 0; i<tmp.size(); ++i){
      out<<(int)tmp[i]<<" ";
   }
   out<<" (Freq = "<<modFreq<<")"<<endl;
   printMelodyNames(out);

}

void Melody::printMelody(ostream& out) const
{
   out<<"Melody "<<m_sName.c_str()<<endl;
   for(uint i = 0; i<m_pitches.size(); ++i){
      out<<(int)m_pitches[i]<<"("<<m_durations[i]<<") ";
   }
   out<<endl;
}

uint Melody::getBarCount() const
{
   uint count = m_timeLine.size();
   if(m_timeLine.size() == 0) return 0;
   return m_timeLine[m_timeLine.size()-1].getBar()- m_timeLine[0].getBar() + (m_timeLine[0].getBar()!=0 ? 1 : 0);
}
uint Melody::getNonEmptyBarCount() const
{
   uint count = m_timeLine.size();
   if(count == 0) return 0;
   uint sum = 1;
   uint oldBar = m_timeLine[0].getBar();
   for(uint i = 1; i<m_timeLine.size(); ++i){
      if(oldBar != m_timeLine[i].getBar()){
         oldBar = m_timeLine[i].getBar();
         sum++;
      }
   } 
   return sum;
}

void Melody::insert(uint idx, MIDIPitch pitch, uint dur)
{
   durations::iterator dIt = m_durations.begin();
   PitchSequence::iterator mIt = m_pitches.begin();

   if(idx>= m_pitches.size()){
      append(pitch, dur);
      return;
   }
   for(uint i = 0; i<idx; ++i){
      ++dIt;
      ++mIt;
   }
   m_pitches.insert(mIt, pitch);

   m_durations.insert(dIt, dur);
   
}

void Melody::append(MIDIPitch pitch, uint dur, uint count)

{
   for(uint i = 0; i< count; ++i){
      m_pitches.push_back(pitch);
      m_durations.push_back(dur);
   }
}


void Melody::zeroPad(uint newLen, uint dur)
{
   int diff  = (int)newLen - getSize();
   if(diff <= 0) return;
   append(0, dur, diff);
}

double Melody::absNorm(const data& pitches1,const data& pitches2)
{
   uint len1 = pitches1.size(), len2 = pitches2.size();
   double sum = 0, diff = 0;
   if(len1 != len2 || len1 == 0) return 0;
   for(uint i = 0; i<len1; ++i){
      diff = fabs(pitches1[i]-pitches2[i]);
      sum += diff;
   }
   return sum;
}

double Melody::absNormNorm(const data& pitches1,const data& pitches2, const uint mode)const
{
   uint i = 0;
   uint len1 = pitches1.size(), len2 = pitches2.size();
   double sum = 0, diff = 0, max = 0, norm = 0.0, norm2 = 0.0, max2 = 0.0;
   double mean = 0;
   if(len1 != len2 || len1 == 0) return 0;
   for(i = 0; i<len1; ++i){
      diff = fabs(pitches1[i]-pitches2[i]);
      sum  += diff;
   }
   mean = sum/len1;
   sum = 0;

   for(i = 0; i<len1; ++i){
      diff = fabs(pitches1[i]-pitches2[i]);
      sum  += diff;
      norm += (diff-mean)*(diff-mean);
      norm2 = fabs(pitches1[i])+fabs(pitches2[i]);
      if(diff>max) max = diff;
      if(norm2>max2) max2 = norm2;
   }
   //dumpData("Intervals1", pitches1);
   //dumpData("Intervals2", pitches2);
   //out()<<"Absnormnorm: sum="<<sum<<" var ="<<sqrt(norm)<<" max= "<<max<<" mean= "<<mean<<" len="<<len1<<endl;
   if(mode & EXPONENTIAL) return exp(-sum/len1);

   if(mode & MAXNORM2 ) 
      if(max2) return 1-sum/(len1*max2);
      else return 0;
   norm /= len1;
   norm = max+sqrt(norm);

   if(norm)   return 1-sum/(len1*norm);
   return 1;
}
double Melody::getFuzzyVal(double val, INTERVALFEATURE ivalF) const
{
   double x1 = 0.0, x2 = 0.0, y1 = 0.0, y2 = 0.0;

   if(val<ivalF.left || val>ivalF.right){
      return 0;
   }

   if(ivalF.left<=val && val<= ivalF.mid1){
      y1 = 0;
      y2 = 1;
      x1 = ivalF.left;
      x2 = ivalF.mid1;
   }
   if(ivalF.mid1<=val && val<=ivalF.mid2){
      return 1.0;

   }

   if(ivalF.mid2<=val && val<= ivalF.right){
      y1 = 1;
      y2 = 0;
      x1 = ivalF.mid1;
      x2 = ivalF.right;
    }
   return (y2-y1)/(x2-x1)*(val - x1) +y1;
}

double Melody::getFuzzyFeatureCorr(const data& feat1, const data& feat2)const
{
   return correlation(feat1, feat2, 0);
}

data Melody::getFuzzyFeature(const double interval) const
{
   static INTERVALFEATURE ivalFs[] ={
      {-4, -1000, -24, -8,-7},
      {-3, -8, -7, -5,-4},
      {-2, -5, -4, -3,-2},
      {-1, -3, -2, -1,-0},
      {0, 0, 0, 0, 0},
      {1, 0, 1, 2, 3},

      {2, 2, 3, 4, 5},
      {3, 4, 5, 7, 8},
      {4, 7, 8, 24, 1000}
   };
/*
   static INTERVALFEATURE ivalFs1[] ={
      {-4, -1000, -24, -8,-6},
      {-3, -9, -7, -5,-3},
      {-2, -6, -4, -3,-1},
      {-1, -4, -2, -1,-0},
      {0, 0, 0, 0, 0},

      {1, 0, 1, 2, 4},
      {2, 1, 3, 4, 6},
      {3, 3, 5, 7, 9},
      {4, 6, 8, 24, 1000}
   };
   static INTERVALFEATURE ivalFs2[] ={
      {-11,-9,-11, -12, -1000},
      {-10,-8,-10,-10,-12},
      {-9, -7,-9, -9, -11},
      {-8, -6,-8, -8, -10},
      {-7, -5, -7, -7,-9},
      {-6, -4, -6, -6,-8},
      {-5, -3, -5, -5,-7},
      {-4, -2, -4, -4,-6},
      {-3, -1, -3, -3,-5},
      {-2,  0, -2, -2,-4},
      {-1,  1, -1, -1,-3},
      { 0, -2, 0, 0,   2},
      { 1, -1, 1, 1,   3},
      { 2,  0, 2, 2,   4},
      { 3,  1, 3, 3,   5},
      { 4,  2, 4, 4,   6},
      { 5,  3, 5, 5,   7},
      { 6,  4, 6, 6,   8},
      { 7,  5, 7, 7,   9},
      { 8,  6, 8, 8,  10},
      { 9,  7, 9, 9,  11},
      {10, 8,10,10, 12},
      {11, 9,11, 12, 1000},
   };
   */
   data retVec;
   uint i = 0;
   uint len = sizeof(ivalFs)/ sizeof(ivalFs[0]);
   for(i = 0; i<len; ++i){
      double val = getFuzzyVal(interval, ivalFs[i]);
      retVec.push_back(val);
   }
   //out()<<"Intervall = "<<interval<<endl;
   //dumpData("Fuzzy Feature Vector", retVec);
   return retVec;


}

double Melody::fuzzyCorrelation(const data& intervals1,const data& intervals2)const
{
   double sum = 0, sum1 = 0, sum2 = 0;
   uint len1 = intervals1.size(), len2 = intervals2.size();
   if(len1 != len2 ) return -1.0;
   for(uint i = 0; i<len1; ++i){
      data tmp1 = getFuzzyFeature(intervals1[i]);
      data tmp2 = getFuzzyFeature(intervals2[i]);
      
      //dumpData<double>("Fuzzy feature:", tmp1, cout);
      sum  += getFuzzyFeatureCorr(tmp1, tmp2);
      sum1 += getFuzzyFeatureCorr(tmp1, tmp1);
      sum2 += getFuzzyFeatureCorr(tmp2, tmp2);      
   }
   if(sum1 == 0.0 || sum2 == 0.0){
       sum == 0.0 ? sum = 1.0 : sum = 0.0;
   }
   else sum /= sqrt(sum1*sum2);

   return sum;   
}
typedef enum eType_rhythmCategories{ eRC_VeryShort = 0, eRC_Short, eRC_Medium, eRC_Long, eRC_VeryLong}; 
durations Melody::fuzzyRhythm(const durations &dur, const uint refDur)const
{
   durations tmp;
   double m = refDur;
   for(uint i = 0; i<dur.size(); ++i){
      double dDur = (double) dur[i];
      //cout<<dur[i]<<" ";
      if(dDur> 3.3* m){
         tmp.push_back(eRC_VeryLong);
         continue;
      }
      if(dDur> 1.8 * m){
         tmp.push_back(eRC_Long);
         continue;
      }
      if(dDur > 0.9 * m){
         tmp.push_back(eRC_Medium);
         continue;
      }
      if(dDur> 0.45* m){
         tmp.push_back(eRC_Short);
         continue;
      }
      tmp.push_back(eRC_VeryShort);

   }
   //cout<<endl;
   //out()<<"Modus: "<<refDur<<endl;
   //dumpData("Fuzzy Rhythm", tmp, out());
   return tmp;      
}

double Melody::fuzzyRhythmCorr(Melody& mel, double &tmp, const uint mode, const uint submode)const
{

   EditDistance<durations> editDistance;

   double corr = 0, maxCorr = 0, maxCorr2 = 0;
   durations tmpDur1, tmpDur2;
   durations fuzzyDur1, fuzzyDur2;
   durations durTicks1, durTicks2;
   uint dummy = 0;
   vector<uint> modDur1, modDur2; 
   if(submode & TICKS){
      tmpDur1 = getTickDurations();
      tmpDur2 = mel.getTickDurations();
      modDur1.push_back(m_signature.m_nTicks);
      modDur2.push_back(mel.m_signature.m_nTicks);
      //dumpData<uint>("Durations 1", tmpDur1, out());
      //dumpData<uint>("Durations 2", tmpDur2, out());
   }
   else{
      tmpDur1 = getDurations();
      tmpDur2 = mel.getDurations();
      modDur1 = getDurationModus(dummy);
      modDur2 = mel.getDurationModus(dummy);
   }

   uint len1 = tmpDur1.size(), len2 = tmpDur2.size();
   uint maxLen = len1<len2? len2 : len1;
   for(uint i = 0; i<modDur1.size(); ++i){
      for(uint j = 0; j<modDur2.size(); ++j){
         double ed = 0;
         fuzzyDur1 = fuzzyRhythm(tmpDur1, modDur1[i]);
         fuzzyDur2 = fuzzyRhythm(tmpDur2, modDur2[j]);
         if(mode & NGRAM){
            uint nGramLen = 2<submode ? 2 : submode;
            corr = getNGramWeight<uint>(fuzzyDur1, fuzzyDur2, mode, nGramLen);
         } 
         else{
            ed = (double) editDistance(fuzzyDur1, fuzzyDur1.size(),fuzzyDur2, fuzzyDur2.size());
            corr = 1 - ed/maxLen;
         }
         if(corr>maxCorr) {
            //out()<<"Corr: "<<corr<<endl;
            maxCorr = corr;
            tmp = ed;
         }
      }

   }
   return maxCorr;
}
double  Melody::correlation(const data& vec1,const data& vec2, const uint MODE)const
{
   double mean1 = 0, mean2 = 0;
   double sum1 = 0.0, sum2 = 0.0, sum3 = 0.0;
   if(vec1.size()!= vec2.size()) return 0;

   if(MODE & INTERVALS){
      if(MODE & FUZZY)
         return fuzzyCorrelation(vec1, vec2);
      else
         if(MODE & FOURIER){
            data p1, p2;
            for(uint i = 0; i<vec1.size()-1; ++i){
               p1.push_back(vec1[i+1]-vec1[i]);
               p2.push_back(vec2[i+1]-vec2[i]);
            }
            return absNormNorm(p1,p2, MODE);
         }
         else
            return absNormNorm(vec1, vec2, MODE);
   }
   if(MODE & PEARSON_BRAVAIS){
      mean1 = getMean<double>(vec1);
      mean2 = getMean<double>(vec2);
   }
   if(MODE & CROSS){
      mean1 = vec1[0];
      mean2 = vec2[0];
   }
   for(uint i = 0; i<vec1.size() ; ++i){
      double x1 = (vec1[i] - mean1);
      double y1 = (vec2[i]  -mean2);
      sum1 += x1 * y1;
      sum2 += x1 * x1;
      sum3 += y1 * y1;
   }

   if(sum2 == 0.0 || sum3 == 0.0){
       sum1 == 0.0 ? sum1 = 1.0 : sum1 = 0.0;
   }

   else sum1 /= sqrt(sum2*sum3);
   if(MODE & SQUASH) if(sum1 < 0) sum1  = 0;
   return sum1;

}

double Melody::getMaxCorrelation(const data& shortMel, const data& longMel, uint leftOffset, double &tmp, const uint mode, 
                                 const uint  submode)const
{
   data mel1, mel2;
   uint i = 0;
   uint diff = 0;
   uint len1 = shortMel.size(), len2 = longMel.size();
   double corr = 0, max = -1.0;

   if(len1>len2) return -1.0;

   diff = len2-len1;

   tmp = 0;

   if(mode & NGRAM){
      uint nGramLen = submode<2 ? 2 : submode;
      corr = getNGramWeight<double>(shortMel, longMel, mode, nGramLen);
      return corr;
   }
   for(i = 1; i<=leftOffset; ++i){
      mel1 = getSubVector(shortMel, i, len1);
      mel2 = getSubVector(longMel,  0, len1-i);
      if(mode & RHYTHM){
         corr = gaussCrossCorrelation(mel1, mel2, 0.030, 0.0);
      }
      else if( mode & FOURIER){
         data ps1, ps2;
         realFT(mel1, -1, ps1);
         realFT(mel2, -1, ps2);
         corr = correlation(ps1, ps2, mode);
      }
      else{
         corr = correlation(mel1, mel2, mode);
      }
      if(corr>max){
         max = corr;
         tmp = - (double)i;
      }
      //out()<<"i = "<<(int)-i<<" corr = "<<corr<<" max = "<<max<<endl;
   }
   for(i = 0; i<diff+1; ++i){
      mel1 = shortMel;
      mel2 = getSubVector(longMel, i, i+len1);
      if(mode & RHYTHM){
         corr = gaussCrossCorrelation(mel1, mel2, 0.030, 0.0);
      }
      else if( mode & FOURIER){
         data ps1, ps2;
         realFT(mel1, -1, ps1);
         realFT(mel2, -1, ps2);
         corr = correlation(mel1, mel2, mode);
      }
      else{
         corr = correlation(mel1, mel2, mode);
      }
      if(corr>max){
         max = corr;
         tmp = i;
      }
      //out()<<"i = "<<i<<" corr = "<<corr<<" max = "<<max<<endl;
   }
   
   return max;
}

void Melody::getContourPoints(indizes& idx, PitchSequence& extrema,
                              const PitchSequence &pitches, const eContourMode contMode, bool const bUseEndpoints) const
{
   uint i = 0, j = 0;
   PitchSequence tmpM;
   indizes    tmpI;
   uint len = pitches.size();
   idx.clear();
   extrema.clear();
   if(len == 0 ) return;
   
   //out()<<"Entered getContourPoints()"<<endl;
   //dumpData("Melody", pitches, out());
   tmpM.push_back(pitches[0]);
   tmpI.push_back(0);

   // we have to take the last of equal tones, except for the first note
   for(i = 1; i<len; ++i){
      if(pitches[i]!= pitches[i-1]){
         for(j = 0; j<len-i; ++j){
            if(pitches[i+j] != pitches[i])
               break;
         }
         tmpM.push_back(pitches[i+j-1]);
         tmpI.push_back(i+j-1);
         //out()<<"Norm:Added Pitch "<<(uint)pitches[i+j-1]<<" with indx: "<<i+j-1<<endl;
      }
   }
   /*
   for(i = 1; i<len; ){

      while(i<len-1 && pitches[i+1] == pitches[i]){
         out()<<"Skipping "<<(uint)pitches[i]<<" at: "<<i<<endl;
         ++i;
      }
      tmpM.push_back(pitches[i]);
      tmpI.push_back(i);
      out()<<"Norm: Added Pitch "<<(uint)pitches[i]<<" with indx: "<<i<<endl;
      ++i;
   }
   */
   //first   
   if(bUseEndpoints){
      idx.push_back(0);
      extrema.push_back(pitches[0]);
      //out()<<"added ("<<0<<","<<(uint)pitches[0]<<")"<<endl;
   }
   switch(contMode){
      case eCont_classic:
      case eCont_all:
         for(i = 1; i< tmpM.size()-1; ++i){
            MIDIPitch p1 = tmpM[i-1], p2= tmpM[i], p3 = tmpM[i+1];
            //out()<<"Testing: "<<i<<" with "<<(uint)p1<<" "<<(uint)p2<<" "<<(uint)p3<<endl;
            if( (p2>p1 && p2 > p3) || ( p2 < p1 && p2 < p3)){
               if( contMode == eCont_all || p1 != p3){ //Wechselnote
                  extrema.push_back(p2);
                  idx.push_back(tmpI[i]);
                  //out()<<"added ("<<tmpI[i]<<","<<(uint)p2<<")"<<endl;
               }
            }
         }
      break;
      case eCont_steinbeck:
         if(tmpM.size()<=4){
            break;
         }
         
         for(i = 0; i< tmpM.size()-4; ++i){
            MIDIPitch d1 = sgn(tmpM[i+1]-tmpM[i]  );
            MIDIPitch d2 = sgn(tmpM[i+2]-tmpM[i+1]);
            MIDIPitch d3 = sgn(tmpM[i+3]-tmpM[i+2]);
            MIDIPitch d4 = sgn(tmpM[i+4]-tmpM[i+3]);
            //out()<<"Steinbeck:Testing: "<<tmpI[i]<<" with "<<(int)d1<<" "<<(int)d2<<" "<<(int)d3<<" "<<(int)d4<<endl;
            if(d1 == d2){
               if( d3 == d4 && d3 != d2){
                  extrema.push_back(tmpM[i+2]);
                  idx.push_back(tmpI[i+2]);
                  //out()<<"added ("<<tmpI[i+2]<<","<<(uint)tmpM[i+2]<<")"<<endl;
               }
            }
         }         
      break;
      default:

      break;
   }
   //last
   if(bUseEndpoints && len>1){
      idx.push_back(len-1);
      extrema.push_back(pitches[len-1]);
      //out()<<"added ("<<len-1<<","<<(uint)pitches[len-1]<<")"<<endl;
   }

}
void Melody::getContourPoints(indizes& idx, PitchSequence& extrema, 
                              const eContourMode contMode, const bool bUseEndpoints, const bool bRespectPhrases) const
{
   uint i = 0, j = 0;
   PitchSequence tmpM;   
   PitchSequence tmpE;   
   indizes    tmpI;
   idx.clear();
   extrema.clear();

   uint len = m_pitches.size();

   if(len == 0) return;
   if(!bRespectPhrases){
      getContourPoints(idx, extrema, m_pitches, contMode, bUseEndpoints);
      return;
   }   
   while(i<len){
      tmpM.clear();
      tmpI.clear();
      tmpE.clear();
      uint start = i;
      while(i<len-1 && !(m_phrases[i+1] & BEGIN_PHRASE)){     
         tmpM.push_back(m_pitches[i]);
         //out()<<"i = "<<i<<" pitches: "<<(uint)m_pitches[i]<<" phrase = "<<m_phrases[i]<<endl;
         i++;
      } 
      tmpM.push_back(m_pitches[i]);
      //out()<<"i = "<<i<<" pitches: "<<(uint)m_pitches[i]<<" phrase = "<<m_phrases[i]<<endl;
      ++i;
      getContourPoints(tmpI, tmpE, tmpM, contMode, bUseEndpoints);
      
      for(size_t k = 0; k<tmpI.size(); ++k){
         idx.push_back(tmpI[k] + start);
         extrema.push_back(tmpE[k]);
      }
   }
   

}
uint Melody::writeNormalization(uint contourMode)
{
   MIDIPitch low = 0, high = 0;
   if(m_onsets.size()==0) return 0;
   double duration = getTotalDuration();
   double t0 = m_onsets[0];
   char buf[1024] = "";

   sprintf(buf, ".norm.%x", contourMode);
   string fName =  stripExtension(getName().c_str()) +string(buf)  ;
   FILE *fp = 0;
   indizes idx;
   PitchSequence extrema;
   uint c_mode = contourMode & 0x0000000F;
   uint n_mode = contourMode & 0x000000F0;
   switch(c_mode){
      case CONTOUR_MODE_MF: getContourPoints(idx, extrema, eCont_classic); break;
      case CONTOUR_MODE_STEINBECK:getContourPoints(idx, extrema, eCont_steinbeck);break;
      case CONTOUR_MODE_ALL:getContourPoints(idx, extrema, eCont_all);break;
      case CONTOUR_MODE_NO_CONTOUR:
         for(uint i = 0; i<m_pitches.size();++i){
            idx.push_back(i);
         }
         extrema = m_pitches;
         break;
      default: return 0;         
   }
   if(idx.size() == 0) return 0;
   //out()<<"Found "<<idx.size()<<" extrema"<<endl;
   getPitchRange(low, high);
   double meanPitch = getMeanPitch();

   
   fp = fopen(fName.c_str(), "w");
   if(fp == 0) return 0;
   out()<<"Writing: "<<fName.c_str()<<endl;

   for(uint i = 0; i<idx.size();++i){
      double t = (m_onsets[idx[i]]-t0)/duration;
      double p = 0;
      switch(n_mode){
         case NORM_MODE_AMBITUS:
            p = (extrema[i] - low + .0);
            if(low != high){
               p = 2 * p/(high-low + .0)-1;
            }
            else p = 1;
            break;
         case NORM_MODE_MEAN:
            p = (extrema[i] - meanPitch + .0)/12;
      }
      fprintf(fp, "%.4f %.4f\n", t, p);
      //out()<<"idx["<<i<<"] = "<<idx[i]<<" t_orig = "<<m_onsets[idx[i]]<<" t_neu = "<<t;
      //out()<<" p_orig = "<<(int)extrema[i]<<" p_neu= "<<p<<endl;
      if(ferror(fp)) break;
   }
   fclose(fp);
   return idx.size();
}
/****************************************************
   returns code for Contour after D.Huron:
   0: horizontal
   1: horizontal-descending
   2: horizontal-ascending
   3: descending-horizontal
   4: descending
   5: concave
   6: ascending-horizontal
   7: convex
   8: ascending
   The value is coming from a 3-adic code, where horinzontal is mapped to 0, ascending to 2 and descending to 1.
   The combined value is 3*1st interval + 2nd interval
   Aggregated classes:
   0: horizontal (0)
   1: descending (4)
   2: ascending  (8)
   3: concave    (2, 3, 5)
   4: convex     (1, 6, 7)
****************************************************/

uint Melody::getHuronCode(const int i1, const int i2, bool bAggregated) const
{
   uint rel1 = i1>=0 ? (i1 == 0 ? 0 : 2 ) : 1; 
   uint rel2 = i2>=0 ? (i2 == 0 ? 0 : 2 ) : 1;
   uint code = 3* rel1 + rel2;
   if(bAggregated){
      switch(code){
         case 0:
            code = 0;
            break;
         case 1:
         case 6:
         case 7:
            code = 4;
            break;
         case 2:
         case 3:
         case 5:
            code = 3;
            break;
         case 4:
            code = 1;
            break;
         case 8:
            code = 2;
            break;

      }
   }
   return code;
}


uint  Melody::getHuronContour(bool bAggregated) const
{
   if(m_pitches.empty() || m_pitches.size()<3) return 9;

   uint first = m_pitches[0];
   uint last  = m_pitches[m_pitches.size()-1];
   //calc mid-average pitch
   double sum = 0;
   for(uint i = 1; i<m_pitches.size()-1;++i){
      sum +=m_pitches[i];
   }
   sum /=m_pitches.size()-2;

   return getHuronCode(int(floor(sum-first)), int(floor(last-sum)), bAggregated);
}

data Melody::realFT(data& vec, int isign, data &ps)const
{
   uint i = 0;
   uint len = vec.size();
   double *tmp = 0;
   double *tmp2 = 0;
      
   tmp = new double[len];
   if(!tmp) return vec;
   tmp2 = new double[1+(len>>1)];
   if(!tmp2) return vec;

   ps.clear();
   for(i = 0; i<len; ++i){
      tmp[i] = vec[i];
   }

   Fourier::slowRealFour(tmp, len, isign, tmp2);

   for(i = 0; i<len; ++i){
      vec[i] = tmp[i];
      if( i>0 && i<=len>>1){
         ps.push_back(tmp2[i]);
         //out()<<"i = "<<i<<" r = "<<tmp[i]<<" f = "<<tmp2[i]<<endl;
      }
   }
   
   delete[] tmp;
   delete[] tmp2;
   return vec;
}


data Melody::realPS(const data& vec, data& ps, uint start, uint end, double factor)
{
   uint len = vec.size()-1;

   if(end>len/2) end = len/2;
   out()<<"vec size = "<< vec.size()<<" ps size = "<<ps.size()<<endl;

   ps.clear();

   if(start == 0){
       ps.push_back(fabs(vec[0])*factor);
       start = 1;
   }
   
   for(uint i = start; i<=end; ++i){
      double val  = sqrt(vec[2*i]*vec[2*i] + vec[2*i+1]*vec[2*i+1])*factor;
      out()<<"idx: "<<i<<" val = "<<val<<endl;
      ps.push_back(val);

   }

   if(end >= len/2)
      ps.push_back(fabs(vec[1])*factor);

   return ps;
}

/** No descriptions */
data Melody::getRanks() const
{

   bool bFound = false;
   data pitchRanks;
   uint len = m_pitches.size();

   uint i = 0, j = 0;
   list<MIDIPitch> ranks;
   list<MIDIPitch>::iterator it;
   for(i = 0; i<len; ++i){
      bFound = false; 
      for(it = ranks.begin(); it!=ranks.end(); ++it){
         if(*it == m_pitches[i])
            bFound = true;
      }
       if(!bFound)ranks.push_back(m_pitches[i]);     
   }

   ranks.sort();
   /*
   for(it = ranks.begin(); it!=ranks.end(); ++it){
      out()<<(uint)*it<<" ";
   }
   out()<<endl;
   */
   for(i = 0; i<len; ++i){
      for(it = ranks.begin(), j = 0; it!=ranks.end(); ++it, ++j){
         if(m_pitches[i] == *it){
            pitchRanks.push_back(j);
            break;
         }
      }      
   }
   //dumpData("Rank melody", pitchRanks);
   return pitchRanks;
}
/** No descriptions */
data Melody::getIntervals(bool bAbs, uint mode) const
{
   data intervals;
   uint i = 0;

   for(i = 0; i<m_pitches.size()-1; ++i){
      double diff = m_pitches[i+1]-m_pitches[i];
      if(bAbs) diff = fabs(diff);
      if(mode & PARSONS){
         if(diff>0) diff = 1.0;
         else if(diff < 0) diff = -1.0;
      }
      else if(mode & FUZZY){
         uint idx = 0;
         getMax(getFuzzyFeature(diff), idx);
         diff = (double )idx-4;
      }
      intervals.push_back(diff);
   }
   //dumpData("Melody:", getPitchesAsDouble());
   //dumpData("Intervals:", intervals);
   return intervals;
}
data  Melody::getBarTimeAsPhases(const uint mode) const
{
   data ret;
   for(uint i = 0; i<m_timeLine.size(); ++i){
      double phase = m_timeLine[i].getAsPhase();
      double division = 0;
      switch (mode){
         case 0:
           ret.push_back(phase);
           break;
         case 1:
           ret.push_back(cos(phase));
           break;
         case 2:
           ret.push_back(sin(phase));
           break;
         case 3:
           division = m_timeLine[i].getDivisionAsPhase();
           ret.push_back(division);
           break;
         case 4:
           phase = m_timeLine[i].getAsPhase(0.75);
           ret.push_back(phase);
           break;
         default:
            break;
      }
   }
   return ret;
}

data  Melody::getTonalPitchesAsPhases(const uint mode) const
{
   data ret;
   for(uint i = 0; i<m_pitches.size(); ++i){
      Chroma c(m_pitches[i]-m_annotatedKey);
      double phase = 2 * PI * c.toDouble() /12.0;
      ret.push_back(phase);
   }
   return ret;
}

double normRand()
{
   return (rand()+0.0)/RAND_MAX;
}

MIDIPitch randPitch(MIDIPitch min, MIDIPitch max)
{
   return static_cast<MIDIPitch>((max-min)*normRand()) + min;
}
void Melody::setRandom(const uint length, const double density, const MIDIPitch min, const MIDIPitch max)
{
   m_pitches.clear();
   m_durations.clear();
   m_onsets.clear();
   m_pitches = getRandomMelody(length, density, min, max);
   if(m_pitches.size()){
      for(uint i = 0; i<m_pitches.size(); ++i){

         m_durations.push_back(1);
         m_onsets.push_back(.125*i);//120 bpm...
      }
   }
}

PitchSequence Melody::getRandomMelody(const uint length, const double density, const MIDIPitch min, const MIDIPitch max)
{
   PitchSequence retMel;
   MIDIPitch pitch = 0;
   double dens = density;
   
   if(min>max) return retMel;

   if(dens<0) dens = 0;
   if(dens>1.0) dens = 1.0;

   srand(time(0));
   
   pitch = randPitch(min,max);
   retMel.push_back(pitch);
   //out()<<" "<<(uint)pitch;
   for(uint i = 1; i<length; ++i){
      if(normRand()<dens){
         pitch = randPitch(min,max);
      }
      retMel.push_back(pitch);
      //out()<<" "<<(uint)pitch;
   }
   //out()<<endl;
   return retMel;
}

void Melody::printMelodyNames(uint start, uint end, ostream& out) const
{
   char buf[16] ="";
   uint nCount = 0, nWhiteSpace = 0;
   uint i = 0, j = 0;
   if(start> end){
      uint tmp = end;
      end = start;
      start = tmp;
   }
   if(end>m_pitches.size()) end = m_pitches.size();
   for( i = start; i<end; ++i){
      int pitchClass = getTonalityIdx((int)m_pitches[i], C0);
      int octave     = ((int)m_pitches[i]-C0)/12;

      
      sprintf(buf, "%s%d",genToneNames[pitchClass], octave);
      out<<buf;
      nWhiteSpace = 3 - strlen(genToneNames[pitchClass])-(octave<0 ? 1 : 0);
      for(j = 0; j<nWhiteSpace; ++j)
         out<<" ";
      nCount++;

      if(nCount == 16 ){
         out<<endl;
         nCount = 0;
      }
      for(j = 1; j<m_durations[i]; ++j){
         out<<".   ";
         nCount++;
         if(nCount == 16 ){
            out<<endl;
            nCount = 0;
         }
      }
   }
   out<<endl;
}

void Melody::normalizeTimeLine()
{
   m_timeLine.normalize();
}


double Melody::major[12] = {6.33, 2.23, 3.48, 2.33, 4.38, 4.09, 2.52, 5.19, 2.39, 3.66, 2.29, 2.88 };
double Melody::minor[12] = {6.33, 2.68, 3.52, 5.38, 2.60, 3.53, 2.54, 4.75, 3.98, 2.69, 3.34, 3.17 };

int Melody::getTonalityIdx(const int pitchClass, const int trans) const
{
   if(pitchClass>= trans) return (pitchClass-trans)%12;
   return (12+(pitchClass-trans)%12)%12;
}

PitchSequence Melody::getBar(const uint ndx) const
{
   return getBars(ndx, 1);
}

PitchSequence Melody::getBars(const uint iStart, const uint nCount) const
{
   PitchSequence mel;

   if(nCount<1) return mel;
   for( uint i = 0; i<m_timeLine.size(); ++i){
      if(m_timeLine[i].getBar() >= iStart && m_timeLine[i].getBar() < iStart+ nCount ){
         mel.push_back(m_pitches[i]);
         //printMelodyNames(i, i+1);
      }
      
      if(m_timeLine[i].getBar() == iStart+ nCount){
         break;
      }
   }
   return mel;
}

data Melody::getRawTonality(const uint barNdx, const uint mode, const bool dump) const
{
   uint i = 0, j = 0;
   uint pitchClass = 0;
   PitchSequence bar;
   data tonality;
   tonalityVector tonVec;
   if(barNdx >0){
      bar = getBar(barNdx);
      //out()<<"Bar "<<barNdx<<"> ";
   }
   else{
      bar = m_pitches;
   }
   uint ndx = barNdx-1;
   while(bar.empty() && ndx>1){
      bar = getBar(ndx);
      ndx--;
   }
   if(bar.empty()){
      out()<<"Empty Bar!"<<endl;
      return tonality;
   }
   for(i = 0; i<12; ++i){
      double sum = 0;
      double val = 0;
      //out()<<"Testing: "<<genToneNames[i]<<endl;
      
      for(j = 0; j<bar.size();++j){
         pitchClass = getTonalityIdx(bar[j], C0);
         if( mode & MAJOR){
            val = Melody::major[getTonalityIdx(pitchClass, i)];
         }
         else
            val = Melody::minor[getTonalityIdx(pitchClass, i)];
         //out()<<"Tone: "<<(uint)bar[j]<<" Pitchclass: "<<pitchClass<<" tonIdx = "<<getTonalityIdx(pitchClass, i)<<" val = "<<val<<endl;
         sum += val;
      }
      tonality.push_back(sum);
      tonVec.push_back(pair<uint, double>(i, sum));
      /*
      char buf[128] = "";
      sprintf(buf, "%s: %5.2f ", genToneNames[i], sum);
      cout<<buf;
      */
   }
   //cout<<endl;
   //printf("\n");
   if(dump){
      dumpData("Tonality vector:", tonality, out());
      tonVec = sortTonality(tonVec);
      dumpTonVec("Sorted:", tonVec, out());
   }
   return tonality;
}

data Melody::getTonicsKrumhansl() const
{
   data dummy;
   data ret;
   uint nBars = getBarCount();
   uint barIdx = (uint)-1;
   uint pitch = (uint)-1; 
   for(uint i = 0; i<m_timeLine.size(); ++i){
      uint tmpBarIdx = m_timeLine[i].getBar();
      if(tmpBarIdx != barIdx){
         pitch = getRawTonalityKrumhansl(tmpBarIdx, dummy, false); 
         if(pitch>12) pitch -=12;//ignore mode
         barIdx = tmpBarIdx; 
      }
      ret.push_back((double)(pitch + 60));
   }
   return ret;
}

uint Melody::getRawTonalityKrumhansl(const int barNdx, data& vec, const bool dump) const
{
   uint i = 0, j = 0;
   uint pitchClass = 0;
   PitchSequence bar;
   data vec1;
   double max = 0;
   uint maxIdx = 0;
   if(barNdx >= 0){
      bar = getBar(barNdx);
   }
   else{
      bar = m_pitches;
   }
   if(bar.empty()){
      if(dump){
         out()<<"Empty bar.";
      }
      return uint(-1);
   }
   vec.clear();
   if(dump) {
      char buf[128]="";
      if(barNdx>0){
         sprintf(buf, "Bar %d", barNdx); 
         dumpMelody(buf, bar);
      }
      else{
         dumpMelody("Melody", bar);
      }
   }
   for(i = 0; i<12; ++i){
      double sum1 = 0, sum2 = 0;
      double val1 = 0, val2 = 0.;
      if(dump)out()<<"Testing Root: "<<genToneNames[i]<<endl;
      for(j = 0; j<bar.size();++j){
         pitchClass = getTonalityIdx(bar[j], C0);
         val1 = Melody::major[getTonalityIdx(pitchClass, i)];
         val2 = Melody::minor[getTonalityIdx(pitchClass, i)];
         sum1 += val1;
         sum2 += val2;
         if(dump) {
            out()<<"Tone: "<<(uint)bar[j]<<" Pitchclass: "<<genToneNames[pitchClass]<<" tonIdx = "<<getTonalityIdx(pitchClass, i)<<endl;
            out()<<"---> Val(maj)= "<<val1<<" Val(min) ="<<val2<<endl;
         }
      }
      if(dump)out()<<"Sum "<<genToneNames[i]<<"(maj) = "<<sum1<<" sum "<<genToneNames[i]<<"(min) = "<<sum2<<endl;
      vec.push_back(sum1);
      vec1.push_back(sum2);
      if(sum1>max){
         max = sum1;
         maxIdx = i;
      }
      if(sum2>max){
         max = sum2;
         maxIdx = i+12;
      }
   }
   
   vec.insert(vec.end(), vec1.begin(), vec1.end());
   if(dump) {
      if(maxIdx>=12){
         if(barNdx>0){
            out()<<"Winner (Bar "<<barNdx<<") = "<<genToneNames[maxIdx-12]<<" minor"<<endl;
         }
         else{
            out()<<"Winner = "<<genToneNames[maxIdx-12]<<" minor"<<endl;         
         }
      }
      else{
         if(barNdx>0){
            out()<<"Winner (Bar "<<barNdx<<") = "<<genToneNames[maxIdx]<<" major"<<endl;
         }
         else{
            out()<<"Winner = "<<genToneNames[maxIdx]<<" minor"<<endl;         
         }
      }
   }
   return maxIdx;
}






void Melody::scaleTimeLine(double factor)
{
   uint i = 0;
   uint oldNum = m_signature.m_nCount;
   uint num = static_cast<uint>(floor(m_signature.m_nCount*factor));
   uint startBar = m_timeLine[0].getBar();
   if(num == 0) return;
   
   for(i = 0; i<m_timeLine.size(); ++i){
      uint totalBeat = (m_timeLine[i].getBar() - startBar)*oldNum + (m_timeLine[i].getBeat()-1);
      m_timeLine[i].dump();
      m_timeLine[i].setBar(startBar + totalBeat/num);
      m_timeLine[i].setBeat((totalBeat % num) +1);
      m_timeLine[i].setTicks(m_timeLine[i].getTicks());
      out()<<"->"<<endl;
      m_timeLine[i].dump();
   }
   m_signature.m_nCount = num;
   
}


MIDIPitch Melody::getPitch(const BarPosition& pos)const

{
   for(uint i = 0; i<m_timeLine.size(); ++i){
      if(m_timeLine[i] == pos) return m_pitches[i];
   }
   return 0;
}

uint Melody::sharps[6] = { 6,  1, 8, 3, 10, 5 };
uint Melody::flats[6]  = { 10, 3, 8, 1, 6,  11};

void Melody::countSigns(const PitchSequence& melody, vector<uint>& nFlats, vector<uint>& nSharps)const
{
   
   uint pitchClass = 0;
   uint i = 0, j = 0;
   for(j = 0; j<6; ++j){
      nFlats.push_back(0);
      nSharps.push_back(0);
   }
   for(i = 0; i<melody.size(); ++i){
      pitchClass = getTonalityIdx(melody[i], C0);
      //out()<<"Testing: "<<pitchClass<<endl;
      for(j = 0; j<6; ++j){
         if(sharps[j] == pitchClass){
            nSharps[j] += 1;
         }
         if(flats[j] == pitchClass){
            nFlats[j] += 1;            
         }
      }
   }
}

int Melody::getMajorBaseTone(const PitchSequence& melody)

{
   vector<uint> nFlats;
   vector<uint> nSharps;
   countSigns(melody, nFlats, nSharps);
   uint flatSum  = vectorSum<uint>(nFlats);
   uint sharpSum = vectorSum<uint>(nSharps);
   uint flatWSum  = vectorWeightedSum<uint>(nFlats);

   uint sharpWSum = vectorWeightedSum<uint>(nSharps);
   uint pitch = 0, ndx = 0;
   //dumpData<uint>("Flats:  ", nFlats, out());
   //dumpData<uint>("Sharps: ", nSharps, out());
   out()<<"Flats sum:  "<<flatSum<<endl;
   out()<<"Sharps sum: "<<sharpSum<<endl;
   out()<<"Flats Wsum:  "<<flatWSum<<endl;
   out()<<"Sharps Wsum: "<<sharpWSum<<endl;
   if(flatWSum == sharpWSum) {
      if(flatSum == 0)
         return 0;
      
   }
   if(flatWSum>=sharpWSum){
      uint ndx = getLast<uint>(nFlats);
      if(ndx == 0) return 0;
      uint pitch = Melody::flats[ndx-1];      
      return getTonalityIdx(pitch, 5);

   }
   ndx = getLast<uint>(nSharps);
   if(ndx == 0) return 0;
   pitch = Melody::sharps[ndx-1];
   return getTonalityIdx(pitch, +11);

}


double Melody::getHarmCorr(const Melody &mel1, const Melody &mel2, uint bar, uint mode)
{
   double corr = 0, sum1 = 0, sum2 = 0;
   if(bar<1 || bar>mel1.getBarCount() || bar>mel2.getBarCount())
      return 0;
   data ton1 = mel1.getRawTonality(bar, mode);
   data ton2 = mel2.getRawTonality(bar, mode);

   if(ton1.size()!=ton2.size()){
      return 0.0;
   }
   //dumpData("Tonal1", ton1);
   //dumpData("Tonal2", ton2);
   for(uint i = 0; i<ton1.size();++i){
      corr += ton1[i]*ton2[i];
      sum1 += ton1[i]*ton1[i];
      sum2 += ton2[i]*ton2[i];
   }

   //out()<<"Sum1 = "<<sum1<<" Sum2 = "<<sum2<<endl;  
   if(sum1 && sum2) corr /= sqrt(sum1*sum2);
   //out()<<"Corr: ="<<corr<<endl;
   return corr;
}
static string getTonalityName(uint ton)
{
   bool minor = ton>11;
   ton = ton % 12;
   string name = genToneNames[ton];
   if(minor) name += "m";
   return name;
}
static double circleOfFifth[] = { 5, 0, 7, 2, 9, 4, 11, 6, 1, 8, 3, 10, 2.5, 9.5, 4.5, 11.5, 6.5, 1.5, 8.5, 3.5, 10.5, 5.5, 0.5, 7.5};

double Melody::getHarmCorrKrumhansl(const Melody &mel, const uint mode, const uint submode)
{
   EditDistance<data> editDistance;
   data tonV1, tonV2;
   data sumV1, sumV2, sumV12;
   double corr = 0.;
   uint nBar1 = getBarCount(), nBar2 = mel.getBarCount();
   int start = 1, end = -1;
   if(mode & EDIT_DISTANCE){
      end = nBar1<nBar2 ? nBar2:nBar1;
   }
   else{
      end = nBar1<nBar2 ? nBar1:nBar2;
   }
   if(!end) return 0;
   if(submode==ALL_BARS){
      start = -1; end = -1;
   }
   for(int i = start; i<=end;++i){
      uint ton1 = getRawTonalityKrumhansl(i, tonV1);
      uint ton2 = mel.getRawTonalityKrumhansl(i, tonV2);

      //out()<<"ton1 = "<<ton1<<" ton2= "<<ton2<<endl;
      //out()<<"ton1 = "<<getTonalityName(ton1).c_str()<<" ton2= "<<getTonalityName(ton2).c_str()<<endl;
      //dumpData("ton1", tonV1, out());
      //dumpData("ton2", tonV2, out());
      if(mode & CROSS){
         sumV1.push_back(scalarProduct<double>(tonV1, tonV1));  
         sumV2.push_back(scalarProduct<double>(tonV2, tonV2));  
         sumV12.push_back(scalarProduct<double>(tonV1, tonV2));  
         //out()<<"Bar "<<i<<" Corr = "<<VecMult<double>(tonV1, tonV2)<<endl;

      }
      else if( mode & EDIT_DISTANCE){
         if(ton1>0)sumV1.push_back((double)ton1);
         if(ton2>0)sumV2.push_back((double)ton2);
      }
      else{
         double dTon1 = circleOfFifth[ton1];
         double dTon2 = circleOfFifth[ton2];
         //out()<<"ton1 = "<<dTon1<<" ton2= "<<dTon2<<endl;

         double omega = 2*PI/12.0;
         //out()<<"Corr = "<<cos(omega *(dTon1-dTon2))<<endl;
         corr += cos(omega *(dTon1-dTon2));
      }
   }
   if(mode & EDIT_DISTANCE){
      corr = (double)editDistance(sumV1, sumV1.size(),sumV2, sumV2.size());
      uint maxLen = getMaxLen<double>(sumV1, sumV2);
      if(maxLen) corr = 1- corr/(maxLen+.0);
      else corr = 0.0;
   } 
   else if(mode & CROSS){
      double sum1  = sqrt(vectorSum<double>(sumV1));
      double sum2  = sqrt(vectorSum<double>(sumV2));
      double sum12 = vectorSum<double>(sumV12);
      if(sum1 && sum2) corr = sum12/sum1/sum2;
      else corr = 0;
   }
   else{
   //out()<<"Corr: ="<<corr<<endl;
      if(submode!=ALL_BARS)
         corr /= nBar1;
   }
   //out()<<"Corr: ="<<corr<<endl;
   if(corr<0) corr = 0.0;
   return corr;
}

double Melody::harmonicCorrelation(const Melody& mel, const uint mode, const uint submode, double &tmp)
{
   double corr = 0;
   data majCorr, minCorr;
   double sum1 = 0, sum2 = 0;
   
   Melody mel1(*this), mel2(mel);

   if(mode & KRUMHANSL)
      return getHarmCorrKrumhansl(mel, mode, submode);
   
   mel1.normalizeTimeLine();
   mel2.normalizeTimeLine();
   uint barCount = mel1.getBarCount();
   if(mel2.getBarCount()<barCount) barCount = mel2.getBarCount();
   if(!barCount) return 0;
   
   for(uint i = 1; i<=barCount; ++i){
      sum1 +=getHarmCorr(mel1, mel2, i, MAJOR);
      sum2 +=getHarmCorr(mel1, mel2, i, MINOR);

   }
   corr = sum1/barCount;
   sum2 /= barCount;
   tmp = MAJOR;
   if(sum2>corr){
      corr = sum2;
      tmp = MINOR;
   }
   return corr;

}
double Melody::guessOnsetUnit() const
{
   if(m_onsets.size() == 0) return 1;
   double last = m_onsets[1];
   if(last>10)return 1000;
   return 1;
   
}
iData Melody::getPhraseBoundaries() const
{
   iData  phrases;
   uint  i = 0;
   uint dummy = 0;
   iData modi = getDurationModus(dummy);
   uint modus = getMinimum<uint>(modi);

   if(getSize() == 0){
      return phrases;
   }
   for(i = 0; i< getSize(); ++i){
      phrases.push_back(IN_PHRASE);
   }

   phrases[0] = BEGIN_PHRASE;
   phrases[getSize()-1] = END_PHRASE;
   if(phrases.size()<3) return phrases;
   double unit = guessOnsetUnit();
   //out()<<"unit = "<<unit<<endl;
   for(i = 1; i< m_durations.size()-2; ++i){
      double ioi = (m_onsets[i+1]-m_onsets[i])/unit;
      if(ioi>1.5 || m_durations[i]>=3.9*modus){
         phrases[i]   |= END_PHRASE;
         phrases[i+1] |= BEGIN_PHRASE;
      }
   }
   //dumpData("Durations", m_durations, out());
   //out()<<"Modus: "<<modus<<endl;
   //dumpData("Phrases  ", phrases, out());
   //out()<<"Lengths: "<<phrases.size()<<" "<<m_durations.size()<<" "<<m_pitches.size()<<endl;
   return phrases;   
}

Melody  Melody::getPhrase(const uint n) const
{
   Melody mel;
   if(m_phrases.empty()) return mel;

   uint nCount = 0;
   bool inPhrase = false;

   mel.setOStream(m_pOut);
   
   for(uint i = 0; i<m_phrases.size(); ++i){
      if(inPhrase && (m_phrases[i] & BEGIN_PHRASE))         
         break; 
      if(m_phrases[i] & BEGIN_PHRASE){
         if(nCount == n){
            inPhrase = true;
         }
         nCount++;
      }
      if(!inPhrase) continue;
      mel.m_pitches.push_back(m_pitches[i]);
      mel.m_durations.push_back(m_durations[i]);
      mel.m_onsets.push_back(m_onsets[i]);
      mel.m_timeLine.push_back(m_timeLine[i]);
   }
   
   if(mel.m_pitches.empty()) return mel;
   char suffix[1024] = "";
   sprintf(suffix, "_P%d", n);
   string name = getName() + string(suffix);
   mel.setName(name.c_str());
   mel.setSignature(m_signature);
   mel.setPhraseBoundaries();
   mel.normalize();
   //mel.printInfo();
   //mel.printMelody();   
   return mel;
}

MelodyList Melody::getPhrases() const
{
   MelodyList melList;
   for(uint i = 0; ;++i){
      Melody *pMel = new Melody();
      Melody tmp = getPhrase(i);
      if(tmp.empty()) break;
      *pMel = tmp;
      melList.push_back(pMel);   
   }
   return melList;
}

uint Melody::getPhraseCount() const
{
   uint nCount = 0;
   for(uint i = 0; i<m_phrases.size(); ++i){
      if(m_phrases[i] == BEGIN_PHRASE) nCount++;
   }
   return nCount;
}

uint Melody::findPhraseBegin(uint pos) const
{
   if(pos>=m_phrases.size()){
      return 0;
   }
   for(int i = pos; i>0; --i){
      if(m_phrases[i] == BEGIN_PHRASE) return i;
   }
   return 0;
}

uint Melody::findPhraseEnd(uint pos) const
{
   if(pos>=m_phrases.size()){
      return m_phrases.size()-1;
   }
   for(uint i = pos; i<m_phrases.size(); ++i){
      if(m_phrases[i] & END_PHRASE) return i;
   }

   return m_phrases.size()-1;
}

void Melody::setAccent(data & accents, const uint rules, const uint flag, const eAccentPos pos, const double acc)const
{
   if(rules & flag){
      if((size_t)pos <accents.size()){
         accents[pos] = acc;
      }
   }
}
dataVec Melody::getAccents(const uint rules, const double baseAcc, const double acc) const
{
   uint     i = 0, j = 0, k1 = 0, k2 = 0, lastLongPr = 0;
   uint     length = getSize();
   uint     ton = 0, pitch = 0;
   dataVec  ret;
   data     accents;
   data     intervals   = getIntervals(true,0);
   data     tonalities;
   uint     dummy = 0;
   indizes extremaSteinbeck, extremaMF;
   PitchSequence dummyMel;

   if(m_phrases.empty()) return ret;
   if(intervals.empty()) return ret;
   if(checkAccents(rules, ACCENTS_CONTOUR)){ 
      getContourPoints(extremaSteinbeck, dummyMel, eCont_steinbeck,false, true);
      getContourPoints(extremaMF, dummyMel, eCont_classic, false, true);
    }
   iData modi;
   uint modus = 0;
   if(checkAccents(rules, ACCENTS_DURATIONAL)){ 
      modi = getDurationModus(dummy);
      modus = getMinimum<uint>(modi);
   }
   //printInfo();
   //dumpData("Phrases", m_phrases, out());
   //dumpData("Intervals", intervals, out());
   //dumpData("Contour St", extremaSteinbeck, cerr);
   //dumpData("Contour MF", extremaMF, out());
   //dumpData("Durations", getDurations(), out());
   for(i = 0; i<length; ++i){
      //init accent vector
      accents.assign(eAccPos_MAX, baseAcc);      
      //const accents
      setAccent(accents, rules, ACCENT_CONST, eAccPos_CONST, acc);
      //pitch accents
      if(checkAccents(rules, ACCENTS_PITCH)){
         if(i<intervals.size()){
            if(intervals[i]>=3 ){
               setAccent(accents, rules, ACCENT_JUMPBEF3,  eAccPos_JUMPBEF3, acc);
               setAccent(accents, rules, ACCENT_JUMPBEA3,  eAccPos_JUMPBEA3, acc);
            }
            if(intervals[i]>=4){
               setAccent(accents, rules, ACCENT_JUMPBEF4,  eAccPos_JUMPBEF4, acc);
               setAccent(accents, rules, ACCENT_JUMPBEA4,  eAccPos_JUMPBEA4, acc);
            }
            if(intervals[i]>=5){
               setAccent(accents, rules, ACCENT_JUMPBEF5,  eAccPos_JUMPBEF5, acc);
               setAccent(accents, rules, ACCENT_JUMPBEA5,  eAccPos_JUMPBEA5, acc);
            }
         }
         if(i>0){
            if(intervals[i-1]>=3){
               setAccent(accents, rules, ACCENT_JUMPAFT3,  eAccPos_JUMPAFT3, acc);
               setAccent(accents, rules, ACCENT_JUMPBEA3,  eAccPos_JUMPBEA3, acc);
            }
            if(intervals[i-1]>=4){
               setAccent(accents, rules, ACCENT_JUMPAFT4,  eAccPos_JUMPAFT4, acc);
               setAccent(accents, rules, ACCENT_JUMPBEA4,  eAccPos_JUMPBEA4, acc);
            }
            if(intervals[i-1]>=5){
               setAccent(accents, rules, ACCENT_JUMPAFT5,  eAccPos_JUMPAFT5, acc);
               setAccent(accents, rules, ACCENT_JUMPBEA5,  eAccPos_JUMPBEA5, acc);
            }
            if(i>1 && i<intervals.size()){
               //out()<<"i = "<<i <<" "<<intervals[i-2]<<" "<<intervals[i-1]<<" "<<intervals[i]<<endl;
               if(abs(intervals[i-2]-intervals[i-1])>2 && abs(intervals[i-1]-intervals[i])>2){
                  if(!(m_phrases[i-1] & END_PHRASE) && !(m_phrases[i] & END_PHRASE)){
                     setAccent(accents, rules, ACCENT_JUMPLOC,  eAccPos_JUMPLOC, acc);
                  }
               }
            }
         }
      }
      //contour accents
      if(checkAccents(rules, ACCENTS_CONTOUR)){ 
         if(i>0 && i<length-1 && m_phrases[i] == IN_PHRASE){         
            if(m_pitches[i]> m_pitches[i-1] && m_pitches[i]>m_pitches[i+1]){
               setAccent(accents, rules, ACCENT_PEXTREM,  eAccPos_PEXTREM, acc);
            }
            if(m_pitches[i] < m_pitches[i-1] && m_pitches[i]<m_pitches[i+1]){
               setAccent(accents, rules, ACCENT_PEXTREM,  eAccPos_PEXTREM, acc);
            }
         }
         if(!extremaSteinbeck.empty() && k1>0 && m_phrases[i] == IN_PHRASE ){
            if(i == extremaSteinbeck[k1-1]+1){
               setAccent(accents, rules, ACCENT_PEXTRSTA,  eAccPos_PEXTRSTA, acc);
            }
         }
         if(!extremaSteinbeck.empty() && k1<extremaSteinbeck.size() && i == extremaSteinbeck[k1]){
            k1++;
            if(i>0 && i<length-1 && m_phrases[i] == IN_PHRASE){
               setAccent(accents, rules, ACCENT_PEXTRST,  eAccPos_PEXTRST, acc);
            }
         }
         if(!extremaMF.empty() && k2<extremaMF.size() && i == extremaMF[k2]){
            k2++;
            if(i>0 && i<length-1 && m_phrases[i] == IN_PHRASE){
               setAccent(accents, rules, ACCENT_PEXTRMF,  eAccPos_PEXTRMF, acc);
            }
         }
      }
      //Durational accents
      if(checkAccents(rules, ACCENTS_DURATIONAL)){ 
         if(m_durations[i]>modus){
            setAccent(accents, rules, ACCENT_LONGMOD,  eAccPos_LONGMOD, acc);
         }
         if(m_durations[i]>2*modus){
            setAccent(accents, rules, ACCENT_LONG2MOD,  eAccPos_LONG2MOD, acc);
         }

         if(i>0){
            uint dur1 = m_durations[i-1];
            uint dur2 = m_durations[i];

            if(dur1>=2*modus && lastLongPr != i-1){
               setAccent(accents, rules, ACCENT_ENDLOIOI,  eAccPos_ENDLIOI, acc);
            }
            if(dur2>dur1){
               setAccent(accents, rules, ACCENT_LONGPR,  eAccPos_LONGPR, acc);
               lastLongPr = i;
            }
            if(dur2>=2*dur1){
               setAccent(accents, rules, ACCENT_LONG2PR,  eAccPos_LONG2PR, acc);
            }

            if(i<length-2){
               if( !(m_phrases[i] & END_PHRASE) && !(m_phrases[i+1] & END_PHRASE) && !(m_phrases[i+2] & END_PHRASE)){
                  //out()<<"i = "<<i<<" dur1 = "<<m_durations[i]<<" dur2 ="<<m_durations[i+1]<<" "<<endl;
                  if(dur2<dur1 ){
                     setAccent(accents, rules, ACCENT_SHORTPR,  eAccPos_SHORTPR, acc);
                  }
                  if(2*dur2<=dur1){
                     setAccent(accents, rules, ACCENT_SHORT2PR,  eAccPos_SHORT2PR, acc);
                  }
               }
            }
         }
      }
      //meter accents      
      if(checkAccents(rules, ACCENTS_METER)){
         MetricTree2     mTree(m_timeLine[i].getSignature().m_nCount);
         eBeatType bt1 = mTree.getBeatType(m_timeLine[i].getBeat());
         BarPosition barPos = m_timeLine[i];

         if( m_timeLine[i].getTicks() == 0){
            if(bt1 == eBeat_primary){
               setAccent(accents, rules, ACCENT_BEAT13,  eAccPos_BEAT13, acc);
            }
            if(m_timeLine[i].getBeat() == 1){
               setAccent(accents, rules, ACCENT_BEAT1,  eAccPos_BEAT1, acc);
            }
            setAccent(accents, rules, ACCENT_BEATALL,  eAccPos_BEATALL, acc);
         }
         if(i<length-1){
            if(m_timeLine[i].getTicks() != 0 && m_timeLine[i+1] != m_timeLine[i].nextFullBeat()){
               if(m_timeLine[i+1].getBeat()>m_timeLine[i].getBeat()
                  || m_timeLine[i+1].getBar()>m_timeLine[i].getBar()
                  )
               {
                  setAccent(accents, rules, ACCENT_SYNCOP1,  eAccPos_SYNCOP1, acc);
                  eBeatType bt2 = mTree.getBeatType(m_timeLine[i].nextFullBeat().getBeat());
                  if(bt2 == eBeat_primary){
                     setAccent(accents, rules, ACCENT_SYNCOP2,  eAccPos_SYNCOP2, acc);
                  } 
               }
            }
         }
         if(i == length-1){
            if(m_timeLine[i].getTicks() != 0){
               setAccent(accents, rules, ACCENT_SYNCOP1,  eAccPos_SYNCOP1, acc);
               eBeatType bt2 = mTree.getBeatType(m_timeLine[i].nextFullBeat().getBeat());
               if(bt2 == eBeat_primary){
                  setAccent(accents, rules, ACCENT_SYNCOP2,  eAccPos_SYNCOP2, acc);
               } 
            }         
         }
      }
      //phrase accents
      if(checkAccents(rules, ACCENTS_PHRASE)){ 
         if(m_phrases[i] & BEGIN_PHRASE){
            setAccent(accents, rules, ACCENT_PHRASBEG,  eAccPos_PHRASBEG, acc);
         }
         if(m_phrases[i] & END_PHRASE){
            setAccent(accents, rules, ACCENT_PHRASEND,  eAccPos_PHRASEND, acc);
         }
         if(i>0){
            if((m_phrases[i] & END_PHRASE) && (m_phrases[i-1] & BEGIN_PHRASE)){
               setAccent(accents, rules, ACCENT_SHORTPHR, eAccPos_SHORTPHR,acc);               
            }
         }
      }
      //harmonic accents
      if(checkAccents(rules,ACCENT_HARM)){
         ton   = getRawTonalityKrumhansl(m_timeLine[i].getBar(), tonalities);
         pitch = getTonalityIdx(m_pitches[i], C0);
         uint baseTone = ton>12 ? ton-12: ton;

         int interval = pitch>=baseTone ? pitch - baseTone: (12+pitch)-baseTone;
         bool prime = interval == 0;
         bool third = ton>=12 ? interval  == 3 : interval  == 4;
         bool quint = interval == 7;
         
         //out()<<"Tonality: "<<ton<<" pitch = "<<pitch<<" interval = "<<interval<<" Prime = "<<prime<<" Quint = "<<quint<<" third = "<<third<<endl;
         
         if(prime|| third || quint){
            if(m_phrases[i] & END_PHRASE){
               setAccent(accents, rules, ACCENT_TRIAPHEN,  eAccPos_TRIAPHEN, acc);
            }
            if(rules & ACCENT_TRIAD){
               setAccent(accents, rules, ACCENT_TRIAD,  eAccPos_TRIAD, acc);
            }
         }
         
      }
      ret.push_back(accents);

      //out()<<"#"<<i<<endl;
      //if(i==0)dumpAccents("Accents", accents);
   }
   return ret;
}


data Melody::getSyncopes(const int level)const
{
   data syncopes;
   for(uint i = 0; i<m_timeLine.size(); ++i){
      syncopes.push_back(m_timeLine.isSyncope(i, level));
   }
   return syncopes;

}
data Melody::getOptiAccentsLinear(bool bThresholded) const
{
   int optiAcc = ACCENT_CONST|ACCENT_JUMPBEA4|ACCENT_PEXTRMF|ACCENT_SYNCOP1|ACCENT_PEXTREM|ACCENT_LONGPR|ACCENT_BEAT1|ACCENT_JUMPAFT4|ACCENT_JUMPLOC|ACCENT_PHRASEND;
   data weights;
   weights.assign(eAccPos_MAX, 0);
   weights[eAccPos_CONST] = -1.55;
   weights[eAccPos_JUMPBEA4] = -.52;
   weights[eAccPos_PEXTRMF] = .53; 
   weights[eAccPos_SYNCOP1] = .55;
   weights[eAccPos_PEXTREM] = 0.57;
   weights[eAccPos_LONGPR] = .68;
   weights[eAccPos_BEAT1] = .7;
   weights[eAccPos_JUMPAFT4] = .71;
   weights[eAccPos_JUMPLOC] = .94;
   weights[eAccPos_PHRASEND] = 1.44;
   dataVec accents = getAccents(optiAcc, 0, 1);
   data ret = getJointAccents(accents, weights);
   for(uint i = 0; i<ret.size(); i++){
      double val = exp(ret[i]);  
      if(val != -1){
         ret[i] = val/(1+val); 
      }
      if(bThresholded){
         double check = ret[i];
         ret[i] = (ret[i] >= 1/2.9 ? 1 : 0);
         check = ret[i];
      }
   }
   return ret;
}

data Melody::getOptiAccentsTree(bool bThresholded) const
{
   //FeatureDef syncopicity(STATISTIC|TRANS_SYNCOPES0|TRANS_SYNCOPES1|TRANS_SYNCOPES2|TRANS_SYNCOPES3, STAT_AVERAGE);
   int optiAcc       = ACCENT_PHRASEND|ACCENT_BEATALL|ACCENT_SYNCOP1|ACCENT_LONGPR|ACCENT_JUMPBEA5;
   //double sync       = getFeature(syncopicity)[0];
   data thom  = getThomassenAccents(ACCENT_THOMASSEN);

   dataVec accents   = getAccents(optiAcc, 0, 1);
   data ret;
   for(uint i = 0; i<m_timeLine.size(); ++i){
      double val = 0;
      if(accents[i][eAccPos_LONGPR]>0.5){
         if(thom[i]>0.1082){
            if(accents[i][eAccPos_BEATALL]>.5){
               val = .6773;  
            }
            else{
               if(accents[i][eAccPos_PHRASEND]>.5){
                  val = .6631;
               }
               else{
                  val = .4229;
               }
            }
         }
         else{
            val = .3673;
         }
      }
      else{
         if(accents[i][eAccPos_SYNCOP1]>.5){
            if(accents[i][eAccPos_JUMPBEA5]>.5){
               val = .5369;
            }
            else{
               val = .3211;
            }
         }
         else{
            val = .2114;
         }
      }
      if(bThresholded){
         val = val >= 1/2.9 ? 1 : 0;
      }
      ret.push_back(val);
   }

   return ret;
}

double Melody::getAccentSimilarity(const Melody& mel, const uint type, const uint rules, 
                                   const double baseAcc, const double acc)const
{
   double corr = 0.0;
   data  classBorder1, classBorder2, classBorder3;
   DiffCost<double>     diffCost;
   data accents1, accents2; 
   uint maxLen = 1;
   classBorder1.push_back(0);
   classBorder1.push_back(1);

   classBorder2.push_back(0);
   classBorder2.push_back(3);
   
   classBorder3.push_back(0);
   classBorder3.push_back(1);
   classBorder3.push_back(3);

   if(type & ACCENTS_CLASS_OPTI1){
      accents1= getOptiAccentsLinear();
      accents2 = mel.getOptiAccentsLinear();
      
   }
   else if(type & ACCENTS_CLASS_OPTI2){
      accents1= getOptiAccentsTree();
      accents2 = mel.getOptiAccentsTree();
      
   }
   else{      
      dataVec tmp1     = getAccents(rules, baseAcc, acc);
      dataVec tmp2     = mel.getAccents(rules, baseAcc, acc);
      accents1 = getJointAccents(tmp1);   
      accents2 = getJointAccents(tmp2);
      maxLen   = accents1.size()<accents2.size()? accents2.size() : accents1.size();
   }
   //dumpData("Joint1", accMel1);
   //dumpData("Joint2", accMel2);

   if(type & ACCENTS_CLASS1){
      accents1 = getClasses<double>(accents1, classBorder1);
      accents2 = getClasses<double>(accents2, classBorder1);
   }
   else if(type & ACCENTS_CLASS2){
      accents1 = getClasses<double>(accents1, classBorder2);
      accents2 = getClasses<double>(accents2, classBorder2);
   }
   else if(type & ACCENTS_CLASS3){
      accents1 = getClasses<double>(accents1, classBorder3);
      accents2 = getClasses<double>(accents2, classBorder3);
   }
   //dumpData("Joint1", accMel1);
   //dumpData("Joint2", accMel2);

   if(type & EDIT_DISTANCE){
      EditDistance<data>   ed;
      double tmp = (double) ed(accents1, accents1.size(), accents2, accents2.size());
      corr = tmp /maxLen;
   }
   
   if( type & ED_DIST_WEIGHT ){
      EditDistance2<double>   ed(&diffCost);
      double tmp = (double)   ed(accents1, accents2);
      double maxDist = ed.getMaxDist();
      //out()<<"ed = "<<tmp<<" maxDist = "<<maxDist<<endl;
      if(maxDist) corr = tmp/maxDist;
      else corr = 1;  
   }
   return 1-corr;
}

data Melody::getTransformation(const uint mode)const
{   
   data ret;

   if(mode & TRANS_PITCH){
      return getPitchesAsDouble();
   }
   if(mode & TRANS_DURATIONS16){
      return getDurationsAsDoubles();
   }
   data sync0, sync1, sync2, sync3;
   bool bSyncRet = false;

   if(mode & TRANS_SYNCOPES0){
      data sync0 = getSyncopes(1);
      if(ret.empty())ret = sync0;
      bSyncRet = true;
   }
   if(mode & TRANS_SYNCOPES1){
      data sync1 = getSyncopes(0);
      if(ret.empty())ret = sync1;
      else ret = fuzzyOr(ret, sync1);
      bSyncRet = true;
   }
   if(mode & TRANS_SYNCOPES2){
      data sync2 = getSyncopes(-1);
      if(ret.empty())ret = sync2;
      else ret = fuzzyOr(ret, sync2);
      bSyncRet = true;
   }
   if(mode & TRANS_SYNCOPES3){
      data sync3 = getSyncopes(-2);
      if(ret.empty())ret = sync3;
      else ret = fuzzyOr(ret, sync3);
      bSyncRet = true;
   }
   if(bSyncRet){
      //dumpData<double>("ret", ret, cout);
      return ret;
   }
   if(mode & TRANS_ACCENT_OPTI1){
      return getOptiAccentsLinear(false);
   }
   if(mode & TRANS_ACCENT_OPTI2){
      return getOptiAccentsTree(false);
   }
   if(mode & TRANS_ACC_OPTI1_THR){
      return getOptiAccentsLinear(true);
   }
   if(mode & TRANS_ACC_OPTI2_THR){
      return getOptiAccentsTree(true);
   }
   if(mode & TRANS_INTERVALS){
      return getIntervals(false, 0);
   }
   if(mode & TRANS_FUZZYINTERVALS){
      return getIntervals(false, FUZZY);
   }
   if(mode & TRANS_PARSONS){
      return getIntervals(false, PARSONS);
   }

   return ret;
}
vector< pair<double, uint> > Melody::getAsOnsetPitchPairs() const
{
   vector< pair<double, uint> > ret;
   if(m_pitches.empty()) return ret;
   uint l = m_pitches.size();
   for(uint i = 0; i<m_pitches.size(); ++i){
      double a  = m_onsets[i];
      uint b = m_pitches[i];
      pair<double, uint> tmp(m_onsets[i], m_pitches[i]);
      ret.push_back(tmp);
   }
   return ret;
}

data Melody::getFeature(const FeatureDef &feature) const
{
   data ret;
   if(feature.getMode() & TRANS_HURONCONTOUR){
      bool bAggregated = feature.getMode() & SIMPLE ? true: false;
      ret.push_back((double)getHuronContour(bAggregated));

   }
   else if(feature.getMode()& TRANS_BARS){
      if(feature.getMode()& PHASES){
         return getBarTimeAsPhases(0);
      }
      if(feature.getMode()& PHASES_DELTA){
         return getBarTimeAsPhases(4);
      }
      if(feature.getMode()& PHASES_RE){
         return getBarTimeAsPhases(1);
      }
      if(feature.getMode()& PHASES_IM){
         return getBarTimeAsPhases(2);
      }
      if(feature.getMode()& KEY){
         return getTonicsKrumhansl();
      }
      if(feature.getMode()& SUBBEAT_PHASES){
         return getBarTimeAsPhases(3);
      }
      if(feature.getMode()& TONAL_PITCH_PHASES){
         return getTonalPitchesAsPhases(feature.getMode());      
      }
      ret.push_back(getBarCount());
   }
   else if(feature.getMode() & STATISTIC){
      data trans = getTransformation(feature.getMode() );

      if(feature.getSubMode() == STAT_NONE){
         return trans;
      }
      else if(feature.getSubMode() & STAT_COUNT){
         ret.push_back(trans.size());
      }
      else if(feature.getSubMode() & STAT_RANGE){
         double max = 0;
         double min = 0;
         max = getMaximum<double>(trans);
         min = getMinimum<double>(trans);
         ret.push_back(max-min);      
      }
      else if(feature.getSubMode() & STAT_AVERAGE){
         double avg = getMean<double>(trans);
         ret.push_back(avg);
      }
      else if(feature.getSubMode() & STAT_MODUS){
         data modus;
         uint dummy = 0;
         modus= getModus<double>(trans, dummy);
         ret.push_back(modus[0]);
      }
      else if((feature.getSubMode() & STAT_MEDIAN) || (feature.getSubMode() & STAT_RUNTEST)){
         double med = getMedian<double>(trans);
         double val = med;
         if(feature.getSubMode() & STAT_RUNTEST){
            uint N = trans.size();
            if(N<=1){
               val  = 1;
            }
            else{
               uint n1 = 0, n2 = 0;
               uint nRuns= getNumberOfRuns<double>(trans, med, n1,n2);
               double n = double(n1 + n2);
               double mean = 2.0*n1*n2/n + 1.0;
               double var = 2.0 * n1 * n2 *(2* n1*n2 - n)/(n*n*(n-1));  
               double z  = (nRuns - mean)/sqrt(var);  
               val =  fabs(z)>1.96 ? 0 : 1 ;
               cout<<"n = "<<n<<" n1 = "<<n1<<" n2 = "<<n2<<" r = "<<nRuns<<" m = "<<mean<<" var = "<<var<<" sqrt(var)="<<sqrt(var)<<" z = "<<z<<endl;
            }
         }
         ret.push_back(val);
      }
      else if(feature.getSubMode() & STAT_VARIANCE){
         double var = getVariance<double>(trans);
         ret.push_back(var);
      }
      else if(feature.getSubMode() & STAT_MIN){
         double min = getMinimum<double>(trans);
         ret.push_back(min);
      }
      else if(feature.getSubMode() & STAT_MAX){
         double max = getMaximum<double>(trans);
         ret.push_back(max);
      }
      else if((feature.getSubMode() & STAT_ENTROPY) && (feature.getMode() & NGRAM) == 0){
         double entropy = getEntropy<double>(trans);
         ret.push_back(entropy);
      }
      else if(feature.getSubMode() & STAT_SKEWNESS){
         double skew = getSkewness<double>(trans);
         ret.push_back(skew);
      }
      else if(feature.getSubMode() & STAT_KURTOSIS){
         data test;
         test.push_back(1);test.push_back(2);test.push_back(3);
         getKurtosis<double>(test);
         double kurt = getKurtosis<double>(trans);
         ret.push_back(kurt);
      }
      else if(feature.getSubMode() & STAT_JSD || feature.getSubMode() & STAT_JSDMAX){
         double max = 0;
         uint idx = 0;
         data jsd = getJensenShannonDivergence(trans, max, idx);
         if(feature.getSubMode() & STAT_JSDMAX){
            ret.push_back(double(idx));
         }
         else{
            ret = jsd;
         }
      }
   }
   else if(feature.getMode() & NGRAM){
      uint maxLen = 0;
      double val = getNgramComplexity(feature.getMode(), feature.getSubMode(), maxLen); 
      //ret.push_back(getTransformedSize(feature.getMode(), feature.getSubMode()));
      //ret.push_back(maxLen);
      ret.push_back(val);
   }
   else if(feature.getMode() & SPECIAL){
      if(feature.getSubMode() & ACCENT_THOMASSEN){
         ret = getThomassenAccents(feature.getSubMode());
      }
   }
   else if(feature.getMode() & ACCENTS){
      ret = getAccentFeature(feature.getMode(), feature.getSubMode()); 
   }
   else if(feature.getMode() & ACCENTCOHERENCE){
      double val = getAccentCoherence(feature.getSubMode(), 0, 1, feature.getMode());
      ret.push_back(val);
   }
   
   return ret;
}

double Melody::getThomassenWeight( const uint v1, const uint v2, const uint v3,  const uint type, const uint mode) const
{
   double val = 0;
   static double values[9] = { 0, 0, 0, 1, .5, .71, 1, .83, .33}; 
   int i1 = v2-v1;
   int i2 = v3-v2;
   uint code  = getHuronCode(i1, i2, false);
   val = values[code]; 
   if(type == 2 && code != 0) val = 1-val;
   return val;
}
data Melody::getThomassenAccents(const uint mode) const
{
   PitchSequence    pitches = m_pitches;
   data     accents;
   pitches.insert(pitches.begin(), 2, m_pitches[0]);
   pitches.push_back(*pitches.rbegin());
   for(uint i = 2; i<pitches.size()-1; ++i){
      double p1 = getThomassenWeight(pitches[i-2], pitches[i-1], pitches[i], 2, mode); 
      double p2 = getThomassenWeight(pitches[i-1], pitches[i], pitches[i+1], 1, mode); 
      double w = p1 * p2;
      if(mode & THOMASSEN_THRESHOLD){
         w = w>=.5 ? 1 : 0; 
      }
      accents.push_back(w);
   }
   return accents;
}

data  Melody::getJensenShannonDivergence(const data& vec, double& max, uint &idx) const
{
   data ret;
   uint N = vec.size();
   JensenShannonDivergence<double> jsd(vec);
   
   for(uint i = 0; i<N-1; ++i){
      double a = jsd.Gamma(i+1);
      //cout<<"Gamma("<<i<<")="<<a<<endl; 
      double b = jsd.IncGamma(i+1,1);
      //cout<<"IncGamma("<<i<<")="<<b<<endl; 
      double v = jsd.getAt(i);
      if(v>max){
         max = v;
         idx = i;
      }
      ret.push_back(v);
   } 
   return ret;
}

data Melody::getAccentFeature(const uint type, const uint rules, const double baseAcc, const double acc) const
{
   data  classBorder1, classBorder2, classBorder3;

   classBorder1.push_back(0);
   classBorder1.push_back(1);

   classBorder2.push_back(0);
   classBorder2.push_back(3);
   
   classBorder3.push_back(0);
   classBorder3.push_back(1);
   classBorder3.push_back(3);

   
   dataVec accents     = getAccents(rules, baseAcc, acc);
   data    accMel      = getJointAccents(accents);   

   //dumpData("Joint1", accMel, cout);

   if(type & ACCENTS_CLASS1){
      accMel = getClasses<double>(accMel, classBorder1);
   }
   else if(type & ACCENTS_CLASS2){
      accMel = getClasses<double>(accMel, classBorder2);
   }
   else if(type & ACCENTS_CLASS3){
      accMel = getClasses<double>(accMel, classBorder3);
   }
   
   //dumpData("Joint1", accMel1);
   //temporary hack!!
   /*
   dataVec beat13   = getAccents(ACCENT_BEAT13, baseAcc, acc);
   data    accbeat13   = getJointAccents(beat13);   
   //double matchval = match<double>(accMel, accbeat13, true);
   double matchval = cosangle<double>(accMel, accbeat13);
   accMel.clear();
   accMel.push_back(matchval);
   */
   return accMel;
}

data Melody::getJointAccents(const dataVec& accents) const
{
   data ret;
   for(uint i = 0; i<accents.size(); ++i){
      double sum = vectorSum(accents[i]);
      ret.push_back(sum);
   }
   return ret;
}
data Melody::getJointAccents(const dataVec& accents, const data& weights) const
{
   data ret;
   for(uint i = 0; i<accents.size(); ++i){
      double sum = 0;
      if(weights.size() != accents[i].size()) return data();
      for(uint j = 0; j<weights.size(); ++j){
         sum += accents[i][j] * weights[j];
      }
      ret.push_back(sum);
   }
   return ret;
}

indizes Melody::getAccentTypePositions(const uint rules) const
{
   indizes ret;
   for(uint i = 0; i<32; ++i){
      if(GETBIT(rules,i)){
         ret.push_back(i);
      }
   }
   return ret;
   
}
data Melody::getAccentVector(const dataVec& accents, const uint accentTypePos)const
{
   data ret;
   for(uint i = 0; i<accents.size(); ++i){
      double   val = accents[i][accentTypePos];
      ret.push_back(val);
   }
   return ret;
}

double Melody::getMeanAccentCorrelation(dataVec accents, const uint rules, const uint mode)const
{
   indizes idz = getAccentTypePositions(rules);
   double sum = 0.0;
   uint N = idz.size();
   if(N <= 1) return 1;
   uint count = 0;
   for(uint i = 0; i<N; ++i){
      data vec1 = getAccentVector(accents, idz[i]);
      //dumpData("Accents1", vec1);
      sum +=1;//the diagonal...+
      for(uint j = i+1; j<N; ++j){
         data vec2 = getAccentVector(accents, idz[j]);
         //dumpData("Accents2", vec2);
         //out()<<"Corr "<<i<<", "<<j<<" = "<<cosangle<double>(vec1, vec2)<<endl;
         if((mode & 0x0000f00f) == ACCENT_COHERENCE_MEAN_KORR){ 
            sum += cosangle<double>(vec1, vec2);
         }
         else{
            EditDistance<data> ed;
            double tmp = (double)ed(vec1, vec1.size(),vec2, vec2.size());
            if((mode & 0x0000f00f) == ACCENT_COHERENCE_MEAN_ED_KORR)
               tmp = 1-tmp/vec1.size();
             else
               tmp = exp(-tmp/vec1.size());
             sum += tmp;
         }
      }
   }
   return 2*sum/(N*(N+1));
}

double Melody::getAccentCoherence(const uint rules, const double baseAcc, const double acc,  const uint mode)const
{
   double ret = 0.0, mean = .0, var = 0.0, var2 = 0.0, max = 0.0, sumOfSquares = 0;
   uint nAccents = 0;
   dataVec accents  = getAccents(rules, baseAcc, acc);
   data jointAccents;
   switch(mode & 0x0000f00f){
      case   ACCENT_COHERENCE_MAX_VARIANCE:
         jointAccents = getJointAccents(accents);
         //dumpData<double>("Joint1: ", jointAccents, out());
         if(jointAccents.empty()) return -1.0;
         max = getMaximum<double>(jointAccents);
         nAccents = jointAccents.size() - getElementCount<double>(jointAccents, baseAcc, -1 );
         if(max == 0 || nAccents == 0) return -1.0; // should not happen
         sumOfSquares = scalarProduct(jointAccents, jointAccents);
         ret = sumOfSquares / ( nAccents * max * max);
         //mean  =  getMean<double>(jointAccents);
         //var   =  getVariance<double>(jointAccents);
         //var2  =  getVariance2<double>(jointAccents);
         //if(mean!=0) ret = sqrt(var)/mean;
         //out()<<"AccCoh: mean:"<<mean<<" var = "<<var<<" var2 = "<<var2<<" sdv = "<<sqrt(var)<<endl;
      break;
      case   ACCENT_COHERENCE_MEAN_KORR:
         ret = getMeanAccentCorrelation(accents, rules, mode);
      break;
      case   ACCENT_COHERENCE_MEAN_ED_KORR:
      case   ACCENT_COHERENCE_MEAN_ED_KORR_EXP:
      
         ret = getMeanAccentCorrelation(accents, rules, mode);
      break;
      default:
         out()<<"Undefined accent coherence mode!"<<endl;
      break;

   }
   return ret;
}

uint Melody::getTransformedSize(const uint mode, const uint submode)const
{
   if(mode & WEIGHTED){
      Melody mel = clone();
      mel = mel.normalize();
      return getSize();
   }

   if(mode & INTERVALS){
      return getSize()-1;
   }
   if(mode & TRANS_BARS){
      return getBarCount();
   }

   return getSize(); 
}
double getEntropyFromFrequencies(const vector<uint> &freq, bool bNorm = true)
{
   double sum = 0;
   uint len = vectorSum<uint>(freq);
   if(len == 0) return 0;

   for(uint i = 0; i< freq.size(); ++i){
      //cerr<<"Val: "<<dist[i].first<<" Freq="<<dist[i].second<<endl;
      double p = (freq[i] +.0)/ len;
      double entropy = p * log((double)p) / log(2.0);
      sum += entropy;
   }
   double norm = bNorm ? log(freq.size()+.0)/log(2.0) : 1;
   return norm? -sum/norm : -sum;
}

double Melody::getNgramComplexity(const uint mode, const uint maxNgramLen, uint &maxLen) const
{
   Melody mel = clone();
   data dMel;
   if(mode & CONTOUR){
      dMel = mel.contourize(mode);
   }
   else if(mode & INTERVALS){
      dMel = mel.getIntervals(false, mode);
   }
   else if(mode & DURATION16){
      dMel = mel.getDurationsAsDoubles();
   }
   else{
      dMel = mel.getPitchesAsDouble();
   }
   double sum = 0;
   uint N = dMel.size();
   vector<uint> freq, indices;
   maxLen = (maxNgramLen==0 || maxNgramLen>N)  ? N : maxNgramLen; 
   if(mode & ZIPF && maxNgramLen != 0){
      uint count = getAllNGrams<double>(dMel, maxLen, indices, freq);
      double val = getZipfComplexity(freq);
      return val;
   }
   double norm = 0;
   for(uint i = 1; i<=maxLen; ++i){
      uint count = getAllNGrams<double>(dMel, i, indices, freq);
      double val = 0;
      if(count == 1){
         maxLen = i-1;
         break;
      }
      if(mode & STAT_ENTROPY){
         val = getEntropyFromFrequencies(freq);
      }
      else if(mode & ZIPF){
         val = getZipfComplexity(freq);         
      }
      else{
         val = (count-1.0)/(N-i);
      }
      double w = 1;
      if(mode & WEIGHTED) w = 1/i;
      norm +=w;
      //sum += w*(count-1.0)/(N-i);
      sum += w*val;
      if(count == N-i+1 && maxNgramLen==0){ 
         maxLen = i;
         break;//maxim reached.
      }
   }
   //sum = maxLen ? sum / maxLen : 0;
   sum = norm ? sum / norm : 0;
   return sum;
}


double  Melody::getZipfComplexity(const iData& vec)const
{
   size_t n = vec.size();
   if(n <= 1) return 0;

   iData sorted = sort<uint>(vec, 1);
   //dumpData<uint>("Original", vec, out());
   //dumpData<uint>("Sorted", sorted, out());

     
   double sum_xx = 0;
   double sum_xy = 0; 
   double sum_x = 0; 
   double sum_y = 0; 

   for(size_t i = 0; i<sorted.size(); ++i){
      double x_t = log(i + 1.0);
      double y_t = log(sorted[i]+.0); 
      sum_xy += x_t*y_t;
      sum_xx += x_t*x_t;
      sum_x  += x_t;
      sum_y  += y_t;
   }
   if(n*sum_xx==sum_x*sum_x ) return 0;
   double gamma = (n*sum_xy  - sum_x*sum_y)/(n*sum_xx-sum_x*sum_x);
   double zipfC  = pow(2.0,gamma);
   return zipfC;
}

bool tonEntryComp(const tonalityEntry& e1,const tonalityEntry& e2)
{
   if(e1.second<e2.second) return false;
   return true;
} 
tonalityVector sortTonality(const tonalityVector& vec)
{
   tonalityVector ret = vec;
   sort(ret.begin(), ret.end(), tonEntryComp);
   return ret;
}

void  dumpTonVec(const char* header, const tonalityVector& vec, ostream& out)
{
   char buf[128] = "";
   out<<header<<" ";
   for(uint i = 0; i<vec.size(); ++i){
      sprintf(buf, "%-2s: %5.2f|", genToneNames[vec[i].first], vec[i].second);
      out<<buf;
   }
   out<<endl;

}