Alternative way of coding

I'm trying to create an application that, based on a simple assumption of water flow, results in flooded area's on a DTM. The principle works, but the way it is done now the tool is extremely slow. I think there should be an alternative way to speed things up using the matrix functionalities in Matlab a bit better, but I don't know how. Any idea anyone?

dbstop if error

[FileName,PathName] = uigetfile('*.asc;*.txt','Selecteer ASCII gridfile');
fid = fopen([PathName FileName], 'r'); 
Heading = cell(6,1);

for ii = 1:6

      AsciiLine = fgetl(fid);
      Heading{ii,1} = AsciiLine;

end

ii = 1;
AsciiLine = fgetl(fid);

while ischar(AsciiLine)

      CellLine = textscan(AsciiLine,'%f');
      Grid(ii,:) = CellLine{1,1}';
      ii = ii + 1;
      AsciiLine = fgetl(fid);

end

fclose(fid);

TestGrid = Grid;
NoData = -9999;
TestGrid(logical(TestGrid==NoData)) = NaN;

[y x] = ndgrid(1:size(TestGrid,1),1:size(TestGrid,2));

CellSize = char(Heading(5,1));
CellSize = str2double(CellSize(max(regexp(Heading{5,1}, '\s'))+1:end));

Rainfall = 19.8;
InitialDepth = (ones(size(TestGrid))*(Rainfall/1000)).*logical(~isnan(TestGrid));
TotalVolume = sum(sum(InitialDepth))*(CellSize^2);
Waterlevel = TestGrid+InitialDepth;
WaterlevelTemp = TestGrid+InitialDepth;

Stabiel =0;
LoopTeller = 0;

[WetCellsRow WetCellsColumn WetCellsVector] = find(InitialDepth);

