Live Editor in MATLAB

An entirely new type of script has been introduced in MATLAB as of Version 2016a. The script is called a live script and is created using the Live Editor. A live script is much more dynamic than a simple script; it can embed equations, images, and hyperlinks in addition to formatted text. Instead of having graphs in separate windows, the graphics are shown next to the code that created them. It is also possible to put the graphs inline, under the code. The scripts that we have seen so far have been simple scripts, stored in files that have an extension of .m. Live scripts are instead stored using the .mix file format.

An example of a live script is shown in Fig. 1. In this live script named “sintest.mix,” there is text, followed by an equation, then code, more text, another equation, and more code. All of these are in separate sections. The sections are created by clicking on “code,” “text,” “equation” and so forth from the “Insert” section, with “section break” in between each type of section. All output, including error messages if there are any, are shown to the right of the code. In this case, the graphs produced by the code sections are shown to the right of the code. Clicking on the icon above the graphs will move them inline. There are several ways to create a live script. The simplest is to click on New, then Live Script. A simple script can also be transformed into a live script by choosing Save As and then Live Script. Right clicking on a set of commands from the Command History window also pops up an option to save as a Live Script.

a1

 

All of the code from the live script can be executed by choosing the Run All button. Alternatively, individual sections can be executed by clicking on the bar to the left of the section as seen in Fig. 2.

b1
Once a live script has been completed, it can be shared with others as an.mix file, or it can be converted to a PDF or HTML format. To do this, click on the down arrow under “Save,” and then either “Export to PDF” or “Export to HTML.” Live scripts can also be converted to code files with the .m extension by choosing Save As and then choosing MATLAB Code file from the drop down menu for the Format.

Link

Simple MATLAB code for Topographic ICA of images

Topographic ICA

function W=tica(Z,topoxdim,topoydim,neighbourhoodsize)

global convergencecriterion

%------------------------------------------------------------
% input parameters settings
%------------------------------------------------------------
%
% Z : whitened data
%
% topoxdim : number of components or each row of topographic grid
%
% topoydim : number of components or each col of topographic grid
%
% maxiter = 100, say : max number of iterations in algorithm
%
% neighbourhoodsize= 1...5 : neighbourhood size for ISA


%------------------------------------------------------------
% Initialize algorithm
%------------------------------------------------------------

%create matrix where the i,j-th element is one if they are neighbours
%and zero if not
H = neighbourhoodmatrix(topoxdim,topoydim,neighbourhoodsize);

%create random initial value of W, and orthogonalize it
W = orthogonalizerows(randn(topoxdim*topoydim,size(Z,1)));

%read sample size from data matrix
N = size(Z,2);

%Initial step size
mu = 1;

%------------------------------------------------------------
% Start algorithm
%------------------------------------------------------------

writeline('Doing topographic ICA. Iteration count: ')

iter = 0;
notconverged = 1;

while notconverged & (iter<2000) %max 2000 iterations

 iter=iter+1;

 %print iterations left
 writenum(iter);

 %-------------------------------------------------------------
 % Gradient step for topographic ICA
 %-------------------------------------------------------------

 % Compute separately estimates of independent components to speed up
 Y=W*Z;

 %compute local energies K
 K=H*Y.^2;

 % This is nonlinearity corresponding to generalized exponential density
 % (note minus sign)
 epsilon=0.1;
 gK = -(epsilon+K).^(-0.5);

 % Calculate convolution with neighborhood
 F=H*gK;

 % This is the basic gradient
 grad = (Y.*F)*Z'/N;

 % project gradient on tangent space of orthogonal matrices (optional)
 grad=grad-W*grad'*W;

 %store old value
 Wold = W;

 % do gradient step
 W = W + mu*grad;

 % Orthogonalize rows or decorrelate estimated components
 W = orthogonalizerows(W);


 % Adapt step size every tenth step
 if rem(iter,1)==0 | iter==1

 changefactor=4/3;

 % Take different length steps
 Wnew1 = Wold + 1/changefactor*mu*grad;
 Wnew2 = Wold + changefactor*mu*grad;
 Wnew1=orthogonalizerows(Wnew1);
 Wnew2=orthogonalizerows(Wnew2);

 % Compute objective function values
 J1=-sum(sum(sqrt(epsilon+H*(Wnew1*Z).^2)));
 J2=-sum(sum(sqrt(epsilon+H*(Wnew2*Z).^2)));
 J0=-sum(sum(sqrt(epsilon+H*(W*Z).^2)));

 % Compare objective function values, pick step giving minimum
 if J1>J2 & J1>J0
 % Take shorter step because it is best
 mu = 1/changefactor*mu;
 %%fprintf(' [mu=%f] ',mu);
 W=Wnew1;
 elseif J2>J1 & J2>J0
 % Take longer step because it is best
 mu = changefactor*mu;
 %%fprintf(' [mu=%f] ',mu);
 W=Wnew2;
 end
 end

 % Check if we have converged
 if norm(grad,'fro') < convergencecriterion *topoxdim*topoydim;
 notconverged=0; end

