Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
for loop hanging during interp1

Subject: for loop hanging during interp1

From: Jonathan W Smith

Date: 3 May, 2013 13:29:09

Message: 1 of 9

Hello

I'm interpolating Model ozone (Y) on a model pressure 3-dimensional grid (x_perm) in time (time is 4th dimension), to a satellite Pressure grid (xi). I rearranged the dimensions of the arrays to make vertical the first dimension. Since x has to be vector I used the squeeze function on x_perm. With x_perm, we have a vector of 25 pressure levels for each longitude/latitude grid cell at each of the 30 times. The interpolation completes but the loop hangs for an hour and hence I have to terminate it. I know that the loop completes because when I print out the new_o3 function, the values that I am expecting are displayed. How do I make the loop finish. Once I get it to finish how can I make it run faster?

x_perm = permute(wrf_p_final, [3 2 1 4]);

Y = permute(wrf_o3_daily, [3 2 1 4]); %OMI Pressure

xi = permute(omi_p, [3 2 1 4]);

for lon = 1:1:71;

       for lat = 1:1:41;

                for time_h = 1:1:30;

                        x = squeeze(x_perm(:,lon,lat,time_h));

                        new_o3 = interp1(x, Y, xi(:,lon,lat,time_h));

                end
       end
end

Subject: for loop hanging during interp1

From: Steven_Lord

Date: 3 May, 2013 17:59:07

Message: 2 of 9



"Jonathan W Smith" <jwsmith9@gmail.com> wrote in message
news:km0e35$836$1@newscl01ah.mathworks.com...
> Hello
>
> I'm interpolating Model ozone (Y) on a model pressure 3-dimensional grid
> (x_perm) in time (time is 4th dimension), to a satellite Pressure grid
> (xi). I rearranged the dimensions of the arrays to make vertical the
> first dimension. Since x has to be vector I used the squeeze function on
> x_perm. With x_perm, we have a vector of 25 pressure levels for each
> longitude/latitude grid cell at each of the 30 times. The interpolation
> completes but the loop hangs for an hour and hence I have to terminate it.
> I know that the loop completes because when I print out the new_o3
> function, the values that I am expecting are displayed. How do I make the
> loop finish. Once I get it to finish how can I make it run faster?
>
> x_perm = permute(wrf_p_final, [3 2 1 4]);
>
> Y = permute(wrf_o3_daily, [3 2 1 4]); %OMI Pressure
>
> xi = permute(omi_p, [3 2 1 4]);
>
> for lon = 1:1:71;
>
> for lat = 1:1:41;
>
> for time_h = 1:1:30;
>
> x = squeeze(x_perm(:,lon,lat,time_h));
>
> new_o3 = interp1(x, Y, xi(:,lon,lat,time_h));

You're calling INTERP1 71*41*130 = 378430 times and throwing away the
answers from all but the last call. If that's truly your intention, replace
all three loops with:

new_o3 = interp1(x, Y, xi(:,71,141,30));

If not, state the sizes of x_perm, Y, and xi and the size you want new_o3 to
be.

--
Steve Lord
slord@mathworks.com
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

Subject: for loop hanging during interp1

From: Jonathan W Smith

Date: 3 May, 2013 20:01:10

Message: 3 of 9

"Steven_Lord" <slord@mathworks.com> wrote in message <km0tta$18f$1@newscl01ah.mathworks.com>...
>
>
> "Jonathan W Smith" <jwsmith9@gmail.com> wrote in message
> news:km0e35$836$1@newscl01ah.mathworks.com...
> > Hello
> >
> > I'm interpolating Model ozone (Y) on a model pressure 3-dimensional grid
> > (x_perm) in time (time is 4th dimension), to a satellite Pressure grid
> > (xi). I rearranged the dimensions of the arrays to make vertical the
> > first dimension. Since x has to be vector I used the squeeze function on
> > x_perm. With x_perm, we have a vector of 25 pressure levels for each
> > longitude/latitude grid cell at each of the 30 times. The interpolation
> > completes but the loop hangs for an hour and hence I have to terminate it.
> > I know that the loop completes because when I print out the new_o3
> > function, the values that I am expecting are displayed. How do I make the
> > loop finish. Once I get it to finish how can I make it run faster?
> >
> > x_perm = permute(wrf_p_final, [3 2 1 4]);
> >
> > Y = permute(wrf_o3_daily, [3 2 1 4]); %OMI Pressure
> >
> > xi = permute(omi_p, [3 2 1 4]);
> >
> > for lon = 1:1:71;
> >
> > for lat = 1:1:41;
> >
> > for time_h = 1:1:30;
> >
> > x = squeeze(x_perm(:,lon,lat,time_h));
> >
> > new_o3 = interp1(x, Y, xi(:,lon,lat,time_h));
>
> You're calling INTERP1 71*41*130 = 378430 times and throwing away the
> answers from all but the last call. If that's truly your intention, replace
> all three loops with:
>
> new_o3 = interp1(x, Y, xi(:,71,141,30));
>
> If not, state the sizes of x_perm, Y, and xi and the size you want new_o3 to
> be.
>
> --
> Steve Lord
> slord@mathworks.com
> To contact Technical Support use the Contact Us link on
> http://www.mathworks.com


