Browse code

starting som prediction fine-tuned class-performance visualisation

git-svn-id: https://svn.discofish.de/MATLAB/spmtoolbox/SVMCrossVal@112 83ab2cfd-5345-466c-8aeb-2b2739fb922d

Christoph Budziszewski authored on 21/01/2009 16:34:25
Showing 147 changed files
1 1
new file mode 100644
... ...
@@ -0,0 +1,100 @@
1
+function classify(varargin)
2
+
3
+
4
+
5
+switch nargin
6
+    case 1
7
+        paramModel = varargin{1};
8
+        % PROJECT_BASE_PATH = 'D:\Analyze\Stimolos';
9
+        PROJECT_BASE_PATH = 'D:\Analyze\Choice\24pilot';
10
+        PROJECT_RESULT_PATH = 'results\SPM.mat';
11
+    otherwise
12
+        error('spmtoolbox:SVMCrossVal:arginError','Please Specify action and parameter model');
13
+end
14
+
15
+        
16
+        % common params
17
+        calculateParams  = struct;
18
+        calculateParams.smoothed        = getDouble(paramModel.txtSmoothed);
19
+
20
+        calculateParams.frameShiftStart = getDouble(paramModel.txtFrameShiftStart);  % -20;
21
+        calculateParams.frameShiftEnd   = getDouble(paramModel.txtFrameShiftEnd); %15;
22
+        calculateParams.decodeDuration  = getDouble(paramModel.txtFrameShiftDur);
23
+        calculateParams.psthStart       = getDouble(paramModel.txtPSTHStart); % -25;
24
+        calculateParams.psthEnd         = getDouble(paramModel.txtPSTHEnd); % 20;
25
+        calculateParams.baselineStart   = getDouble(paramModel.txtBaselineStart); % -22;
26
+        calculateParams.baselineEnd     = getDouble(paramModel.txtBaselineEnd); % -20;
27
+
28
+        calculateParams.svmargs         = get(paramModel.txtSVMopts,'String');
29
+        calculateParams.sessionList     = 1:3;
30
+
31
+        classStruct = parseClassDef(paramModel);
32
+        
33
+        
34
+        calculateParams.labelMap        = LabelMap(classStruct.labelCells , classStruct.conditionCells, 'auto'); % LabelMap({'<','>','<+<','>+>','<+>','>+<'},{-2,-1,1,2,3,4}); 0 is autolabel
35
+        calculateParams.classList       = getClasses(calculateParams.labelMap);
36
+        calculateParams.eventList       = classStruct.eventMatrix; %[9,11,13; 10,12,14];
37
+%         calculateParams.eventList       = getPSTEventMatrix(calculateParams.labelMap);
38
+        
39
+        subjectSelection = getSubjectIDString(paramModel);
40
+        decode = struct;
41
+        decode.decodePerformance = [];
42
+        decode.rawTimeCourse     = [];
43
+        
44
+        for subjectCell = subjectSelection
45
+            SubjectID = cell2mat(subjectCell);
46
+            namehelper = strcat('s',SubjectID); %Vars can not start with numbers.
47
+
48
+            display('loading SPM.mat ...');
49
+            spm = load(fullfile(PROJECT_BASE_PATH,SubjectID,PROJECT_RESULT_PATH));
50
+            display('... done.');
51
+
52
+            %% calculate
53
+            calculateParams.(namehelper).des             = spm.SPM;
54
+            calculateParams.(namehelper).voxelList       = parseVoxelList(paramModel,SubjectID);
55
+            assignin('base','calculateParams',calculateParams);
56
+
57
+            display(sprintf('calculating cross-validation performance time-shift for Subject %s. Please Wait. ...',SubjectID));
58
+            display('switching off all warnings');
59
+            warning_state               = warning('off','all');
60
+            display('calculating ...');
61
+            decode.(namehelper)         = calculateDecodePerformance(calculateParams,SubjectID);
62
+            display('... done');
63
+            display('restoring warnings');
64
+            warning(warning_state);
65
+            
66
+            decode.decodePerformance    = [decode.decodePerformance decode.(namehelper).decodePerformance];
67
+            decode.rawTimeCourse        = [decode.rawTimeCourse decode.(namehelper).rawTimeCourse];
68
+
69
+            assignin('base','decode',decode);
70
+        end
71
+
72
+        display('Finished calculations.');
73
+        display('Plotting...');
74
+
75
+        plotParams                  = struct;
76
+        plotParams.psthStart        = calculateParams.psthStart;
77
+        plotParams.psthEnd   =  calculateParams.psthEnd;
78
+        plotParams.nClasses  = length(calculateParams.classList);
79
+        
80
+        plotParams.frameShiftStart   = calculateParams.frameShiftStart;
81
+        plotParams.frameShiftEnd     = calculateParams.frameShiftEnd;
82
+        plotParams.decodePerformance = decode.decodePerformance;
83
+        plotParams.rawTimeCourse     = decode.rawTimeCourse;
84
+        
85
+        if numel(subjectSelection) == 1
86
+          plotParams.SubjectID         = SubjectID;
87
+        else
88
+          plotParams.SubjectID         = 'Multiple';
89
+        end
90
+
91
+        plotParams.smoothed          = boolToYesNoString(calculateParams.smoothed);
92
+         
93
+
94
+        assignin('base','plotParams',plotParams);
95
+%         plotDecodePerformance(params.psthStart,params.psthEnd,params.nClasses,decode.decodeTable,params.frameShiftStart,params.frameShiftEnd,decode.rawTimeCourse);
96
+        plotDecodePerformance(plotParams);
97
+            
98
+        display('all done.');
99
+
100
+    end
0 101
\ No newline at end of file
1 102
new file mode 100644
... ...
@@ -0,0 +1,103 @@
1
+function plotDecodePerformance(varargin)
2
+% plotDecodePerformance(timeline,decodePerformance,nClasses,rawData)
3
+
4
+PSTH_AXIS_MIN = -1;
5
+PSTH_AXIS_MAX = 1;
6
+
7
+switch nargin
8
+    
9
+    case 1
10
+        inputStruct       = cell2mat(varargin(1));
11
+
12
+        psthStart         = inputStruct.psthStart;
13
+        psthEnd           = inputStruct.psthEnd;
14
+        nClasses          = inputStruct.nClasses;
15
+        decodePerformance = inputStruct.decodePerformance;
16
+        frameStart        = inputStruct.frameShiftStart;
17
+        frameEnd          = inputStruct.frameShiftEnd;
18
+        psth              = inputStruct.rawTimeCourse;
19
+        SubjectID         = inputStruct.SubjectID;
20
+        smoothed          = inputStruct.smoothed;
21
+
22
+    otherwise
23
+        error('spmtoolbox:SVMCrossVal:plotDecodePerformance:WrongArgument','Wrong Arguments');
24
+end
25
+
26
+    f = figure;
27
+    subplot(2,1,1);
28
+    hold on;
29
+      for voxel = 1:size(psth,2)
30
+          for label = 1:size(psth{voxel},2)
31
+              psthData = [];
32
+              for timepoint = 1:size(psth{voxel}{label},2)
33
+                  psthData = nanmean(psth{voxel}{label});
34
+              end
35
+              plot(psthStart:psthEnd,psthData,[colorChooser(voxel), lineStyleChooser(label)]);
36
+          end
37
+      end
38
+    axis([psthStart psthEnd PSTH_AXIS_MIN PSTH_AXIS_MAX])
39
+    hold off
40
+    
41
+    subplot(2,1,2)    
42
+    hold on;
43
+    
44
+    chanceLevel = 100/nClasses;
45
+    goodPredictionLevel = chanceLevel*1.5;
46
+    plot([psthStart psthEnd],[chanceLevel chanceLevel],'r');
47
+    plot([psthStart psthEnd],[goodPredictionLevel goodPredictionLevel],'g');
48
+    axis([psthStart psthEnd 0 100])
49
+    
50
+    plot(frameStart:frameEnd, mean(decodePerformance,2) ,'b');
51
+    PLOT_STD_ERR = 1;
52
+    PLOT_CLASS_PERFORMANCE = 1;
53
+    if PLOT_STD_ERR 
54
+        se = myStdErr(decodePerformance,2);
55
+        plot(frameStart:frameEnd, mean(decodePerformance,2)+se ,'b:');
56
+        plot(frameStart:frameEnd, mean(decodePerformance,2)-se ,'b:');
57
+    end
58
+    if PLOT_CLASS_PERFORMANCE
59
+        for c = 1:nClasses
60
+            plot(frameStart:frameEnd, decodePerformance() ,[colorChooser(c+2) '-']);
61
+        end
62
+    end
63
+    
64
+    
65
+    hold off;
66
+
67
+    title = sprintf('Subject %s, over %g voxel, smoothed %s',SubjectID,size(psth,2),smoothed);
68
+    set(f,'Name',title);
69
+    display(sprintf('%s',title));
70
+
71
+
72
+
73
+end
74
+
75
+
76
+function color = colorChooser(n)
77
+    switch (mod(n,8))
78
+    case 0
79
+        color = 'y';
80
+    case 1
81
+        color = 'r';
82
+    case 2
83
+        color = 'b';
84
+    case 3
85
+        color = 'g';
86
+    otherwise
87
+        color = 'k';
88
+    end
89
+end
90
+
91
+function style = lineStyleChooser(n)
92
+switch(mod(n,4))
93
+    case 0
94
+      style = '--';
95
+    case 1
96
+        style = '-';
97
+    case 2 
98
+        style = ':';
99
+    case 3
100
+        style = '-.';
101
+end
102
+end
103
+
... ...
@@ -43,14 +43,28 @@ end
43 43
     
