Alternate search functions to speed up code

I have tried profiling my code and apparently it is very slow to the use of the desarchn algorithm.

The whole program intital takes around 400 seconds to run with this one function shown below being the bottle neck taking 350 seconds. In particular, the dsearchn function takes a very long time. What can I do to make it run faster?

Other things I have tried

1. I tried implementing the desarchn function but, the code took signficiantly longer to run (even) 1000 seconds the function had to finish exectuing)
2. I tried using parfor loops but, MATLAB gives me an error saying that the code cannot use parfor due to the loop structures.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function connectivity_coords = mapcoordinates(num_voxels, sz, vox_x_dim, vox_y_dim, vox_z_dim, node_list_matrix,index, corner_coords, counter, connectivity_coords)
for ind = 1:num_voxels
    % Convert voxel index to subscripts
    [x, y, z] = ind2sub(sz, ind);
    %reset the origin to zero 
    x= (x-1)*vox_x_dim;
    y =(y-1)*vox_y_dim;
    z=(z-1)*vox_z_dim;
    % Calculate the coordinates of the cube corners
    corners = [
        x, y, z %4;
        x+vox_x_dim, y, z %3;
        x+vox_x_dim, y+vox_y_dim, z %2
        x, y+vox_y_dim, z %1; 
        x, y, z+vox_z_dim %8
        x+vox_x_dim, y, z+vox_z_dim %7; 
        x+vox_x_dim, y+vox_y_dim, z+vox_z_dim %6;
        x, y+vox_y_dim, z+vox_z_dim %5;
        ];
    smallindex = dsearchn(node_list_matrix(:,2:4), corners(1:end,1:3));
    index = [index;smallindex];
    start = ind*8-7;
    % Assign the element set 
    connectivity  = [
        index(start);
        index(start+1);
        index(start+2);
        index(start+3);
        index(start+4);
        index(start+5);
        index(start+6);
        index(start+7); 
        ];
    % Store the corner coordinates in the array
    corner_coords((ind-1)*8 + 1 : ind*8, :) = corners;
    % Store the connectivity information
    if counter<= num_voxels
        connectivity_coords(counter,1:8) = connectivity;
    end 
    counter=counter+1;
end
end

[Les Beckham]

If you provide all of the variable values and the command that you used to call the function so that people can test, you will be more likely to get an answer.

load('variables.mat');

whos

  Name                          Size               Bytes  Class     Attributes

  ans                           1x38                  76  char                
  cmdout                        1x33                  66  char                
  connectivity_coords       10584x8               677376  double              
  corner_coords                 8x3                  192  double              
  counter                       1x1                    8  double              
  no_vox_x                      1x1                    8  double              
  no_vox_y                      1x1                    8  double              
  no_vox_z                      1x1                    8  double              
  node_list_matrix         551905x4             17660960  double              
  sz                            1x2                   16  double              

coords = mapcoordinates(num_voxels, sz, vox_x_dim, vox_y_dim, vox_z_dim, node_list_matrix,index, corner_coords, counter, connectivity_coords);

Unrecognized function or variable 'num_voxels'.

function connectivity_coords = mapcoordinates(num_voxels, sz, vox_x_dim, vox_y_dim, vox_z_dim, node_list_matrix,index, corner_coords, counter, connectivity_coords)
for ind = 1:num_voxels
    % Convert voxel index to subscripts
    [x, y, z] = ind2sub(sz, ind);
    %reset the origin to zero 
    x= (x-1)*vox_x_dim;
    y =(y-1)*vox_y_dim;
    z=(z-1)*vox_z_dim;
    % Calculate the coordinates of the cube corners
    corners = [
        x, y, z %4;
        x+vox_x_dim, y, z %3;
        x+vox_x_dim, y+vox_y_dim, z %2
        x, y+vox_y_dim, z %1; 
        x, y, z+vox_z_dim %8
        x+vox_x_dim, y, z+vox_z_dim %7; 
        x+vox_x_dim, y+vox_y_dim, z+vox_z_dim %6;
        x, y+vox_y_dim, z+vox_z_dim %5;
        ];
    smallindex = dsearchn(node_list_matrix(:,2:4), corners(1:end,1:3));
    index = [index;smallindex];
    start = ind*8-7;
    % Assign the element set 
    connectivity  = [
        index(start);
        index(start+1);
        index(start+2);
        index(start+3);
        index(start+4);
        index(start+5);
        index(start+6);
        index(start+7); 
        ];
    % Store the corner coordinates in the array
    corner_coords((ind-1)*8 + 1 : ind*8, :) = corners;
    % Store the connectivity information
    if counter<= num_voxels
        connectivity_coords(counter,1:8) = connectivity;
    end 
    counter=counter+1;
end
end

Answer 1

[Bruno Luong]

The below code takes

% Elapsed time is 48.221007 seconds.

on my PC.

Improvement are:

• use delaunayTriangulation instead of dsearchn
• vectorize the for loop in order to call nearestNeighbor only once
• agregate vortex corners in case there are duplication

load('variables.mat')

sz = [floor(no_vox_x), floor(no_vox_y), floor(no_vox_z)];
num_voxels = prod(sz);

% Fix missing parameters not given by OP
connectivity_coords = [];
dxyz=max(node_list_matrix(:,2:4),[],1)./sz;
vox_x_dim = dxyz(1);
vox_y_dim = dxyz(2);
vox_z_dim = dxyz(3);

% Code starts here
ind = 1:num_voxels;
[x, y, z] = ind2sub(sz, ind);

x = (x-1)*vox_x_dim;
y = (y-1)*vox_y_dim;
z = (z-1)*vox_z_dim;

[xc,yc,zc] = ndgrid([0 vox_x_dim], [0 vox_y_dim], [0 vox_z_dim]);
cube = [xc(:) yc(:) zc(:)];
cube = cube([1 2 4 3 5 6 8 7],:);

tic
cube  = reshape(cube, [8 1 3]);
xyz   = reshape([x(:) y(:) z(:)], 1, [], 3);
corners = xyz + cube;
corners = reshape(corners, [], 3);
[ucorners,~,J] = uniquetol(corners, 'ByRows', 1);
DT = delaunayTriangulation(node_list_matrix(:,2:4));
ID = nearestNeighbor(DT,ucorners);
smallindex = ID(J);
smallindex = reshape(smallindex, [8 num_voxels]);
toc

%check matching of the 1000th data
smallindex1000 = dsearchn(node_list_matrix(:,2:4), corners(999*8+(1:8),:));
smallindex1000
smallindex(:,1000)
isequal(smallindex1000, smallindex(:,1000))

index = smallindex.';
h = size(connectivity_coords,1);
connectivity_coords = [connectivity_coords; index(1:min(num_voxels-h,end),:)];

[Maggie Chong]

Hi I have a follow up question, I was wondering if there is a way to use parallel computing with the delaunayTriangulation(node_list_matrix(:,2:4))?

My code is currently running but it requires nearly 98% of my CPU resources (on a single node in a high performance computing system).

[Bruno Luong]

IMO there is no obvious way to run delaunayTriangulation in parallel.

But you can ask this question to another thread so that other person can see. if they know how to accelerate it