Chris@19: 
Chris@19: module em_onecolumn;
Chris@19: 
Chris@19: mm = load may.mathmisc;
Chris@19: vec = load may.vector;
Chris@19: mat = load may.matrix;
Chris@19: 
Chris@19: inRange ranges instrument note =
Chris@19:     note >= ranges[instrument].lowest and note <= ranges[instrument].highest;
Chris@19: 
Chris@19: initialise ranges templates notes =
Chris@19:    (instruments = keys ranges;
Chris@19:     {
Chris@19:         pitches = // z in the original. 1 per note
Chris@19:             vec.randoms notes,
Chris@19:         sources = // u in the original. 1 per note-instrument
Chris@19:             mapIntoHash id
Chris@19:                 do instrument:
Chris@19:                     vec.fromList
Chris@19:                        (map do note:
Chris@19:                             if inRange ranges instrument note then 1 else 0 fi
Chris@19:                         done [0..notes-1])
Chris@19:                 done instruments,
Chris@19:         instruments,
Chris@19:         templates,
Chris@19:         ranges,
Chris@19:         lowest = head
Chris@19:            (sort (map do i: ranges[i].lowest done instruments)),
Chris@19:         highest = head (reverse
Chris@19:            (sort (map do i: ranges[i].highest done instruments))),
Chris@19:     });
Chris@19: 
Chris@19: epsilon = 1e-16;
Chris@19: 
Chris@19: select predicate = concatMap do v: if predicate v then [v] else [] fi done;
Chris@19: 
Chris@19: distributionsFor data instrument note =
Chris@19:     {
Chris@19:         w = mat.getColumn note data.templates[instrument],
Chris@19:         p = vec.at data.pitches note,
Chris@19:         s = vec.at data.sources[instrument] note,
Chris@19:     };
Chris@19: 
Chris@19: performExpectation data column =
Chris@19:    (estimate = 
Chris@19:         fold do acc instrument:
Chris@19:             fold do acc note:
Chris@19:                 { w, p, s } = distributionsFor data instrument note;
Chris@19:                 vec.add [acc, vec.scaled (p * s) w];
Chris@19:             done acc [data.ranges[instrument].lowest .. 
Chris@19:                       data.ranges[instrument].highest]
Chris@19:         done (vec.consts epsilon (vec.length column)) data.instruments;
Chris@19:     vec.divide column estimate);
Chris@19: 
Chris@19: performMaximisation data column q =
Chris@19:    (pitchTop = vec.fromList
Chris@19:        (map do note:
Chris@19:             fold do acc instrument:
Chris@19:                 { w, p, s } = distributionsFor data instrument note;
Chris@19:                 fold do acc bin:
Chris@19:                     acc + s * (vec.at w bin) * (vec.at column bin);
Chris@19:                 done acc [0..vec.length column - 1]
Chris@19:             done epsilon data.instruments;
Chris@19:         done [data.lowest .. data.highest]);
Chris@19:     pitchBottom = vec.sum pitchTop;
Chris@19: 
Chris@19: /*    sourceTop = vec.fromList
Chris@19:        (map do instrument:
Chris@19:             fold do acc note:
Chris@19:                 { w, p, s } = distributionsFor data instrument note;
Chris@19:                 fold do acc bin:
Chris@19:                     acc + p * (vec.at w bin) * (vec.at column bin);
Chris@19:                 done acc [0..vec.length column - 1]
Chris@19:             done epsilon [data.lowest .. data.highest]
Chris@19:         done data.instruments);
Chris@19:     sourceBottom = vec.sum sourceTop;
Chris@19: */
Chris@19:     data with { 
Chris@19:         pitches = vec.divideBy pitchBottom pitchTop,
Chris@19: //        sources = vec.divideBy sourceBottom sourceTop
Chris@19:     });
Chris@19: 
Chris@19: {
Chris@19:     initialise,
Chris@19:     performExpectation,
Chris@19:     performMaximisation,
Chris@19: }
Chris@19: 
Chris@19: