jazzparser.misc.raphsto

1 """Raphael and Stoddard's model for chord labelling from midi data. 2 3 An implementation of the models described in Function Harmonic Analysis Using 4 Probabilistic Models, Raphael and Stoddard, 2004. 5 6 In the future, this should be fitted into the tagger framework. For now it's 7 just a rough prototype, so I'm not bothering to work out how it will implement 8 the tagger interface, etc. 9 10 """ 11 """ 12 ============================== License ======================================== 13 Copyright (C) 2008, 2010-12 University of Edinburgh, Mark Granroth-Wilding 14 15 This file is part of The Jazz Parser. 16 17 The Jazz Parser is free software: you can redistribute it and/or modify 18 it under the terms of the GNU General Public License as published by 19 the Free Software Foundation, either version 3 of the License, or 20 (at your option) any later version. 21 22 The Jazz Parser is distributed in the hope that it will be useful, 23 but WITHOUT ANY WARRANTY; without even the implied warranty of 24 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 25 GNU General Public License for more details. 26 27 You should have received a copy of the GNU General Public License 28 along with The Jazz Parser. If not, see <http://www.gnu.org/licenses/>. 29 30 ============================ End license ====================================== 31 32 """ 33 __author__ = "Mark Granroth-Wilding <mark.granroth-wilding@ed.ac.uk>" 34 35 import numpy, os, re 36 from numpy import ones, float64, sum as array_sum, zeros 37 import cPickle as pickle 38 from datetime import datetime 39 40 from jazzparser.utils.nltk.ngram import NgramModel 41 from jazzparser.utils.nltk.probability import mle_estimator, logprob, add_logs, \ 42 sum_logs, prob_dist_to_dictionary_prob_dist, \ 43 cond_prob_dist_to_dictionary_cond_prob_dist 44 from jazzparser import settings 45 from . import constants 46 from .midi import MidiHandler 47 48 from nltk.probability import ConditionalProbDist, FreqDist, \ 49 ConditionalFreqDist, DictionaryProbDist, \ 50 DictionaryConditionalProbDist, MutableProbDist 51 52 FILE_EXTENSION = "mdl"

53 54 -def format_state(state, time=None):

55 """ 56 Formats a state label tuple as a string for output. 57 58 """ 59 tonic,mode,chord = state 60 tonic_str = constants.TONIC_NAMES[tonic] 61 mode_str = constants.MODE_NAMES[mode] 62 chord_str = constants.CHORD_NAMES[chord] 63 return "%s %s: %s" % (tonic_str, mode_str, chord_str)

64

65 -def format_state_as_chord(state, time=None):

66 """ 67 Formats a state label tuple, like L{format_state}, but produces a chord 68 name. This contains less information than the result of L{format_state}, 69 but is easier to read. 70 71 """ 72 tonic,mode,chord = state 73 scale_chord_root = constants.CHORD_NOTES[mode][chord][0] 74 chord_root = (tonic+scale_chord_root) % 12 75 # Work out the triad type of this chord 76 triad_name = constants.TRIAD_TYPE_SYMBOLS[constants.SCALE_TRIADS[mode][chord]] 77 return "%s%s" % (constants.TONIC_NAMES[chord_root], triad_name)

78

79 -def format_state_as_raphsto(state, time=None):

80 """ 81 State representation modeled on the original R&S model's labels that are 82 inserted into the example MIDI files. 83 84 """ 85 tonic,mode,chord = state 86 # First part is the roman numeral chord name 87 roman_chord = constants.CHORD_NAMES[chord].ljust(4) 88 # Then the key name 89 tonic_name = constants.TONIC_NAMES[tonic].lower().ljust(2) 90 # including mode 91 mode_name = constants.MODE_SHORT_NAMES[mode] 92 # Put the time if available 93 if time is None: 94 time = "" 95 else: 96 time = " ms = %d" % time 97 # Work out the indentation on the basis of a line of fifths, starting at C 98 line_of_fifths = [ 0, 7, 2, 9, 4, 11, 6, 1, 8, 3, 10, 5 ] 99 pos = line_of_fifths.index(tonic) 100 padding = " " * (pos*2) 101 return "%s%s(%s%s)%s" % (padding, roman_chord, tonic_name, mode_name, time)

102

103 -def states_to_key_transition(state, previous_state):

104 """ 105 Takes a state and previous state representation and returns the 106 (state, previous state) pair that is used by the model. This mapping 107 effectively implements the parameter tying of the first part of the 108 transition distribution. 109 110 """ 111 tonic, mode, chord = state 112 ptonic, pmode, pchord = previous_state 113 return ( ((tonic-ptonic)%12, mode), pmode )

114

115 -def raphsto_d(pitch_class, state):