while Stabiel == 0

      LoopTeller = LoopTeller + 1;
      WaterlevelChange = zeros(size(Waterlevel));
      tic

      for ii = 1:size(WetCellsVector,1)

          if ~isnan(TestGrid(WetCellsRow(ii),WetCellsColumn(ii))) && (Waterlevel(WetCellsRow(ii),WetCellsColumn(ii))-TestGrid(WetCellsRow(ii),WetCellsColumn(ii))) > 0

              out = min(max(abs(bsxfun(@minus,x,reshape(WetCellsColumn(ii),1,1,[]))),abs(bsxfun(@minus,y,reshape(WetCellsRow(ii),1,1,[])))),[],3);
              dist = min(hypot(bsxfun(@minus,y,permute(WetCellsRow(ii),[3 2 1])),bsxfun(@minus,x,permute(WetCellsColumn(ii),[3 2 1]))),[],3);

              MinSurroundLevel = Waterlevel(logical(out==1));
              MinSurroundLevel = MinSurroundLevel(logical(((Waterlevel(out==1)-Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)))./dist(out==1))==min(((Waterlevel(out==1)-Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)))./dist(out==1)))));

              if Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)) >= MinSurroundLevel(1,1)

                  MeanLevel = (sum(Waterlevel((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)) + Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)))/(size(Waterlevel((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1),1)+1);

                  if ((Waterlevel(WetCellsRow(ii),WetCellsColumn(ii))-TestGrid(WetCellsRow(ii),WetCellsColumn(ii)))/sum(sum(sum(logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1)))))+(Waterlevel((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)) > MeanLevel;

                      WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii)) = WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii))+(MeanLevel-Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)));
                      WaterlevelChange((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1) = WaterlevelChange(((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1))+(MeanLevel-Waterlevel(((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)));

                  else

                      WaterlevelChange((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1) = WaterlevelChange((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)+((Waterlevel(WetCellsRow(ii),WetCellsColumn(ii))-TestGrid(WetCellsRow(ii),WetCellsColumn(ii)))/sum(sum(sum(logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1)))));
                      WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii)) = WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii))+(TestGrid(WetCellsRow(ii),WetCellsColumn(ii)) - Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)));

                  end

              end

          end

      end

      Waterlevel = Waterlevel + WaterlevelChange;
      Waterdepth = Waterlevel-TestGrid;
      Waterdepth(isnan(Waterdepth)) = 0;
      [WetCellsRow WetCellsColumn WetCellsVector] = find(Waterdepth);
      TimerResult(LoopTeller) = toc;

      if round(WaterlevelTemp(~isnan(Waterlevel))*1000) == round(Waterlevel(~isnan(Waterlevel))*1000)

          Stabiel = 1;
          CheckVolume = sum(sum(Waterdepth(~isnan(Waterdepth))))*(CellSize^2);

          Waterlevel(isnan(Waterlevel)) = NoData;
          Waterdepth(isnan(Waterdepth)) = NoData;

          Files = ['Waterlevel' 'Waterdepth'];

          for ii = 1:2

              fid = fopen([PathName Files(ii) datestr(now,'yyyymmdd')], 'wt');

              for jj = 1:6

                  fprintf(fid,'%s\r',char(Heading{jj,1}));

              end

              GridSize = size(Grid);
              fprintf(fid,[repmat('%g\t',1,GridSize(2)-1) '%g\n'],Grid.');

              fclose(fid);

          end

          fid =  fopen([PathName 'TimerResult'], 'wt');

          for ii = 1:size(TimerResult,1)

              fprintf(fid,'%s\r',TimerResult(ii,1));

          end        

          close(fid);

      else

          WaterlevelTemp = Waterlevel;

      end

end

Answer 1

[Jan]

I'd start with cleaning the code. Examples:

• You can omit the "logical" in logical(out==1), because "out==1" replies a logical already. But instead of creating this value again and again, it would be more efficient and much nicer to do this once only and store the rsult in a new variable.
• The repeated calculation of "~isnan(Waterlevel)" wastes time also.
• Files = ['Waterlevel' 'Waterdepth']; fid = fopen([PathName Files(ii)... will most likely not do what you expect, because "Files" is one string and "Files(ii)" is the ii.th character. You need Files = {'Waterlevel' 'Waterdepth'} with curly braces.
• It is not a good idea to ask a forum to improve lines like these, when there is not even a single comment about the meaning:

dist = min(hypot(bsxfun(@minus,y,permute(WetCellsRow(ii),[3 2 1])), ...
           bsxfun(@minus,x,permute(WetCellsColumn(ii),[3 2 1]))),[],3);

• Finally I strongly suggest to use the PROFILEr to find out, which lines consumes the most time.

Answer 2

[M.J. Glasbergen]

Thnx Jan,

Unfortunately, most of the time is lost in the amount of repetitions and not in a single line. I know it helps trying to speed up the code by speeding up the individual lines, but in this case I was more or less lokking for a more efficient apporoach. It works it's way through a grid one cell at a time until it reaches a steady state. I think there should be a way to do that part more efficient. I'll see if I can find my m-file with comments for each line (I should have one, might just have to translate the comments). I've posted this question in March already (before I created a commented version) without recieving any replies and was now just trying to see if it would help if I changed the title of the question :-)

Answer 3

[M.J. Glasbergen]

This time with comments

if true
% Code stops on error while debugging
dbstop if error

%Select file [FileName,PathName] = uigetfile('*.asc;*.txt','Select ASCII gridfile'); fid = fopen([PathName FileName], 'r');

%Blanc array for headings Heading = cell(6,1);

%Get headings from file and fill heading array for ii = 1:6

    AsciiLine = fgetl(fid);
    Heading{ii,1} = AsciiLine;

end

ii = 1; AsciiLine = fgetl(fid);

%Converting lines with data to a 'grid' while ischar(AsciiLine)

    CellLine = textscan(AsciiLine,'%f');
    Grid(ii,:) = CellLine{1,1}';
    ii = ii + 1;
    AsciiLine = fgetl(fid);

end

%Close file fclose(fid);

%Define value for NoData in 'grid' and convert to NaN TestGrid = Grid; NoData = -9999; TestGrid(logical(TestGrid==NoData)) = NaN;

%Code for the replacement of the bwdist function from the image processing toolbox. %Thnx to Andrei Bobrov via MatlabCentral: %[i1,j1] = find(yourarray) %[y x] = ndgrid(1:size(yourarray,1),1:size(yourarray,2)) %out = min(max(abs(bsxfun(@minus,x,reshape(j1,1,1,[]))),abs(bsxfun(@minus,y,reshape(i1,1,1,[])))),[],3) [y x] = ndgrid(1:size(TestGrid,1),1:size(TestGrid,2));

%Determining the grid cell size from the headings CellSize = char(Heading(5,1)); CellSize = str2double(CellSize(max(regexp(Heading{5,1}, '\s'))+1:end));

%Rainfall Rainfall = 19.8; %mm

%Waterdepth on the grid when defined rain falls on the grid InitialDepth = (ones(size(TestGrid))*(Rainfall/1000)).*logical(~isnan(TestGrid));

%Watervolume (for water balance check) TotalVolume = sum(sum(InitialDepth))*(CellSize^2);

%Waterlevel in m + AD (grid level + waterdepth) Waterlevel = TestGrid+InitialDepth; WaterlevelTemp = TestGrid+InitialDepth;

Stabiel =0; LoopTeller = 0;

%Only look at the 'wet' cells [WetCellsRow WetCellsColumn WetCellsVector] = find(InitialDepth);

%As long as there is no stable answer, run the loop while Stabiel == 0

    %Increase counter
    LoopTeller = LoopTeller + 1;

    %Blanc array for water level change
    WaterlevelChange = zeros(size(Waterlevel));

    %Start timer (for debugging only)
    tic

    %For all wet cells (each wet cell is evaluated individually):
    for ii = 1:size(WetCellsVector,1)

        %If there is a value in the current cell and there is water on the cell (waterlevel > gridlevel)
        if ~isnan(TestGrid(WetCellsRow(ii),WetCellsColumn(ii))) && (Waterlevel(WetCellsRow(ii),WetCellsColumn(ii))-TestGrid(WetCellsRow(ii),WetCellsColumn(ii))) > 0

            %Determine the distance from the current cell to the surrounding cells. (max) 8 surrounding cells het the value 1, (max) 16 around them the value 2 etc. Bepalen van de afstand van de huidige cel naar de omliggende
            out = min(max(abs(bsxfun(@minus,x,reshape(WetCellsColumn(ii),1,1,[]))),abs(bsxfun(@minus,y,reshape(WetCellsRow(ii),1,1,[])))),[],3);

            %Determine the distance from the middle of the current cell to the middle of the surrounding cells. (max) 4 directly attached cells (left, right, top and bottom) get the value 1, (max) 4 diagonally attached cells get the value SQRT(2) etc
            dist = min(hypot(bsxfun(@minus,y,permute(WetCellsRow(ii),[3 2 1])),bsxfun(@minus,x,permute(WetCellsColumn(ii),[3 2 1]))),[],3);

            %Determine the surrounding water levels and the largest hydraulic head between the current cell and the surrounding cells. Water will flow from the current cell to the one with the largest hydraulic head
            MinSurroundLevel = Waterlevel(logical(out==1));
            MinSurroundLevel = MinSurroundLevel(logical(((Waterlevel(out==1)-Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)))./dist(out==1))==min(((Waterlevel(out==1)-Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)))./dist(out==1)))));

            %Water will only flow if one of the surrounding cells has a lower value
            if Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)) >= MinSurroundLevel(1,1)

                %Determine the average water level from the current cell and the lowest surrounding cell(s)
                MeanLevel = (sum(Waterlevel((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)) + Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)))/(size(Waterlevel((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1),1)+1);

                %If the resulting water level in the surrounding cells reaches a value that is higher than the average level if the water from the current cell is added...
                if ((Waterlevel(WetCellsRow(ii),WetCellsColumn(ii))-TestGrid(WetCellsRow(ii),WetCellsColumn(ii)))/sum(sum(sum(logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1)))))+(Waterlevel((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)) > MeanLevel;

                    ...then only the amount of water resulting in the average water level will be transported from the current cell. The movement of water is put in a temporary array
                    WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii)) = WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii))+(MeanLevel-Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)));

                    ...then only the amount of water resulting in the average water level will be transported to the surrounding cell(s). The movement of water is put in a temporary array
                    WaterlevelChange((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1) = WaterlevelChange(((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1))+(MeanLevel-Waterlevel(((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)));

                %Otherwise
                else

                    ...all the water from the current cell will go to the surrounding cell. The movement of water is put in a temporary array
                    WaterlevelChange((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1) = WaterlevelChange((logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1))==1)+((Waterlevel(WetCellsRow(ii),WetCellsColumn(ii))-TestGrid(WetCellsRow(ii),WetCellsColumn(ii)))/sum(sum(sum(logical(Waterlevel == MinSurroundLevel(1,1)).*logical(out==1)))));
                    WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii)) = WaterlevelChange(WetCellsRow(ii),WetCellsColumn(ii))+(TestGrid(WetCellsRow(ii),WetCellsColumn(ii)) - Waterlevel(WetCellsRow(ii),WetCellsColumn(ii)));

                end

            end

        end

    end

    %The new levels will be calculated when all water movements have been calculated
    Waterlevel = Waterlevel + WaterlevelChange;

    %The new waterdepths are calculated. NaN values will become 0
    Waterdepth = Waterlevel-TestGrid;
    Waterdepth(isnan(Waterdepth)) = 0;

    %Determine the wet cells again
    [WetCellsRow WetCellsColumn WetCellsVector] = find(Waterdepth);

    %Stop the Timer and save the result to monitor the progress.
    TimerResult(LoopTeller) = toc;

    %If the resulting water level array is the same as the one from the previous timestep, the result is stable and new Ascii grids will be created
    if round(WaterlevelTemp(~isnan(Waterlevel))*1000) == round(Waterlevel(~isnan(Waterlevel))*1000)

        Stabiel = 1;
        CheckVolume = sum(sum(Waterdepth(~isnan(Waterdepth))))*(CellSize^2);

        Waterlevel(isnan(Waterlevel)) = NoData;
        Waterdepth(isnan(Waterdepth)) = NoData;

        Files = ['Waterlevel' 'Waterdepth'];

        for ii = 1:2

            fid = fopen([PathName Files(ii) datestr(now,'yyyymmdd')], 'wt');

            for jj = 1:6

                fprintf(fid,'%s\r',char(Heading{jj,1}));

            end

            GridSize = size(Grid);
            fprintf(fid,[repmat('%g\t',1,GridSize(2)-1) '%g\n'],Grid.');

            fclose(fid);

        end

        fid =  fopen([PathName 'TimerResult'], 'wt');

        for ii = 1:size(TimerResult,1)

            fprintf(fid,'%s\r',TimerResult(ii,1));

        end        

        close(fid);

    else

        WaterlevelTemp = Waterlevel;
%         imagesc (Waterlevel); figure(gcf)

    end

end end