Alternate search functions to speed up code

Question

Maggie Chong 2023-9-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2022757-alternate-search-functions-to-speed-up-code

评论： Bruno Luong 2023-10-9

采纳的回答： Bruno Luong

variables.mat

在 MATLAB Online 中打开

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

I tried implementing the desarchn function but, the code took signficiantly longer to run (even) 1000 seconds the function had to finish exectuing)
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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Les Beckham 2023-9-18

在 MATLAB Online 中打开

variables.mat

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 2023-9-19

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2022757-alternate-search-functions-to-speed-up-code#answer_1313207

编辑：Bruno Luong 2023-9-19

在 MATLAB Online 中打开

variables.mat

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),:)];

2 个评论
显示无隐藏无

Maggie Chong 2023-10-9

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 2023-10-9

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

请先登录，再进行评论。

Alternate search functions to speed up code

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

2 个评论
显示无隐藏无

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Alternate search functions to speed up code

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

采纳的回答

2 个评论 显示 无隐藏 无

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

2 个评论
显示无隐藏无