end %of gradient iterations loop

gathers patches from the grey-scale images

function X = sampleimages(samples, winsize)

% gathers patches from the grey-scale images, no preprocessing done yet
%
% INPUT variables:
% samples total number of patches to take
% winsize patch width in pixels
%
% OUTPUT variables:
% X the image patches as column vectors

%IMPORTANT: IMAGE DATA IS SUPPOSED TO BE IN DIRECTORY ./DATA/
%CHANGE THIS ON LINE 34 IF YOU HAVE MOVED THE DATA ELSEWHERE

%----------------------------------------------------------------------
% Gather rectangular image patches
%----------------------------------------------------------------------

% We have a total of 13 images.
dataNum = 13;

% This is how many patches to take per image
getsample = floor(samples/dataNum);

% Initialize the matrix to hold the patches
X = zeros(winsize^2,samples);

sampleNum = 1;
for i=(1:dataNum)

 % Even things out (take enough from last image)
 if i==dataNum, getsample = samples-sampleNum+1; end

 % Load the image. Change the path here if needed.
 I = imread(['data/' num2str(i) '.tiff']);

 % Transform to double precision
 I = double(I);

 % Normalize to zero mean and unit variance (optional)
 I = I-mean(mean(I));
 I = I/sqrt(mean(mean(I.^2)));

 % Sample patches in random locations
 sizex = size(I,2);
 sizey = size(I,1);
 posx = floor(rand(1,getsample)*(sizex-winsize-2))+1;
 posy = floor(rand(1,getsample)*(sizey-winsize-1))+1;
 for j=1:getsample
 X(:,sampleNum) = reshape( I(posy(1,j):posy(1,j)+winsize-1, ...
 posx(1,j):posx(1,j)+winsize-1),[winsize^2 1]);
 sampleNum=sampleNum+1;
 end

end



Removes DC component from image patches

Removes DC component from image patches
Data given as a matrix where each patch is one column vectors
That is, the patches are vectorized.

function Y=removeDC(X);

% Subtract local mean gray-scale value from each patch in X to give output Y

Y = X-ones(size(X,1),1)*mean(X);

return;

 

display receptive field(s) or basis vector(s) for image patches

function plotrf( A, cols ,filename )

% display receptive field(s) or basis vector(s) for image patches
%
% A the basis, with patches as column vectors
% cols number of columns (x-dimension of grid)
% filename where to save the image, If [], don't save but plot
% In some versions of matlab, must use '' instead of []

global figurepath

%set colormap
colormap(gray(256));

%normalize each patch
A=A./(ones(size(A,1),1)*max(abs(A)));

% This is the side of the window
dim = sqrt(size(A,1));

% Helpful quantities
dimp = dim+1;
rows = floor(size(A,2)/cols); %take floor just in case cols is badly specified

% Initialization of the image
I = ones(dim*rows+rows-1,dim*cols+cols-1);

%Transfer features to this image matrix
for i=0:rows-1
 for j=0:cols-1
 % This sets the patch
 I(i*dimp+1:i*dimp+dim,j*dimp+1:j*dimp+dim) = ...
 reshape(A(:,i*cols+j+1),[dim dim]);
 end
end

%Save of plot results
imagesc(I);
axis equal
axis off
print('-dps',[figurepath,filename,'.eps'])




MATLAB code for PCA on image patches

PCA ,or Principal Component Analysis, is defined as the following in wikipedia[1]:

A statistical procedure that uses an orthogonal transformation to convert a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables called principal components.

In other words, it is a technique that reduces an data to its most important components by removing correlated characteristics. Personally, I like to think of it as reducing an data to its “essence.”

PCA Introduction
Principal component analysis, or what I will throughout the rest of this article refer to as PCA, is considered the go-to tool in the machine learning arsenal. It has applications in computer vision, big data analysis, signal processing, speech recognition, and more. Many articles, professors, and textbooks tend to shroud the method behind a wall of text or equations. Let me tear down that shroud!

function [V,E,D] = pca(X)

% do PCA on image patches
%
% INPUT variables:
% X matrix with image patches as columns
%
% OUTPUT variables:
% V whitening matrix
% E principal component transformation (orthogonal)
% D variances of the principal components


% Calculate the eigenvalues and eigenvectors of the new covariance matrix.
covarianceMatrix = X*X'/size(X,2);
[E, D] = eig(covarianceMatrix);

