Faraday Complexity: A How to Guide

In this Jupyter notebook, we'll give an overview of how our scripts work and how to utilize them for constructing and training neural networks to identify complex Faraday sources.

Generating data with possum and simulate

The modules possum and simulate are used to generate data. By itself, the functions with possum are more focused on creating a spectrum for a single source, while the inherited functions from simulate allow us to generate a large number of sources with a few function calls.

Let's begin by importing the possum class. Since possum inherits simulate, it will automatically load as well.

There are a few frequency bands that are already part of the possum class:

ASKAP 12 (700 - 1300 MHz, 1500 - 1800 MHz)
ASKAP 36 (1130 - 1430 MHZ)
Westerbork Synthesis Radio Telescope (WRST, 310 - 380 MHz)

These can be created by calling the functions

_createASKAP12()
_createASKAP36()
_createWSRT()

which store the frequencies (in Hz) in the variable

nu_

If you need another frequency range, the function _createFrequency will be of good use. It's function call is:

_createFrequency(numin, numax, nchan, store=False)

It takes four parameters: the minimum (numin) and maximum (numax) frequencies, as well as the number of channels (nchan) and a store parameter. If you would like to store the results in self.nu_, then set store=True; otherwise the frequency array will be returned.

As an example of how the functions within this class work, we can create the ASKAP 12 coverage by running the following commands and verify that our range is over the proper frequency range:

Polarization Spectra

Simple and complex polarization spectra can both be created using the _createNspec function, which takes in the parameters flux, depth, chi, and sig, where flux is the strength of the source, depth is the Faraday depth, chi is the phase $\chi$, and sig(ma) is the standard deviation of the noise $\sigma$. The function call is as follows

_createNspec(flux, depth, chi, sig)

For complex sources, one passes in the flux, depth, and chi as an array or list. The polarization spectra is stored in the variable

polarization_

As an illustration, let's create a polarization spectrum of a simple and a complex source and plot the real ($Q$) and imaginary ($U$) components. We'll create a custom frequency range to illustrate the _createFrequency function.

Faraday Rotation

The Faraday rotation spectra can be created using the function _createFaradaySpectrum, which takes in two parameters: the minimum (philo) and maximum (phihi) Faraday depths:

_createFaradaySpectrum(philo, phihi)

and stores the rotation measure $\phi$ and the source spectrum in the following variables:

phi_
faraday_

In the script below, we generate the Faraday spectrum of the polarization spectrum show above and plot the real and imaginary components, along with their magnitude. The faraday depths for the two sources are indicated by vertical dashed lines.

Generating Parameters for a Large Number of Spectra

That's everything you need to generate a simulated polarization and faraday spectra using the possum class. To generate a large number of sources, we turn ourselves to the functions inside the simulate class.

The function _generateParams will generate parameters for $N$ sources, with the phase $\chi$, depth $\phi$, and flux values stored as a tuple for the complex sources. The parameters are randomly sampled over a uniform distributed between the chosen minimum and maximum values. The function call is of the form

_generateParams(N, depthMin=-50, depthMax=50, 
        fluxMin=0.01, fluxMax=1, chiMin=0, chiMax=np.pi,
        noiseMin=0.01, noiseMax=1.0, pcomplex=0.35, seed=8595)

where N is the number of spectral parameters to generate and pcomplex is the probability that the source is complex. seed sets the random seed number, while the other parameters set the range of values to sample over. The flux of the first component is set to 1, while the parameters fluxMin and fluxMax refer only to the strength of the second component.

Note that noiseM-- refers to the standard deviation of the noise; the notation "noise" is a bit inconsistence with the remaining parts of the class, which all reference the noise via the parameter sig.

The sampled values are stored in the arrays

chi_
depth_
flux_
sig_

while the label of the source (simple = 0, complex = 1) is stored in the array

label_

As an example, we can generate parameters for 5 spectra with a complex probabliity of 50% by the following commands:

and output the parameters by printing the appropriate variables. Note that the complex source parameters are stored as tuples.

