//=============================================================================
//  MuseScore
//  Music Composition & Notation
//
//  Copyright (C) 2007-2011 Werner Schweer
//
//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License version 2
//  as published by the Free Software Foundation and appearing in
//  the file LICENCE.GPL
//=============================================================================


//=============================================================================
//  Keyfinder.cpp uses code from the "Melisma Music Analyzer"
//  project:
//
//  Copyright (C) 2000 Daniel Sleator and David Temperley
//  See http://www.link.cs.cmu.edu/music-analysis
//  for information about commercial use of this system
//=============================================================================

/* _Important comment about how pitches and keys are represented_

  Notes may be inputted in either Note format (with just a pitch) or
  TPCNote format (with a TPC). The user may also specify
  "npc_or_tpc_profile": 0 for npc, 1 for tpc. (A TPC profile is like
  an NPC profile except that a user-settable default value is used for
  TPC's outside the range of b2 to #4.) We assume that Note format
  will be used with an npc profile; if not, a fatal error is
  reported. TPCNote format could be used with either; if you use
  TPCNote input with an npc profile, the notes are simply mapped on to
  one cycle of the LOF.  (Scoring mode 0 or 3 requires an npc profile;
  otherwise a fatal error is reported.)

  Pitches (e.g. note[].tpc values) are always represented in
  line-of-fifths order, with C = 14. (This is done when the pitches
  are first read in). (If there's Note input, or TPCNote input and
  npc_or_tpc_profile = 0, all pitches are shifted into one cycle of
  the line of fifths, from 9 to 20 inclusive; again, this is done when
  the input is read.)

  Keys are also represented in line-of-fifths order. Major keys are 0
  to 27, C=14; minor keys are 28-55, C minor = 42. If
  npc_or_tpc_profile = 0, the search is nominally done on all keys,
  but only keys from 9 to 20 (inclusive) and 37 to 48 (inclusive) are
  looked at; others are given large negative values.

  Although key-profiles are read in in pitch-height order, they are
  then adjusted (in generate_tpc/npc_profiles) to line-of-fifths
  order, with the tonic as 5. */

#include <string.h>
#include <math.h>
#include "keyfinder.h"
#include "sig.h"
#include "pitchspelling.h"
#include "synthesizer/event.h"