116 """ 117 The function described in the paper as d. Maps a pitch class to an 118 abstract representation that depends on the state. 119 120 This gets called a bajillion times, so needs to be as efficient as 121 possible. 122 123 """ 124 tonic, mode, chord = state 125 # Get the pitch class relative to the tonic 126 key_pitch_class = (pitch_class - tonic) % 12 127 128 # Find where the note is in the chord denoted by the state 129 try: 130 chord_pos = constants.CHORD_NOTES[mode][chord].index(key_pitch_class) 131 except ValueError: 132 # Not in the chord 133 # Try the scale 134 scale = constants.SCALES[mode] 135 if key_pitch_class in scale: 136 # Not in chord, but in scale 137 return 3 138 else: 139 # Not in chord or scale 140 return 4 141 else: 142 if chord_pos == 0: 143 # Root of chord 144 return 0 145 elif chord_pos == 1: 146 # Third of chord 147 return 1 148 else: 149 # Fifth of chord, or beyond (e.g. seventh) 150 return 2

151

152 153 -class RaphstoHmm(NgramModel):

154 """ 155 Hidden Markov Model that implements the model described in the paper. 156 157 States are in the form of a tuple (tonic,mode,chord) where tonic is in 158 \{0, ..., 11\}; mode is one of \{L{constants.MODE_MAJOR}, 159 L{constants.MODE_MINOR}\}; chord is one of \{L{constants.CHORD_I}, 160 ..., L{constants.CHORD_VII}\}. 161 162 Emissions are in the form of a list of pairs (pc,r), where pc is a pitch 163 class (like C{tonic} above) and r is an onset time abstraction and is 164 one of \{0, ...,3\}. 165 166 Unlike with NgramModel, the emission domain is the domain of values 167 from which each element of an emission is selected. In other words, 168 the actually domain of emissions is the powerset of C{emission_dom}. 169 170 In the description of the model, r is described as a condition of the 171 emission distribution. Although the model is truly replicated here, the 172 interface suggests otherwise, since we treat the rythmic markers as if 173 they're part of the emissions. From a conceptual point of view, this makes 174 more sense and I think it's rather odd that the model doesn't treat them 175 this way. 176 177 As for prior distributions (start state distribution), we ignore 178 the tonic of the first state - it doesn't make any sense to look 179 at it since the model is pitch-invariant throughout. We then 180 just use our marginalized chord distribution and assume that the 181 mode distribution is uniform (there are only two and it probably 182 won't make much difference). 183 184 @note: B{mutable distributions}: if you use mutable distributions for 185 transition or emission distributions, make sure you invalidate the cache 186 by calling L{clear_cache} after updating the distributions. Various 187 caches are used to speed up retreival of probabilities. If you fail to 188 do this, you'll end up getting some values unpredictably from the old 189 distributions 190 191 """ 192 V = { 193 0 : 1, 194 1 : 1, 195 2 : 1, 196 3 : 4, 197 4 : 5 198 } 199 """ This is the function (mapping) described in the model as V. """ 200 LABEL_DOM = None 201

202 - def __init__(self, key_transition_dist, chord_transition_dist, \ 203 emission_dist, chord_dist, model_name="default", 204 history="", description="", chord_set="scale+dom7"):

205 # We override the init completely because we need a different set of 206 # distributions 207 208 # HMM model: order 2 ngram 209 self.order = 2 210 211 self.chord_set = chord_set 212 213 # We define the domain for the distributions here, because they're fixed 214 self.key_transition_dom = [(tonic,mode) for tonic in range(12) for mode in constants.MODES] 215 self.chord_transition_dom = constants.CHORD_SETS[chord_set] 216 self.emission_dom = [(pc,rhythm) for pc in range(12) for rhythm in range(4)] 217 # The actual domain of the stored emission distribution is the domain 218 # of the d-function (i.e. 0-4) 219 self.emission_dist_dom = range(5) 220 self.beat_dom = range(4) 221 # This label_dom is not used for the distributions, since the 222 # transition distribution is split up into two and has tied parameters 223 # We still need to know all possible labels for decoding, though 224 self.label_dom = self.get_label_dom(chord_set=chord_set) 225 226 self.num_key_transitions = len(self.key_transition_dom) 227 self.num_chord_transitions = len(self.chord_transition_dom) 228 self.num_labels = len(self.label_dom) 229 self.num_emissions = len(self.emission_dom) 230 231 self.key_transition_dist = key_transition_dist 232 self.chord_transition_dist = chord_transition_dist 233 self.emission_dist = emission_dist 234 self.chord_dist = chord_dist 235 236 # We don't use backoff with this kind of model, so this is always None 237 self.backoff_model = None 238 239 # Name to save the model under 240 self.model_name = model_name 241 242 # Store a string with information about training, etc 243 self.history = history 244 # Store another string with an editable description 245 self.description = description 246 247 # Initialize the various caches 248 # These will be filled as we access probabilities 249 self.clear_cache()

