IDL to Python Translation Guide

This part of the guide is going to be more personal and contains lessons from my own (Joel Dunstan, SkyWa7ch3r) approx. 3 years experience of translating the majority of FHD from IDL to Python to create pyfhd.

Introduction

Hello there weary traveler, you wish to translate from IDL to Python or have been tasked with doing so, first thing my sincerest apologies for your future sanity and also to your colleagues as you likely vent to them about the annoyances of doing the translation. Here I’ll try to put all the discoveries, roadbumps and potholes that will slow you down along the way. Be prepared to laugh, sigh, be filled with rage and maybe shed a tear if you find out a problem you’ve been facing for a week is because of some dumb interesting behaviour.

Getting IDL

First thing, you’ll need access to IDL, this is easier said than done as you require a license to use the language if you’re lucky your instituion has an IDL license or the HPC you use has an IDL license. If you plan on installing IDL to your local machine you will need to contact NV5 Geospatial directly to get a license and get the download link through their locked download links. If you’re institution has a license server you should be able to get a license from the server and an install package from your institution. There are open source versions of IDL, GDL and Fawlty, however in my experience they didn’t support running all the functions needed for FHD or produced different results to what’s expected even in simple scenarios, hopefully that has changed and you can try to use them, feel free to let me know how that went.

Now that you have IDL, you’ll need to follow the FHD install steps. You’ll notice that all the functions get added to the IDL path, IDL only really has one namespace for functions, do keep that in mind unlike Python which has packages/libraries, modules/subpackages which divide up the name spaces and require explicit imports.

FHD Dependency Equivalents

Now in terms of matching up some of the FHD dependencies and pyfhd dependencies, anything from the FHD, fhdps_utils, and pipeline_scripts repositories will require straight translation line by line. For the IDLAstro library from NASA, a mix of astropy, numpy and scipy should allow you to either drectly replace functions when they appear in the FHD code or allow much easier translation/rewriting of the function. For the IDLAstro library all the pro files are publicly available so if there is a need to directly translate the function into Python that is possible. coyote does much of the graphical heavy lifting in FHD, for any plots inside of pyfhd, the plotting/imaging libraries matplotlib, seaborn, plotly, Pillow and OpenCV should get you what you need to make almost any visualization you wish to make (plotly is excellent when you need interactivity). For HEALPix you have the healpy library maintained by the same organization, so functions should have similar names and similar purposes without any translation needed.

The Right Tools for the Right Job

In terms of tools you can use to translate IDL to Python, after all this time, unfortunately for you, I’ve decided the best tool is you. None of the tools out there like idl2py successfully translate to Python in a way that makes it faster than you translating the function line by line, they even occasionally miss things like FLTARR which translate perfectly into np.zeros. To be fair on these tools its an almost impossible task to automate because in IDL keywords in functions can entirely change how the function works thus any code you make to automatically translate will be full of exceptions and will still have problems translating pointer arrays due to the shapes of said arrays not always being easy to know ahead of time. Unfortunately, not even ChatGPT or LLMs can completely help you either, although to be honest they do a somewhat decent job, feel free to experiment with them and let me know how it goes, although I suspect it will miss some of the interesting behaviour that I will point out later. If you’re using VSCode, NV5 recently brought out an official VSCode extension for IDL, check it out here. IDL for VSCode will give you proper syntax highlighting for pro files, and hopefully they will continue to increase it’s usefulness in the future.

Another pair of tools that are worth talking about is the IDL to Python Bridge and Python to IDL Bridge that is provided as a part of the IDL package. The Python to IDL Bridge can provide some nice to have features such as running IDL code in a Jupyter notebook which can be particularly useful for debugging some IDL to Python translations because you can rerun just parts of the code you’re having issues with translating. The IDL to Python Bridge does also make it possible to use Python within IDL, while the Python to IDL Bridge makes it possible to use IDL in Python. At the time of writing this its now possible to pass objects to and from IDL in the Python to IDL Bridge, however I cannot confirm if you can still only run the built-in IDL functions in the Python to IDL Bridge. When I previously used the bridge I wasn’t able to run functions made inside of FHD such as weight_invert (I suspected this was due to the IDL_PATH being setup everytime you re-ran IDL.run for every command). Hopefully if you try the bridges on the latest version of IDL, you may have better luck meaning you could have more power for debugging than I did because you can directly run the IDL code side by side and directly compare the IDL and Python outputs without the need for sav files. Feel free to let me know how this changes or to update this section yourself.

