Rewrite readWDM.py to read by data group & block

[aufdenkampe]

From Jack in Feb 5 email:

I have been investigating bugs related to WDM data with variable timesteps. My review of the readWDM code that loops through WDM data Groups (after comment #Get timeseries data) with calls to 'getFloats' has led me to believe that 'getFloats' should be replaced by a new routine called something like 'getTimeseriesDataGroup'.

The new routine would do the following:

1. check to see that the Group date word matches the expected date.
2. loop through WDM Blocks within Group

• if Block timeunit/timestep does not match HDF5 timeunit/timestep then adjust number of values in Block for use filling array to be written to HDF5. Then process into HDF5 data array like compressed WDM data
• if Block is compressed, add multiple identical values to HDF5 data array.
• if Block is uncompressed, march through WDM record copying data values into HDF5 data array.

There need to be robust optional debug messages written to a debug file as there can be subtle issues related to midnight, new year, leap year and end of WDM records.

Code from atcWDMvb routine TimeseriesDataFromWdm can be used as an example of how the new routine might work.

Paul D. thinks: "This vb code seems to handle your complex WDM data correctly, so if we could just translate this vb to python we should be golden."

Last, @ptomasula created this EFDC fortran binary reader in Python, which could be worth following as an example designed to efficiently create Pandas data frames. https://github.com/LimnoTech/LTPy/blob/master/ltpy/ltpy/model_tools/efdc/binary_reader.py

[aufdenkampe]

Connected to #20 for info, #22 for testing, 23 for reading into Pandas data frames and writing outputs to Parquet.

[aufdenkampe]

Update: LimnoTech & Respec met Feb. 10 and outlined to rewrite the getfloats function in readWDM.py by mimicking the TimeseriesDataFromWdm in the BASINS atcWDMfile.vb script, which does not rely on a Fortran DDL file.

@jlkittle pointed out that this WDM debugger tool would be very helpful: https://github.com/respec/FORTRAN/blob/master/lib3.0/BIN/WDMRX.EXE.

@PaulDudaRESPEC shared the scanned "WDMProgrammersGuide.pdf", which I'll post somewhere in this repo.

[ptomasula]

Background

An initial rewrite of getfloats was implemented as the new function process_groups. The process_groups function is based off the logic in TimeseriesDataFromWdm from RESPEC's BASINS code, which processes timeseries as groups and blocks. Each group contains one or more blocks, which are structured data records. Each block begins with a control word (a 32bit integer) that contains information about the block's structure, and then continues with one more 32bit float values based on the number of values parameters specified in the block control word. The logic was implemented in commit c5b2a1cdd6a55eccc0db67d7840ec3eaf904dcec and was able to successfully read and convert two testing WDM files: rpo772.wdm and the test10b file TEST.WDM. However additional testing with RPO_SWMM48LINKS2017_wCBOD_June2019.wdm revealed two issues discussed further below.

Numpy Array IndexError

The irregular timestep means that the exact dimensions required for a numpy array to hold the converted timeseries could not precisely determined prior to processing the timeseries groups. Allocating a numpy array requires a fixed size, and misjudging that size results in an IndexError (which also isn't handled by the current code) and the file processing to fail. The current solution, which was committed as bad3771d4d6c06e0b2aefbd322273c990a025021, allocates the numpy array in chucks and expands the array when processing the next block would store values that exceed the current boundary of the array. However this solution is not ideal. Arrays are not really extendable. Reallocating the array, even in chunks, means copy the previous array to a larger array in a new block of memory. This drastically impacts performance and the impact of this approach will be felt even more with larger files.

Next Steps:

1. Check with @PaulDudaRESPEC and @jlkittle to see if there is a method to determine the number of data values in a timeseries prior to processing the groups.
2. If it is not possible to determine number of values prior to processing, consider replacing data_array & value_array with a more suitable data structure (e.g. python List)

Blocks Ending on 511th Element of a Record

Timeseries 134 of RPO_SWMM48LINKS2017_wCBOD_June2019.wdm contains a block which end on the 511th element of a record. The previous logic tried to parse the 512th element of that record as a block control word, however the 512th element was not a valid block control word. This would cause issues in processing the remaining blocks in the timeseries. Working on the assumption that blocks are always at least two words long I implement logic in commit 129a228922539c524fcfe7c90ec65b9b19291091, which moves to the next record in the timeseries chain of records when a block ends either on 511th or 512th element of the record.

Next Steps:

1. Check with @PaulDudaRESPEC and @jlkittle about the significance of the 512th element when a block end on the 511th element of a record.
2. Modify the code if the above assumption proved to be inaccurate.

[aufdenkampe]

Looping in @bcous for notifications.

[htaolimno]

I've made some good progress based on Paul's previous work and VB code. I added an alternative implementation of process_groups() function (commit 0f5a2d1). The two major changes I made from Paul's version are:

1. used lists to replace numpy matrixes; Lists seem to run faster on reading rpo772.wdm file;
2. added a loop to iterate through each group and used ending date as the ending condition. VB code was implemented this way.

All three of these WDM files can be read successfully now.

• rpo772.wdm
• test.wdm
• 2_RPO_SWMM48LINKS2017_CBOD.wdm

[aufdenkampe]

Moving excerpts of conversation from Teams on Mar. 1, 2021:

[Monday 10:10 AM] Paul Tomasula

Thanks Hua, excellent work on this! I reviewed the code and it looks solid to me. Have you done any comparisons between the converted time series with those from the BASINS code?

Ideally we automate this comparison as a test in the future, but for now spot checking a few time series by hand should do. I'd mostly be interested in making sure some of the time series in the RPO_SWMM48LINKS2017_wCBOD_June2019.wdm match since that file has the funky block ending on the 511th element of a record.

[Monday 10:24 AM] Hua Tao

Now that I remember that dates in the VB generated csv file is Julian Year. You must convert them back to normal date/time format. I can point you to those. I placed the two folders of csv here: S:\DOESBIR2\HSP2-Dev\WDM-TestFiles

[Yesterday 11:12 AM] Paul Tomasula

Hey Hua, I did some spot checking on times series and values are near identical. The run time is considerable though. I'm going to take a crack at some refactoring. I'll primarily focus on getting numba implemented and we'll see how it looks from there.

[aufdenkampe]

Moving excerpts of conversation from Teams on Mar. 3, 2021:

[9:03 AM] Paul Tomasula

Refactoring and Reader Performance

Hua and Anthony, I did some general refactoring of the WDM reader yesterday to tidy up the code (remove deprecated functions and irrelevant comments, clarify variables through better naming, raise exceptions in place of print statement for errors, etc). (Commits b600ebf, 2561869, 2e205f1). I also tried to implement numba, but was ultimately unsuccessful.

Numba implementation is best accomplished through the @njit decorator, which performs just in time compiling of a function without using python structures. When done correctly, this will allow python code to approach the speeds C or Fortran with a little extra runtime for compiling. However numba also has rather restrictive data type allowances and is particularly restrictive with datetimes. The only supported structures are numpy's datatime64 and timedelta64 objects, and even then it's not a full implementation. For example creation on a datetime64 object in not presently possible within a numba function. Numpy's datetime object is also extremely difficult to use if you have a variable unit for timedetla (i.e. some in years, others in day, minutes, etc) like the WDM files do. Years and months are nonlinear so there's no easy way to convert those units to a timedelta in seconds. I did start to think about a work around to manually convert years and months. The implementation will be pretty messy though and would reduce code readability and maintainability. Also the complexity of the workaround doesn't guarantee a runtime boast because it needs to run a bunch of extra calculations to convert years and months to days.

I guess my question is, what do you all think out about this? The runtime is still considerable. It took my machine over 20 minutes to convert the RPO_SWMM48LINKS2017_wCBOD_June2019.wdm to an HDF file. Should we leave it at this? Any other thoughts to bump up performance? Parallelization will help a little bit but we'd still be looking at a considerable runtime.

[9:32 AM] Hua Tao

Paul, I never used Numba so I don't have a lot insight about it. The code that we developed to read WDM files are pretty simple. The only part that might be time consuming is probably those date/time calculation in Pandas. The rest of the code is very primitive calculations and loops. I don't know whether and how much Numba is able to recompile and speed up those. I suspect that the writing a HD5 file (hundreds Meg) to disk is probably the most time consuming part. If we only keep the WDM data in the RAM, it is probably not too bad. We can strip the HD5 out from the reader and make it a separate ExporttoHD5 method() and only call it when needed.

[9:35 AM] Anthony Aufdenkampe

Paul, I'm wondering if we follow Hua's suggestion to skip use of Numpy until the very end of the ingestion process, or at the very least, don't convert the datetime from Fortran (isn't it just a float?) until after all the looping is done. That way you can use Numba.

[9:51 AM] Paul Tomasula

I couple of quick response though. I isolated the process_groups method for testing and just running that I average about 3 seconds per time series. The RPO...June2019.wdm file has 432 time series so its 21ish minutes for executing that step. This definitely appears to be at least one of the bottle necks. The writing could also be pretty substantial though I haven't tested that.

The issue with isolating the datetime conversions is that's currently integral to the processing loop. We use the datetime to indicate when the loop needs to stop executing. Even putting that piece aside for now. The only datetime information stored is in the block control word and it specifies the timedelta for each value. This means that each timestep is dependent on the previous timestep in its calculation. So we could maybe store the timedelta information to have a numba compatible loop but I'd imagine looping through a second time to convert the timedelta to a timeseries is probably slower than our current approach.

[9:57 AM] Hua Tao

Paul, then I would think parallelization would help improve speed. There is no dependencies between time series. This approach is a low hanging fruit, I would think. What do you think ?

[10:00 AM] Paul Tomasula

Yeah I agree. We definitely want parallelization added in as well regardless. Maybe as a first step we implement parallelization and see what performance looks like? Then revisit refining the processing method only if necessary.

[11:43 AM] Anthony Aufdenkampe

​​Paul, that's great that you were able to profile the typical speed for each time series. I agree that 3 s per series sounds fast, but perhaps you could speed it up to 0.3 s or faster by using Numba. Let's try get the time stamp conversion outside the loop.

Hua, what is the total run time of the VB code to run that same RPO...June2019.wdm file? Knowing that would be very helpful.

Paul, when its time to parallelize, I would really like to go with Dask for this and for the rest of the rep. It looks relatively simple according to this tutorial: http://gallery.pangeo.io/repos/pangeo-data/pangeo-tutorial-gallery/dask.html. Also, Dask also provides Dask Dataframe, which is very similar to Pandas Dataframe but has built-in parallelization and also allows for very large datasets that are bigger than RAM. We might want to use these instead of Pandas dataframes throughout the repo.

[12:02 PM] Hua Tao

VB.NET uses Just-In-Time compiler so performance wise it will be worse but comparable to FORTRAN. I haven't profiled performance for this particular reader. I would think it's 10 to 20 times faster than Python based on my past experience. It's great that you suggested going with Dask. There are many options. This saves us time.

Anthony Aufdenkampe

Hua, great! Could you try timing the VB.net code today or tomorrow?

​[12:22 PM] Paul Tomasula

Thanks Hua, Sounds good! I need to change gears for this afternoon but knowing the runtime for VB.NET will be helpful to see the level of performance we're aiming for. In the meantime I'll start thinking about if it is possible to pull the time conversion out into a separate loop or some other work around. Our call with Jack on Friday might also offer up some clarity. If we can figure out how to decouple datetime from the loop control (using an indexing approach or something similar) then it might be possible to perform those conversions in a later step.

[12:42 PM] Hua Tao

VB.NET: half second on average for each time series. Yes, only 1/2 second each time series. This is still in debug mode. Release mode will be even faster.

​[1:56 PM] Anthony Aufdenkampe

Hua, thanks for timing it! So ~6x faster in debug mode, and faster once compiled. That's our goal! I suspect we can get close to that fast if we can use numba, and faster once we apply Dask.

Paul & Hua, let's chat sometime before Friday to figure this out. Remind me how datetime is stored in the WDM again? Maybe we emulate that.