Matlab: Vector indexing;

This is related to this post, and this post.

Q:
Suppose I have a vector like: v = [ 3 2 4 1 3 2...]
I want to obtain a new vector like: w = [3 3 3 2 2 4 4 4 4 1 3 3 3 2 2...]
basically every value repeates it's own number of times.

A1 (one-liner):

w = cell2mat(arrayfun(@(x) x+zeros(1,x), v, 'uni', 0))

A2:

v(v == 0) = [] % Need to get rid of zeros in the general case.
c = cumsum(v)
w = zeros(1,c(end))
w([1,c(1:end-1)+1]) = 1
w = v(cumsum(w))

Speed test:

% v = randi(800,1,90000)+1;
Method 1:
Elapsed
time is 1.997060 seconds.
Method 2:
Elapsed time is 1.630775 seconds.

Two methods are about same speed.

Matlab: vector indexing;

OP here.

Q:

I have a reference vector r, and a target vector v, I want to find, for every element in v, the index of first element that is greater than or equal to it in r.

With a for loop it shall look like this:

for k = 1:length(r), idx(k) = find(r >= v(k), 1, 'first'); end

A1 (the 'standard' vectorization):

[dummy,idx] = max(bsxfun(@ge,r(:),v(:).'));

A2 (interesting):

idx = ceil(interp1(r,1:length(r),v));

A3 (same logic as A2 but better):

[dummy, idx] = histc(v,r);
idx = idx+1;

Speed test:

% r = cumsum(rand(100000,1));
% r = r./r(end);
% v = rand(20000,1);

Method 1: % actually slower than the for loop ... memory.
Elapsed time is 6.582422 seconds.

Method 2:
Elapsed time is 0.015725 seconds.

Method 3:
Elapsed time is 0.006501 seconds.

The For loop:
Elapsed time is 4.209633 seconds.

Matlab: indexing; find first value in each row

OP here.

Q:
How can I find the first index of certain value in each row without using a for loop?

A1:

m = round(rand(9,4)*3); % the data;
[r,c] = size(m);
n = 3; % number to find.

idx = mod(sum(cumprod(double(m ~= n),2),2)+1,c+1)

A2:

idx = zeros(size(m,1),1);
tmp = m == n;
[r,c] = find(tmp & cumsum(tmp,2) == 1);
idx(r) = c;

A3 (fastest):

[v,idx] = max(m == n,[],2)
idx = idx.*(v==1);
% or idx(v==0) = 0;

Speed test:

% m = round(rand(1600000,4)*3);
A1:
Elapsed time is 0.293340 seconds.
A2:
Elapsed time is 0.407424 seconds.
A3:
Elapsed time is 0.110964 seconds.

Matlab: matrix indexing; find duplicate entries in a matrix and replace them with average

OP here.

Q:
I've got a matrix of integer values with three columns (X,Y,Z) and about 10k rows.
X and Y are in range 0..1000 and Z is in the range of 0..200.

In this matrix, there are some pairs of (X,Y) that occur more than once, but with different Z values. Now I want to remove these duplicates and replace all of them with (X,Y,mean of Z).

A:

[u,id1,id2] = unique(A(:,1:2),'rows');
B = [u,accumarray(id2,A(:,3))./accumarray(id2,1)];
% or accumarray(id2,A(:,3),[],@mean) but this is much slower.

Speed test:

% A = [round(rand(10000,2)*1000),round(rand(10000,1))*200];
Elapsed time is 0.005959 seconds.

Matlab: Block randomization

Randomized block design is a common experiment design paradigm. Suppose there are 3 conditions, it seeks to create a random sequence like this:

[3 1 2 2 1 3 1 3 2 1 3 2 1 2 3 2 3 1]

Basically random perm is performed within every size-3 block.

Here's a tip to do this in matlab in one line:

[dump,idx] = sort(rand(nCondition,nTrial));

Matlab: Vector operation

OP here.

Q:

Is there a vectorized way for this code?

for i=1:size(m,2)
m(:,i) = m(:,i).*v(i);
end

A1:

mv = bsxfun(@times,m,v');

A2:

mv = m*diag(v);

all these solutions are much better than using repmat etc.

Matlab: Matrix Indexing

OP here.


Q:

Suppose I have a 3D matrix A:

A(:,:,1) =
[0.1, 0.2, 0;
0.2, 0, 0;
0.6, 0.4, 0]

A(:,:,2) =
[0.2, 0.7, 0;
0.3, 0.8, 0;
0.1, 0, 0]

I want to find, for each row, the element before the first occurrence of zero and subtract it from 1. The outcome should look like this:

B(:,:,1) =
[0.1, 0.8, 0;
0.8, 0, 0;
0.4, 0.6, 0]

B(:,:,2) =
[0.2,0.3,0;
0.3, 0.2, 0;
0.9, 0, 0]

A1:

% location of first occurrence of zero in each row.
z = cumsum(A == 0,2) == 1;
% the element before it.
y = circshift(z,[0,-1]);
B = A;
B(y) = 1-B(y);

A2 (a beautiful solution):

% convolution along columns.
y = convn(logical(A),[-1 1],'same') > 0;
B = A;
B(y) = 1-B(y);

EEG: Topographic contour on 3d head model.

Took me a couple days to program it. But I'm happy.

Matlab: Vectorization.

OP here.

Q:

I have

c = cell(n,1); % Each cell is an array of integers.

How do you vectorize this loop? it takes very long.

for i=1:n-1
for j=i+1:n
lic(i,j)=length(intersect(c{i},c{j}));
end
end

A1:

Generate some data for testing:

the data:

n=200; % number of sets
m=100; % element value are from 0 to m-1
meannumel=50; % average number of elements in lists
s = 20; % standard deviation of number of elements
c=cell(1,n);

for k=1:n
nelk = meannumel+round(s*randn);
nelk = min(max(nelk,1),m);
c{k} = floor(m*rand(1,nelk));
end

The engine:

n = length(c);
c = cellfun(@(x) unique(x(:)), c, 'uni', false);

s = cumsum([1 cellfun(@length,c)]).'; % boundary of cell, in vector x;
x = cat(1,c{:}); % all elements, in one long vector;
[x,p] = sort(x); % x is sorted now;

u = cumsum([1;diff(x) ~= 0]);
% before cumsum: boundary of unique (different) elements in x; cumsum: to get a vector u, so that u = 11111223333, # of index equal to count of value;

t = cumsum(accumarray(s,1));
% accumarray: here equal to tmp = zeros(s(end),1); tmp(s) = 1; cumsum: to get a vector t, so that t = 1111222223333 # of index equal to count of that value. here s is boundary of cell, so this t is actually which cell each element belongs;

A = sparse(u,t(p),1,u(end),n);
% sparse: only store non-zero elements, save memory space; syntax: p = sparce(i,j,s,m,n) will create a m x n sparse matrix, with elements p(i(k),j(k)) = s(k); here the intention is to add 1 at, the unique element index, vs. the cell index that each elements belongs.
% A is in the format of: unique item index by which cell, here A is essentially a logical matrix (valued 0, 1 only)

A = triu(full(A'*A),1);
% A'*A is the trick to count intersected values; colums and rows in A are only 1 or 0. matrix muptiplication sums the ones.

A2:
A2 is slightly improved over A1 in that the 'unique' is only called once.

% c must be a list of *row* vectors
[dummy dummy cmap] = unique([c{:}]);
% length(dummy) is the number of unique values in c. cmap is in the form of [1 2 3 1 1 2 2 2 4 ...], which specifies unique values indeces in [c{:}]. length(cmap) = numel([c{:}];
nele = cellfun(@length,c(:));
setid = cumsum(accumarray(cumsum([1; nele]),1));
A = sparse(setid(1:end-1),cmap(:),1,length(c),length(dummy));
A = spones(A); % SPONES Replace nonzero sparse matrix elements with ones.
A = triu(full(A*A.'),1);

Speed comparison:

The for loop: 3.55 sec.
A1: 0.022 sec.
A2: 0.015 sec.

Matlab: Vector manipulation.

OP here.


Q:

Suppose I have

p = [ 0.81 0.90 0.12 0.91 0.63 0.09 0.27 0.54 0.95 0.96 0.15 0.97 0.95 0.48 0.80 0.14 0.42 0.21 0.54 0.14];

v = [ 1 0 0 1 0 1 0 1 1 1 1 1 0 0 0 0 1 1 0 1];

Where the 1 defines the boundary of blocks. So the first block is 1 0 0, second 1 0, third 1 0, fourth 1, fifth 1, ... and so on.

I want to sum values in p according to blocks defined by v. The output in this case should look like:

r = [1.83 1.54 0.36 0.54 0.95 0.96 0.15 3.34 0.42 0.75 0.14];

A:

r = accumarray(cumsum(v(:)),p(:))

Matlab: Vector indexing;

Problem:

Suppose I have a vector:

x = [1 1 1 1 2 2 2 8 8 8 8 8 8 0 0 9 9 7 5 5 2 2 2 2 2]

I want to obtain:

b = [1 2 3 4 1 2 3 1 2 3 4 5 6 1 2 1 2 1 1 2 1 2 3 4 5]

Basically each block of integer is replaced by continuous index of the count.


Solution:

Some initial processing to detect block size and locations:

y = find([true;diff(x(:)) ~= 0;true]);
c = diff(y);
v = x(y(1:end-1));

method 1:

b = cell2mat(arrayfun(@(z)(1:z)',c,'uniformout',false));

method 2:

b = ones(size(x));
b(y(2:end-1)) = 1-c(1:end-1);
b = cumsum(b);

method 3:

g = (1:max(c))'*ones(1,numel(c));
b = g(bsxfun(@ge,c,1:max(c))');

Speed comparison:

% The testing data;
x = cell2mat(arrayfun(@(y,z)(y+zeros(1,z)), ...
randi(100,10000,1),randi(2000,10000,1),'uniformout',false));

% Method 1; 0.25 sec
% Method 2; 0.17 sec
% Method 3; 1.05 sec

Method 2 is better.

Matlab: Matrix Indexing

OP here.


Q:

How do you replace the first nonzero element in each row with -1 ?


A = [
50 75 0;
50 0 100;
0 75 100;
75 100 0;
0 75 100;
0 75 100];

to
-1 75 0
-1 0 100
0 -1 100
-1 100 0
0 -1 100
0 -1 100

A1:

[I, J] = ind2sub(size(A), find(A ~= 0));
[b, c] = unique(I, 'first');
A(sub2ind(size(A), b, J(c))) = -1

A2:

B=cumsum(A~=0,2)>0
B=[false(size(B,1),1) B]
A(logical(diff(B,1,2))) = -1

A3:

A((cumsum(~~A, 2) == 1) & (A ~= 0)) = -1
% notice ~~A is equal to A ~= 0

DSP: How to understand convolution.

OP here.

Definition from wikipedia:

The convolution of ƒ and g is written ƒ∗g. It is defined as the integral of the product of the two functions after one is reversed and shifted.

Let the unit response of a filter be denoted by u, and let x denote the input, y denote the output. Think of the input as the sum of many different delayed and scaled unit inputs:

t = [1, 2, 3, 4, ... n]

x = [x(1), 0, 0, 0, 0 ... 0]
+ [0, x(2), 0, 0, 0 ... 0]
+ [0, 0, x(3), 0, 0 ... 0]
+ ...
+ [0 ... 0, 0, 0, 0, x(n)]

then the output could be seen as the sum of responses to these unit inputs:

y = [u(1)*x(1), u(2)*x(1), u(3)*x(1) ... u(n)*x(1)]
+ [0, u(1)*x(2), u(2)*x(2), u(3)*x(2) ... u(n-1)*x(2)]
+ [0, 0, u(1)*x(3), u(2)*x(3), u(3)*x(3) ... u(n-2)*x(3)]
+ ...
+ [0, 0, ... 0, 0, u(1)*x(n)]

so at particular time n, the output takes the form of:

y(n) = u(n)*x(1) + u(n-1)*x(2) + u(n-2)*x(3) + ... + u(1)*x(n)

which explains why "flipping" is needed in convolution.

Matlab: Matrix manipulation.

OP here.

Q:

I have the following matrix:

[1 5
4 9]

What's the fastest way, and without a for loop, to expand it to:

[1 1 5 5
1 1 5 5
4 4 9 9
4 4 9 9]

A1:

A = [1 4; 5 9] ;
kron(A, ones(3)) ;

A2:

idx = cumsum(ones(2),2 )
B = A(idx, idx)

A3:

[m, n] = size(A);
mdup = 4;
ndup = 3;
vx = ceil((1:m*mdup)/mdup);
vy = ceil((1:n*ndup)/ndup);

B = A(vx, vy)

A4:

[m,n] = size(A);
mdup = 4;
ndup = 3;
AX = repmat(A, [1 1 mdup ndup]);
AX = permute(AX, [3 1 4 2]);
AX = reshape(AX, m*mdup, n*ndup);

Matlab: Vector manipulation.

OP here.


Q:

Suppose

a = [0 0 0 1 1 0 0 0 1 0 0 0 0 1 0 0 1 0 1 ];

I'd like to organize this vector into a matrix with the following criteria: every time I see a 1,
I want to copy that 1 and 2 other members before it into a new line, so my resulting matrix will look like this:

b = [
0 0 1
0 1 1
0 0 1
0 0 1
0 0 1
1 0 1];

Is there a vectorized solution?

A:

apad=[0 0 a];
i1 = find(apad == 1);
b = apad(bsxfun(@plus, i1(:), (-2:0)))

Matlab: Split a vector into continuous chunks.

Original post here.


Q:

Suppose I have [1 2 3 4 8 9 16 17 18], how do I get [1 2 3 4], [8 9], [16 17 18] in a vectorized way?

A:

a = [1 2 3 4 8 9 16 17 18];

b=a-(1:length(a));
b = [true; diff(b(:)) ~= 0; true];
split = mat2cell(a(:).', 1, diff(find(b)));

split{:}

Matlab: Automatic conversion to element-by-element operation.

Original post here.

Q:

Is there a function which automatically puts points in front of / , * or ^ in a formula?

e.g.
from
f = 3*sin(2*x^2)

to
f = 3.*sin(2.*x.^2)

A:

help vectorize

Matlab: Including function definitions when Publishing

Q:

I'm using MATLAB to generate a report from a script file, and was wondering if it's possible to somehow include function definitions which will be automatically be added to my published document. Instead of having to copy/paste between seperate published documents into a final one.

A:

help dbtype

Original post here.

Matlab: Product Support 1109 - Code Vectroization Guide

Original post here.

This section might be useful:

Suppose you want to multiply each column of a very large matrix by the same vector. There is a way to do this using sparse matrices that saves space and can also be faster than matrix construction...

F = rand(1024,1024);
X = rand(1024,1);

Y = F * diag(sparse(X));
Y = diag(sparse(X)) * F;

On the other hand, most of times bsxfun would be a more efficient solution.

Matlab: search via command line

lookfor Search all M-files for keyword.

docsearch Search HTML documentation in the Help browser.



Original post here.

Q:
Is there a general command in command window to display all conversion functions used in matlab, i.e. num2str, int2str, setstr, hex2num, ... etc?

A:
lookfor convert;
docsearch convert;
docsearch conversion;