The Hidden Gems in IDL to make you tear your hair out

There are several gotchas hidden thoughout IDL, usually most of them won’t appear to an unexperienced IDL person who hasn’t been burned by them before until they try to run the Python equivalent, for example, me and hopefully not you. I will try to list as many as I can here, lot of these can be found in IDL documentation as extra notes or footnotes (rather than a big notice saying “Hey this is important, if you use this function this way, this function is entirely different” which would be a little nicer).

Column Based vs Row-Based Indexing

IDL is column based while Python is row based, meaning you will likely need to change the shape of certain arrays for the multiplication of arrays and matrices to work as you expect. It’s the reason why the visibility arrays go by polarization then frequency then baselines in pyfhd compared to polarization then baselines then frequencies in FHD as it allows us to translate the indexing as it’s done in FHD without having to constantly swap indexes around during translation which gets confusing.

Single Precision wins everytime

IDL is single precision even if it kills itself, in many cases even if you specify double precision using the /DOUBLE keyword, there’s bound to be some function somewhere that is restricted to single precision. This is slowly changing.

Input or Output…why not both? (Function Parameters)

Function parameters in IDL can be input or output! If during your translation a variable exists seemingly without being initialized somewhere, check the function calls, parameters store the results, this is the case with many IDL functions,like HISTOGRAM where the REVERSE_INDICES=ri part of the call actually returns the reverse indices into an ri variable. This is also used extensively for many FHD functions and can get confusing, in Python doing returns this way is possible but for the sake of being consistent, any and all returns are done with an explicit return statement and put into the docstrings as a return.

Python Subsetting vs IDL subsetting

When dealing with how IDL deals with subsetting, it’s important to remember than IDL is inclusive of both the start and the end index, i.e. if we index an array between the 0th value and the 3rd value, we will get the 0th value and all the values upto and including the 3rd value. In Python the subsetting is only inclusive of the start but not the end index, i.e. if we index an array between the 0th value and the 3rd value, we will get the 0th value and all the values upto the 3rd value, but not including the third value.

IDL
```idl
IDL> test = [1,2,3,4,5]
IDL> test[0:3]
       1       2       3       4
```

Python
```python
>>> test = np.array([1,2,3,4,5])
>>> test[0:3]
array([1, 2, 3])
```

The inclusion of both ends in IDL also works the same for loops of any kind,
so when translating you should be careful with how you index the array at the
end primarily or the end conditions for your loop. You will find many users
of IDL have already taken this into account and loop `n` times going from
`0` to `n - 1`, so in python using the `range(n)` function should do the trick.

What the /EVEN (IDL Median)

The median function in IDL doesn’t do the average between the two middle numbers by default when dealing with an array of even size unless you use the /EVEN keyword, it will take the maximum of the two numbers when dealing with an even sized array be default. This has often made a difference when dealing with close to single precision values, in my experience you’re usually better off using np.median in most cases and accepting the differences, there has only been one notable case where we needed the IDL median over np.median and that’s during the calibration process (as the comparison to conv_thresh would fail due to the precision differences).

MATRIX_MULTIPLY: Is it a dot product or outer product, you’ll find out at run time!

Matrix Multiply does the dot product for matrices, and because IDL is column based while Python is row based, you will need to switch around the order of multiplication in Python compared to IDL. Furthermore, matrix multiply does the outer product when you make on the provided matrices/arrays 1D array. So everytime you see matrix multiply you’ll likely have to check at runtime during testing if the multiply is meant to be the dot product or outer product.

NaN Shenanigans (IDL total and mean)

In the case of an array containing only NaNs there is a slight difference between IDL’s total and mean, where total(x,/nan) will give you 0 while mean(x, /nan) will give you a NaN

COMPLEX ATAN

In IDL ATAN when used with complex numbers and the /PHASE keyword gives you the phase of the imaginary number. In Python the equivalent function is np.angle. Do not be tempted to compute it using np.atan – it can give you NaNs unexpectedly.