250

251 - def label(self, handler):

252 """ 253 Produces labels for the midi data using the model. 254 255 Input is given in the form of a 256 L{jazzparser.misc.raphsto.midi.MidiHandler} instance. 257 258 Uses Viterbi on the model to get a state sequence and returns a list 259 containing a (state,time) pair for each state change. 260 261 """ 262 emissions,times = handler.get_emission_stream() 263 # Decode using viterbi to get a list of states 264 states = self.viterbi_decode(emissions) 265 266 # Remove states that don't change from the previous 267 last_state = None 268 state_changes = [] 269 for state,time in zip(states,times): 270 if state != last_state: 271 state_changes.append((state,time)) 272 last_state = state 273 return state_changes

274 275 @classmethod

276 - def get_label_dom(cls, chord_set="scale+dom7"):

277 if cls.LABEL_DOM is not None: 278 # Use the predefined label domain (used by subclasses) 279 return cls.LABEL_DOM 280 else: 281 # Define the label domain from the component domains 282 return [(tonic,mode,chord) for tonic in range(12) \ 283 for mode in constants.MODES \ 284 for chord in constants.CHORD_SETS[chord_set]]

285 286 @staticmethod

287 - def get_trainer():

288 from .train import RaphstoBaumWelchTrainer 289 return RaphstoBaumWelchTrainer

290

291 - def clear_cache(self):

292 """ 293 Initializes or empties probability distribution caches. 294 295 Make sure to call this if you change or update the distributions. 296 297 """ 298 # Whole emission identity 299 self._emission_cache = {} 300 # Whole transition identity 301 self._transition_cache = {} 302 # Single note observation 303 self._note_emission_cache = {}

304

305 - def add_history(self, string):

306 """ 307 Adds a line to the end of this model's history string. 308 309 """ 310 self.history += "%s: %s\n" % (datetime.now().isoformat(' '), string)

311

312 - def sequence_to_ngram(self, seq):

313 if len(seq) > 1: 314 raise RaphstoHmmError, "sequence_to_ngram() got a sequence with "\ 315 "more than one value. This shouldn't happen on an HMM" 316 elif len(seq) == 0: 317 return None 318 else: 319 return seq[0]

320

321 - def ngram_to_sequence(self, ngram):

322 if ngram is None: 323 return [] 324 else: 325 return [ngram]

326

327 - def last_label_in_ngram(self, ngram):

328 return ngram

329

330 - def backoff_ngram(self, ngram):

331 raise RaphstoHmmError, "backoff_ngram() should never be called on "\ 332 "a RaphstoHmm, since we don't use backoff models"

333 334 @staticmethod

335 - def train(*args, **kwargs):

336 """ 337 We don't train a RaphstoHmm using the C{train} method, since our 338 training procedure is not the same as the superclass, so this would 339 be confusing, as this method would require completely different 340 input. 341 342 We train our models by initializing in some way (usually hand-setting 343 of parameters), then using Baum-Welch on unlabelled data. 344 345 This method will just raise a RaphstoHmmError. 346 347 """ 348 raise RaphstoHmmError, "don't use train() to train a RaphstoHmm. "\ 349 "Instead initialize and then train unsupervisedly"

350 351 @classmethod

352 - def initialize_chord_types(cls, probs, model_name="default", chord_set="scale+dom7"):