Steve
Thanks for your response.

The sizes of x_perm, Y, and x are 7 by 71 by 41 x 30.

I want to keep the arrays from the last call. The 71 and 41 are for longitude and latitude respectively. The interpolation has to be done for 30 71x 41 grids.
Jonathan

Subject: for loop hanging during interp1

From: dpb

Date: 3 May, 2013 20:23:42

Message: 4 of 9

On 5/3/2013 3:01 PM, Jonathan W Smith wrote:
> "Steven_Lord" <slord@mathworks.com> wrote in message
> <km0tta$18f$1@newscl01ah.mathworks.com>...
>> "Jonathan W Smith" <jwsmith9@gmail.com> wrote in message
>> news:km0e35$836$1@newscl01ah.mathworks.com...
...
...

>> >
>> > x_perm = permute(wrf_p_final, [3 2 1 4]);
>> > Y = permute(wrf_o3_daily, [3 2 1 4]); %OMI Pressure
>> > xi = permute(omi_p, [3 2 1 4]);
>> > for lon = 1:1:71;
>> > for lat = 1:1:41;
>> > for time_h = 1:1:30;
>> > x = squeeze(x_perm(:,lon,lat,time_h));
>> > new_o3 = interp1(x, Y, xi(:,lon,lat,time_h));
>>
>> You're calling INTERP1 71*41*130 = 378430 times and throwing away the
>> answers from all but the last call. If that's truly your intention,
>> replace all three loops with:
>>
>> new_o3 = interp1(x, Y, xi(:,71,141,30));
>>
>> If not, state the sizes of x_perm, Y, and xi and the size you want
>> new_o3 to be.

...

> The sizes of x_perm, Y, and x are 7 by 71 by 41 x 30.
>
> I want to keep the arrays from the last call. The 71 and 41 are for
> longitude and latitude respectively. The interpolation has to be done
> for 30 71x 41 grids.

What you didn't say is what the end result dimension of new_o3 is
intended to be.

I'd also ask if the xi are variable across the various dimensions or is
it a fixed spacing by any chance? That could cut down just a tiny
fraction of the overhead...

--

Subject: for loop hanging during interp1

From: Jonathan W Smith

Date: 3 May, 2013 20:38:08

Message: 5 of 9

dpb <none@non.net> wrote in message <km16bk$6pv$1@speranza.aioe.org>...
> On 5/3/2013 3:01 PM, Jonathan W Smith wrote:
> > "Steven_Lord" <slord@mathworks.com> wrote in message
> > <km0tta$18f$1@newscl01ah.mathworks.com>...
> >> "Jonathan W Smith" <jwsmith9@gmail.com> wrote in message
> >> news:km0e35$836$1@newscl01ah.mathworks.com...
> ...
> ...
>
> >> >
> >> > x_perm = permute(wrf_p_final, [3 2 1 4]);
> >> > Y = permute(wrf_o3_daily, [3 2 1 4]); %OMI Pressure
> >> > xi = permute(omi_p, [3 2 1 4]);
> >> > for lon = 1:1:71;
> >> > for lat = 1:1:41;
> >> > for time_h = 1:1:30;
> >> > x = squeeze(x_perm(:,lon,lat,time_h));
> >> > new_o3 = interp1(x, Y, xi(:,lon,lat,time_h));
> >>
> >> You're calling INTERP1 71*41*130 = 378430 times and throwing away the
> >> answers from all but the last call. If that's truly your intention,
> >> replace all three loops with:
> >>
> >> new_o3 = interp1(x, Y, xi(:,71,141,30));
> >>
> >> If not, state the sizes of x_perm, Y, and xi and the size you want
> >> new_o3 to be.
>
> ...
>
> > The sizes of x_perm, Y, and x are 7 by 71 by 41 x 30.
> >
> > I want to keep the arrays from the last call. The 71 and 41 are for
> > longitude and latitude respectively. The interpolation has to be done
> > for 30 71x 41 grids.
>
> What you didn't say is what the end result dimension of new_o3 is
> intended to be.
>
> I'd also ask if the xi are variable across the various dimensions or is
> it a fixed spacing by any chance? That could cut down just a tiny
> fraction of the overhead...
>
> --

new_o3 should is size 7 x 71 x 41 x 30. (7 vertical levels, 71 lon's, 41 lat's, and 30 days). And the first dimension of the array (size 7) is variable across the 71 x 41 grid cells and the 30 days.

