Merge branch 'master' of https://github.com/albertoandreottiATgmail/DeepLearnToolbox

Author: albertoandreottiATgmail

CNN/cnnavg.m

@@ -0,0 +1,15 @@
+function avgnet = cnnavg(avgnet, net)
+    for l = 2 : numel(net.layers)
+        if strcmp(net.layers{l}.type, 'c')
+            for j = 1 : numel(net.layers{l}.a)
+                for ii = 1 : numel(net.layers{l - 1}.a)
+                    avgnet.layers{l}.k{ii}{j} = (avgnet.layers{l}.k{ii}{j} + net.layers{l}.k{ii}{j})/2;
+                end
+                avgnet.layers{l}.b{j} = avgnet.layers{l}.b{j} + net.layers{l}.b{j};
+            end
+        end
+    end
+
+    avgnet.ffW = (avgnet.ffW + net.ffW)/2;
+    avgnet.ffb = (avgnet.ffb + net.ffb)/2 ;
+end

CNN/cnnbp.m

@@ -41,7 +41,7 @@
         if strcmp(net.layers{l}.type, 'c')
             for j = 1 : numel(net.layers{l}.a)
                 for i = 1 : numel(net.layers{l - 1}.a)
-                    net.layers{l}.dk{i}{j} = convn(flipall(net.layers{l - 1}.a{i}), net.layers{l}.d{j}, 'valid') / size(net.layers{l}.d{j}, 3);
+                    net.layers{l}.dk{i}{j} = custom_convn(flipall(net.layers{l - 1}.a{i}), net.layers{l}.d{j}, 'valid') / size(net.layers{l}.d{j}, 3);
                 end
                 net.layers{l}.db{j} = sum(net.layers{l}.d{j}(:)) / size(net.layers{l}.d{j}, 3);
             end

CNN/cnncopy.m

@@ -0,0 +1,18 @@
+function cpnet = cnncopy(cpnet, net)
+
+
+
+    for l = 2 : numel(net.layers)
+        if strcmp(net.layers{l}.type, 'c')
+            for j = 1 : numel(net.layers{l}.a)
+                for ii = 1 : numel(net.layers{l - 1}.a)
+                    cpnet.layers{l}.k{ii}{j} = net.layers{l}.k{ii}{j};
+                end
+                cpnet.layers{l}.b{j} = net.layers{l}.b{j};
+            end
+        end
+    end
+
+    cpnet.ffW = net.ffW;
+    cpnet.ffb = net.ffb;
+end

CNN/cnntrain.m

@@ -1,27 +1,29 @@
+
 function net = cnntrain(net, x, y, opts)
     m = size(x, 3);
     numbatches = m / opts.batchsize;
     if rem(numbatches, 1) ~= 0
         error('numbatches not integer');
     end
+
     net.rL = [];
     for i = 1 : opts.numepochs
         disp(['epoch ' num2str(i) '/' num2str(opts.numepochs)]);
         tic;
         kk = randperm(m);
-        for l = 1 : numbatches
-            batch_x = x(:, :, kk((l - 1) * opts.batchsize + 1 : l * opts.batchsize));
-            batch_y = y(:,    kk((l - 1) * opts.batchsize + 1 : l * opts.batchsize));
+        %how many processes?
+        numWorkers = 4
+		pids = 0:numWorkers-1;
+		starts = pids * numbatches / numWorkers
+	
+		%process starts
+		turn = 0;
 
-            net = cnnff(net, batch_x);
-            net = cnnbp(net, batch_y);
-            net = cnnapplygrads(net, opts);
-            if isempty(net.rL)
-                net.rL(1) = net.L;
-            end
-            net.rL(end + 1) = 0.99 * net.rL(end) + 0.01 * net.L;
-        end
+        pararrayfun(numWorkers,
+                    @(starts, pids)process_batch(x, y, kk, net, turn, starts, (numbatches/numWorkers),  pids, numWorkers, opts),
+                    starts,
+					pids,
+					"ErrorHandler" , @eh);
         toc;
     end
-    
 end

CNN/convadd.m

@@ -0,0 +1,8 @@
+	%a is the image, k the kernel, o the output, i is the index.
+function result = convadd(a, k, m, pid, chunkSize)
+	init = pid*chunkSize + 1;
+	result = zeros(m, m);
+	for i=init:init+chunkSize -1
+	    result += conv2D(a(:,:, i), k(:,:, i));
+	end
+end

CNN/convn_valid.m