353 """ 354 Creates a new model with the distributions initialized naively to 355 favour simple chord-types, as R&S do in the paper. They don't say 356 what values they use for C{probs}, except that they're high, medium 357 and low respectively. 358 359 The transition distribution is initialized so that everything is 360 equiprobable. 361 362 @type probs: 3-tuple of floats 363 @param probs: probability mass to assign to (0.) chord notes, (1.) 364 scale notes and (2.) other notes. The three values should sum to 365 1.0 (but will be normalized to if they don't) 366 367 """ 368 prob_sum = sum(probs) 369 probs = [p/prob_sum for p in probs] 370 371 # Create a probability distribution for the emission 372 # distribution 373 dists = {} 374 # Create the distribution for each possible r-value 375 for r in range(4): 376 probabilities = {} 377 for d in [0,1,2]: 378 probabilities[d] = probs[0]/3.0 379 probabilities[3] = probs[1] 380 probabilities[4] = probs[2] 381 dists[r] = DictionaryProbDist(probabilities) 382 emission_dist = DictionaryConditionalProbDist(dists) 383 384 # These distributions will make everything equiprobable 385 key_transition_counts = ConditionalFreqDist() 386 chord_transition_counts = ConditionalFreqDist() 387 chord_counts = {} 388 # Get all possible labels 389 label_dom = cls.get_label_dom(chord_set=chord_set) 390 391 for label0 in label_dom: 392 for label1 in label_dom: 393 key,pkey = states_to_key_transition(label1, label0) 394 # Give one count to the key transition corresponding to this state transition 395 key_transition_counts[pkey].inc(key) 396 # And one to the chord transition corresponding to this state transition 397 if label0[0] == label1[0] and label0[1] == label1[1]: 398 # tonic = tonic', mode = mode' 399 chord_transition_counts[label0[2]].inc(label1[2]) 400 else: 401 chord_counts.setdefault(label1[2], 0) 402 chord_counts[label1[2]] += 1 403 404 # Estimate distributions from these frequency distributions 405 key_dist = ConditionalProbDist(key_transition_counts, mle_estimator, None) 406 chord_trans_dist = ConditionalProbDist(chord_transition_counts, mle_estimator, None) 407 chord_dist = DictionaryProbDist(chord_counts) 408 # Sample these to get dictionary prob dists 409 key_dist = cond_prob_dist_to_dictionary_cond_prob_dist(key_dist) 410 chord_trans_dist = cond_prob_dist_to_dictionary_cond_prob_dist(chord_trans_dist) 411 chord_dist = prob_dist_to_dictionary_prob_dist(chord_dist) 412 413 model = cls(key_dist, \ 414 chord_trans_dist, \ 415 emission_dist, \ 416 chord_dist, \ 417 model_name=model_name, 418 chord_set=chord_set) 419 model.add_history(\ 420 "Initialized model '%s' to chord type probabilities, using "\ 421 "parameters: %s, %s, %s" % (model_name, probs[0], probs[1], probs[2])) 422 return model

423 424 @classmethod

425 - def initialize_existing_model(cls, old_model_name, model_name="default"):

426 """ 427 Initializes a model using parameters from an already trained model. 428 429 """ 430 # Load the trained model 431 model = cls.load_model(old_model_name) 432 # Change the model name and save it to a new file 433 model.model_name = model_name 434 # Continue the history from the old model and note the change of name 435 model.add_history("Parameters from '%s' used as initialization for "\ 436 "model '%s'" % (old_model_name,model_name)) 437 model.save() 438 # Now continue to use this model under its new name 439 return model

440

441 - def set_chord_transition_probabilities(self, spec):

442 """ 443 Sets the parameters of the chord transition distribution. This is used 444 in initialization. The parameters are extracted from a string: this is 445 so that it can be specified in a script option. 446 447 The required format of the string is a comma-separated list of 448 parameters given as C0->C1-P, where C0 and C1 are chords (I, II, etc) 449 that are in the model's distribution and P is a float probability. 450 Parameters not specified will be evenly distributed the remaining 451 probability mass. 452 453 """ 454 params = {} 455 param_re = re.compile(r'(?P<chord0>.+)->(?P<chord1>.+)-(?P<prob>.+)') 456 chord_ids = dict((name,num) for (num,name) in constants.CHORD_NAMES.items()) 457 458 def _chord_id(name): 459 # Get the id for the named chord 460 if name not in chord_ids: 461 raise RaphstoHmmParameterError, "unrecognised chord name '%s' "\ 462 "in parameter spec: %s" % (name,spec) 463 cid = chord_ids[name] 464 if cid not in self.chord_transition_dom: 465 raise RaphstoHmmParameterError, "chord %s is not used with this "\ 466 "model (in parameter spec: %s)" % (name,spec) 467 return cid

468 469 for param_str in spec.split(","): 470 # Pull out the bits of the parameter specification 471 match = param_re.match(param_str.strip()) 472 if not match: 473 raise RaphstoHmmParameterError, "could not parse parameter "\ 474 "spec: %s (in: %s)" % (param_str, spec) 475 parts = match.groupdict() 476 chord0 = _chord_id(parts['chord0']) 477 chord1 = _chord_id(parts['chord1']) 478 try: 479 prob = float(parts['prob']) 480 except ValueError: 481 raise RaphstoHmmParameterError, "not a valid probability: %s "\ 482 "(in %s)" % (parts['prob'], spec) 483 # Store the parameter value 484 params.setdefault(chord0, {})[chord1] = prob 485 486 # Set the values in the transition distribution 487 dists = {} 488 for chord0 in self.chord_transition_dom: 489 dist_params = {} 490 if chord0 not in params: 491 # Not given in the spec: uniform distribution 492 uniform_mass = 1.0 / len(self.chord_transition_dom) 493 for chord1 in self.chord_transition_dom: 494 dist_params[chord1] = uniform_mass 495 else: 496 # Work out the prob mass to be distributed among unspecified parameters 497 not_given = len(self.chord_transition_dom) - len(params[chord0]) 498 if not_given > 0: 499 given_mass = sum(params[chord0].values(), 0.0) 500 uniform_mass = (1.0 - given_mass) / not_given 501 else: 502 uniform_mass = 0.0 503 # Calculate the whole distribution 504 for chord1 in self.chord_transition_dom: 505 if chord1 in params[chord0]: 506 dist_params[chord1] = params[chord0][chord1] 507 else: 508 dist_params[chord1] = uniform_mass 509 dists[chord0] = DictionaryProbDist(dist_params) 510 # Use this distribution instead of what's already there 511 self.chord_transition_dist = DictionaryConditionalProbDist(dists) 512 513 self.add_history("Set chord transition distribution using "\ 514 "parameters: %s" % spec)