Jonathan

Subject: for loop hanging during interp1

From: dpb

Date: 3 May, 2013 21:38:16

Message: 6 of 9

On 5/3/2013 3:38 PM, Jonathan W Smith wrote:
> dpb <none@non.net> wrote in message <km16bk$6pv$1@speranza.aioe.org>...
>> On 5/3/2013 3:01 PM, Jonathan W Smith wrote:
>> > "Steven_Lord" <slord@mathworks.com> wrote in message
>> > <km0tta$18f$1@newscl01ah.mathworks.com>...
>> >> "Jonathan W Smith" <jwsmith9@gmail.com> wrote in message
>> >> news:km0e35$836$1@newscl01ah.mathworks.com...
>> ...
>> ...
>>
>> >> >
>> >> > x_perm = permute(wrf_p_final, [3 2 1 4]);
>> >> > Y = permute(wrf_o3_daily, [3 2 1 4]); %OMI Pressure
>> >> > xi = permute(omi_p, [3 2 1 4]);
>> >> > for lon = 1:1:71;
>> >> > for lat = 1:1:41;
>> >> > for time_h = 1:1:30;
>> >> > x = squeeze(x_perm(:,lon,lat,time_h));
>> >> > new_o3 = interp1(x, Y, xi(:,lon,lat,time_h));
>> >>
>> >> You're calling INTERP1 71*41*130 = 378430 times and throwing away the
>> >> answers from all but the last call. If that's truly your intention,
>> >> replace all three loops with:
>> >>
>> >> new_o3 = interp1(x, Y, xi(:,71,141,30));
>> >>
>> >> If not, state the sizes of x_perm, Y, and xi and the size you want
>> >> new_o3 to be.
>>
>> ...
>>
>> > The sizes of x_perm, Y, and x are 7 by 71 by 41 x 30.
>> >
>> > I want to keep the arrays from the last call. The 71 and 41 are for
>> > longitude and latitude respectively. The interpolation has to be done
>> > for 30 71x 41 grids.
>>
>> What you didn't say is what the end result dimension of new_o3 is
>> intended to be.
>>
>> I'd also ask if the xi are variable across the various dimensions or
>> is it a fixed spacing by any chance? That could cut down just a tiny
>> fraction of the overhead...
>>
...

> new_o3 should is size 7 x 71 x 41 x 30. (7 vertical levels, 71 lon's, 41
> lat's, and 30 days). And the first dimension of the array (size 7) is
> variable across the 71 x 41 grid cells and the 30 days.

Then you need to preallocate it and index into it to avoid just
overwriting the same return vector on every pass.

I'm tied up on other things right now, but when you say you're sure it
quits but "hangs", have you

a) just cut the upper limits of the loops down and time it to do
sanity-checks on the code, and

b) if you interrupt it after some point in time what are the values of
lot, lan, time_h? I'm betting there not all at the end yet if what
you've shown is all the code.

--

Subject: for loop hanging during interp1

From: Jonathan W Smith

Date: 3 May, 2013 21:56:09

Message: 7 of 9

dpb <none@non.net> wrote in message <km1anc$j6g$1@speranza.aioe.org>...
> On 5/3/2013 3:38 PM, Jonathan W Smith wrote:
> > dpb <none@non.net> wrote in message <km16bk$6pv$1@speranza.aioe.org>...
> >> On 5/3/2013 3:01 PM, Jonathan W Smith wrote:
> >> > "Steven_Lord" <slord@mathworks.com> wrote in message
> >> > <km0tta$18f$1@newscl01ah.mathworks.com>...
> >> >> "Jonathan W Smith" <jwsmith9@gmail.com> wrote in message
> >> >> news:km0e35$836$1@newscl01ah.mathworks.com...
> >> ...
> >> ...
> >>
> >> >> >
> >> >> > x_perm = permute(wrf_p_final, [3 2 1 4]);
> >> >> > Y = permute(wrf_o3_daily, [3 2 1 4]); %OMI Pressure
> >> >> > xi = permute(omi_p, [3 2 1 4]);
> >> >> > for lon = 1:1:71;
> >> >> > for lat = 1:1:41;
> >> >> > for time_h = 1:1:30;
> >> >> > x = squeeze(x_perm(:,lon,lat,time_h));
> >> >> > new_o3 = interp1(x, Y, xi(:,lon,lat,time_h));
> >> >>
> >> >> You're calling INTERP1 71*41*130 = 378430 times and throwing away the
> >> >> answers from all but the last call. If that's truly your intention,
> >> >> replace all three loops with:
> >> >>
> >> >> new_o3 = interp1(x, Y, xi(:,71,141,30));
> >> >>
> >> >> If not, state the sizes of x_perm, Y, and xi and the size you want
> >> >> new_o3 to be.
> >>
> >> ...
> >>
> >> > The sizes of x_perm, Y, and x are 7 by 71 by 41 x 30.
> >> >
> >> > I want to keep the arrays from the last call. The 71 and 41 are for
> >> > longitude and latitude respectively. The interpolation has to be done
> >> > for 30 71x 41 grids.
> >>
> >> What you didn't say is what the end result dimension of new_o3 is
> >> intended to be.
> >>
> >> I'd also ask if the xi are variable across the various dimensions or
> >> is it a fixed spacing by any chance? That could cut down just a tiny
> >> fraction of the overhead...
> >>
> ...
>
> > new_o3 should is size 7 x 71 x 41 x 30. (7 vertical levels, 71 lon's, 41
> > lat's, and 30 days). And the first dimension of the array (size 7) is
> > variable across the 71 x 41 grid cells and the 30 days.
>
> Then you need to preallocate it and index into it to avoid just
> overwriting the same return vector on every pass.
>
> I'm tied up on other things right now, but when you say you're sure it
> quits but "hangs", have you
>
> a) just cut the upper limits of the loops down and time it to do
> sanity-checks on the code, and
>
> b) if you interrupt it after some point in time what are the values of
> lot, lan, time_h? I'm betting there not all at the end yet if what
> you've shown is all the code.
>
> --

