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