Loop matrix

Share icon

Hi everyone,

I have written a procedure where the waves are sorted by time of day and then the statistics is calculated for each hour. There are several waves and I would like to do a loop in the Function CalcConcStats_BC_Test(sortedTimeOfDay, sortedBC1, sortedBC2, averagingIntervalInMinutes) to avoid writing code for each wave (I mean: BC1Stats, BC2Stats,...).
I have tried to do a wavelist, but it does not work.


Thanks in advance, I would really appreciate any help you could provide!

Please find attached a .txt file with some data (Prueba_Stats_BC.txt).

Here is the procedure:

#pragma rtGlobals=3     // Use modern global access method and strict wave access.


//------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
//-----------------------------------------------------------------------------------------  SELECT BC waves to be sorted  -----------------------------------------------------------------------------------------
//------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Function WavesSorted_BC_Test()
    String date_utc_time        //Name of the Date&Time wave that we want to analyze
    string BC1
    string BC2
    prompt date_utc_time, "Date&Time wave to be sorted", popup WaveList("*", ";", "")
    prompt BC1, "BC1 wave to be sorted", popup WaveList("*", ";", "")
    prompt BC2, "BC2 wave to be sorted", popup WaveList("*", ";", "")
    DoPrompt "Waves sorted", date_utc_time, BC1, BC2
    SortedByTimeOfDay_BC_Test($date_utc_time, $BC1, $BC2)
End


//--------------------------------------------------------------------------------------------------------------------------------------------------------------------------
//--------------------  Function to sort input Date&Time wave by the time of day. NaNs (blanks) go to the end  --------------------
//--------------------------------------------------------------------------------------------------------------------------------------------------------------------------

Function/S SortedByTimeOfDay_BC_Test(date_utc_time, BC1, BC2)
    wave date_utc_time                          //Date&Time wave to be sorted
    wave BC1, BC2       //BC waves to be sorted
    variable i
   
    duplicate/O date_utc_time, TimeOfDay_Sorted
    duplicate/O BC1,  BC1_Sorted
    duplicate/O BC2,  BC2_Sorted
   
    TimeOfDay_Sorted = mod(date_utc_time, 24*60*60)     // TimeOfDay now contains just time of day (remove date)
    Sort {TimeOfDay_Sorted,date_utc_time}, TimeOfDay_Sorted, BC1_Sorted, BC2_Sorted     // Sort by TimeOfDay, then by date


    CalcConcStats_BC_Test(TimeOfDay_Sorted, BC1_Sorted, BC2_Sorted, 60)
   
End


//----------------------------------------------------------------------------------------------------------------
//--------------------  Function to estimate the statistics of the waves  --------------------
//----------------------------------------------------------------------------------------------------------------