DIVIDE BY ZERO: 1/0 + 10 = 11 Folks

When doing a divide by 0 in IDL, it does produce an error % Program caused arithmetic error: Integer divide by 0 however it doesn’t actually break the program, IDL continues to work, watch out for it! Don’t ask me why they’ve done it…and if you don’t believe me here is an example in IDL 8.8.0:

IDL> test = [[1.,2.],[3.,4.]]
IDL> test[0] = 1/0 + 10
% Program caused arithmetic error: Integer divide by 0
IDL> test
11.000000       2.0000000
3.0000000       4.0000000

Indexing beyond the size of the array will work

Indexing arrays using other arrays can also lead to dumb interesting behaviour, if we follow our test array again of 4 values and we try to access the 5th element of the array we expectedly get an error. idl     IDL> test = [[1.,2.],[3.,4.]]     IDL> test[5]     % Attempt to subscript TEST with <INT      (       5)> is out of range.     % Execution halted at: $MAIN$     However if we put that index into another array, let’s call it idx and then access the test array with the idx array we get the following: idl     IDL> test = [[1.,2.],[3.,4.]]     IDL> idx = [5]     IDL> test[idx]        4.0000000     Yep that’s the last value of the array, and to really get the ball into the endzone: idl     IDL> idx = indgen(10) + 1000     IDL> idx       1000    1001    1002    1003    1004    1005    1006    1007    1008    1009     IDL> test[idx]        4.0000000       4.0000000       4.0000000       4.0000000       4.0000000       4.0000000       4.0000000        4.0000000       4.0000000       4.0000000     In summary, watch out for this behaviour, it usually comes up in C when you use pointers beyond their designated memory addresses, resulting in overflows, however for whatever reason IDL does their own version of overflow. In Python, this functionality of indexing beyond the array, does not work, and is expected to never work (and it never should).

Poly fit/val

The POLY_FIT function in IDL has the ability to do both polyfit and polyval, so watch out for the use of the YFIT keyword. In general, poly_fit(a, b, 2) in IDL is np.polynomial.polynomial.Polynomial.fit(a, b, deg=2).convert().coef in Python. If you see the YFIT keyword, then it’s a case of calling np.polynomial.polynomial.polyval on the calculated coefficients i.e. phase_fit = np.polynomial.polynomial.polyval(a, np.polynomial.polynomial.Polynomial.fit(a, b, deg=2).convert().coef).

Standard Linear Algebra Shenanigans (LA_LEAST_SQUARES)

LA_LEAST_SQUARES in IDL and np.linalg.lstsq in Python (and the SciPy least square functions as well) will produce different results due to the solvers in use, and the assumption that the design matrix given to the solvers is full rank in LA_LEAST_SQUARES while NumPy and SciPy assuming that the design matrix is not full rank. Even using the same solvers doesn’t give the same results, so it might come down to how each of the functions were compiled, either way, keep in consideration they just won’t get the same results in most cirumstances.

REBIN

The REBIN has been replicated in pyfhd, however it’s use should be used sparingly as its doing interpolation when increasing an array in size, and averaging when decreasing the size of an array i.e. it’s upscaling or downscaling and it’s treating your array or matrix as an image rather than any ordinary array. In many cases it could be best to use np.tile, np.pad or np.repeat to do the same task more consistently.

That’s a smooth move (IDL SMOOTH)

The SMOOTH function in IDL is a boxcar averaging function so it should be replaced with scipy.ndimage.uniform_filter and the /edge keywords will dictate the mode you need to use for the uniform_filter function.

Where is COMPLEX? (IDL WHERE and COMPLEX numbers)

When dealing with complex numbers in IDL, it’s possible to get different results when using WHERE in IDL and np.where as WHERE in IDL applies to the absolute values of the numbers, while np.where looks at just the real numbers for example:

IDL

IDL> test = COMPLEX([1,2,0],[1,2,3])
IDL> WHERE(test gt 1)
         0           1           2

Python

>>> test = np.array([1 + 1j, 2 + 2j, 0 + 3j])
>>> np.where(test > 1)
(array([0, 1]),)
>>> np.where(np.abs(test) > 1)
(array([0, 1, 2]),)

It has only affected my translation once, but something to be aware of.