@@ -0,0 +1,23 @@
+%Convolution for 3 dimensional vectors using conv2
+%equivalent to convn(A,B, 'valid')
+
+function result = convn_valid(A, B)
+
+    m = size(A, 1) - size(B, 1) + 1;
+    numWorkers = 2;
+	
+    function retcode = eh(error)
+        a = error
+        retcode = zeros(25, 1);	
+    end
+	
+	%each worker will write its output to specific part of the output
+	chunkSize = size(A,3)/numWorkers;
+	result = pararrayfun(numWorkers,  @(i)convadd(A, B, m, i, chunkSize), 0:numWorkers-1, "ErrorHandler" , @eh);
+	
+    for j=m:numWorkers*m:m  
+	    result(:,1:m) += result(:,j+1:j+m);
+	end
+	
+	result = result(1:m,1:m);
+end

CNN/convnfft.m

@@ -0,0 +1,221 @@
+function A = convnfft(A, B, shape, dims, options)
+% CONVNFFT  FFT-BASED N-dimensional convolution.
+%   C = CONVNFFT(A, B) performs the N-dimensional convolution of
+%   matrices A and B. If nak = size(A,k) and nbk = size(B,k), then
+%   size(C,k) = max([nak+nbk-1,nak,nbk]);
+% 
+%   C = CONVNFFT(A, B, SHAPE) controls the size of the answer C:
+%       'full'   - (default) returns the full N-D convolution
+%       'same'   - returns the central part of the convolution that
+%                  is the same size as A.
+%       'valid'  - returns only the part of the result that can be
+%                  computed without assuming zero-padded arrays.
+%                  size(C,k) = max([nak-max(0,nbk-1)],0).
+%
+%   C = CONVNFFT(..., SHAPE, DIMS) with DIMS is vector of dimensions where
+%       the convolution will be carried out. By default DIMS is
+%       [1:max(ndims(A),ndims(B))] (all dimensions). A and B must have the
+%       same lengths on other dimensions.
+%   C = CONVNFFT(..., SHAPE, DIMS, GPU)
+%       GPU is boolean flag, see next
+%
+%   C = CONVNFFT(..., SHAPE, DIMS, OPTIONS)
+%
+%   OPTIONS is structure with following optional fields
+%       - 'GPU', boolean. If GPU is TRUE Jacket/GPU FFT engine will be used
+%       By default GPU is FALSE.
+%       - 'Power2Flag', boolean. If it is TRUE, use FFT with length rounded
+%       to the next power-two. It is faster but requires more memory.
+%       Default value is TRUE.
+%
+% Class support for inputs A,B:
+% float: double, single
+%
+% METHOD: CONVNFFT uses Fourier transform (FT) convolution theorem, i.e.
+%         FT of the convolution is equal to the product of the FTs of the
+%         input functions.
+%         In 1-D, the complexity is O((na+nb)*log(na+nb)), where na/nb are
+%         respectively the lengths of A and B.
+%
+% Usage recommendation:
+%         In 1D, this function is faster than CONV for nA, nB > 1000.
+%         In 2D, this function is faster than CONV2 for nA, nB > 20.
+%         In 3D, this function is faster than CONVN for nA, nB > 5.
+% 
+% See also conv, conv2, convn.
+% 
+%   Author: Bruno Luong <[email protected]>
+%   History:
+%       Original: 21-Jun-2009
+%       23-Jun-2009: correct bug when ndims(A)<ndims(B)
+%       02-Sep-2009: GPU/JACKET option
+%       04-Sep-2009: options structure
+%       16-Sep-2009: inplace product
+
+if nargin<3 || isempty(shape)
+    shape = 'full';
+end
+
+if nargin<5 || isempty(options)
+    options = struct();
+elseif ~isstruct(options) % GPU options
+    options = struct('GPU', options);
+end
+
+nd = max(ndims(A),ndims(B));
+% work on all dimensions by default
+if nargin<4 || isempty(dims)
+    dims = 1:nd;
+end
+dims = reshape(dims, 1, []); % row (needed for for-loop index)
+
+% GPU enable flag
+GPU = getoption(options, 'GPU', false);
+% Check if Jacket is installed
+GPU = GPU && ~isempty(which('ginfo'));
+
+% IFUN function will be used later to truncate the result
+% M and N are respectively the length of A and B in some dimension
+switch lower(shape)
+    case 'full',
+        ifun = @(m,n) 1:m+n-1;
+    case 'same',
+        ifun = @(m,n) ceil((n-1)/2)+(1:m);
+    case 'valid',
+        ifun = @(m,n) n:m;
+    otherwise
+        error('convnfft: unknown shape %s', shape);
+end
+
+classA = class(A);
+classB = class(B);
+ABreal = isreal(A) && isreal(B);
+
+% Special case, empty convolution, try to follow MATLAB CONVN convention
+if any(size(A)==0) || any(size(B)==0)
+    szA = zeros(1,nd); szA(1:ndims(A))=size(A);
+    szB = zeros(1,nd); szB(1:ndims(B))=size(B);
+    % Matlab wants these:
+    szA = max(szA,1); szB = max(szB,1);
+    szC = szA;
+    for dim=dims
+        szC(dim) = length(ifun(szA(dim),szB(dim)));
+    end
+    A = zeros(szC,classA); % empty -> return zeros
+    return
+end
+
+power2flag = getoption(options, 'Power2Flag', true);
+if power2flag
+    % faster FFT if the dimension is power of 2
+    lfftfun = @(l) 2^nextpow2(l);
+else
+    % slower, but smaller temporary arrays
+    lfftfun = @(l) l;
+end
+
+if GPU % GPU/Jacket FFT
+    if strcmp(classA,'single')
+        A = gsingle(A);
+    else
+        A = gdouble(A);
+    end
+    if strcmp(classB,'single')
+        B = gsingle(B);
+    else
+        B = gdouble(B);
+    end
+    % Do the FFT
+    subs(1:ndims(A)) = {':'};
+    for dim=dims
+        m = size(A,dim);
+        n = size(B,dim);
+        % compute the FFT length
+        l = lfftfun(m+n-1);
+        % We need to swap dimensions because GPU FFT works along the
+        % first dimension
+        if dim~=1 % do the work when only required
+            swap = 1:nd;
+            swap([1 dim]) = swap([dim 1]);
+            A = permute(A, swap);
+            B = permute(B, swap);
+        end
+        A = fft(A,l);
+        B = fft(B,l);
+        subs{dim} = ifun(m,n);
+    end
+else % Matlab FFT
+    % Do the FFT
+    subs(1:ndims(A)) = {':'};
+    for dim=dims
+        m = size(A,dim);
+        n = size(B,dim);
+        % compute the FFT length
+        l = lfftfun(m+n-1);
+        A = fft(A,l,dim);
+        B = fft(B,l,dim);
+        subs{dim} = ifun(m,n);
+    end
+end
+ 
+if GPU
+    A = A.*B;
+    clear B
+else
+    % inplace product to save 1/3 of the memory
+	% Modified by Alberto Andreotti([email protected])
+    %inplaceprod(A,B);
+	A(:) = A(:).*B(:);
+end
+
+% Back to the non-Fourier space
+if GPU % GPU/Jacket FFT
+    for dim=dims(end:-1:1) % reverse loop
+        A = ifft(A,[]);
+        % Swap back the dimensions
+        if dim~=1 % do the work when only required
+            swap = 1:nd;
+            swap([1 dim]) = swap([dim 1]);
+            A = permute(A, swap);
+        end        
+    end   
+else % Matlab IFFT  
+    for dim=dims
+        A = ifft(A,[],dim);
+    end
+end
+
+% Truncate the results
+if ABreal
+    % Make sure the result is real
+    A = real(A(subs{:}));
+else
+    A = A(subs{:});
+end
+
+% GPU/Jacket
+if GPU
+    % Cast the type back
+    if strcmp(class(A),'gsingle')
+        A = single(A);
+    else
+        A = double(A);
+    end
+end
+
+end % convnfft
+
+
+%% Get defaut option
+function value = getoption(options, name, defaultvalue)
+% function value = getoption(options, name, defaultvalue)
+    value = defaultvalue;
+    fields = fieldnames(options);
+    found = strcmpi(name,fields);
+    if any(found)
+        i = find(found,1,'first');
+        if ~isempty(options.(fields{i}))
+            value = options.(fields{i});
+        end
+    end
+end