Function CalcConcStats_BC_Test(sortedTimeOfDay, sortedBC1, sortedBC2, averagingIntervalInMinutes)
    wave sortedTimeOfDay        // Time-of-day wave
    wave sortedBC1, sortedBC2       // Concentrations wave
    variable averagingIntervalInMinutes     // Interval over which to average in minutes; e.g., 60 for 1 hour
   
    variable numBins = (24 * 60) / averagingIntervalInMinutes
   
    variable averagingIntervalInSeconds = averagingIntervalInMinutes * 60
   
    //Number of columns of Statistical tables
    variable numStatsColumns = 12
   
    // Variables to find the first and last point of each bin
    variable numInputPoints = numpnts(sortedTimeOfDay)
    variable i, bin, timeOfDay
   
    Variable BC1, BC2
   
    // To find the median and standard deviation for each bin
    variable first
    variable last
    variable count
    variable median

    //-------------------------------------------------------------------------------------------//
    //--------------------  Make output waves of BC waves  --------------------//
    //-------------------------------------------------------------------------------------------//
   
    //--------------------------------------------------
    //--------------------  BC1  --------------------
    //--------------------------------------------------
   
    Make /O /D /N=(numBins,numStatsColumns) BC1Stats = 0        // 2D wave
    SetDimLabel 1, 0, Count, BC1Stats           // Column 0 is count of concentrations in bin
    SetDimLabel 1, 1, First, BC1Stats           // Column 1 is first input point in bin
    SetDimLabel 1, 2, Last, BC1Stats            // Column 2 is last input point in bin
    SetDimLabel 1, 3, Min, BC1Stats         // Column 3 is median concentration for bin
    SetDimLabel 1, 4, Max, BC1Stats         // Column 4 is median concentration for bin
    SetDimLabel 1, 5, Mean, BC1Stats            // Column 5 is mean concentration for bin
    SetDimLabel 1, 6, StdDev, BC1Stats      // Column 6 is standard deviation for bin
    SetDimLabel 1, 7, Median, BC1Stats      // Column 7 is median concentration for bin
    SetDimLabel 1, 8, Q25, BC1Stats     // Column 8 is Quartil25 for bin
    SetDimLabel 1, 9, Q75, BC1Stats     // Column 9 is Quartil75 for bin
    SetDimLabel 1, 10, upperInnerFence1, BC1Stats       // Column 10 is  for bin
    SetDimLabel 1, 11, lowerInnerFence1, BC1Stats       // Column 11 is  for bin
   
   
    // Set X scaling such that the statistics for each bin fall in the middle of bin when graphed
    SetScale/P x, averagingIntervalInSeconds-3600, averagingIntervalInSeconds, "dat", BC1Stats
   
    // Preset columns of output
    BC1Stats[][%first] = -1
    BC1Stats[][%last] = -1

    // Find the first and last point of each bin, and count
    for (i=0; i<numInputPoints; i+=1)           // Bin and sum up the concentrations
        timeOfDay = sortedTimeOfDay[i]
        bin = trunc(timeOfDay / averagingIntervalInSeconds)
        if (BC1Stats[bin][%first] == -1)
            BC1Stats[bin][%first] = i
        endif
        BC1Stats[bin][%last] = i
        if (numtype(BC1) == 0)              // Exclude NaNs and INFs
            BC1Stats[bin][%count] += 1          // Increment bin count
        endif
        count = BC1Stats[bin][%count]
    endfor 
   
    // Find the minimum for each bin
    BC1Stats[][%min] = WaveMin(sortedBC1, BC1Stats[p][%first], BC1Stats[p][%last])

    // Find the maximum for each bin
    BC1Stats[][%max] = WaveMax(sortedBC1, BC1Stats[p][%first], BC1Stats[p][%last])

    // Find the mean, median and standard deviation for each bin
    for(bin=0; bin<numBins; bin+=1)
        first = BC1Stats[bin][%first]       // First input point for this bin
        last = BC1Stats[bin][%last]     // Last input point for this bin
        Duplicate /R=[first,last] /FREE sortedBC1, temp1        // Needed because StatsMedian works on the whole wave

        WaveTransform zapNaNs, temp1        // Needed because StatsMedian does not support NaNs
       
        median = StatsMedian(temp1)
        BC1Stats[bin][%median] = median

//      WaveStats /M=2 /Q temp1
        WaveStats/Q temp1
        BC1Stats[bin][%mean] = V_avg
        BC1Stats[bin][%stddev] = V_sdev

        StatsQuantiles/Q temp1
        BC1Stats[bin][%Q25] = V_Q25
        BC1Stats[bin][%Q75] = V_Q75
       
        Variable H1=V_Q75 - V_Q25
        Variable step1=1.5*H1
        Variable upperInnerFence1=V_Q75 + step1
        Variable upperOuterFence1=V_Q75 + 2*step1
        Variable lowerInnerFence1=V_Q25 - step1
        Variable lowerOuterFence1=V_Q25 - 2*step1
       
        BC1Stats[bin][%upperInnerFence1] = upperInnerFence1
        BC1Stats[bin][%lowerInnerFence1] = lowerInnerFence1

    endfor
   

    //--------------------------------------------------
    //--------------------  BC2  --------------------
    //--------------------------------------------------
   
    Make /O /D /N=(numBins,numStatsColumns) BC2Stats = 0        // 2D wave
    SetDimLabel 1, 0, Count, BC2Stats           // Column 0 is count of concentrations in bin
    SetDimLabel 1, 1, First, BC2Stats           // Column 1 is first input point in bin
    SetDimLabel 1, 2, Last, BC2Stats            // Column 2 is last input point in bin
    SetDimLabel 1, 3, Min, BC2Stats         // Column 3 is median concentration for bin
    SetDimLabel 1, 4, Max, BC2Stats         // Column 4 is median concentration for bin
    SetDimLabel 1, 5, Mean, BC2Stats            // Column 5 is mean concentration for bin
    SetDimLabel 1, 6, StdDev, BC2Stats      // Column 6 is standard deviation for bin
    SetDimLabel 1, 7, Median, BC2Stats      // Column 7 is median concentration for bin
    SetDimLabel 1, 8, Q25, BC2Stats     // Column 8 is Quartil25 for bin
    SetDimLabel 1, 9, Q75, BC2Stats     // Column 9 is Quartil75 for bin
    SetDimLabel 1, 10, upperInnerFence2, BC2Stats       // Column 10 is  for bin
    SetDimLabel 1, 11, lowerInnerFence2, BC2Stats       // Column 10 is  for bin
   
    // Set X scaling such that the statistics for each bin fall in the middle of bin when graphed
    SetScale/P x, averagingIntervalInSeconds-3600, averagingIntervalInSeconds, "dat", BC2Stats
   
    // Preset columns of output
    BC2Stats[][%first] = -1
    BC2Stats[][%last] = -1

    // Find the first and last point of each bin
    for (i=0; i<numInputPoints; i+=1)           // Bin and sum up the concentrations
        timeOfDay = sortedTimeOfDay[i]
        bin = trunc(timeOfDay / averagingIntervalInSeconds)
        if (BC2Stats[bin][%first] == -1)
            BC2Stats[bin][%first] = i
        endif
        BC2Stats[bin][%last] = i
        if (numtype(BC2) == 0)              // Exclude NaNs and INFs
            BC2Stats[bin][%count] += 1          // Increment bin count
        endif
        count = BC2Stats[bin][%count]
    endfor 
   
    // Find the minimum for each bin
    BC2Stats[][%min] = WaveMin(sortedBC2, BC2Stats[p][%first], BC2Stats[p][%last])

    // Find the maximum for each bin
    BC2Stats[][%max] = WaveMax(sortedBC2, BC2Stats[p][%first], BC2Stats[p][%last])

    // Find the median and standard deviation for each bin
    for(bin=0; bin<numBins; bin+=1)
        first = BC2Stats[bin][%first]       // First input point for this bin
        last = BC2Stats[bin][%last]     // Last input point for this bin
        Duplicate /R=[first,last] /FREE sortedBC2, temp2        // Needed because StatsMedian works on the whole wave

        WaveTransform zapNaNs, temp2        // Needed because StatsMedian does not support NaNs
       
        median = StatsMedian(temp2)
        BC2Stats[bin][%median] = median

