view yeti/silvet.yeti @ 20:982aa1197a7e

Getting there, slowly, sort of, with EM

program silvet;

{ prepareTimeFrequency } = load timefreq;
{ loadTemplates, extractRanges } = load templates;

em = load em;
em1 = load em_onecolumn;

mat = load may.matrix;
vec = load may.vector;
plot = load may.plot;

templates = loadTemplates ();

ranges = extractRanges templates;

eprintln "\nWe have \(length (keys templates)) instruments:";
for (sort (keys templates)) do k:
    eprintln " * \(k) \(mat.size templates[k]) range \(ranges[k].lowest) -> \(ranges[k].highest)";
done;
eprintln "";

columns = prepareTimeFrequency "test.wav";

/*
height = if empty? columns then 0 else vec.length (head columns) fi;

chunkSize = { rows = height, columns = 100 };

emdata = em.initialise ranges templates 88 chunkSize;

eprintln "initialised EM data: overall pitch range \(emdata.lowest) -> \(emdata.highest)";

chunkify cols = 
    if empty? cols then []
    else
       (mat.resizedTo chunkSize
           (mat.fromColumns (take chunkSize.columns cols)))
        :. \(chunkify (drop chunkSize.columns cols));
    fi;

chunks = chunkify columns;

eprintln "we have \(length chunks) chunks of size \(mat.size (head chunks))";

eprintln "attempting one expectation phase...";

error = em.performExpectation emdata (head chunks);

eprintln "done, result has dimension \(mat.size error)";

eprintln "attempting one maximisation phase...";

newP = em.performMaximisation emdata (head chunks) error;

eprintln "done";

\() (plot.plot [ Grid (head chunks) ]);
\() (plot.plot [ Grid error ]);

\() (plot.plot [ Grid newP ]);
*/

em1data = em1.initialise ranges templates 88;

col = head (drop 35 columns);

\() (plot.plot [ Caption "Source frequency distribution", Vector col ]);

sourceGrid d =
    mat.fromColumns (map do k: d.sources[k] done (sort (keys d.sources)));

\() (plot.plot [ Caption "Source distribution beforehand", Grid (sourceGrid em1data)]);

oneIteration em1data col n =
   (q = em1.performExpectation em1data col;
//    \() (plot.plot [ Caption "Frequency distribution and output of E-step for iteration \(n)", Vector col, Vector q ]);
    newdata = em1.performMaximisation em1data col q;
    if (n % 2 == 0) then
        \() (plot.plot [ Caption "Pitch distribution before and after M-step update for iteration \(n)", Vector (em1data.pitches), Vector (newdata.pitches) ]);
        \() (plot.plot [ Caption "Source distribution after M-step update for iteration \(n)", Grid (sourceGrid newdata) ]);
    fi;
    newdata);

iterations = 6;

var d = em1data;

for [1..iterations] do i:
    d := oneIteration d col i;
done;    

();