515

516 - def retrain_unsupervised(self, *args, **kwargs):

517 """ 518 Unsupervised training. Passes straight over to the 519 L{train.RaphstoBaumWelchTrainer}. 520 521 @see: jazzparser.misc.raphsto.train.RaphstoBaumWelchTrainer 522 523 """ 524 from .train import RaphstoBaumWelchTrainer 525 trainer = RaphstoBaumWelchTrainer(self) 526 trainer.train(*args, **kwargs)

527

528 - def transition_log_probability(self, state, previous_state):

529 # If we're transitioning to None, the sequence is over 530 # We only ask this probability when we know the sequence has finished, 531 # so the probability of the end state is 1.0. 532 if state is None: 533 return 0.0 534 cache_key = (state, previous_state) 535 if cache_key not in self._transition_cache: 536 tonic, mode, chord = state 537 if previous_state is None: 538 # Start state: prior state distribution 539 # The tonic doesn't affect this: we ignore it and distribute 540 # prob mass uniformly (divide by 12) 541 # We also ignore the mode: we could train a mode 542 # distribution, but at the moment we don't bother (divide by 2) 543 chord_transition_prob = self.chord_dist.logprob(chord) 544 self._transition_cache[cache_key] = chord_transition_prob - 24.0 545 else: 546 ptonic, pmode, pchord = previous_state 547 # The first part is the key transition 548 # p(t'-t, m' | m) 549 key,pkey = states_to_key_transition(state, previous_state) 550 key_transition_prob = self.key_transition_dist[pkey].logprob(key) 551 # The second part is the chord transition 552 if tonic == ptonic and mode == pmode: 553 # p(c' | c) 554 chord_transition_prob = self.chord_transition_dist[pchord].logprob(chord) 555 else: 556 # p(c') 557 chord_transition_prob = self.chord_dist.logprob(chord) 558 self._transition_cache[cache_key] = \ 559 key_transition_prob + chord_transition_prob 560 return self._transition_cache[cache_key]

561

562 - def emission_log_probability(self, emission, state):

563 if len(emission) == 0: 564 # No notes emitted 565 # Not clear from paper what the probability of this should be 566 # Give it zero for now 567 return float('-inf') 568 cache_key = (tuple(sorted(emission)), state) 569 if cache_key not in self._emission_cache: 570 prob = 0.0 571 # Take the product of the probability for each note in the set 572 for pc,beat in emission: 573 # Get the value of the d function for this note 574 d = raphsto_d(pc, state) 575 576 note_cache_key = (beat, d) 577 if note_cache_key not in self._note_emission_cache: 578 # Now get the probability of this d value, conditioned on the r value 579 # We divide by V, to spread the probability over other notes 580 # that share the same d value 581 self._note_emission_cache[note_cache_key] = \ 582 self.emission_dist[beat].logprob(d) - RaphstoHmm.V[d] 583 prob += self._note_emission_cache[note_cache_key] 584 self._emission_cache[cache_key] = prob 585 return self._emission_cache[cache_key]

586 587

588 - def forward_log_probabilities(self, sequence, normalize=True):

589 """We override this to provide a faster implementation. 590 591 It might also be possible to speed up the superclass' implementation 592 using numpy, but it's easier here because we know we're using an 593 HMM, not a higher-order ngram. 594 595 This is based on the fwd prob calculation in NLTK's HMM implementation. 596 597 Returns a numpy 2d array instead of a list of lists. 598 599 """ 600 T = len(sequence) 601 N = len(self.label_dom) 602 alpha = numpy.zeros((T, N), numpy.float64) 603 604 # Prepare the first column of the matrix: probs of all states in the 605 # first timestep 606 for i,state in enumerate(self.label_dom): 607 alpha[0,i] = self.transition_log_probability(state, None) + \ 608 self.emission_log_probability(sequence[0], state) 609 610 # Iterate over the other timesteps 611 for t in range(1, T): 612 for j,sj in enumerate(self.label_dom): 613 # Multiply each previous state's prob by the transition prob 614 # to this state and sum them all together 615 log_probs = [ 616 alpha[t-1, i] + self.transition_log_probability(sj, si) \ 617 for i,si in enumerate(self.label_dom)] 618 # Also multiply this by the emission probability 619 alpha[t, j] = sum_logs(log_probs) + \ 620 self.emission_log_probability(sequence[t], sj) 621 # Normalize by dividing all values by the total probability 622 if normalize: 623 for t in range(T): 624 total = sum_logs(alpha[t,:]) 625 for j in range(N): 626 alpha[t,j] -= total 627 return alpha