Generating a Large Number of Spectra

The function call of _simulateNspec is as follows:

_simulateNspec(N=5, pcomplex=0.35, width=100, center=False,
                seed=8008, save=False, dir='./', timeit=False)

The function first checks to see if the function _generateParams has been previously run by checking if the parameter variable self.label_ has been stored. If not, then the function runs _generateParams using the defaut parameter ranges, but inserting pcomplex and seed into the function call.

With the generated parameters, the function then iterates through each parameter space and generates both the polarization and Faraday spectra, storing the results in the 3D tensors

X_
S_

where X_ is the Faraday spectra and S_ is the polarization spectra. The first dimension of each tensor is the oberservation instance (that is, it has length $N$). The second dimension contains the spectra, while the third dimension is for the real (0) and imaginary (1) components. The $Q$ and $U$ values can be accessed via

Q = S_[:,:,0]
U = S_[:,:,1]

The width parameter is used for cropping the Faraday spectrum. If center is set to true, then the cut is +/- width from the peak in the spectrum (otherwise at $\phi=0$). If the spectra is noisy, centering may cause some issues.

The parameter save is a boolean asking whether to save the spectra (and parameters) as numpy arrays, and dir is a path to the directory that you want the arrays to be save to. If the directory doesn't exist, then the function creates the directory.

The final parameter timeit is a boolean to keep track of how long it takes to run. It will output the total time in minutes once the script is complete.

In the following cell, 6 complex spectra are created using the default ranges in _generateParams:

The real, imaginary, and absolute values of the Faraday spectrum are shown in the plot below:

This concludes our discussion of the data generation process. In the next session, we will look at creating a 1D convolutional neural network (CNN) using the inception class, as well as how to utilize some of the analysis and plot functions that come with it.

Building CNNs with inception

Recently, there has been a shift away from traditional CNNs to an "inception" developed by the GoogLeNet team (Szegedy et al. 2014). In the inception model a series of convolutions is done in parallel, each with a different kernel size. A key feature of the inception model, however, is the 1x1 convolution (which replaces MaxPooling from earlier CNN models) that serves to collapse the data at intervals. To build an inception model, we can use the functions within the inception module.

We can create a simple inception layer using the _inception function, whose call is of the form

_inception(convl=[3,5], pool=[3], act='relu', strides=2, padding='same')

The parameters of the _inception function allows us to simply and succinctly create a multitude of inception variations. For instance, the arguments [3,5,23], pool=[3] create four "branches;" one with a kernel of size 3, one with a kernel of size 5, one with a kernel of size 23, and one max-pooling branch. Then, we want to add another inception cell in sequence, we can simply call _inception([a,b,...,k], pool=[p]).

Before using this function, however, we first have to load in a training dataset so that the module can extract the input dimensions. The following functions can be used to load in a dataset and are part of the loadData module, inherited by inception

_loadTrain(dir)
_loadTest(dir)
_loadValid(dir)

Each function takes in the directory for the corresponding dataset and opens the necessary files. Since we have a training and validation stored in the directories "data/train" and "data/valid", let's load these datasets

For the training dataset, the following variables are created

trainX_
trainY_
trainLabel_
trainChi_
trainDepth_
trainFlux_
trainSig_

The variable trainX_ contains the feature tensor for the $N$ faraday spectra, while trainY_ contains the labels in a format readily usable by Keras. The original labels are stored in the variable trainLabel_, while the remaining four variables hold the parameters that generated the data.

With the dataset loaded, we can now construct an CNN with an initial inception layer. Let's use the convolution kernels 3, 5, and 23 (the FWHM)

If we want more layers, we can just call the inception function more times. Note, however, that the function _inception is only a basic block and doesn't implement batch normalization, and consequently is only able to produce workable models that contain up to 2 inception layers. The function _inception2 implements batch normalization along with an activation function, so this function should be used when constructing deeper networks.

To finish constructing a 1-layer inception model, we need to add the final few layers. We first need to flatten our tensor, which can be done using the function _flatten

