Matrix Indexing

I have two index files:

On times:

 t_on = [2 10 20 ...];

Off times:

 t_off = [5 15 23 ...];

I want to get a vector with ones from 2..5, 10...15, 20...23

 [0 1 1 1 1 0 0 0 0 1 1 1 1 1 1 0 ...]

I am looking for a way to do so without using loops.

[Martijn]

Why do you want to do this without a loop? With a loop it is so simple, easy and readable:

A = zeros(1,t_off(end));

for i=1:length(t_on)

A(t_on(i):t_off(i)) = 1;

end

[Joerg]

Because the vectors are large (1e6) and I assue working with loops is timeconsuming.

[Doug Hull]

I would not assume. Run the profiler.

Answer 1

[Matt Fig]

Another alternative:

A = zeros(1,t_off(end)+1);
A(t_on) = 1;
A(t_off+1) = -1;
A = cumsum(A(1:t_off(end)));

[Sean de Wolski]

This is a better answer; it'll be much faster.

[Matt Fig]

I want to point out that the above method depends on there being a zero between each run of ones. If this is not the case, use this instead:

L = length(t_on);

A = t_on(2:L)-t_off(1:L-1) > 1;

t_on = t_on([true A]);

t_off = t_off([A true]);

A = zeros(1,t_off(end));

A(t_on) = 1;

A(t_off(1:end-1)+1) = -1;

A = cumsum(A);

This doesn't add much time. Note that a well designed FOR loop is not much slower on my machine. For example, this:

R = zeros(1,t_off(end));

for ii = 1:length(t_off)

R(t_on(ii):t_off(ii)) = 1;

end

when called with this data:

t_on = 2:8:10000000;

t_off = t_on+ceil(rand(1,length(t_on))*8)+1;

runs in 1 second, whereas the vectorized solution above runs in .5 seconds. The CELLFUN method runs in 13 seconds.

Answer 2

[Matt Fig]

Bruno's MCOLON function may be what you need.

http://www.mathworks.com/matlabcentral/fileexchange/29854-multiple-colon

Also, see the lengthy discussion about this here:

http://www.mathworks.com/matlabcentral/newsreader/view_thread/299764#809227

Answer 3

[Sean de Wolski]

This will work, may not be faster than a loop.

 t_on = [2 10 20];
 t_off = [5 15 23];
 pts = cellfun(@(x)(x(1):x(2)),num2cell([t_on.',t_off.'],2),'uni',false);
 t = false(1,t_off(end));
 t([pts{:}]) = true;

Answer 4

[James Tursa]

As an exercise in speed I wrote a mex routine to do this. It is faster then the previously posted m-code, but not by a huge amount. Here it is:

/* Program: onesx
 * Programmer: James Tursa
 *
 * Syntax:  C = onesx( A , B [,largest_index] )
 *              A = real double vector of indexes to turn on 1's
 *              B = real double vector of indexes to turn off 1's (end inclusive)
 *              largest_index = optional largest index for result
 *
 * Example: C = onesx( [3 7 11] , [4 9 14] )
 *          C =  0   0   1   1   0   0   1   1   1   0   1   1   1   1
 *
 *
 * Example: C = onesx( [3 7 11] , [4 9 14] , 16 )
 *          C =  0   0   1   1   0   0   1   1   1   0   1   1   1   1  0  0
 * 
 */

#include "mex.h"

#ifndef MWSIZE_MAX
#define  mwIndex        int
#define  mwSignedIndex  int
#define  mwSize         int
#endif

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    mwSize i, k, k1, k2, n1, n2, largest_index;
  double *pr1, *pr2, *ones;
  double x, z;

  if( nrhs < 2 || nrhs > 3 ) {
    mexErrMsgTxt("Expected two double vector inputs (and optional third max index).");
  }
  if( nlhs > 1 ) {
    mexErrMsgTxt("Too many outputs.");
  }
  n1 = mxGetNumberOfElements(prhs[0]);
  n2 = mxGetNumberOfElements(prhs[1]);
  if( n1 == 0 || !mxIsDouble(prhs[0]) || mxGetNumberOfDimensions(prhs[0]) != 2 ||
    mxIsComplex(prhs[0]) || (mxGetM(prhs[0]) != 1 && mxGetN(prhs[0]) != 1) ||
      n2 == 0 || !mxIsDouble(prhs[1]) || mxGetNumberOfDimensions(prhs[1]) != 2 ||
    mxIsComplex(prhs[1]) || (mxGetM(prhs[1]) != 1 && mxGetN(prhs[1]) != 1) ) {
    mexErrMsgTxt("Inputs must be non-empty real double vectors.");
  }
  if( n1 != n2 ) {
    mexErrMsgTxt("Inputs must have the same number of elements.");
  }
  if( nrhs == 3 ) {
    largest_index = z = mxGetScalar(prhs[2]);
    if( largest_index != z || largest_index < 0 ) {
      mexErrMsgTxt("3rd argument largest index must be positive integral value and not too large.");
    }
  }
  pr1 = mxGetPr(prhs[0]);
  pr2 = mxGetPr(prhs[1]);
  x = 0.0;
  for( i=0; i<n1; i++ ) {
    if( ((mwSize) pr1[i]) != pr1[i] ) {
      mexPrintf("1st(%d) = %f\n",i+1,pr1[i]);
      mexErrMsgTxt("Invalid input element. Must be positive integral value and not too large.");
    }
    if( ((mwSize) pr2[i]) != pr2[i] ) {
      mexPrintf("2nd(%d) = %f\n",i+1,pr2[i]);
      mexErrMsgTxt("Invalid input element. Must be positive integral value and not too large.");
    }
    if( pr1[i] > pr2[i] || pr1[i] <= 0 || pr2[i] <= 0 ) {
      mexPrintf("1st(%d) = %d , 2nd(%d) = %d\n",i+1,(mwSize)pr1[i],i+1,(mwSize)pr2[i]);
      mexErrMsgTxt("Invalid input element pair. Must have 0 < 1st <= 2nd.");
    }
    if( pr2[i] > x ) x = pr2[i];
    if( nrhs == 3 && x > z ) {
      mexPrintf("2nd(%d) = %d > 3rd = %d\n",i+1,(mwSize)x,largest_index);
      mexErrMsgTxt("Input element is larger than 3rd argument largest index.");
    }
  }
  if( nrhs == 3 ) x = z;
  plhs[0] = mxCreateDoubleMatrix( 1, x, mxREAL);
  ones = mxGetPr(plhs[0]);
  for( i=0; i<n1; i++ ) {
    k1 = ((mwSize) pr1[i]) - 1;
    k2 = ((mwSize) pr2[i]) - 1;
    for( k=k1; k<=k2; k++ ) {
      ones[k] = 1.0;
    }
  }
}

[James Tursa]

GOOD GRIEF! Apologies for the code formatting above. This was a simple copy & paste from the editor. I think I need to just post this stuff on the regular newsgroup or FEX and direct people there ... posting code here is useless!

[Ned Gulley]

Hi James:

If you add a space to the beginning of each line, it will be formatted as monospace text. Another way to do this is to select all of the code and then select the "code" button in the toolbar above the text input. I added the necessary spaces to your answer to make it look better.

[James Tursa]

Aha! Thanks ... (I'm still getting used to this Answers section)