44 44
     chanceLevel = 100/nClasses;
45 45
     goodPredictionLevel = chanceLevel*1.5;
46
-    plot([psthStart psthEnd],[chanceLevel chanceLevel],'r');
47
-    plot([psthStart psthEnd],[goodPredictionLevel goodPredictionLevel],'g');
46
+    plot([psthStart psthEnd],[chanceLevel chanceLevel],'k:');
47
+    plot([psthStart psthEnd],[goodPredictionLevel goodPredictionLevel],'k:');
48 48
     axis([psthStart psthEnd 0 100])
49 49
     
50
-    plot(frameStart:frameEnd, mean(decodePerformance,2) ,'b');
51
-    se = myStdErr(decodePerformance,2);
52
-    plot(frameStart:frameEnd, mean(decodePerformance,2)+se ,'b:');
53
-    plot(frameStart:frameEnd, mean(decodePerformance,2)-se ,'b:');
50
+    
51
+    PLOT_EXTRAS = 'class performance'
52
+    
53
+    switch PLOT_EXTRAS
54
+        case 'stderr'
55
+            plot(frameStart:frameEnd, mean(decodePerformance,2) ,'b');
56
+
57
+            se = myStdErr(decodePerformance,2);
58
+            plot(frameStart:frameEnd, mean(decodePerformance,2)+se ,'b:');
59
+            plot(frameStart:frameEnd, mean(decodePerformance,2)-se ,'b:');
60
+        case 'class performance'
61
+            for c = 1:size(decodePerformance,2)
62
+                plot(frameStart:frameEnd, decodePerformance(:,c) ,[colorChooser(mod(c,nClasses)+3) '-']);
63
+            end
64
+
65
+            plot(frameStart:frameEnd, mean(decodePerformance,2) ,'b','LineWidth',2);
66
+
67
+    end
54 68
     
55 69
     
56 70
     hold off;
... ...
@@ -67,13 +81,17 @@ end
67 81
 function color = colorChooser(n)
68 82
     switch (mod(n,8))
69 83
     case 0
70
-        color = 'y';
71
-    case 1
72 84
         color = 'r';
85
+    case 1
86
+        color = 'g';
73 87
     case 2
74 88
         color = 'b';
75 89
     case 3
76
-        color = 'g';
90
+        color = 'c';
91
+    case 4
92
+        color = 'm';
93
+    case 5
94
+        color = 'y';
77 95
     otherwise
78 96
         color = 'k';
79 97
     end
... ...
@@ -59,33 +59,38 @@ maxPerformance = -inf;
59 59
                 svmlabel  = svmlabel(rndindex);
60 60
             end
61 61
 
62
-            SVM_METHOD = 2;
62
+            SVM_METHOD = 'class performance'
63 63
             switch SVM_METHOD;
64
-                case 1
64
+                case 'libsvm crossval'
65 65
                     performance  = svmtrain(svmlabel, svmdata, svmargs);
66 66
 
67 67
                     minPerformance = min(minPerformance,performance);
68 68
                     maxPerformance = max(maxPerformance,performance);
69 69
 
70 70
                     decodePerformance = [decodePerformance; performance];
71
-                case 2
71
+                case 'class performance'
72 72
                     newsvmopt = killCrossvalOpt(svmargs);
73 73
                     
74 74
                     model = svmtrain(svmlabel,svmdata,newsvmopt);
75 75
                     classperformance = [];
76 76
                     for class = unique(svmlabel)';
77
-%                         assignin('base','uniquelabel',unique(svmlabel));
78
-%                         assignin('base','class',class);
79
-%                         assignin('base','svmlabel',svmlabel);
77
+
80 78
                         filterindex = find(class == svmlabel);
81
-                        testing_label = svmlabel(filterindex)
82
-                        testing_data  = svmdata(filterindex)
83
-                        [plabel accuracy dvalue] = svmpredict(testing_label,testing_data,model,'')
84
-%                         assignin('base','accuracy',accuracy);
79
+                        testing_label = svmlabel(filterindex);
80
+                        testing_data  = svmdata(filterindex);
81
+                        [plabel accuracy dvalue] = svmpredict(testing_label,testing_data,model,'');
82
+
85 83
                         classperformance = [classperformance accuracy(1)];
86 84
                     end
87 85
                     decodePerformance = [decodePerformance; classperformance];
88 86
                     
87
+                case 'som training'
88
+                    display('SOM TRAINING');
89
+                    addpath 'somtoolbox2';
90
+                    sD = som_data_struct(svmdata);
91
+                    assignin('base','sD',sD);
92
+                    sM = som_make(sD,'msize', [3 4],'lattice', 'rect');
93
+                    
89 94
             end
90 95
             
91 96
         end