//      WaveStats /M=2 /Q temp2
        WaveStats/Q temp2
        BC2Stats[bin][%mean] = V_avg
        BC2Stats[bin][%stddev] = V_sdev
       
        StatsQuantiles/Q temp1
        BC2Stats[bin][%Q25] = V_Q25
        BC2Stats[bin][%Q75] = V_Q75
       
        Variable H2=V_Q75 - V_Q25
        Variable step2=1.5*H2
        Variable upperInnerFence2=V_Q75 + step2
        Variable upperOuterFence2=V_Q75 + 2*step2
        Variable lowerInnerFence2=V_Q25 - step2
        Variable lowerOuterFence2=V_Q25 - 2*step2
       
        BC2Stats[bin][%upperInnerFence2] = upperInnerFence2
        BC2Stats[bin][%lowerInnerFence2] = lowerInnerFence2

    endfor 
   
End
 Prueba_Stats_BC.txt (5.67 KB)
So you want to avoid writing the same code over and over again? Is there a specific meaning for CalcConcStats_BC_Test to take two waves BC1 and BC2 (I guess not from the code)? If not just rewrite the code to use dynamic wave names instead of hard-coded ones (i.e., BC1Stats etc. ... never a good idea but for to simplest cases), and then use a function to loop though a wave list repeatedly calling your stats function in the process. Here this should give you some idea (note the use of the '$', and you should probably read more about it by executing DisplayHelpTopic "Converting a String into a Reference Using $"):
Function DynamicExampe(inwave)
    wave inwave
    String statsname = nameofwave(inwave)+"_stats"
    Make/D/O/N=(round(DimSize(inwave,0)/2)) $statsname
    SetScale/P x,DimOffset(inwave,0),DimDelta(inwave,0)*2, $statsname
    Wave workwave = $statsname
    workwave = inwave[p*2]
End


Btw. this function downsamples the incoming data to halve the points in the process.

July 2, 2015 at 04:04 am - Permalink

Just to clear things up: Of course the function has no real world application besides demonstrating the use of '$'. I thought it might be nice to have at least some effect beyond creating an empty wave. Btw I use Interpolate2 for 'resampling', because I rather dial in the output points.

July 2, 2015 at 05:51 am - Permalink

Thank you chozo and thomas_braun for your help.

The waves come from an Aethalometer instrument, so I have seven black carbon waves (BC1, BC2, ..., BC7). The function CalcConcStats_BC_Test() calculates the statistics of each wave.
I am going to try what you wrote me, I hope it works!

Thanks again!

July 3, 2015 at 12:19 am - Permalink