We can then add a couple fully connected (Dense) layers using the function _dense,

_dense(z, act='relu', drop=0.5, nlayers=1)

where z is the number of neurons in the hidden layer, act is the activation function, drop is the dropout rate, and nlayers is the number of hidden layers to create. So to add two hidden layers with 512 neurons / layer and passing through the ELU activation function, we would just run the following

To complete our model, we can compile it using the function _compile, whose function call is

_compile(classes=2, act='softmax', optimizer='adadelta', loss='binary_crossentropy', 
            metrics=['binary_accuracy'], weights=None)

The parameter classes is the number of output classes that we are trying to predict (which is 2 in our case: simple or complex), while optimizer is the optimization method we would like to use.

In the plots module, there is a function

_plotCNN(to_file)

which will plot our model and save it to the file "to_file". Let's save a plot of our CNN and display the image

As another example, let's build an even deeper network

Analyzing an Imbalanced Dataset using plots and analysis

In this section we'll examine the performance of a trained model on an imbalanced dataset with 100,000 sources (~3% complex) using the plots and analysis modules.

The model we will be using is a 2-layer inception model, which was trained on 30,000 sources and validated on 10,000 sources (each ~50% complex) with GPUs using the Argon cluster at the University of Iowa. The model architecture is displayed below.

The logfile is stored at "model/log/i2layer_final.log". We can load in the logfile using the _loadLog function and display the accuracy and loss using the fuctions _plotAcc and _plotLoss. The _loadLog functions takes in the filename as an input parameter:

_loadLog(logfile)

while the calls for the plotting functions are

_plotLoss(logfile=None, save=False, to_file='loss_vs_epoch.pdf', fontsize=20)
_plotAcc(logfile=None, save=False, to_file='acc_vs_epoch.pdf', fontsize=20)

save is a boolean to indicate whether to save the image to to_file; fontsize is used for adjusting the size of the axes labels. If logfile is not equal to "None", then the function attempts to load in the logfile; otherwise it will attempt to access a pre-loaded logfile.

In this case the model was subjected to Early Stopping by monitoring the validation loss with a patience of 25 epochs. The lowest validation loss occured for epoch = 51 (index starting at 0), although there is a slight downward trend, suggesting we may be able to get slightly better results by running for a larger number of epochs. However, running more epochs was found to produce only a marginal improvement at best.

The best-fit model is saved in "models/model/i2layer_final.h5" and we can load it using the _loadModel function in inception:

To evaluate its performance, we can load in the test set using the _loadTest function and passing in the appropriate directory, which in this case is "data/test"

We can then use our trained model to predict the type of faraday spectrum for each instance using a 50% threshold for being classified as a complex source

and check the confusion matrix, which has already been loaded from sklearn

So most of the sources that were predicted to be complex aren't complex due to the large class imbalance. Another way to evaluate the model performace is to plot an ROC curve, which is conveniently included in the plots module as the function _plotROC. The area under the curve (auc) is computed and stored in the variable self.auc_, so can can check this value as well.

We can separate the complex and simple sources and store their parameters in a dataframe by calling the functions _getComplexParams and _getSimpleParams, which create the pandas dataframes self.dfComplex_ and self.dfSimple_ respectively. These can be created at the same time using the function

_getParams(prob)

which takes the threshold probability prob as an input parameter for labelling a source as simple or complex. These are often created automatically when calling plotting functions within the plots module.

We can do some simple comparisons to see how the various parameters (chi, depth, flux, sig) are related to the probability outputs for the complex sources using the _plotParamProb function, whose function call is as follows:

_plotParamProb(param, kind='kde', gridsize=50, 
           save=False, to_file="FluxProb.pdf", fontscale=1.25)

param is a string that identifies the column name in the dataframe ("chi", "depth", "flux", "sig"). Let's see how well the flux of the second component affects our model's prediction:

The model has a steep drop once the fluxes get below 0.2, but is doing quite well otherwise. Let's check how the model performs on the (absolute) difference in the two Faraday components.