dpb

How do you preallocate?

The code I present here is part of a larger script. for instance, I started the large script at about 10:30 am ET this morning. It took 5 min to reach this part of the script and its been sitting there ever since. I'll try cutting down the upper limits of the loops. Thanks for those suggestions.

Jonathan

Subject: for loop hanging during interp1

From: dpb

Date: 4 May, 2013 02:17:12

Message: 8 of 9

On 5/3/2013 4:56 PM, Jonathan W Smith wrote:
...

> How do you preallocate?

new0=zeros(7, 71, 41, 30];

Inside loop use proper indices to store the vector. I didn't have time
to try to figure out the actual directions/sizes of the return vector at
the moment, sorry.


> The code I present here is part of a larger script. for instance, I
> started the large script at about 10:30 am ET this morning. It took 5
> min to reach this part of the script and its been sitting there ever
> since. I'll try cutting down the upper limits of the loops. Thanks for
> those suggestions.

I have a feeling given the squeeze() and the sizes you're not accessing
sequential elements of the array (remember Matlab storage is in
column-major order) and cache misses are probably killing this loop in
building that vector.

Look into reorganizing the data arrays to be able to address the proper
vectors in order sequentially.

It's also possible you might be able to move up to interp2|3 and let
Matlab deal w/ it internally instead of for each individual vector as
you're doing here.

--

Subject: for loop hanging during interp1

From: Bruno Luong

Date: 4 May, 2013 08:34:10

Message: 9 of 9

"Jonathan W Smith" wrote in message <km0e35$836$1@newscl01ah.mathworks.com>...
> Hello
>
> I'm interpolating Model ozone (Y) on a model pressure 3-dimensional grid (x_perm) in time (time is 4th dimension), to a satellite Pressure grid (xi). I rearranged the dimensions of the arrays to make vertical the first dimension. Since x has to be vector I used the squeeze function on x_perm. With x_perm, we have a vector of 25 pressure levels for each longitude/latitude grid cell at each of the 30 times. The interpolation completes but the loop hangs for an hour and hence I have to terminate it. I know that the loop completes because when I print out the new_o3 function, the values that I am expecting are displayed. How do I make the loop finish. Once I get it to finish how can I make it run faster?
>
> x_perm = permute(wrf_p_final, [3 2 1 4]);
>
> Y = permute(wrf_o3_daily, [3 2 1 4]); %OMI Pressure
>
> xi = permute(omi_p, [3 2 1 4]);
>
> for lon = 1:1:71;
>
> for lat = 1:1:41;
>
> for time_h = 1:1:30;
>
> x = squeeze(x_perm(:,lon,lat,time_h));
>
> new_o3 = interp1(x, Y, xi(:,lon,lat,time_h));
>
> end
> end
> end

It seems very odd that you use the whole array Y to interpolate the O3 at one spatial-temporal point.

I bet you rather want to do something like this:

% Fake data
nx = 50;
nlon = 71;
nlat = 41;
nt = 30;
nxi = 7;
x_perm = sort(rand(nx,nlon,nlat,nt),1);
Y = rand(nx,nlon,nlat,nt);
xi = sort(rand(nxi,nlon,nlat,nt),1);

yi = zeros(nxi,nlon,nlat,nt);
for lon = 1:71
    for lat = 1:41
        for time_h = 1:30
            x = x_perm(:,lon,lat,time_h);
            y = Y(:,lon,lat,time_h); % Interpolate at 1 point
            yi(:,lon,lat,time_h) = interp1(x, y, xi(:,lon,lat,time_h));
        end
    end
end

% Bruno

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us