628

629 - def backward_log_probabilities(self, sequence, normalize=True):

630 """We override this to provide a faster implementation. 631 632 @see: forward_log_probability 633 634 Returns a numpy 2d array instead of a list of lists. 635 636 """ 637 T = len(sequence) 638 N = len(self.label_dom) 639 beta = numpy.zeros((T, N), numpy.float64) 640 641 # Initialize 642 beta[T-1, :] = numpy.log2(1.0/N) 643 644 # Iterate backwards over the other timesteps 645 for t in range(T-2, -1, -1): 646 for i,si in enumerate(self.label_dom): 647 # Multiply each next state's prob by the transition prob 648 # from this state to that and the emission prob in that next 649 # state 650 log_probs = [ 651 beta[t+1, j] + self.transition_log_probability(sj, si) + \ 652 self.emission_log_probability(sequence[t+1], sj) \ 653 for j,sj in enumerate(self.label_dom)] 654 beta[t, i] = sum_logs(log_probs) 655 # Normalize by dividing all values by the total probability 656 if normalize: 657 total = sum_logs(beta[t,:]) 658 for j in range(N): 659 beta[t,j] -= total 660 return beta

661 662

663 - def normal_forward_probabilities(self, sequence):

664 """If you want the normalized matrix of forward probabilities, it's 665 ok to use normal (non-log) probabilities and these can be computed 666 more quickly, since you don't need to sum logs (which is time 667 consuming). 668 669 Returns the matrix, and also the vector of values that each timestep 670 was divided by to normalize (i.e. total probability of each timestep 671 over all states). 672 Also returns the total log probability of the sequence. 673 674 @return: (matrix,normalizing vector,log prob) 675 676 """ 677 T = len(sequence) 678 N = len(self.label_dom) 679 alpha = numpy.zeros((T, N), numpy.float64) 680 scale = numpy.zeros(T, numpy.float64) 681 682 # Prepare the first column of the matrix: probs of all states in the 683 # first timestep 684 for i,state in enumerate(self.label_dom): 685 alpha[0,i] = self.transition_probability(state, None) * \ 686 self.emission_probability(sequence[0], state) 687 # Normalize by dividing all values by the total probability 688 total = sum(alpha[0,:]) 689 for i in range(N): 690 alpha[0,i] /= total 691 scale[0] = total 692 693 # Iterate over the other timesteps 694 for t in range(1, T): 695 for j,sj in enumerate(self.label_dom): 696 # Multiply each previous state's prob by the transition prob 697 # to this state and sum them all together 698 log_prob = sum( 699 (alpha[t-1, i] * self.transition_probability(sj, si) \ 700 for i,si in enumerate(self.label_dom)), 0.0) 701 # Also multiply this by the emission probability 702 alpha[t, j] = log_prob * \ 703 self.emission_probability(sequence[t], sj) 704 # Normalize by dividing all values by the total probability 705 total = sum(alpha[t,:]) 706 for j in range(N): 707 alpha[t,j] /= total 708 scale[t] = total 709 710 # Multiply together the probability of each timestep to get the whole 711 # probability of the sequence 712 log_prob = sum((logprob(total) for total in scale), 0.0) 713 return alpha,scale,log_prob

714

715 - def normal_backward_probabilities(self, sequence):

716 """ 717 @see: normal_forward_probabilities 718 719 (except that this doesn't return the logprob) 720 721 """ 722 T = len(sequence) 723 N = len(self.label_dom) 724 beta = numpy.zeros((T, N), numpy.float64) 725 scale = numpy.zeros(T, numpy.float64) 726 727 # Initialize 728 beta[T-1, :] = 1.0/N 729 scale[T-1] = 1.0 730 731 # Iterate backwards over the other timesteps 732 for t in range(T-2, -1, -1): 733 for i,si in enumerate(self.label_dom): 734 # Multiply each next state's prob by the transition prob 735 # from this state to that and the emission prob in that next 736 # state 737 beta[t, i] = sum( 738 (beta[t+1, j] * self.transition_probability(sj, si) * \ 739 self.emission_probability(sequence[t+1], sj) \ 740 for j,sj in enumerate(self.label_dom)), 0.0) 741 # Normalize by dividing all values by the total probability 742 total = sum(beta[t,:]) 743 for j in range(N): 744 beta[t,j] /= total 745 scale[t] = total 746 return beta,scale

