Indexing a variable using the slice of another matrix

Hello, I'm struggling to wrap my head around what I think it quite a simple problem My digi-bog model is taking an awfully long time to complete (5hours), I checked the profiler and over 3 hours of that runtime was spent on this line

%% convert Q to water stage
for i = 1:numel(river_discharge_per_channel) % 4.4e6 values (very long daily time series)
    [~,index] = min(abs(Q_list-river_discharge_per_channel(i))); %Q_list changing length 2000-7000
    water_depth(i) = d(index);
end 

This step was being perfomed for every 20 interation of the main loop (12000 iterations, with 365 sub steps). I realised that this was not necassary and that I only needed to calculate the water stage for the next 20 years so I rewrote the above code to:

%% convert Q to water stage
        for i = 1:20    % not nessary to go through the entire daily discharge list, we will only do it for 20years  
            river_discharge_per_channel_year = river_discharge_per_channel(lookup_year==t-1+i); %lookup_year is the year corresponding to the daily values (I made a lookup table for this) t = the main yearly loop
            for j = 1:numel(river_discharge_per_channel_year)
                [~,index] = min(abs(Q_list-river_discharge_per_channel_year(j)));
                water_depth(river_discharge_per_channel(lookup_year==t-1+i(j))) = d(index); % this is where I'm going wrong, I want to reference the day within the yearly slice I made earlier.
            end    
        end

It's when I try to assign the water_depth values (d) back into the daily list of water_depths, I'm not managing the index it correctly. I attempt this with this line:

water_depth(river_discharge_per_channel(lookup_year==t-1+i(j))) = d(index);

I'm not sure this is possible the way I'm trying to accomplish it, suggestions most welcome. I'm not sure the issue will be understood easily by readers, not least of becuase of my dyslexic brain, so questions most welcome.

Thanks in advance

Alex

[Jan]

Is water_depth pre-allocated.

[Alexander James]

@Jan yes, its calculated for the entire time series before the main loop, however, since I'm simulating the growth of a bog within a flood plain, I have to re-estimate the water_depth due to the changing hydraulic conductivity of the flood plain due to the process that I'm modelling, this is done every twenty years.

[Matt J]

I don't see how the expression i(j) in this line

water_depth(river_discharge_per_channel(lookup_year==t-1+i(j))) = d(index);

is not producing an error. i is a scalar and cannot be indexed by j>1.

[Alexander James]

it is producing an error, it's why I created this post. I'm looking for alternative ways of accomplishing what I'm attempting in that line.

[Matt J]

Well, I think I've already given it to you in my Answer below. Once you know how to rewrite the first version of your loop, it is straightforward to rewrite your double for-loop in a similar way (however, I don't see any benefit to this):

%% convert Q to water stage
     
          [Q_list, isort]=sort(Q_list); 
           d=d(isort);
  
       for i = 1:20    % not nessary to go through the entire daily discharge list, we will only do it for 20years  
           
            index=(lookup_year==t-1+i);
 
            water_depth(index) = interp1(Q_list,d, river_discharge_per_channel(index),  'nearest');

        end

Answer 1

[Matt J]

Your original loop,

for i = 1:numel(river_discharge_per_channel) % 4.4e6 values (very long daily time series)
    [~,index] = min(abs(Q_list-river_discharge_per_channel(i))); %Q_list changing length 2000-7000
    water_depth(i) = d(index);
end 

can be replaced with the following,

[Q_list, isort]=sort(Q_list); 
d=d(isort);

water_depth = interp1(Q_list,d, river_discharge_per_channel,  'nearest');

You then showed us a modified version with a double for-loop. Adapting the above to this is straightforward,

for i = 1:20    % not nessary to go through the entire daily discharge list, we will only do it for 20years  
   
    index=(lookup_year==t-1+i);
    
    water_depth(index) = interp1(Q_list,d, river_discharge_per_channel(index),  'nearest');

end

but might be unnecessary and sub-optimal. I don't see why you wouldn't use the original version, which requires no loops.

[Alexander James]

Hi Matt, thanks for the suggestion, I don't want to replace the entire water_depth matrix with those values, just add the value's I caclulate into the correct postion within that list.

I was not aware of the interp() function though, and honestly, I don't completly follow the process you're suggesting here, wouldn't sort(Q) result in a distorted matix?

[Matt J]

wouldn't sort(Q) result in a distorted matix?

Yes, it would sort it. But nothing in your original loop depends on the order of Q_list (only that it is ordered consistent with d). If other parts of your code not shown do depend on it, you are free of course to make a back-up copy.

A simple test shows that your original loop and my proposal produce the same results,

[river_discharge_per_channel,Q_list,d, water_depth]=deal(rand(1,10));

%% convert Q to water stage
for i = 1:numel(river_discharge_per_channel) % 4.4e6 values (very long daily time series)
    [~,index] = min(abs(Q_list-river_discharge_per_channel(i))); %Q_list changing length 2000-7000
    water_depth(i) = d(index);
end 

water_depth

water_depth = 1×10

0.4557 0.8646 0.3858 0.6567 0.4446 0.8195 0.2802 0.5644 0.0991 0.1318

[Q_list, isort]=sort(Q_list); 
indices=1:numel(Q_list);
d=d(isort);

water_depth2 = interp1(Q_list,d, river_discharge_per_channel,  'nearest')

water_depth2 = 1×10

0.4557 0.8646 0.3858 0.6567 0.4446 0.8195 0.2802 0.5644 0.0991 0.1318

[Alexander James]

Thanks, I will have to try this and see if it works for my application, indeed water_depth is used by other functions in the loop and my gut feeling is this would not work. I will try it, and if it does work, I'll accept the solution you suggested, regardless, thank you for your suggestion.

[Alexander James]

Having checked the results, you are correct, this was a better, and faster way of achieving the same results. It took me a while to wrap my head around how interp1 'nearest' achieved the same results as min(abs(Qlist-river_discharge_per_channel, but now I understand. Thank you very much for taking the time to look into the problem and I am grateful for your solution and explanation.

Kind Regards

Alex

[Alexander James]

Just one final enquiry, you create a variable

indices=1:numel(Q_list);

Then you don't use this later, why include this line?

[Matt J]

You can discard that. It's a remnant of an earlier version of my answer which turned out to be unnecessary.

[Alexander James]

Great, I thought it was something like that. many thanks. FYI your suggestion helped to improve the runtime of my code from 5 hours to 45minutes.