CNN/custom_convn.m

@@ -0,0 +1,29 @@
+% This function is used to choose between different implementations of the convn function
+% according to the platform.
+% This is mainly useful for overcoming an Octave's bug in convn(), https://savannah.gnu.org/bugs/?39314
+% parameters: x is the chunk of the image/samples, k is the kernel.
+
+function result = custom_convn(x, k, convmode) 
+    
+    %FFT works with batchsize>50, but still really slow, about 10000 seconds. It should produce an error near 0.18.
+	%this convnfft is taken from http://www.mathworks.com/matlabcentral/fileexchange/24504-fft-based-convolution.
+    if exist('convmode') && strcmp(convmode, 'fft')
+       result = convnfft(x,k, 'valid');
+       return
+	end
+	
+	%the 'valid' version of convolution has problems in Octave, use 'same' instead.
+	if isOctave()
+	   %Alternative to convnftt, use for small batch size ~ 5, will give 2676.56 seconds, otherwise it will explode(too long running time).
+       start = size(x,1) - size(k,1);
+       fin = 2*start;
+       %note: if x and k have not the same size in the third dimension, middle could be a range.
+       middle = floor(size(x,3)/2) + 1;
+       result = convn (x, k, "same")(start:fin, start:fin, middle);
+	else
+	   %we're running matlab
+	   result = convn(x, k, 'valid');
+	end
+
+end
+

CNN/eh.m

@@ -0,0 +1,5 @@
+%here goes what to put in the output when the function fails.
+function retcode = eh(err)
+    a = err
+    retcode = zeros(26,1).+255;	
+end