1function [grid D] = mvn_som_skl(cm, n, global_training, local_training)
2%MVN_SOM_SKL  Computes a NxN SOM using the MVN models and the
3%   Symmetric Kullback-Leibler divergence.
4%
5%The algorithm was published in 2010:
6%   ``Islands of Gaussians: The Self Organizing Map and Gaussian Music Similarity Features'', D. Schnitzer, A. Flexer, G. Widmer and M. Gasser, Proceedings of the 11th International Society for Music Information Retrieval Conference, 2010.
7%
8%   [grid D] = MVN_SOM_SKL(models, n) Select the models to compute
9%   the NxN SOM.
10%
11%   [grid D] = MVN_SOM_SKL(models, n, global_training, local_training)
12%   Select the models to compute the NxN SOM; global/local_training sets the global/local training iterations of the SOM
13
14%   MVN is (c) 2010-2011, Dominik Schnitzer, <dominik.schnitzer@ofai.at>
15%   <http://www.ofai.at/~dominik.schnitzer/mvn>
16%
17%   This file is part of the MVN Octave/Matlab Toolbox
18%
19%   MVN is free software: you can redistribute it and/or modify
20%   it under the terms of the GNU General Public License as published by
21%   the Free Software Foundation, either version 3 of the License, or
22%   (at your option) any later version.
23%
24%   MVN is distributed in the hope that it will be useful,
25%   but WITHOUT ANY WARRANTY; without even the implied warranty of
26%   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
27%   GNU General Public License for more details.
28%
29%   You should have received a copy of the GNU General Public License
30%   along with MVN.  If not, see <http://www.gnu.org/licenses/>.
31
32    cmsel = 1:length(cm);
33
34    % initialization
35    grid = [];
36    for i = 0:n-1
37        for j = 0:n-1
38            grid(i*n + (j+1)).x = i;
39            grid(i*n + (j+1)).y = j;
40            grid(i*n + (j+1)).n = [];
41        end
42    end
43
44    cmlen = length(cmsel);
45
46    mlen = n*n;
47
48    % RANDOM INIT
49    mperm = randperm(cmlen);
50    m = cm(mperm(1:mlen));
51
52    if (nargin <= 2)
53        global_training = 100;
54        local_training = 200;
55    end
56
57    for learniter = 1:global_training+local_training
58
59        if (learniter < global_training)
60
61            % select appropriate alpha
62            alpha = 0.9 * (1 - learniter/global_training);
63
64            % relatively large radius, decay to 1 during global learning
65            sigma = 1 + 2/3*n*(1 - learniter/global_training);
66            fprintf(1, 'global training; som iteration=%d, alpha=%f, sigma=%f\n', learniter,...
67                alpha, sigma);
68        else
69            local_learniter = learniter - global_training;
70
71            % attain alpha values in the range of 0.06 - 0.01
72            alpha = 0.01 + 0.05*(1 - local_learniter/local_training);
73
74            % only modify neighbours
75            sigma = 1.5;
76            fprintf(1, 'local training; som iteration=%d, alpha=%f, sigma=%f\n',...
77                local_learniter, alpha, sigma);
78        end
79
80        % select a random input
81        cm_now = randi(cmlen);
82
83        % compute the winner unit (c)
84        dc = realmax();
85        c = 0;
86        for i = 1:mlen
87            d = mvn_div_skl(cm(cmsel(cm_now)), m(i));
88            if (d < dc)
89                c = i;
90                dc = d;
91            end
92        end
93
94        for i = 1:mlen
95
96            e2dist = (grid(i).x - grid(c).x)^2 + (grid(i).y - grid(c).y)^2;
97            nkernel = alpha * exp(-e2dist/(2*sigma*sigma));
98
99
100            try
101                r = bregman_combine_r(m(i), cm(cmsel(cm_now)), nkernel);
102                l = bregman_combine_l(m(i), cm(cmsel(cm_now)), nkernel);
103
104                lrm = [mvn_new(l.cov, l.m), mvn_new(r.cov, r.m)];
105                sym = mvn_bregmancentroid_skl(lrm);
106                m(i).cov = sym.cov;
107                m(i).icov = sym.icov;
108                m(i).m = sym.m;
109
110            catch ex
111                fprintf(1, 'Skip this, this is dangerous (mainloop, i=%d)\n', i);
112                fprintf(1, 'Exception: %s\n', ex.message);
113                break;
114            end
115        end
116
117    end
118
119    % build the final grid
120    D = zeros(cmlen, mlen);
121    for i = 1:cmlen
122        fprintf(1, 'assignment iteration = %d\n', i);
123
124        dnn = realmax();
125        nn = 0;
126        for j = 1:mlen
127            D(i, j) = mvn_div_skl(cm(cmsel(i)), m(j));
128            if (D(i, j) < dnn)
129                dnn = D(i, j);
130                nn = j;
131            end
132        end
133        grid(nn).n = [grid(nn).n cmsel(i)];
134    end
135end
136
137
138function r = bregman_combine_r(x, y, ytheta)
139% compute right centroid mixture
140
141    r.m = (1-ytheta) * (x.icov * x.m) + (ytheta)*(y.icov * y.m);
142    r.cov = (1-ytheta) * (0.5 * x.icov) + (ytheta)*(0.5 * y.icov);
143    r.icov = 2 * r.cov;
144    r.cov = inv(r.icov);
145    r.m = r.cov * r.m;
146end
147
148function l = bregman_combine_l(x, y, ytheta)
149% compute left centroid mixture
150
151    l.m = (1-ytheta)*(x.m) + (ytheta)*(y.m);
152    l.cov = (1-ytheta)* (-(x.cov + x.m*x.m')) +...
153       (ytheta)*(-(y.cov + y.m*y.m'));
154    l.cov = -(l.cov + l.m*l.m');
155    l.icov = inv(l.cov);
156end
157