92 97
new file mode 100644
... ...
@@ -0,0 +1,196 @@
1
+% SOM Toolbox
2
+% Version 2.0beta, May 30 2002
3
+% 
4
+% Copyright 1997-2000 by
5
+% Esa Alhoniemi, Johan Himberg, Juha Parhankangas and Juha Vesanto
6
+% Contributed files may contain copyrights of their own.
7
+% 
8
+% SOM Toolbox comes with ABSOLUTELY NO WARRANTY; for details
9
+% see License.txt in the program package. This is free software,
10
+% and you are welcome to redistribute it under certain conditions;
11
+% see License.txt for details.
12
+% 
13
+% 
14
+% Demos
15
+% 
16
+%            som_demo1   SOM Toolbox demo 1: basic properties
17
+%            som_demo2   SOM Toolbox demo 2: basic usage
18
+%            som_demo3   SOM Toolbox demo 3: visualization
19
+%            som_demo4   SOM Toolbox demo 4: data analysis
20
+% 
21
+% Creation of structs
22
+% 
23
+%              som_set   create & set (& check) values to structs
24
+%             som_info   print out information on a given struct  
25
+%      som_data_struct   create & initialize a data struct 
26
+%       som_map_struct   create & initialize a map struct 
27
+%     som_topol_struct   create & initialize a topology struct 
28
+%     som_train_struct   create & initialize a train struct 
29
+%         som_clstruct   create a cluster struct
30
+%            som_clset   set properties in a cluster struct
31
+%            som_clget   get stuff from a cluster struct
32
+% 
33
+% Struct conversion and file I/O
34
+% 
35
+%           som_vs1to2   converts a version 1.0 struct to version 2.0 struct
36
+%           som_vs2to1   converts a version 2.0 struct to version 1.0 struct
37
+%        som_read_data   reads a (SOM_PAK format) ASCII data file
38
+%       som_write_data   writes a SOM_PAK format codebook file
39
+%        som_write_cod   writes a SOM_PAK format data file
40
+%         som_read_cod   reads a SOM_PAK format codebook file
41
+% 
42
+% Data preprocessing
43
+% 
44
+%        som_normalize   normalize data set
45
+%      som_denormalize   denormalize data set 
46
+%    som_norm_variable   (de)normalize one variable
47
+%           preprocess   preprocessing GUI
48
+% 
49
+% Initialization and training functions
50
+% 
51
+%             som_make   create, initialize and train a SOM
52
+%         som_randinit   random initialization algorithm
53
+%          som_lininit   linear initialization algorithm
54
+%         som_seqtrain   sequential training algorithm
55
+%       som_batchtrain   batch training algorithm
56
+%              som_gui   SOM initialization and training GUI
57
+%       som_prototrain   a simple version of sequential training: easy to modify
58
+% 
59
+% Clustering algorithms
60
+% 
61
+%           som_kmeans   k-means algorithm (was earlier kmeans)
62
+%      kmeans_clusters   try and evaluate several k-means clusterings
63
+%           neural_gas   neural gas vector quantization algorithm
64
+%          som_linkage   hierarchical clustering algorithms
65
+%        som_cllinkage   hierarchical clustering of SOM
66
+%       som_dmatminima   local minima from distance (or U-) matrix
67
+%     som_dmatclusters   distance (or U-) matrix based clustering
68
+%         som_clspread   spreads clusters to unassinged map units
69
+%           som_cldist   calculate distances between clusters
70
+%         som_gapindex   gap validity index of clustering
71
+%             db_index   Davies-Bouldin validity index of clustering  
72
+% 
73
+% Supervised/classification algorithms
74
+% 
75
+%       som_supervised   supervised SOM algorithm
76
+%                 lvq1   LVQ1 algorithm
77
+%                 lvq3   LVQ3 algorithm
78
+%                  knn   k-NN classification algorithm 
79
+%              knn_old   k-NN classification algorithm (old version)
80
+% 
81
+% SOM error measures
82
+% 
83
+%          som_quality   quantization and topographic error of SOM
84
+%       som_distortion   SOM distortion measure
85
+%      som_distortion3   elements of the SOM distortion measure
86
+% 
87
+% Auxiliary functions
88
+% 
89
+%             som_bmus   calculates BMUs for given data vectors
90
+%         som_eucdist2   pairwise squared euclidian distances between vectors
91
+%            som_mdist   calculates pairwise distances between vectors 
92
+%           som_divide   extract subsets of data based on map
93
+%            som_label   give labels to map units
94
+%        som_label2num   rcodes string data labels to interger class labels 
95
+%        som_autolabel   automatically labels the SOM based on given data
96
+%      som_unit_coords   calculates coordinates in output space for map units
97
+%       som_unit_dists   distances in output space between map units
98
+%      som_unit_neighs   units in 1-neighborhood for each map unit
99
+%     som_neighborhood   calculates neighborhood matrix for the given map
100
+%        som_neighbors   calculates different kinds of neighborhoods 
101
+%           som_neighf   calculates neighborhood function values
102
+%           som_select   GUI for manual selection of map units
103
+%     som_estimate_gmm   create Gaussian mixture model on top of SOM
104
+%  som_probability_gmm   evaluate Gaussian mixture model
105
+%          som_ind2sub   from linear index to subscript index 
106
+%          som_sub2ind   from subscript index to linear index
107
+%          som_ind2cod   from linear index to SOM_PAK linear index 
108
+%          som_cod2ind   from SOM_linear index to SOM_PAK linear index 
109
+%             nanstats   mean, std and median which ignore NaNs
110
+%   som_modify_dataset   add, remove, or extract samples and components
111
+%         som_fillnans   fill NaNs in a data set based on given SOM
112
+%            som_stats   statistics of a data set
113
+%           som_drmake   calculate descriptive rules for a cluster
114
+%           som_dreval   evaluate descriptive rules for a cluster
115
+%         som_drsignif   rule significance measures
116
+% 
117
+% Using SOM_PAK from Matlab
118
+% 
119
+%      som_sompaktrain   uses SOM_PAK to train a map
120
+%           sompak_gui   GUI for using SOM_PAK from Matlab
121
+%          sompak_init   call SOM_PAK's initialization programs from Matlab
122
+%      sompak_init_gui   GUI for using SOM_PAK's initialization from Matlab
123
+%    sompak_rb_control   an auxiliary function for sompak_*_gui functions.
124
+%        sompak_sammon   call SOM_PAK's Sammon program from Matlab
125
+%    sompak_sammon_gui   GUI for using SOM_PAK's Sammon program from Matlab
126
+%         sompak_train   call SOM_PAK's training program from Matlab
127
+%     sompak_train_gui   GUI for using SOM_PAK's training program from Matlab 
128
+% 
129
+% Visualization
130
+% 
131
+%             som_show   basic visualization
132
+%         som_show_add   add labels, hits and trajectories
133
+%       som_show_clear   remove extra markers
134
+%       som_recolorbar   refresh/reconfigure colorbars
135
+%         som_show_gui   GUI for using som_show and associated functions
136
+%             som_grid   visualization of SOM grid
137
+%           som_cplane   component planes and U-matrices
138
+%         som_barplane   bar chart visualization of map
139
+%         som_pieplane   pie chart visualization of map
140
+%        som_plotplane   plot chart visualization of map
141
+%       som_trajectory   launches a GUI for presenting comet-trajectories 
142
+%       som_dendrogram   visualization of clustering tree
143
+%       som_plotmatrix   pairwise scatter plots and histograms
144
+%    som_order_cplanes   order and visualize the component planes
145
+%           som_clplot   plots of clusters (based on cluster struct)
146
+% som_projections_plot   projections plots (see som_projections)
147
+%       som_stats_plot   plots of statistics (see som_stats)
148
+% 
149
+% Auxiliary functions for visualization
150
+% 
151
+%                 hits   calculates hits, or sum of values for each map unit
152
+%             som_hits   calculates the response of data on the map
153
+%             som_umat   calculates the U-matrix
154
+%                  cca   curvilinear component analysis projection algorithm
155
+%              pcaproj   principal component projection algorithm
156
+%               sammon   Sammon's mapping projection algorithm
157
+%       som_connection   connection matrix for map 
158
+%       som_vis_coords   map unit coordinates used in visualizations
159
+%        som_colorcode   create color coding for map/2D data
160
+%         som_bmucolor   colors of the BMUs from a given map color code
161
+%        som_normcolor   simulate indexed colormap
162
+%     som_clustercolor   color coding which depends on clustering structure
163
+%      som_kmeanscolor   color coding according to k-means clustering
164
+%     som_kmeanscolor2   a newer version of the som_kmeanscolor function
165
+%       som_fuzzycolor   a fuzzy color coding 
166
+%         som_coloring   a SOM-based color coding 
167
+%      som_projections   calculates a default set of projections
168
+% 
169
+% Report generation stuff
170
+% 
171
+%     som_table_struct   creates a table struct
172
+%     som_table_modify   modifies a table struct
173
+%      som_table_print   print a table in various formats
174
+%            rep_utils   various utilities for printing report elements
175
+%      som_stats_table   a table of data set statistics
176
+%     som_stats_report   report on data set statistics
177
+% 
178
+% Low level routines used by visualization functions
179
+% 
180
+%            vis_patch   defines hexagonal and rectangular patches
181
+%    vis_som_show_data   returns UserData and subplot handles stored by som_show.m
182
+%        vis_valuetype   used for type checks 
183
+%         vis_footnote   adds a movable text to the current figure 
184
+%          vis_trajgui   the actual GUI started by som_trajectory.m 
185
+% vis_PlaneAxisProperties   set axis properties in visualization functions
186
+% vis_footnoteButtonDownFcn   callback function for vis_footnote.m
187
+%     vis_planeGetArgs   converts topol struct to lattice, msize argument pair
188
+%    vis_show_gui_comp   internal function used by som_show_gui.m
189
+%    vis_show_gui_tool   internal function used by som_show_gui.m
190
+% 
191
+% Other
192
+% 
193
+%           somtoolbox   this file
194
+%            iris.data   IRIS data set (used in demos)
195
+%          License.txt   GNU General Public License 
196
+%        Copyright.txt   Copyright notice
0 197
new file mode 100644
... ...
@@ -0,0 +1,55 @@
1
+COPYRIGHT
2
+---------
3
+
4
+SOM Toolbox 2.0, a software library for Matlab 5 implementing the
5
+Self-Organizing Map algorithm is Copyright (C) 1999 by Esa Alhoniemi,
6
+Johan Himberg, Jukka Parviainen and Juha Vesanto.
7
+
8
+This package is free software; you can redistribute it and/or
9
+modify it under the terms of the GNU General Public License
10
+as published by the Free Software Foundation; either version 2
11
+of the License, or any later version.
12
+
13
+Note: only part of the files distributed in the package belong to the
14
+SOM Toolbox. The package also contains contributed files, which may
15
+have their own copyright notices. If not, the GNU General Public
16
+License holds for them, too, but so that the author(s) of the file
17
+have the Copyright.
18
+
19
+This package is distributed in the hope that it will be useful,
20
+but WITHOUT ANY WARRANTY; without even the implied warranty of
21
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
22
+GNU General Public License for more details.
23
+
24
+You should have received a copy of the GNU General Public License
25
+along with this package (file License.txt); if not, write to the Free
26
+Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
27
+02111-1307, USA.
28
+
29
+
30
+LICENSING THE LIBRARY FOR PROPRIETARY PROGRAMS
31
+----------------------------------------------
32
+
33
+As stated in the GNU General Public License (see License.txt) it is
34
+not possible to include this software library in a commercial
35
+proprietary program without written permission from the owners of the
36
+copyright. If you wish to obtain such permission, you can reach us by
37
+paper mail:
38
+
39
+  SOM Toolbox team
40
+  Laboratory of Computer and Information Science
41
+  P.O.Box 5400
42
+  FIN-02015 HUT 
43
+  Finland
44
+  Europe 
45
+
46
+and by email:
47
+
48
+  somtlbx@mail.cis.hut.fi
49
+
50
+
51
+
52
+
53
+
54
+
55
+
0 56
new file mode 100644
... ...
@@ -0,0 +1,386 @@
1
+GNU GENERAL PUBLIC LICENSE
2
+
3
+Version 2, June 1991
4
+
5
+Copyright (C) 1989, 1991 Free Software Foundation, Inc.  
6
+59 Temple Place - Suite 330, Boston, MA  02111-1307, USA
7
+
8
+Everyone is permitted to copy and distribute verbatim copies
9
+of this license document, but changing it is not allowed.
10
+
11
+------------------------------------------------------------------------
12
+
13
+Preamble
14
+
15
+  The licenses for most software are designed to take away your
16
+freedom to share and change it.  By contrast, the GNU General Public
17
+License is intended to guarantee your freedom to share and change free
18
+software--to make sure the software is free for all its users.  This
19
+General Public License applies to most of the Free Software
20
+Foundation's software and to any other program whose authors commit to
21
+using it.  (Some other Free Software Foundation software is covered by
22
+the GNU Library General Public License instead.)  You can apply it to
23
+your programs, too.
24
+
25
+  When we speak of free software, we are referring to freedom, not
26
+price.  Our General Public Licenses are designed to make sure that you
27
+have the freedom to distribute copies of free software (and charge for
28
+this service if you wish), that you receive source code or can get it
29
+if you want it, that you can change the software or use pieces of it
30
+in new free programs; and that you know you can do these things.
31
+
32
+  To protect your rights, we need to make restrictions that forbid
33
+anyone to deny you these rights or to ask you to surrender the rights.
34
+These restrictions translate to certain responsibilities for you if you
35
+distribute copies of the software, or if you modify it.
36
+
37
+  For example, if you distribute copies of such a program, whether
38
+gratis or for a fee, you must give the recipients all the rights that
39
+you have.  You must make sure that they, too, receive or can get the
40
+source code.  And you must show them these terms so they know their
41
+rights.
42
+
43
+  We protect your rights with two steps: (1) copyright the software, and
44
+(2) offer you this license which gives you legal permission to copy,
45
+distribute and/or modify the software.
46
+
47
+  Also, for each author's protection and ours, we want to make certain
48
+that everyone understands that there is no warranty for this free
49
+software.  If the software is modified by someone else and passed on, we
50
+want its recipients to know that what they have is not the original, so
51
+that any problems introduced by others will not reflect on the original
52
+authors' reputations.
53
+
54
+  Finally, any free program is threatened constantly by software
55
+patents.  We wish to avoid the danger that redistributors of a free
56
+program will individually obtain patent licenses, in effect making the
57
+program proprietary.  To prevent this, we have made it clear that any
58
+patent must be licensed for everyone's free use or not licensed at all.
59
+
60
+  The precise terms and conditions for copying, distribution and
61
+modification follow.
62
+
63
+---------------------------------------------------------------------
64
+
65
+TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION
66
+
67
+0.
68
+
69
+This License applies to any program or other work which contains a
70
+notice placed by the copyright holder saying it may be distributed
71
+under the terms of this General Public License.  The "Program", below,
72
+refers to any such program or work, and a "work based on the Program"
73
+means either the Program or any derivative work under copyright law:
74
+that is to say, a work containing the Program or a portion of it,
75
+either verbatim or with modifications and/or translated into another
76
+language.  (Hereinafter, translation is included without limitation in
77
+the term "modification".)  Each licensee is addressed as "you".
78
+
79
+Activities other than copying, distribution and modification are not
80
+covered by this License; they are outside its scope.  The act of
81
+running the Program is not restricted, and the output from the Program
82
+is covered only if its contents constitute a work based on the
83
+Program (independent of having been made by running the Program).
84
+Whether that is true depends on what the Program does.
85
+
86
+1.
87
+
88
+You may copy and distribute verbatim copies of the Program's source
89
+code as you receive it, in any medium, provided that you conspicuously
90
+and appropriately publish on each copy an appropriate copyright notice
91
+and disclaimer of warranty; keep intact all the notices that refer to
92
+this License and to the absence of any warranty; and give any other
93
+recipients of the Program a copy of this License along with the
94
+Program.
95
+
96
+You may charge a fee for the physical act of transferring a copy, and
97
+you may at your option offer warranty protection in exchange for a fee.
98
+
99
+2.
100
+
101
+You may modify your copy or copies of the Program or any portion of
102
+it, thus forming a work based on the Program, and copy and distribute
103
+such modifications or work under the terms of Section 1 above,
104
+provided that you also meet all of these conditions:
105
+
106
+a) You must cause the modified files to carry prominent notices
107
+   stating that you changed the files and the date of any change.
108
+
109
+b) You must cause any work that you distribute or publish, that in
110
+   whole or in part contains or is derived from the Program or any part
111
+   thereof, to be licensed as a whole at no charge to all third parties
112
+   under the terms of this License.
113
+
114
+c) If the modified program normally reads commands interactively when
115
+   run, you must cause it, when started running for such interactive use
116
+   in the most ordinary way, to print or display an announcement
117
+   including an appropriate copyright notice and a notice that there is
118
+   no warranty (or else, saying that you provide a warranty) and that
119
+   users may redistribute the program under these conditions, and telling
120
+   the user how to view a copy of this License.  (Exception: if the
121
+   Program itself is interactive but does not normally print such an
122
+   announcement, your work based on the Program is not required to print
123
+   an announcement.)
124
+
125
+These requirements apply to the modified work as a whole.  If
126
+identifiable sections of that work are not derived from the Program,
127
+and can be reasonably considered independent and separate works in
128
+themselves, then this License, and its terms, do not apply to those
129
+sections when you distribute them as separate works.  But when you
130
+distribute the same sections as part of a whole which is a work based
131
+on the Program, the distribution of the whole must be on the terms of
132
+this License, whose permissions for other licensees extend to the
133
+entire whole, and thus to each and every part regardless of who wrote
134
+it.
135
+
136
+Thus, it is not the intent of this section to claim rights or contest
137
+your rights to work written entirely by you; rather, the intent is to
138
+exercise the right to control the distribution of derivative or
139
+collective works based on the Program.
140
+
141
+In addition, mere aggregation of another work not based on the Program
142
+with the Program (or with a work based on the Program) on a volume of
143
+a storage or distribution medium does not bring the other work under
144
+the scope of this License.
145
+
146
+3.
147
+
148
+You may copy and distribute the Program (or a work based on it, under
149
+Section 2) in object code or executable form under the terms of
150
+Sections 1 and 2 above provided that you also do one of the following:
151
+
152
+a) Accompany it with the complete corresponding machine-readable
153
+   source code, which must be distributed under the terms of Sections 1
154
+   and 2 above on a medium customarily used for software interchange; or,
155
+
156
+b) Accompany it with a written offer, valid for at least three years,
157
+   to give any third party, for a charge no more than your cost of
158
+   physically performing source distribution, a complete machine-readable
159
+   copy of the corresponding source code, to be distributed under the
160
+   terms of Sections 1 and 2 above on a medium customarily used for
161
+   software interchange; or,
162
+
163
+c) Accompany it with the information you received as to the offer to
164
+   distribute corresponding source code.  (This alternative is allowed
165
+   only for noncommercial distribution and only if you received the
166
+   program in object code or executable form with such an offer, in
167
+   accord with Subsection b above.)
168
+
169
+The source code for a work means the preferred form of the work for
170
+making modifications to it.  For an executable work, complete source
171
+code means all the source code for all modules it contains, plus any
172
+associated interface definition files, plus the scripts used to
173
+control compilation and installation of the executable.  However, as a
174
+special exception, the source code distributed need not include
175
+anything that is normally distributed (in either source or binary
176
+form) with the major components (compiler, kernel, and so on) of the
177
+operating system on which the executable runs, unless that component
178
+itself accompanies the executable.
179
+
180
+If distribution of executable or object code is made by offering
181
+access to copy from a designated place, then offering equivalent
182
+access to copy the source code from the same place counts as
183
+distribution of the source code, even though third parties are not
184
+compelled to copy the source along with the object code.
185
+
186
+4.
187
+
188
+You may not copy, modify, sublicense, or distribute the Program except
189
+as expressly provided under this License.  Any attempt otherwise to
190
+copy, modify, sublicense or distribute the Program is void, and will
191
+automatically terminate your rights under this License.  However,
192
+parties who have received copies, or rights, from you under this
193
+License will not have their licenses terminated so long as such
194
+parties remain in full compliance.
195
+
196
+5.
197
+
198
+You are not required to accept this License, since you have not signed
199
+it.  However, nothing else grants you permission to modify or
200
+distribute the Program or its derivative works.  These actions are
201
+prohibited by law if you do not accept this License.  Therefore, by
202
+modifying or distributing the Program (or any work based on the
203
+Program), you indicate your acceptance of this License to do so, and
204
+all its terms and conditions for copying, distributing or modifying
205
+the Program or works based on it.
206
+
207
+6.
208
+
209
+Each time you redistribute the Program (or any work based on the
210
+Program), the recipient automatically receives a license from the
211
+original licensor to copy, distribute or modify the Program subject to
212
+these terms and conditions.  You may not impose any further
213
+restrictions on the recipients' exercise of the rights granted herein.
214
+You are not responsible for enforcing compliance by third parties to
215
+this License.
216
+
217
+7.
218
+
219
+If, as a consequence of a court judgment or allegation of patent
220
+infringement or for any other reason (not limited to patent issues),
221
+conditions are imposed on you (whether by court order, agreement or
222
+otherwise) that contradict the conditions of this License, they do not
223
+excuse you from the conditions of this License.  If you cannot
224
+distribute so as to satisfy simultaneously your obligations under this
225
+License and any other pertinent obligations, then as a consequence you
226
+may not distribute the Program at all.  For example, if a patent
227
+license would not permit royalty-free redistribution of the Program by
228
+all those who receive copies directly or indirectly through you, then
229
+the only way you could satisfy both it and this License would be to
230
+refrain entirely from distribution of the Program.
231
+
232
+If any portion of this section is held invalid or unenforceable under
233
+any particular circumstance, the balance of the section is intended to
234
+apply and the section as a whole is intended to apply in other
235
+circumstances.
236
+
237
+It is not the purpose of this section to induce you to infringe any
238
+patents or other property right claims or to contest validity of any
239
+such claims; this section has the sole purpose of protecting the
240
+integrity of the free software distribution system, which is
241
+implemented by public license practices.  Many people have made
242
+generous contributions to the wide range of software distributed
243
+through that system in reliance on consistent application of that
244
+system; it is up to the author/donor to decide if he or she is willing
245
+to distribute software through any other system and a licensee cannot
246
+impose that choice.
247
+
248
+This section is intended to make thoroughly clear what is believed to
249
+be a consequence of the rest of this License.
250
+
251
+8.
252
+
253
+If the distribution and/or use of the Program is restricted in certain
254
+countries either by patents or by copyrighted interfaces, the original
255
+copyright holder who places the Program under this License may add an
256
+explicit geographical distribution limitation excluding those
257
+countries, so that distribution is permitted only in or among
258
+countries not thus excluded.  In such case, this License incorporates
259
+the limitation as if written in the body of this License.
260
+
261
+9.
262
+
263
+The Free Software Foundation may publish revised and/or new versions
264
+of the General Public License from time to time.  Such new versions
265
+will be similar in spirit to the present version, but may differ in
266
+detail to address new problems or concerns.
267
+
268
+Each version is given a distinguishing version number.  If the Program
269
+specifies a version number of this License which applies to it and
270
+"any later version", you have the option of following the terms and
271
+conditions either of that version or of any later version published by
272
+the Free Software Foundation.  If the Program does not specify a
273
+version number of this License, you may choose any version ever
274
+published by the Free Software Foundation.
275
+
276
+10.
277
+
278
+If you wish to incorporate parts of the Program into other free
279
+programs whose distribution conditions are different, write to the
280
+author to ask for permission.  For software which is copyrighted by
281
+the Free Software Foundation, write to the Free Software Foundation;
282
+we sometimes make exceptions for this.  Our decision will be guided by
283
+the two goals of preserving the free status of all derivatives of our
284
+free software and of promoting the sharing and reuse of software
285
+generally.
286
+
287
+NO WARRANTY
288
+
289
+11.
290
+ 
291
+BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY
292
+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW.  EXCEPT
293
+WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER
294
+PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND,
295
+EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE
296
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
297
+PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE
298
+PROGRAM IS WITH YOU.  SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME
299
+THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
300
+
301
+12.
302
+
303
+IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
304
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR
305
+REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR
306
+DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL
307
+DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE PROGRAM
308
+(INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING RENDERED
309
+INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A FAILURE OF
310
+THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS), EVEN IF SUCH HOLDER
311
+OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH DAMAGES.
312
+
313
+END OF TERMS AND CONDITIONS
314
+
315
+-----------------------------------------------------------------------
316
+
317
+How to Apply These Terms to Your New Programs 
318
+
319
+If you develop a new program, and you want it to be of the greatest
320
+possible use to the public, the best way to achieve this is to make it
321
+free software which everyone can redistribute and change under these
322
+terms.
323
+
324
+To do so, attach the following notices to the program.  It is safest
325
+to attach them to the start of each source file to most effectively
326
+convey the exclusion of warranty; and each file should have at least
327
+the "copyright" line and a pointer to where the full notice is found.
328
+
329
+one line to give the program's name and an idea of what it does.
330
+Copyright (C) 20<yy>  <name of author>
331
+
332
+This program is free software; you can redistribute it and/or modify
333
+it under the terms of the GNU General Public License as published by
334
+the Free Software Foundation; either version 2 of the License, or (at
335
+your option) any later version.
336
+
337
+This program is distributed in the hope that it will be useful, but
338
+WITHOUT ANY WARRANTY; without even the implied warranty of
339
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
340
+General Public License for more details.
341
+
342
+You should have received a copy of the GNU General Public License
343
+along with this program; if not, write to the Free Software
344
+Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
345
+USA. 
346
+
347
+Also add information on how to contact you by electronic and paper mail.
348
+
349
+If the program is interactive, make it output a short notice like this
350
+when it starts in an interactive mode:
351
+
352
+ Gnomovision version 69, Copyright (C) 20<yy> <name of author> 
353
+ Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show
354
+ w'.  This is free software, and you are welcome to redistribute it
355
+ under certain conditions; type `show c' for details.  
356
+
357
+The hypothetical commands `show w' and `show c' should show the
358
+appropriate parts of the General Public License.  Of course, the
359
+commands you use may be called something other than `show w' and `show
360
+c'; they could even be mouse-clicks or menu items--whatever suits your
361
+program.
362
+
363
+You should also get your employer (if you work as a programmer) or
364
+your school, if any, to sign a "copyright disclaimer" for the program,
365
+if necessary.  Here is a sample; alter the names:
366
+
367
+ Yoyodyne, Inc., hereby disclaims all copyright
368
+ interest in the program `Gnomovision'
369
+ (which makes passes at compilers) written 
370
+ by James Hacker.
371
+
372
+ <signature of Ty Coon>, 1 April 1989
373
+ Ty Coon, President of Vice
374
+
375
+This General Public License does not permit incorporating your program
376
+into proprietary programs.  If your program is a subroutine library,
377
+you may consider it more useful to permit linking proprietary
378
+applications with the library.  If this is what you want to do, use
379
+the GNU Library General Public License instead of this License.
380
+
381
+FSF & GNU inquiries & questions to gnu@prep.ai.mit.edu.
382
+
383
+Copyright notice above.
384
+Free Software Foundation, Inc.,
385
+59 Temple Place - Suite 330, Boston, MA  02111,  USA
386
+
0 387
new file mode 100644
... ...
@@ -0,0 +1,282 @@
1
+function [P] = cca(D, P, epochs, Mdist, alpha0, lambda0)
2
+
3
+%CCA Projects data vectors using Curvilinear Component Analysis.
4
+%
5
+% P = cca(D, P, epochs, [Dist], [alpha0], [lambda0])
6
+%
7
+%  P = cca(D,2,10);           % projects the given data to a plane
8
+%  P = cca(D,pcaproj(D,2),5); % same, but with PCA initialization
9
+%  P = cca(D, 2, 10, Dist);   % same, but the given distance matrix is used
10
+%  
11
+%  Input and output arguments ([]'s are optional):
12
+%   D          (matrix) the data matrix, size dlen x dim
13
+%              (struct) data or map struct            
14
+%   P          (scalar) output dimension
15
+%              (matrix) size dlen x odim, the initial projection
16
+%   epochs     (scalar) training length
17
+%   [Dist]     (matrix) pairwise distance matrix, size dlen x dlen.
18
+%                       If the distances in the input space should
19
+%                       be calculated otherwise than as euclidian
20
+%                       distances, the distance from each vector
21
+%                       to each other vector can be given here,
22
+%                       size dlen x dlen. For example PDIST
23
+%                       function can be used to calculate the
24
+%                       distances: Dist = squareform(pdist(D,'mahal'));
25
+%   [alpha0]   (scalar) initial step size, 0.5 by default
26
+%   [lambda0]  (scalar) initial radius of influence, 3*max(std(D)) by default
27
+%  
28
+%   P          (matrix) size dlen x odim, the projections
29
+%
30
+% Unknown values (NaN's) in the data: projections of vectors with
31
+% unknown components tend to drift towards the center of the
32
+% projection distribution. Projections of totally unknown vectors are
33
+% set to unknown (NaN).
34
+%
35
+% See also SAMMON, PCAPROJ. 
36
+
37
+% Reference: Demartines, P., Herault, J., "Curvilinear Component
38
+%   Analysis: a Self-Organizing Neural Network for Nonlinear
39
+%   Mapping of Data Sets", IEEE Transactions on Neural Networks, 
40
+%   vol 8, no 1, 1997, pp. 148-154.
41
+
42
+% Contributed to SOM Toolbox 2.0, February 2nd, 2000 by Juha Vesanto
43
+% Copyright (c) by Juha Vesanto
44
+% http://www.cis.hut.fi/projects/somtoolbox/
45
+
46
+% juuso 171297 040100
47
+
48
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49
+%% Check arguments 
50
+
51
+error(nargchk(3, 6, nargin)); % check the number of input arguments
52
+
53
+% input data
54
+if isstruct(D), 
55
+  if strcmp(D.type,'som_map'), D = D.codebook; else D = D.data; end
56
+end
57
+[noc dim] = size(D);
58
+noc_x_1  = ones(noc, 1); % used frequently
59
+me = zeros(1,dim); st = zeros(1,dim);
60
+for i=1:dim,
61
+  me(i) = mean(D(find(isfinite(D(:,i))),i));
62
+  st(i) = std(D(find(isfinite(D(:,i))),i));
63
+end
64
+
65
+% initial projection
66
+if prod(size(P))==1, 
67
+  P = (2*rand(noc,P)-1).*st(noc_x_1,1:P) + me(noc_x_1,1:P); 
68
+else
69
+  % replace unknown projections with known values
70
+  inds = find(isnan(P)); P(inds) = rand(size(inds));
71
+end
72
+[dummy odim] = size(P);
73
+odim_x_1  = ones(odim, 1); % this is used frequently
74
+
75
+% training length
76
+train_len = epochs*noc;
77
+
78
+% random sample order
79
+rand('state',sum(100*clock));
80
+sample_inds = ceil(noc*rand(train_len,1));
81
+
82
+% mutual distances
83
+if nargin<4 | isempty(Mdist) | all(isnan(Mdist(:))),
84
+  fprintf(2, 'computing mutual distances\r');
85
+  dim_x_1 = ones(dim,1);
86
+  for i = 1:noc,
87
+    x = D(i,:); 
88
+    Diff = D - x(noc_x_1,:);
89
+    N = isnan(Diff);
90
+    Diff(find(N)) = 0; 
91
+    Mdist(:,i) = sqrt((Diff.^2)*dim_x_1);
92
+    N = find(sum(N')==dim); %mutual distance unknown
93
+    if ~isempty(N), Mdist(N,i) = NaN; end
94
+  end
95
+else
96
+  % if the distance matrix is output from PDIST function
97
+  if size(Mdist,1)==1, Mdist = squareform(Mdist); end
98
+  if size(Mdist,1)~=noc, 
99
+    error('Mutual distance matrix size and data set size do not match'); 
100
+  end
101
+end
102
+
103
+% alpha and lambda
104
+if nargin<5 | isempty(alpha0) | isnan(alpha0), alpha0 = 0.5; end
105
+alpha = potency_curve(alpha0,alpha0/100,train_len);
106
+
107
+if nargin<6 | isempty(lambda0) | isnan(lambda0), lambda0 = max(st)*3; end
108
+lambda = potency_curve(lambda0,0.01,train_len);
109
+
110
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
111
+%% Action
112
+
113
+k=0; fprintf(2, 'iterating: %d / %d epochs\r',k,epochs);
114
+
115
+for i=1:train_len, 
116
+  
117
+  ind = sample_inds(i);     % sample index
118
+  dx = Mdist(:,ind);        % mutual distances in input space
119
+  known = find(~isnan(dx)); % known distances
120
+
121
+  if ~isempty(known),
122
+    % sample vector's projection
123
+    y = P(ind,:);                 
124
+
125
+    % distances in output space
126
+    Dy = P(known,:) - y(noc_x_1(known),:); 
127
+    dy = sqrt((Dy.^2)*odim_x_1);           
128
+  
129
+    % relative effect
130
+    dy(find(dy==0)) = 1;        % to get rid of div-by-zero's
131
+    fy = exp(-dy/lambda(i)) .* (dx(known) ./ dy - 1);
132
+
133
+    % Note that the function F here is e^(-dy/lambda)) 
134
+    % instead of the bubble function 1(lambda-dy) used in the 
135
+    % paper.
136
+    
137
+    % Note that here a simplification has been made: the derivatives of the
138
+    % F function have been ignored in calculating the gradient of error
139
+    % function w.r.t. to changes in dy.
140
+    
141
+    % update
142
+    P(known,:) = P(known,:) + alpha(i)*fy(:,odim_x_1).*Dy;
143
+  end
144
+
145
+  % track
146
+  if rem(i,noc)==0, 
147
+    k=k+1; fprintf(2, 'iterating: %d / %d epochs\r',k,epochs);
148
+  end
149
+
150
+end
151
+
152
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
153
+%% clear up
154
+
155
+% calculate error
156
+error = cca_error(P,Mdist,lambda(train_len));
157
+fprintf(2,'%d iterations, error %f          \n', epochs, error);
158
+
159
+% set projections of totally unknown vectors as unknown
160
+unknown = find(sum(isnan(D)')==dim);
161
+P(unknown,:) = NaN;
162
+
163
+return;
164
+
165
+
166
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
167
+%% tips
168
+
169
+% to plot the results, use the code below
170
+
171
+%subplot(2,1,1), 
172
+%switch(odim), 
173
+%  case 1, plot(P(:,1),ones(dlen,1),'x')
174
+%  case 2, plot(P(:,1),P(:,2),'x'); 
175
+%  otherwise, plot3(P(:,1),P(:,2),P(:,3),'x'); rotate3d on
176
+%end
177
+%subplot(2,1,2), dydxplot(P,Mdist);
178
+
179
+% to a project a new point x in the input space to the output space
180
+% do the following:
181
+
182
+% Diff = D - x(noc_x_1,:); Diff(find(isnan(Diff))) = 0; 
183
+% dx = sqrt((Diff.^2)*dim_x_1);
184
+% p = project_point(P,x,dx); % this function can be found from below
185
+% tlen = size(p,1);
186
+% plot(P(:,1),P(:,2),'bx',p(tlen,1),p(tlen,2),'ro',p(:,1),p(:,2),'r-')
187
+
188
+% similar trick can be made to the other direction
189
+
190
+
191
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
192
+%% subfunctions
193
+
194
+function vals = potency_curve(v0,vn,l)
195
+
196
+  % curve that decreases from v0 to vn with a rate that is 
197
+  % somewhere between linear and 1/t
198
+  vals = v0 * (vn/v0).^([0:(l-1)]/(l-1));
199
+
200
+
201
+function error = cca_error(P,Mdist,lambda)
202
+
203
+  [noc odim] = size(P);
204
+  noc_x_1 = ones(noc,1);
205
+  odim_x_1 = ones(odim,1);
206
+
207
+  error = 0;
208
+  for i=1:noc,
209
+    known = find(~isnan(Mdist(:,i)));
210
+    if ~isempty(known),   
211
+      y = P(i,:);                 
212
+      Dy = P(known,:) - y(noc_x_1(known),:);
213
+      dy = sqrt((Dy.^2)*odim_x_1);
214
+      fy = exp(-dy/lambda);
215
+      error = error + sum(((Mdist(known,i) - dy).^2).*fy);
216
+    end
217
+  end
218
+  error = error/2;
219
+
220
+
221
+function [] = dydxplot(P,Mdist)
222
+
223
+  [noc odim] = size(P);
224
+  noc_x_1 = ones(noc,1);
225
+  odim_x_1 = ones(odim,1);
226
+  Pdist = zeros(noc,noc);
227
+    
228
+  for i=1:noc,
229
+    y = P(i,:);                 
230
+    Dy = P - y(noc_x_1,:);
231
+    Pdist(:,i) = sqrt((Dy.^2)*odim_x_1);
232
+  end
233
+
234
+  Pdist = tril(Pdist,-1); 
235
+  inds = find(Pdist > 0); 
236
+  n = length(inds);
237
+  plot(Pdist(inds),Mdist(inds),'.');
238
+  xlabel('dy'), ylabel('dx')
239
+
240
+
241
+function p = project_point(P,x,dx)
242
+
243
+  [noc odim] = size(P);
244
+  noc_x_1 = ones(noc,1);
245
+  odim_x_1 = ones(odim,1);
246
+
247
+  % initial projection
248
+  [dummy,i] = min(dx);
249
+  y = P(i,:)+rand(1,odim)*norm(P(i,:))/20;
250
+ 
251
+  % lambda 
252
+  lambda = norm(std(P));
253
+
254
+  % termination
255
+  eps = 1e-3; i_max = noc*10;
256
+  
257
+  i=1; p(i,:) = y; 
258
+  ready = 0;
259
+  while ~ready,
260
+
261
+    % mutual distances
262
+    Dy = P - y(noc_x_1,:);        % differences in output space
263
+    dy = sqrt((Dy.^2)*odim_x_1);  % distances in output space
264
+    f = exp(-dy/lambda);
265
+  
266
+    fprintf(2,'iteration %d, error %g \r',i,sum(((dx - dy).^2).*f));
267
+
268
+    % all the other vectors push the projected one
269
+    fy = f .* (dx ./ dy - 1) / sum(f);
270
+  
271
+    % update    
272
+    step = - sum(fy(:,odim_x_1).*Dy);
273
+    y = y + step;
274
+  
275
+    i=i+1;
276
+    p(i,:) = y;   
277
+    ready = (norm(step)/norm(y) < eps | i > i_max);
278
+
279
+  end
280
+  fprintf(2,'\n');
281
+     
282
+  
0 283
\ No newline at end of file
1 284
new file mode 100644
... ...
@@ -0,0 +1,83 @@
1
+function [t,r] = db_index(D, cl, C, p, q)
2
+ 
3
+% DB_INDEX Davies-Bouldin clustering evaluation index.
4
+%
5
+% [t,r] = db_index(D, cl, C, p, q)
6
+%
7
+%  Input and output arguments ([]'s are optional):  
8
+%    D     (matrix) data (n x dim)
9
+%          (struct) map or data struct
10
+%    cl    (vector) cluster numbers corresponding to data samples (n x 1)
11
+%    [C]   (matrix) prototype vectors (c x dim) (default = cluster means)
12
+%    [p]   (scalar) norm used in the computation (default == 2)
13
+%    [q]   (scalar) moment used to calculate cluster dispersions (default = 2)
14
+% 
15
+%    t     (scalar) Davies-Bouldin index for the clustering (=mean(r))
16
+%    r     (vector) maximum DB index for each cluster (size c x 1)    
17
+% 
18
+% See also  KMEANS, KMEANS_CLUSTERS, SOM_GAPINDEX.
19
+
20
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
21
+%% input arguments
22
+
23
+if isstruct(D), 
24
+    switch D.type,
25
+    case 'som_map', D = D.codebook; 
26
+    case 'som_data', D = D.data; 
27
+    end
28
+end
29
+
30
+% cluster centroids
31
+[l dim] = size(D);
32
+u = unique(cl); 
33
+c = length(u); 
34
+if nargin <3, 
35
+  C = zeros(c,dim); 
36
+  for i=1:c, 
37
+      me = nanstats(D(find(cl==u(i)),:));
38
+      C(i,:) = me';
39
+  end 
40
+end
41
+
42
+u2i = zeros(max(u),1); u2i(u) = 1:c; 
43
+D = som_fillnans(D,C,u2i(cl)); % replace NaN's with cluster centroid values
44
+
45
+if nargin <4, p = 2; end % euclidian distance between cluster centers
46
+if nargin <5, q = 2; end % dispersion = standard deviation
47
+ 
48
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49
+%% action
50
+
51
+% dispersion in each cluster 
52
+for i = 1:c
53
+  ind = find(cl==u(i)); % points in this cluster
54
+  l   = length(ind);
55
+  if l > 0
56
+    S(i) = (mean(sqrt(sum((D(ind,:) - ones(l,1) * C(i,:)).^2,2)).^q))^(1/q);
57
+  else
58
+    S(i) = NaN;
59
+  end
60
+end
61
+ 
62
+% distances between clusters
63
+%for i = 1:c
64
+%  for j = i+1:c
65
+%    M(i,j) = sum(abs(C(i,:) - C(j,:)).^p)^(1/p);
66
+%  end
67
+%end
68
+M = som_mdist(C,p); 
69
+
70
+% Davies-Bouldin index
71
+R = NaN * zeros(c);
72
+r = NaN * zeros(c,1);
73
+for i = 1:c
74
+  for j = i+1:c
75
+    R(i,j) = (S(i) + S(j))/M(i,j);
76
+  end
77
+  r(i) = max(R(i,:));
78
+end
79
+ 
80
+t = mean(r(isfinite(r)));
81
+ 
82
+return;                                                                                                     
83
+
0 84
new file mode 100644
... ...
@@ -0,0 +1,53 @@
1
+function [hits,ninvalid] = hits(bmus, mmax, values)
2
+
3
+%HITS Calculate number of occurances of each value.
4
+%
5
+% hits = hits(bmus,[mmax],[values])
6
+%
7
+%   h = hits(bmus);
8
+%   h = hits(bmus,length(sM.codebook)); 
9
+%
10
+%  Input and output arguments ([]'s are optional): 
11
+%   bmus     (vector) BMU indeces (or other similar) 
12
+%   [mmax]   (scalar) maximum index, default value max(bmus)
13
+%            (struct) map or topology struct from where the maximum
14
+%                     index is acquired
15
+%   [values] (vector) values associated with the data (default = 1)
16
+%
17
+%   hits     (vector) the number of occurances of each index
18
+%                     (or if values are given, their sum for each index)
19
+%   ninvalid (scalar) number of invalid indeces (NaN, Inf or 
20
+%                     <=0 or > mmax)
21
+%
22
+% See also SOM_HITS, SOM_BMUS.    
23
+
24
+% Copyright (c) 2002 by the SOM toolbox programming team.
25
+% Contributed to SOM Toolbox by Juha Vesanto, April 24th, 2002
26
+% http://www.cis.hut.fi/projects/somtoolbox/
27
+
28
+% Version 2.0beta juuso 240402
29
+
30
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
31
+
32
+if nargin<2 | isempty(mmax), 
33
+  mmax = max(bmus);
34
+elseif isstruct(mmax), 
35
+  switch mmax.type, 
36
+   case 'som_map',   mmax = prod(mmax.topol.msize);
37
+   case 'som_topol', mmax = prod(mmax.msize);
38
+   otherwise, 
39
+    error('Illegal struct for 2nd argument.')