The Point (IDL Pointers)

Pointers are used in IDL a lot, especially in FHD, when translating pointer arrays my advice is to try and find out the full shape of the array and set the dtype appropriately, rather than creating an object array or list to represent the pointer array. object arrays in NumPy usually lead to more problems than they’re worth, and make using the built in vectorized functions a real pain.

I feel sorry for the person that had to make this function (scipy.io.readsav)

There is a readsav function in scipy.io that gives you the ability to read in sav files, which is great that it exists. Unfortunately though, readsav is usually quite slow due to IDL’s differences to other languages down to the byte level, as such it’s having to loop through all the bytes to get the necessary bits out (poor thing, imagine having to read a 44GB file several bytes at a time in a python loop). readsav does read the sav file into a python dictionary which sounds good, except that it turns structures that were saved into the sav files into np.recarrays or NumPy record arrays. NumPy record arrays are often slower than dictionaries to access and in the case of pyfhd offer no benefits over dictionaries (in fact I’d argue a dictionary is better is almost all circumstances). readsav also usually makes a mess of pointer arrays too by turning pointer arrays into object arrays which each index contains a numpy.ndarray, meaning to access an array you might have to do something like sav_file['array'][0]array[0][0] to access the first numpy array in the array of arrays. As such I developed the pyfhd.io.recarray_to_dict function which can take in a np.recarray or dict and turn any record arrays into a dict. recarray_to_dict works recursively too, so any sub-record arrays will also turn into a dict, furthermore, additional formatting of arrays will take place to turn any object arrays into a proper NumPy array with a numpy data type (dtype). recarray_to_dict also deals with the values that are scalar values when loaded in from the sav file as well, readsav usually turns scalar values into arrays with a single vaue which are inconvenient to access. recarray_to_dict is always a work in progress, happy to try and edit the function as necessary.

• If you’re dealing with large pointer arrays that can’t convert to a proper dtype array due their size, e.g. the beam_ptr complex array was an example of this with a size of 2*384*8128*51*51*2916*16 bytes or ~757.5 TB, then ideally you need to convert these into HDF5 (with h5py) files which can allow you to chunk large datasets. It’s also worth checking that there aren’t pointers referencing other pointers, this was the case with beam_ptr as such the actual size ended up being 2*384*51*51*2916*16 bytes or ~93.2 GB instead. It’s also possible to use other frameworks like Dask to achieve what we’re doing with h5py, but we can use h5py without having to worry about compatiblity with NumPy functions.

Concatenating Arrays

In IDL concatenating arrays can be done in many ways, one of the ways you may see is like so:

IDL> test = [1,2,3]
IDL> test = [test, test]
IDL> test
       1       2       3       1       2       3

So be on the look out for that sort of behaviour.

It’s the minimum and maximum I can do (IDL < and > operators)

The > and < operators in IDL do the same behaviour as np.maximum and np.minimum respectively, do not get them confused with the GT (greater than) and LT (less than) operators in IDL. Furthermore, the > and < operators in IDL can be chained together, I’ll provide many examples below: idl     IDL> test = [-1,2,3,-4,5]     IDL> test > 0        0       2       3       0       5     IDL> test < 0       -1       0       0      -4       0     IDL> test < 0 > 3        3       3       3       3       3     IDL> test < 1 > 0        0       1       1       0       1    

The bottom example is the most common example you'll see in `FHD` as it
turns all numbers above 1 into 1, and the numbers below 0 into 0 allowing to
quickly make a flagged array. To show the same examples in Python look below:

```python
>>> test = np.array([-1,2,3,-4,5])
>>> np.maximum(test, 0)
array([0, 2, 3, 0, 5])
>>> np.minimum(test, 0)
array([-1,  0,  0, -4,  0])
>>> np.maximum(np.minimum(test, 0), 3)
array([3, 3, 3, 3, 3])
>>> np.maximum(np.minimum(test, 1), 0)
array([0, 1, 1, 0, 1])
```

Thanks for Reading, and Ye Be Warned

I hope this contribution guide helps you in your translation efforts with much less pain than me, and all of this helps you get your function into Python as quick as possible. Translation and testing of IDL to Python code can be a frustrating task even at the best of times. So…