747

748 - def compute_gamma(self, sequence, forward=None, backward=None):

749 """ 750 Computes the gamma matrix used in Baum-Welch. This is the matrix 751 of state occupation probabilities for each timestep. It is computed 752 from the forward and backward matrices. 753 754 These can be passed in as 755 arguments to avoid recomputing if you need to reuse them, but will 756 be computed from the model if not given. They are assumed to be 757 the matrices computed by L{normal_forward_probabilities} and 758 L{normal_backward_probabilities} (i.e. normalized, non-log 759 probabilities). 760 761 """ 762 if forward is None: 763 forward = self.normal_forward_probabilities(sequence) 764 if backward is None: 765 backward = self.normal_backward_probabilities(sequence) 766 # T is the number of timesteps 767 # N is the number of states 768 T,N = forward.shape 769 770 gamma = zeros((T,N), float64) 771 for t in range(T): 772 for i in range(N): 773 gamma[t][i] = forward[t][i]*backward[t][i] 774 denominator = array_sum(gamma[t]) 775 # Normalize 776 for i in range(N): 777 gamma[t][i] /= denominator 778 return gamma

779

780 - def compute_xi(self, sequence, forward=None, backward=None):

781 """ 782 Computes the xi matrix used by Baum-Welch. It is the matrix of joint 783 probabilities of occupation of pairs of conecutive states: 784 P(i_t, j_{t+1} | O). 785 786 As with L{compute_gamma} forward and backward matrices can optionally 787 be passed in to avoid recomputing. 788 789 """ 790 if forward is None: 791 forward = self.normal_forward_probabilities(sequence) 792 if backward is None: 793 backward = self.normal_backward_probabilities(sequence) 794 # T is the number of timesteps 795 # N is the number of states 796 T,N = forward.shape 797 798 xi = zeros((T-1,N,N), float64) 799 for t in range(T-1): 800 total = 0.0 801 # For each first state 802 for i,statei in enumerate(self.label_dom): 803 # and each second state 804 for j,statej in enumerate(self.label_dom): 805 # get the probability of this pair of states given the emissions 806 prob = forward[t][i] * backward[t+1][j] * \ 807 self.transition_probability(statej, statei) * \ 808 self.emission_probability(sequence[t+1], statej) 809 xi[t][i][j] = prob 810 total += prob 811 # Normalize all these probabilities 812 for i in range(N): 813 for j in range(N): 814 xi[t][i][j] /= total 815 return xi

816

817 - def to_picklable_dict(self):

818 """ 819 Produces a picklable representation of model as a dict. 820 You can't just pickle the object directly because some of the 821 NLTK classes can't be pickled. You can pickle this dict and 822 reconstruct the model using NgramModel.from_picklable_dict(dict). 823 824 """ 825 from jazzparser.utils.nltk.storage import object_to_dict 826 return { 827 'key_transition_dist' : object_to_dict(self.key_transition_dist), 828 'chord_transition_dist' : object_to_dict(self.chord_transition_dist), 829 'emission_dist' : object_to_dict(self.emission_dist), 830 'chord_dist' : object_to_dict(self.chord_dist), 831 'history' : self.history, 832 'description' : self.description, 833 'chord_set' : self.chord_set, 834 }

835 836 @classmethod

837 - def from_picklable_dict(cls, data, model_name="default"):

838 """ 839 Reproduces an n-gram model that was converted to a picklable 840 form using to_picklable_dict. 841 842 """ 843 from jazzparser.utils.nltk.storage import dict_to_object 844 return cls(dict_to_object(data['key_transition_dist']), 845 dict_to_object(data['chord_transition_dist']), 846 dict_to_object(data['emission_dist']), 847 dict_to_object(data['chord_dist']), 848 model_name=model_name, 849 history=data.get('history', ''), 850 description=data.get('description', ''), 851 chord_set=data.get('chord_set', 'scale+dom7'))

852 853 @classmethod

854 - def _get_model_dir(cls):

855 return os.path.join(settings.MODEL_DATA_DIR, "raphsto")

856 857 @classmethod

858 - def _get_filename(cls, model_name):

859 return os.path.join(cls._get_model_dir(), "%s.%s" % (model_name, FILE_EXTENSION))

860 - def _get_my_filename(self):

861 return type(self)._get_filename(self.model_name)

862 _filename = property(_get_my_filename) 863 864 @classmethod