% Sort the eigenvalues and recompute matrices
[dummy,order] = sort(diag(-D));
E = E(:,order);
d = diag(D);
dsqrtinv = real(d.^(-0.5));
Dsqrtinv = diag(dsqrtinv(order));
D = diag(d(order));
V = Dsqrtinv*E';

Create neighbourhood structure for TICA

Time-structure Independent Components Analysis

Create neighbourhood structure for TICA, expressed in a single matrix

 

function H=neighbourhoodmatrix(xdim,ydim,size)

% n is dimension of data
% size is longest distance of neighbours (i.e. neighbourhood size) e.g. 1 or 2

% This will hold the neighborhood function entries
H = zeros(xdim*ydim,xdim*ydim);

% Step through nodes one at a time to build the matrix
ind=0;
for y=1:ydim

 for x=1:xdim

 ind=ind+1;


 % Rectangular neighbors
 [xn,yn] = meshgrid( (x-size):(x+size), (y-size):(y+size) );
 xn = reshape(xn,[1 (size*2+1)^2]);
 yn = reshape(yn,[1 (size*2+1)^2]);

 % Cycle round to create torus
 i = find(yn<1); yn(i)=yn(i)+ydim;
 i = find(yn>ydim); yn(i)=yn(i)-ydim;
 i = find(xn<1); xn(i)=xn(i)+xdim;
 i = find(xn>xdim); xn(i)=xn(i)-xdim;

 % Set neighborhood
 H( ind, (yn-1)*xdim + xn )=1;

 end
end


%Normalize to that row norm = 1
H=H./(sqrt(sum(H.^2))'*ones(1,xdim*ydim));

Initialize random number generator in MATLAB

random number MATLAB rng

This is done at the beginning of each section so that each section can be run separately
Seed values 0 and 0 are used in the book

 

function initializerandomseeds()

%This is mainly used in selecting initial points for algorithms
randn('state',0);

%This one is used in random sampling of patches from images
rand('state',0);

return;

Simple MATLAB code for ICA of images

independent components analysis MATLAB code

%Simple code for ICA of images
%Aapo Hyvنrinen, for the book Natural Image Statistics

function W=ica(Z,n)

global convergencecriterion

%------------------------------------------------------------
% input parameters settings
%------------------------------------------------------------
%
% Z : whitened image patch data
%
% n = 1..windowsize^2-1 : number of independent components to be estimated


%------------------------------------------------------------
% Initialize algorithm
%------------------------------------------------------------

%create random initial value of W, and orthogonalize it
W = orthogonalizerows(randn(n,size(Z,1)));

%read sample size from data matrix
N=size(Z,2);

%------------------------------------------------------------
% Start algorithm
%------------------------------------------------------------

writeline('Doing FastICA. Iteration count: ')

iter = 0;
notconverged = 1;

while notconverged & (iter<2000) %maximum of 2000 iterations

 iter=iter+1;

 %print iteration count
 writenum(iter);

 % Store old value
 Wold=W;

 %-------------------------------------------------------------
 % FastICA step
 %-------------------------------------------------------------

 % Compute estimates of independent components
 Y=W*Z;

 % Use tanh non-linearity
 gY = tanh(Y);

 % This is the fixed-point step.
 % Note that 1-(tanh y)^2 is the derivative of the function tanh y
 W = gY*Z'/N - (mean(1-gY'.^2)'*ones(1,size(W,2))).*W;

 % Orthogonalize rows or decorrelate estimated components
 W = orthogonalizerows(W);

 % Check if converged by comparing change in matrix with small number
 % which is scaled with the dimensions of the data
 if norm(abs(W*Wold')-eye(n),'fro') < convergencecriterion * n;
 notconverged=0; end

end %of fixed-point iterations loop








subroutine for creating gabors MATLAB code

Gabor_filter

%this routine is not really used much:
%it is *not* used in analysis of RF's but only in creating
%a single gabor illustration of in Chapter 1.

function gaborf=gabor(windowsize,xo,yo,sigma,freq,type);

%size of window is [-5,5] x [-5,5]
%y-axis is now horizontal
%"type" controls whether is it is odd- or even-symmetric
%"windowsize" is size of window in pixels

g2d=zeros(windowsize);
m=(windowsize+1)/2;
fact=5/(windowsize-m);

for x=1:windowsize
for y=1:windowsize

if type=='o';
g2d(y,x)=exp(-(((x-m)*fact-xo)/sigma)^2-(((y-m)*fact-yo)/sigma)^2)*sin(((y-m)*fact-yo)*freq);
else
g2d(y,x)=exp(-(((x-m)*fact-xo)/sigma)^2-(((y-m)*fact-yo)/sigma)^2)*cos(((y-m)*fact-yo)*freq);
end

end
end

gaborf=reshape( g2d,[windowsize^2 1]);