ISO a faster way to process 150 years of daily data

From: Bridget Thrasher <bthrasher_at_nyahnyahspammersnyahnyah>
Date: Tue Feb 23 2010 - 13:02:54 MST

I've got a script that is on track to take over 10 weeks(!!!) to complete,
but I can't figure out how to eliminate any of the loops to speed up the
process. The bottleneck portion is as follows. I have 150 years worth of
daily data. For each day, I need to compare the field to those within 45
days over a different data set of 50 years. For example, if I'm processing
Feb 23, 2010, I need to compare that day to Jan 9 through Apr 9 in each of
the 50 years 1950-1999 by calculating the RMSD between the day of interest
and each of the 4550 (91*50) comparison days. I then need to store the
comparison fields and dates that correspond to the 30 smallest RMSDs. So
I've got a loop of 91 iterations inside a loop of 50 iterations inside a
loop of 365 iterations inside a loop of 150 iterations. Crazy! Any ideas for
speeding this up would be greatly appreciated. See loop below.

nyrs = 150
ndays = 365
naggobsyrs = 50
ns = 45
; Some code here to read files and create arrays
do yr = 0,nyrs-1

    do d = 0,ndays-1
        zanalogs(0,:,:) = 1.
        day_data_20c =
bcfilesall[0:(nyrs20c/nyrsperf)-1]->bc_data(d::ndays,:,:)
        gclimo = dim_avg_n(day_data_20c,0)
        delete(day_data_20c)
        day_data_all = bcfilesall[:]->bc_data(d::ndays,:,:)
        day_data = day_data_all(yr,:,:) - gclimo
        day_data@_FillValue = fill
        delete(day_data_all)
        cna = fspan(10000,10000+na-1,na)

        do obs = 0,naggobsyrs-1
            curryr = aggobsyrsall[obs]->$var$ - aggobsdaily
            if (obs.gt.0) then prevyr = aggobsyrsall[obs-1]->$var$ -
aggobsdaily end if
            if (obs.lt.naggobsyrs-1) then nextyr =
aggobsyrsall[obs+1]->$var$ - aggobsdaily end if

            do buff = d-ns,d+ns
                if (buff.ge.0.and.buff.lt.ndays) then
                    reanal = curryr(buff,:,:)
                else
                    if (buff.lt.0) then
                        if (obs.gt.0) then
                            reanal = prevyr(buff+ndays,:,:)
                        else
                            continue
                        end if
                    else ;if (buff.ge.ndays)
                        if (obs.lt.naggobsyrs-1) then
                            reanal = nextyr(buff-ndays,:,:)
                        else
                            continue
                        end if
                    end if
                end if
                reanal@_FillValue = fill
                if (.not.all(ismissing(reanal))) then
                    c = dim_rmsd_n(reanal,day_data,(/0,1/))
                    if (c.lt.max(cna)) then
                        ii = ind(cna.eq.max(cna))
                        cna(ii(0)) = c
                        zanalogs(ii(0)+1,:,:) = reanal
                        if (buff.lt.0) then
                            ana_list(ii(0),0) = obs-1
                            ana_list(ii(0),1) = buff + ndays
                        else
                            if (buff.ge.ndays) then
                                ana_list(ii(0),0) = obs + 1
                                ana_list(ii(0),1) = buff-ndays
                            else
                                ana_list(ii(0),0) = obs
                                ana_list(ii(0),1) = buff
                            end if
                        end if
                        delete(ii)
                    end if
                    delete(c)
                end if
            end do ; end buff loop (+/- 45 days)

        end do ; end obs loop (50 years of obs)

; Some code here to do additional, relatively fast, computations

    end do ; end d loop (365 days)

; Some more code here to write out yearly files

end do ; end yr loop (150 years)

-- 
Bridget Thrasher, PhD
Postdoctoral Researcher
Climate Central
www.climatecentral.org

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Feb 23 13:03:01 2010

This archive was generated by hypermail 2.1.8 : Fri Mar 12 2010 - 09:11:56 MST