865 - def list_models(cls):

866 """ 867 Returns a list of the names of available models. 868 869 """ 870 model_dir = cls._get_model_dir() 871 if not os.path.exists(model_dir): 872 return [] 873 model_ext = ".%s" % FILE_EXTENSION 874 names = [name.rpartition(model_ext) for name in os.listdir(model_dir)] 875 return [name for name,ext,right in names if ext == model_ext and len(right) == 0]

876

877 - def save(self):

878 """ 879 Saves the model data to a file. 880 881 """ 882 data = pickle.dumps(self.to_picklable_dict(), 2) 883 filename = self._filename 884 # Check the directory exists 885 filedir = os.path.dirname(filename) 886 if not os.path.exists(filedir): 887 os.mkdir(filedir) 888 f = open(filename, 'w') 889 f.write(data) 890 f.close()

891

892 - def delete(self):

893 """ 894 Removes all the model's data. 895 896 """ 897 fn = self._filename 898 if os.path.exists(fn): 899 os.remove(fn)

900 901 @classmethod

902 - def load_model(cls, model_name):

903 filename = cls._get_filename(model_name) 904 # Load the model from a file 905 if os.path.exists(filename): 906 f = open(filename, 'r') 907 model_data = f.read() 908 model_data = pickle.loads(model_data) 909 f.close() 910 else: 911 raise RaphstoModelLoadError, "the model '%s' has not been trained" % model_name 912 return cls.from_picklable_dict(model_data, model_name=model_name)

913

914 915 916 -class RaphstoHmmThreeChord(RaphstoHmm):

917 """ 918 Modification of the Raphsto algorithm that allows it only to assign one 919 of three chords: I, IV and V. They say that secondary dominants are 920 treated using modulation, but in fact their model uses IIs, VIs, etc. 921 This version forces these things to be handled using modulation. 922 923 """

924 - def __init__(self, *args, **kwargs):

925 kwargs['chord_set'] = 'three-chord' 926 RaphstoHmm.__init__(self, *args, **kwargs)

927 928 @classmethod

929 - def _get_model_dir(cls):

930 # Use a different directory for these models 931 return os.path.join(settings.MODEL_DATA_DIR, "raphsto3c")

932

933 934 935 -class RaphstoHmmFourChord(RaphstoHmm):

936 """ 937 Like L{RaphstoHmmThreeChord}, but also allows a minor secondary dominant 938 chord (II). 939 940 """

941 - def __init__(self, *args, **kwargs):

942 kwargs['chord_set'] = 'four-chord' 943 RaphstoHmm.__init__(self, *args, **kwargs)

944 945 @classmethod

946 - def _get_model_dir(cls):

947 # Use a different directory for these models 948 return os.path.join(settings.MODEL_DATA_DIR, "raphsto4c")

949

950 -class RaphstoHmmUnigram(RaphstoHmm):

951 """ 952 Like L{RaphstoHmm}, but always gives all state transitions equal 953 probability, so that it is effectively a unigram model. 954 955 Train with L{jazzparser.misc.raphsto.train.RaphstoBaumWelchUnigramTrainer}. 956 957 """

958 - def __init__(self, *args, **kwargs):

959 super(RaphstoHmmUnigram, self).__init__(*args, **kwargs) 960 # Precompute the uniform transition probability 961 self._uniform_transition = 1.0 / len(self.label_dom) 962 self._uniform_transition_log = - logprob(len(self.label_dom))

963

964 - def transition_log_probability(self, state, previous_state):

965 # If we're transitioning to None, the sequence is over 966 if state is None: 967 return 0.0 968 else: 969 return self._uniform_transition_log

970

971 - def transition_probability(self, state, previous_state):

972 if state is None: 973 return 1.0 974 else: 975 return self._uniform_transition

976 977 @staticmethod

978 - def get_trainer():

979 from .train import RaphstoBaumWelchUnigramTrainer 980 return RaphstoBaumWelchUnigramTrainer

981 982 @classmethod

983 - def _get_model_dir(cls):

984 # Use a different directory for these models 985 return os.path.join(settings.MODEL_DATA_DIR, "raphstouni")

986

987 -class RaphstoHmmError(Exception):

988 pass

989

990 -class RaphstoHmmParameterError(Exception):

991 pass

992

993 -class RaphstoModelLoadError(Exception):

994 pass

995

996 -class RaphstoModelSaveError(Exception):

997 pass

998 999 MODEL_TYPES = { 1000 'standard' : RaphstoHmm, 1001 'three-chord' : RaphstoHmmThreeChord, 1002 'four-chord' : RaphstoHmmFourChord, 1003 'unigram' : RaphstoHmmUnigram, 1004 } 1005

Source Code for Package jazzparser.misc.raphsto