namespace Ms {

#if 0
//---------------------------------------------------------
//   SBeat
//---------------------------------------------------------

struct SBeat {
      int time;
      };

//---------------------------------------------------------
//   MidiSegment
//---------------------------------------------------------

struct MidiSegment {
      int start;
      int end;
      QList<Event> snote;
      int numnotes;           // number of notes in the segment
      qreal average_dur;     // average input vector value (needed for K-S algorithm)
      };

#define CHANGE_PENALTY        12

static float change_penalty   = CHANGE_PENALTY;
static int npc_or_tpc_profile = 1;
static int scoring_mode       = 1;
static int verbosity          = 1;

static qreal major_profile[12] = {5.0, 2.0, 3.5, 2.0, 4.5, 4.0, 2.0, 4.5, 2.0, 3.5, 1.5, 4.0};
static qreal minor_profile[12] = {5.0, 2.0, 3.5, 4.5, 2.0, 4.0, 2.0, 4.5, 3.5, 2.0, 1.5, 4.0};
static qreal default_profile_value=1.5;

/*

CBMS profiles:

major_profile = 5.0 2.0 3.5 2.0 4.5 4.0 2.0 4.5 2.0 3.5 1.5 4.0
minor_profile = 5.0 2.0 3.5 4.5 2.0 4.0 2.0 4.5 3.5 2.0 1.5 4.0

Bayesian profiles (based on frequencies in Kostka-Payne corpus):

major_profile =  0.748  0.060  0.488  0.082  0.670  0.460  0.096  0.715  0.104  0.366  0.057  0.400
minor_profile =  0.712  0.084  0.474  0.618  0.049  0.460  0.105  0.747  0.404  0.067  0.133  0.330

Krumhansl's profiles:

major_profile = 6.35 2.23 3.48 2.33 4.38 4.09 2.52 5.19 2.39 3.66 2.29 2.88
minor_profile = 6.33 2.68 3.52 5.38 2.60 3.53 2.54 4.75 3.98 2.69 3.34 3.17

Krumhansl's minor, normalized: 5.94 2.51 3.30 5.05 2.44 3.31 2.38 4.46 3.73 2.52 3.13 2.97

*/

static int firstbeat;
static QList<Event> note;
static QList<MidiSegment> segment;  // An array storing the notes in each segment
static int segtotal;                // total number of segments - 1

static qreal seglength;

static QList<int> seg_prof[28];
static QList<qreal> keyscore()[56];
static QList<SBeat> sbeat;
static QList<qreal> analysis[56][56];

static QList<int> best[56];
static QList<int> final;

// static int numnotes, numchords, num_sbeats;
static int num_sbeats;

static QList<int> pc_tally;
static QList<qreal> finalscore();

static qreal key_profile[56][28];
static int final_timepoint;

//---------------------------------------------------------
//   print_keyname
//---------------------------------------------------------

static void print_keyname(int f)
      {
      static const char letter[] = "CGDAEBF";

      int mf = f % 27;
      mf -= 14;

      qDebug("(%d,%d)%c", f, mf, letter[f % 7]);
      if (f < 6 || (f >= 28 && f < 34))
            qDebug("-");
      if ((f >= 6 && f < 13) || (f >= 34 && f < 41))
            qDebug("b");
      if ((f >= 20 && f < 27) || (f >=48 && f < 55))
            qDebug("#");
      if (f == 27 || f == 55)
            qDebug("x");
      if (f >= 28)
            qDebug("m");
      if (f < 28)
            qDebug(" ");
      qDebug(" ");
      }

//---------------------------------------------------------
//   create_segments
//    Each segment starts at a sbeat and ends at the
//    following sbeat
//---------------------------------------------------------

static void create_segments()
      {
      segment.append(MidiSegment());

      segment[0].start = firstbeat; // Always start a segment at the very beginning of the piece (the first beat)
      for (int b = 0; b < num_sbeats; b++) {
            if (b == 0 && (sbeat[0].time - firstbeat) < ((sbeat[1].time-firstbeat) - (sbeat[0].time-firstbeat))/2)
                  continue;
            /* If it's the first beat of the piece, and the upbeat is
               less than half of the first beat interval, don't start a segment
            */
            else {
                  MidiSegment seg;
                  seg.start          = sbeat[b].time;
                  segment.back().end = sbeat[b].time;
                  segment.append(seg);
                  }
            }
      int s = segment.size() - 1;

    /* If final segment starts at or after final timepoint of piece, ignore it,
       decrementing number of segments by 1; if not, set that segment's ending
       to final timepoint of piece
       */
      if (segment[s].start >= final_timepoint) {
            s--;
            }
      else {
            segment[s].end = final_timepoint;
            /* qDebug("Final segment ends at %d", segment[s].end); */
            }
      segtotal = s; // index of final segment
      }

//---------------------------------------------------------
//   fill_segments
//---------------------------------------------------------

static void fill_segments()
      {
      for (int s = 0; s < segment.size(); ++s) {
            foreach (const Event& n, note) {
                  int ontime  = n.ontime();
                  int offtime = n.offtime();
                  int start   = segment[s].start;
                  int end     = segment[s].end;

                  if (ontime >= start && ontime < end && offtime <= end) {
                        // note begins and ends in segment
                        Event sn(ME_NOTE);
                        sn.setDataA(n.dataA());
                        sn.setTpc(n.tpc());
                        sn.setDuration(n.duration());
                        segment[s].snote.append(sn);
                        }
                  if (ontime >= start && ontime < end && offtime > end) {
                        // note begins, doesn't end in segment
                        Event sn(ME_NOTE);
                        sn.setDataA(n.dataA());
                        sn.setTpc(n.tpc());
                        sn.setDuration(end - ontime);
                        segment[s].snote.append(sn);
                        }
                  if (ontime < start && offtime > start && offtime <= end) {
                        // note ends, doesn't begin in segment
                        Event sn(ME_NOTE);
                        sn.setDataA(n.dataA());
                        sn.setTpc(n.tpc());
                        sn.setDuration(offtime - start);
                        segment[s].snote.append(sn);
                        }
                  if (ontime < start && offtime > end) {
                        // note doesn't begin or end in segment
                        Event sn(ME_NOTE);
                        sn.setDataA(n.dataA());
                        sn.setTpc(n.tpc());
                        sn.setDuration(end - start);
                        segment[s].snote.append(sn);
                        }
                  }
            segment[s].numnotes = segment[s].snote.size();
            // qDebug("fillSegments %d: %d-%d  %d", s, segment[s].start,
            //   segment[s].end, segment[s].numnotes);
            }
      }

//---------------------------------------------------------
//   count_segment_notes
//    In each segment, tally up the notes of each TPC
//---------------------------------------------------------

static void count_segment_notes()
      {
      for (int s = 0; s <= segtotal; ++s) {
            pc_tally[s] = 0;

            for (int y = 0; y < 28; ++y)  // cycle through the pc's, make sure all the seg_prof values are zero
                  seg_prof[y].append(0);

            qreal total_dur = 0;

            for (int n = 0; n < segment[s].numnotes; ++n) {
                  if (scoring_mode == 0)
                        total_dur += segment[s].snote[n].duration();
                  for (int y=0; y<28; ++y) {
                        if (segment[s].snote[n].tpc() == y) {
                              if(seg_prof[y][s]==0)
                                    pc_tally[s]++;
                              /* This keeps track of how many different pc's the segment contains. This counts TPCs, not NPCs! */
                              /* If scoring_mode is > 1, set array value to 1. If 0, add the note's duration to the
                                 array value (as in the K-S algorithm) */
                              if (scoring_mode > 0)
                                    seg_prof[y][s] = 1;
                              else {
                                    seg_prof[y][s] += segment[s].snote[n].duration();
                                    }
                              }
                        }
                  }

            if(scoring_mode == 0) {
                  if(pc_tally[s]==0)
                        segment[s].average_dur = 0.0;
                  segment[s].average_dur = total_dur / 12.0;
                  /* qDebug("Segment %d total dur = %6.3f, average dur = %6.3f", s, total_dur, segment[s].average_dur); */
                  }

            if (verbosity>=2) {
                  qDebug("Segment %d: ", s);
                  for (int y=0; y<28; ++y) {
                        if(npc_or_tpc_profile == 0 && (y<9 || y>20))
                              continue;
                        qDebug("%d ", seg_prof[y][s]);
                        }
                  }
            /* qDebug("pc_tally = %d", pc_tally[s]); */
            }
      }

//---------------------------------------------------------
//   prepare_profiles
//    We're only here if scoring_mode is 0 (the K-S algorithm).
//    Sum all the profile values, take the mean,
//    and subtract that from each value
//---------------------------------------------------------

static void prepare_profiles()
      {
      qreal total = 0.0;
      for (int i = 0; i < 12; i++) {
            total += major_profile[i];
            }
      qreal average = total / 12.0;
      for (int i = 0; i < 12; i++)
            major_profile[i]=major_profile[i] - average;

      total = 0;
      for (int i = 0; i < 12; i++) {
            total += minor_profile[i];
            }
      average = total / 12.0;
      for (int i = 0; i < 12; i++)
            minor_profile[i] = minor_profile[i] - average;

      if (verbosity > 2) {
            qDebug("Adjusted major profile: ");
            for(int i = 0; i < 12; i++)
                  qDebug("%6.3f ", major_profile[i]);
            qDebug("Adjusted minor profile: ");
            for (int i = 0; i < 12; i++)
                  qDebug("%6.3f ", minor_profile[i]);
            }
      }

/* Here we generate the key profiles. (This is assuming tpc input.) Key_profile[key] numbers correspond to
   the line of fifths, with C = 14.  Major keys are 0-27, minor keys are 28-55. PCs are also numbered
   according to the line of fifths. The major_step_profile has the tonic in step 5. For a given key, the
   profile value for a given tpc is equal to the line of fifths difference between the tpc and the key, plus
   5. */

//---------------------------------------------------------
//   generate_tpc_profiles
//---------------------------------------------------------

static void generate_tpc_profiles()
      {
      int key, shift, tpc, i;
      float majp[12];
      float minp[12];

      /* First we rearrange the key profile values (inputted in pitch height order) into lof order, C = 5 */
      for(i=0; i<12; i++) {
            majp[((((i * 7) % 12) + 5) % 12)] = major_profile[i];
            minp[((((i * 7) % 12) + 5) % 12)] = minor_profile[i];
            }

      for (key=0, shift=0; key<28; ++key, ++shift) {
            for (tpc=0; tpc<28; ++tpc) {
                  if (tpc-shift >= -5 && tpc-shift <= 6) {
                        key_profile[key][tpc] = majp[5 + (tpc-shift)];    /* For example: for key 14 (C major) and tpc 17 (A),
                                                                             use profile step 5 + (17-14) = 8 */
                        }
                  if (tpc-shift < -5 || tpc-shift > 6) {
                        key_profile[key][tpc] = default_profile_value;
                        }
                  }
            }
      for (key=28, shift=0; key<56; ++key, ++shift) {
            for (tpc=0; tpc<28; ++tpc) {
                  if (tpc-shift >= -5 && tpc-shift <= 6 ) {
                        key_profile[key][tpc] = minp[5 + (tpc-shift)];
                        }
                  if (tpc-shift < -5 || tpc-shift > 6) {
                        key_profile[key][tpc] = default_profile_value;
                        }
                  }
            }
/*
      This routine just prints out the key profiles
      for(key=0; key<56; ++key) {
            for(tpc=0; tpc<28; ++tpc) {
                  qDebug("%1.2f ", key_profile[key][tpc]);
                  }
            }
*/
      }

//-----------------------------------------------------------------------------
//   generate_npc_profiles
//    This is an alternative function for generating profiles given NPC
//    profile. It's similar to the TPC version, except that only keys on
//    a certain range of the line are assigned non-zero values. (This is
//    really an unnecessary step, as steps outside the range will be
//    disqualified in "match_profiles" in any case.) Also, for those keys
//    considered, only key profile slots within the 9-to-20 range are
//    assigned non-zero values. (This is necessary, since the input profile
//    is within this range as well.) We begin by assigning default values
//    of ZERO to everything.
//-----------------------------------------------------------------------------

static void generate_npc_profiles()
      {
      int key, shift, tpc, tpc_to_use, i;
      float majp[12];
      float minp[12];

      /* First we rearrange the key profile values (inputted in pitch height order) into lof order, C = 5 */
      for(i=0; i<12; i++) {
            majp[((((i * 7) % 12) + 5) % 12)] = major_profile[i];
            minp[((((i * 7) % 12) + 5) % 12)] = minor_profile[i];
            }
      for (key=0; key<56; ++key) {
            for(tpc=0; tpc<28; ++tpc) {
                  key_profile[key][tpc]=0;
                  }
            }

      for (key=9, shift=9; key<21; ++key, ++shift) {
            for (tpc=0; tpc<28; ++tpc) {
                  /* tpc_to_use is the profile slot to use for a given value of the key-profile. In this way we keep
                     all profile slots within the 9-to-20 range
                  */
                  if(tpc<9)
                        tpc_to_use=tpc+12;
                  else if(tpc>20)
                        tpc_to_use=tpc-12;
                  else
                        tpc_to_use = tpc;
                  if (tpc-shift >= -5 && tpc-shift <= 6) {
                        /* For example: for key 14 (C major) and tpc 17 (A),
                        read from profile step 5 + (17-14) = 8. For degree
                        6 of B major (key 19), degree 6 (22) is outside the
                        9-to-20 range, so tpc_to_use is 22-12=10; still
                        read from profile step 5 + (22-19) = 8. */

                        key_profile[key][tpc_to_use] = majp[5 + (tpc-shift)];
                        }
                  }
            }
      for (key = 37, shift = 9; key < 49; ++key, ++shift) {
            for (tpc = 0; tpc < 28; ++tpc) {
                  if (tpc<9)
                        tpc_to_use = tpc + 12;
                  else if (tpc > 20)
                        tpc_to_use = tpc - 12;
                  else
                        tpc_to_use = tpc;
                  if (tpc - shift >= -5 && tpc - shift <= 6 ) {
                        key_profile[key][tpc_to_use] = minp[5 + (tpc-shift)];
                        }
                  }
            }

/*
      for(key = 0; key < 56; ++key) {
            for(tpc=0; tpc<28; ++tpc) {
                  qDebug("%1.2f ", key_profile[key][tpc]);
                  }
            }
*/
      }

//---------------------------------------------------------
//   match_profiles
//    Here we generate the "key scores" - the local score for each key
//
//    Notice that we generate profiles for all 56 keys and all 28
//    tpc's, even in the case where npc profiles are being used. In
//    this case, though, all keys outside the allowable NPC range are
//    given large negative values.  And within both the key profiles
//    and the segment profiles, only TPC's within the range 9 to 20
//    have been given nonzero values.
//---------------------------------------------------------

static void match_profiles()
      {
      int key, tpc, s, best_key, i;
      qreal major_sumsq, minor_sumsq, input_sumsq;
      qreal kprob[56];

      for (key = 0; key < 56; ++key) {
            for (s = 0; s <= segtotal; ++s)
                  keyscore()[key].append(0.0);
            }

      if (scoring_mode==0) {
            major_sumsq = 0.0;
            minor_sumsq = 0.0;
            for(i=0; i<12; i++)
                  major_sumsq += major_profile[i]*major_profile[i];
            for(i=0; i<12; i++)
                  minor_sumsq += minor_profile[i]*minor_profile[i];
            if (verbosity==3)
                  qDebug("major_sumsq = %6.3f, minor_sumsq = %6.3f", major_sumsq, minor_sumsq);
            }

      qreal total_prob[segtotal + 1];

      for (s = 0; s <= segtotal; ++s) {
            if (scoring_mode==0) {
                  input_sumsq = 0.0;
                  for (i = 9; i <= 20; i++) {
                        input_sumsq += pow((seg_prof[i][s]-segment[s].average_dur), 2.0);
                        /* qDebug("%d X %6.3f squared is %6.3f", seg_prof[i][s], segment[s].average_dur, pow((seg_prof[i][s]-segment[s].average_dur), 2.0)); */
                        }
                  if (verbosity==3)
                        qDebug("For segment %d: average_dur = %6.3f; input_sumsq = %6.3f", s, segment[s].average_dur, input_sumsq);
                  }
            best_key=0;

            for (key=0; key<56; ++key) {
                  kprob[key] = 0.0;
                  keyscore()[key][s] = -1000000.0;
                  if (npc_or_tpc_profile==0 && (key<9 || (key>20 && key<37) || key>48))
                        continue;
                  kprob[key] = 1.0;
                  keyscore()[key][s] = 0.0;
                  for (tpc=0; tpc<28; ++tpc) {

          /*
             If scoring mode is 0, this is the K-S algorithm (this works for npc mode only). Segment
             profile values represent total duration of each pc (in all other cases, they're just 1
             for present pc's and 0 for absent ones). Key-profiles have been normalized linearly
             around the average key-profile value. We normalize the input values similarly by taking
             (seg_prof[tpc][s]-segment[s].average_dur). Then we multiply each normalized KP value by
             the normalized input value, and sum these products; this gives us the numerator of the
             correlation expression (as commented below). We've summed the squares of the normalized
             key-profile value (major_sumsq and minor_sumsq above) and the normalized input values
             (input_sumsq above), so this allows us to calculate the denominator also.

             If scoring_mode is 1, the key score is the sum of key-profile values for all pc's present
             (this is the algorithm used in CBMS)

             If scoring_mode is 2, calculate key scores as above, but divide each one by the number
             of pc's in the segment

             If scoring_mode is 3: for each key, add the log of the key-profile value for all present pc's;
             subtract values for all absent pc's. (This is the Bayesian approach; assume key-profiles
             represent pc distribution's in a corpus, i.e. the number of segments containing each scale
             degree)
          */

                        if(scoring_mode == 0) {
                              if(tpc<9 || tpc>20)
                                    continue;
                              /* calculate numerator */
                              keyscore()[key][s] += key_profile[key][tpc] * (seg_prof[tpc][s]-segment[s].average_dur);
                              /* qDebug("x-X=%6.3f, y-Y=%6.3f, product=%6.3f, new total=%6.3f", key_profile[key][tpc], seg_prof[tpc][s]-segment[s].average_dur, key_profile[key][tpc] * (seg_prof[tpc][s]-segment[s].average_dur), keyscore()[key][s]); */
                              }

                        if(scoring_mode==1 || scoring_mode==2)
                              keyscore()[key][s] += (key_profile[key][tpc] * seg_prof[tpc][s]);

                        if(scoring_mode == 3) {
                              /* if(tpc>11) continue; */
                              /* if(tpc<9 || tpc>20) continue; */

                              if(seg_prof[tpc][s]==0) {
                                    keyscore()[key][s] += log(1.000 - key_profile[key][tpc]);
                              /* qDebug("kp value = %6.3f: log(1-p) = %6.3f: score = %6.3f", key_profile[key][tpc], log(1.000 - key_profile[key][tpc]), keyscore()[key][s]); */
                                    if(tpc>=9 && tpc<=20)
                                          kprob[key] *= (1.000 - key_profile[key][tpc]);
                                    }
                              else {
                                    keyscore()[key][s] += log(key_profile[key][tpc]);
                                    if(tpc>=9 && tpc<=20)
                                          kprob[key] *= key_profile[key][tpc];
                                    }

                              /* qDebug("kp value = %6.3f: log(p) = %6.3f: score = %6.3f", key_profile[key][tpc], log(key_profile[key][tpc]), keyscore()[key][s]); */
                              }
                        }

                  if(scoring_mode == 0) {
                        /* qDebug("sqrt(major_sumsq * input_sumsq) = %6.3f", sqrt(major_sumsq * input_sumsq)); */
                        /* calculate denominator */
                        if(key<28)
                              keyscore()[key][s] = keyscore()[key][s] / sqrt(major_sumsq * input_sumsq);
                        else
                              keyscore()[key][s] = keyscore()[key][s] / sqrt(minor_sumsq * input_sumsq);
                        }
                  if(scoring_mode == 2) {
                        if(pc_tally[s] == 0)
                              keyscore()[key][s] = 0;
                        else
                              keyscore()[key][s] = keyscore()[key][s] / pc_tally[s];
                        }

                  /* if(s==0) qDebug("local score for key %d on segment %d: %6.3f", key, s, keyscore()[key][s]); */
                  if (keyscore()[key][s] > keyscore()[best_key][s])
                        best_key = key;
                  }

            if(verbosity>=2) {
                  qDebug("The best local key for segment %d at time %d is ", s, segment[s].start);
                  print_keyname(best_key);
                  qDebug("with score %6.3f", keyscore()[best_key][s]);
                  }

            if(scoring_mode==3) {
                  total_prob[s]=0.0;
                  for(key=0; key<56; key++) {
                        total_prob[s] += kprob[key] / 24.0;
                        /* qDebug("  Prob of segment %d given key %d: %6.8f", s, key, kprob[key]);  */
                        }

                  /* Now total_prob[s] is the total probability of the segment occurring: its probability given
                     a key, summed over all keys. But suppose we want to know the probability of ANY major triad
                     occurring? Then we have to multiply this by 12. (But not for something like a diminished
                     seventh which is symmetrical! */

                  if (verbosity>=3) {
                        qDebug("Best key for segment %d = %d, score = %6.8f", s, best_key, kprob[best_key]);
                        qDebug("Total (local) probability of segment %d: %6.8f", s, total_prob[s]);
                        }
                  }
            }
      }

//---------------------------------------------------------
//   choose_best_i
//---------------------------------------------------------

static void choose_best_i(int seg)
      {
      for (int j = 0; j < 56; ++j) {
            int k = 0;
            for (int i = 0; i < 56; ++i) {
                  if (analysis[i][j][seg] > analysis[k][j][seg])
                        k=i;
                  }
            /* For a given segment seg, and the key j at that segment,
               the best previous key is k
            */
            int size = best[j].size();
            for (int i = size; i < seg+1; ++i)
                  best[j].append(0);
            best[j][seg] = k;
            /* qDebug("For segment-%d-key %d, best segment-%d-key is %d, with score %d", seg, j, seg-1, k, analysis[k][j][seg]); */
            }
      }

//---------------------------------------------------------
//   make_first_table
//---------------------------------------------------------

static void make_first_table(int seg)
      {
      int i, j, s;
      qreal seg_factor, mod_factor, nomod_factor;

      if (scoring_mode==3) {
            mod_factor = log(change_penalty);
            nomod_factor = log(1.0 - change_penalty);
            seg_factor = 1.0;
            }
      else {
            mod_factor = -change_penalty;
            nomod_factor = 0.0;
            seg_factor = seglength;
            }

      for(s = 0; s <= segtotal; s++) {
            for(i = 0; i < 56; ++i) {
                  for(j = 0; j < 56; ++j)
                        analysis[i][j].append(-1000.0);
                  }
            }

      for(i = 0; i < 56; ++i) {
            for(j = 0; j < 56; ++j) {
                  if (j != i)
                        analysis[i][j][1] = ((keyscore()[i][0] + keyscore()[j][1]) * seg_factor) + mod_factor;
                  else
                        analysis[i][j][1] = ((keyscore()[i][0] + keyscore()[j][1]) * seg_factor) + nomod_factor;
                  }
            }
      choose_best_i(seg);
      }

//---------------------------------------------------------
//   make_tables
//---------------------------------------------------------

static void make_tables()
      {
      qreal seg_factor, mod_factor, nomod_factor;

      /* When scoring_mode = 3, the change_penalty represents the probability of changing key. So raising
         the penalty actually _increases_ the likelihood of modulations. */

      if (scoring_mode == 3) {
            mod_factor = log(change_penalty / 23.0);
            nomod_factor = log(1.0 - change_penalty);
            seg_factor = 1.0;
            }
      else {
            mod_factor = -change_penalty;
            nomod_factor = 0.0;
            seg_factor = seglength;
            }

      for (int seg = 2; seg <= segtotal; ++seg) {
            /* qDebug("mod_factor = %6.6f; ; nomod_factor = %6.6f", mod_factor, nomod_factor);  */
            for(int j = 0; j < 56; ++j) {
                  for(int i = 0; i < 56; ++i) {
                        int n = best[i][seg-1];
                        if (j != i)
                              analysis[i][j][seg] = analysis[n][i][seg-1] + (keyscore()[j][seg] * seg_factor) + mod_factor;
                        else
                              analysis[i][j][seg] = analysis[n][i][seg-1] + (keyscore()[j][seg] * seg_factor) + nomod_factor;
                        }
                  }
            choose_best_i(seg);
            }
      }

//---------------------------------------------------------
//   best_key_analysis
//---------------------------------------------------------

static void best_key_analysis()
      {
      int n, m, f, tie1=-1, tie2=-1;
      int s = segtotal;
      int k = 0;
      for(int j = 0; j < 56; ++j) {
            n = best[j][s];
            m = best[k][s];

            if (analysis[n][j][s] < analysis[m][k][s] + .001 && analysis[n][j][s] > analysis[m][k][s] - .001 && j!=k) {
                  tie1=j;
                  tie2=k;
                  }

            if (verbosity>1 && !(npc_or_tpc_profile == 0 && (j<9 || (j>20 && j<37) || j>48))) {
                  qDebug("Final score for ");
                  print_keyname(j);
                  /* qDebug("is %6.3f", analysis[n][j][s] * 1000 / (segment[segtotal].end - segment[0].start));   */
                  qDebug("is %6.3f", analysis[n][j][s]);
                  }
            if (analysis[n][j][s] > analysis[m][k][s] + .000001) {
                  /* The .000001 is to fix a strange bug: sometimes it thinks the conditional is satisfied in the case of ties */
                  k = j;   /* compute best key of final segment */
                  }
            }

      // To force a key choice at the final segment, insert key number here

      final[s] = k;
      if (verbosity > 1)
            if (k==tie1 || k==tie2)
                  qDebug("Tie at the end between %d and %d", tie1, tie2);

      // Here's where we take the best key choices and put them into final[s]

      for(s = segtotal; s >= 1; --s) {
            final[s-1] = best[k][s];
            k = final[s-1];
            }

      if (verbosity >= 2) {
            qDebug("Segment 0: key choice %d; total score %6.3f; segment score %6.3f", final[0], keyscore()[final[0]][0],
               keyscore()[final[0]][0]);
            qDebug("Segment 1: key choice %d; total score %6.3f; segment score %6.3f", final[1], analysis[final[0]][final[1]][1], analysis[final[0]][final[1]][1] - keyscore()[final[0]][0]);
            for(s = 2; s <= segtotal; s++) {
                  qDebug("Segment %d: key choice %d; total score %6.3f; segment score %6.3f", s, final[s], analysis[final[s-1]][final[s]][s], analysis[final[s-1]][final[s]][s] - analysis[final[s-2]][final[s-1]][s-1]);
                  }
            }

      if (verbosity > 1)
            qDebug("'Key-fit' scores for preferred analysis:");

      /* This routine calculates the key-fit scores for the final analysis. These are really per-second scores.
         Key-profile scores are not multiplied by seglength (as they would be in actually computing the analyses);
         change penalties are divided by seglength. */

      for (s = 0; s <= segtotal; ++s) {
            f = final[s];
            if (s > 0 && final[s] != final[s-1])
                  finalscore()[s] = (keyscore()[f][s]) - (change_penalty / seglength);
            else
                  finalscore()[s]=keyscore()[f][s];
            if (verbosity > 1) {
                  qDebug(" segment %d: %6.3f", s, finalscore()[s]);
                  }
            }
      }

//---------------------------------------------------------
//   findKey
//---------------------------------------------------------

int findKey(MidiTrack* mt, TimeSigMap* sigmap)
      {
      int tpc_found, npc_found;

      if ((scoring_mode == 0 || scoring_mode == 3) && npc_or_tpc_profile == 1) {
            qDebug("Error: scoring mode %d requires an npc profile", scoring_mode);
            exit(1);
            }

      final_timepoint=0;

      tpc_found = 0;
      npc_found = 0;

      int lastTick = 0;
      const EventList el = mt->events();

      foreach (Event e, el) {
            if (e.type() != ME_NOTE)
                  continue;
            if (e.offtime() > lastTick)
                  lastTick = e.offtime();
            // For note input, generate TPC labels within the 9-to-20 range
            e.setTpc((((((e.pitch() % 12) * 7) % 12) + 5) % 12) + 9);
            note.append(e);
            }
      spell(note, 0);
      npc_found = 1;

      // create one segment for every measure
      for (int i = 0;; ++i) {
            int tick = sigmap->bar2tick(i, 0, 0);
            SBeat b;
            b.time = tick;
            sbeat.append(b);
            if (tick > lastTick)
                  break;
            }

      firstbeat       = sbeat.first().time;
      final_timepoint = sbeat.last().time;
      num_sbeats      = sbeat.size();

      numnotes   = note.size();
      numchords  = 0;

      if (note.empty()) {
            qDebug("Error: No notes in input.");
            return 0;
            }

      seglength  = (sbeat[1].time - sbeat[0].time) / 1000.0; /* define segment length as the length of the first segment */
      if (verbosity > 1)
            qDebug("seglength = %3.3f", seglength);

      create_segments();
      for (int i = 0; i < segtotal+1; ++i) {
            final.append(0);
            pc_tally.append(0);
            finalscore().append(0.0);
            }

      fill_segments();
      count_segment_notes();
      if (scoring_mode==0)
            prepare_profiles();
      if (npc_or_tpc_profile == 1)
            generate_tpc_profiles();
      else
            generate_npc_profiles();
      match_profiles();
      if (segtotal > 0) {
            make_first_table(1);
            make_tables();
            best_key_analysis();
            }

      QList<int> keys;
      for (int i = 0; i < 27; ++i)
            keys.append(0);
      for (int i = 0; i <= segtotal; ++i) {
            keys[final[i] % 27]++;        // fold major/minor
//            qDebug("key %d: %d  %d", i, final[i], (final[i] % 27) - 14);
            }
      int xkey   = 0;
      int xcount = 0;
      for (int i = 0; i < 27; ++i) {
            if (keys[i] > xcount) {
                  xcount = keys[i];
                  xkey   = i;
                  }
            }

      xkey -= 14;
      if (xkey < -7 || xkey > 7) {
            qDebug("findKey(): illegal key %d found", xkey);
            xkey = 0;
            }

      spell(note, xkey);      // spell again with found key

      //
      // clear all arrays
      //
      note.clear();
      segment.clear();
      sbeat.clear();
      for (int i = 0; i < 28; ++i)
            seg_prof[i].clear();
      for (int i = 0; i < 56; ++i) {
            keyscore()[i].clear();
            best[i].clear();
            for (int k = 0; k < 56; ++k)
                  analysis[i][k].clear();
            }
      final.clear();
      pc_tally.clear();
      finalscore().clear();

      return xkey;
      }
#endif

}