Study the Universe with Python tutorial, part 2 -- power spectrum

In the first part of this series we discussed how to download the galaxy catalogue of the Baryon Oscillation Spectroscopic Survey (BOSS). We also made some plots to show the distribution of galaxies. In this blog post, we will calculate a summary statistic of the dataset, which is the first step to use the data to constrain cosmological parameters, like dark matter and dark energy.

As we saw in the last blog post, the BOSS dataset contains about 1 million galaxies and their distribution in the Universe. The position of these galaxies is what carries the cosmological information. Imagine that there would be more dark matter in the Universe. The additional matter would have additional gravitational force, which would clump together the material. On the other hand, with less dark matter, the galaxies would be more spread out.

We can now compare these distributions with our data from the BOSS survey, and depending on which of the two distributions looks more like the data, we can determine how much dark matter there is in the Universe.

In practice, we do not want to compare the actual distribution of galaxies, but instead we compare a summary statistic, meaning a compression of these galaxies, which carries all the information we need. There are many choices for such a summary statistic, but here we will use the power spectrum.

To calculate the power spectrum, we first need to download two more catalogues

wget -N https://data.sdss.org/sas/dr12/boss/lss/random0_DR12v5_CMASSLOWZTOT_North.fits.gz -P path/to/folder/
wget -N https://data.sdss.org/sas/dr12/boss/lss/random0_DR12v5_CMASSLOWZTOT_South.fits.gz -P path/to/folder/

which contain a random distribution of points. These distributions are needed to calibrate the data catalogues.

To read this file, you can use the read_data() function we defined earlier (the file is significantly larger and it can take a while to read it into memory). Now we can combine the data catalogues and random catalogue, assign the point distributions to a 3D grid and calculate the power spectrum, all using triumvirate.

from astropy import units as u
from astropy.coordinates import Distance
from astropy.cosmology import Planck15

from triumvirate.parameters import fetch_paramset_template
from triumvirate.catalogue import ParticleCatalogue
from triumvirate.twopt import compute_powspec
def measure_pk(cat_data, cat_random, tag=''):
    """
    Calculate the power spectrum monopole P0(k) using Triumvirate.
    Assumes cat_data and cat_random were read with read_data() that:
      - converts columns to native-endian
      - promotes RA/DEC/Z/NZ/WEIGHT_* to float64
    """
    # Parameters
    paramset = fetch_paramset_template('dict')
    paramset['boxsize']['x'] = 5000.0
    paramset['boxsize']['y'] = 10000.0
    paramset['boxsize']['z'] = 5000.0
    paramset['ngrid']['x']   = 512
    paramset['ngrid']['y']   = 1024
    paramset['ngrid']['z']   = 512 
    paramset['tags']['output']  = tag
    paramset['degrees']['ELL']  = 0  
    paramset.update({
        'range': [0.005, 0.295],
        'num_bins': 29,
        'norm_convention': 'particle',
        'statistic_type': 'powspec',
    })

    # Data catalogue: spherical -> Cartesian (Mpc)
    r1_mpc = Planck15.comoving_distance(cat_data['Z']).to_value(u.Mpc)
    r1_mpc_h = r1_mpc * Planck15.h  # now in Mpc/h
    ra1  = np.deg2rad(cat_data['RA'])
    dec1 = np.deg2rad(cat_data['DEC'])

    x1 = r1_mpc_h * np.cos(dec1) * np.cos(ra1)
    y1 = r1_mpc_h * np.cos(dec1) * np.sin(ra1)
    z1 = r1_mpc_h * np.sin(dec1)

    ws1 = cat_data['WEIGHT_SYSTOT'] * (cat_data['WEIGHT_NOZ'] + cat_data['WEIGHT_CP'] - 1.0)
    wc1 = cat_data['WEIGHT_FKP']
    nz1 = cat_data['NZ']

    catalogue_data = ParticleCatalogue(x1, y1, z1, nz=nz1, ws=ws1, wc=wc1)

    # Random catalogue
    r2_mpc = Planck15.comoving_distance(cat_random['Z']).to_value(u.Mpc)
    r2_mpc_h = r2_mpc * Planck15.h  # now in Mpc/h
    ra2  = np.deg2rad(cat_random['RA'])
    dec2 = np.deg2rad(cat_random['DEC'])

    x2 = r2_mpc_h * np.cos(dec2) * np.cos(ra2)
    y2 = r2_mpc_h * np.cos(dec2) * np.sin(ra2)
    z2 = r2_mpc_h * np.sin(dec2)

    wc2 = cat_random['WEIGHT_FKP']
    nz2 = cat_random['NZ']

    catalogue_rand = ParticleCatalogue(x2, y2, z2, nz=nz2, wc=wc2)

    # Compute power spectrum (monopole)
    results = compute_powspec(
        catalogue_data,
        catalogue_rand,
        degree=0,
        paramset=paramset,
        save='.txt'  # or a filename/path you prefer
    )
    return results

Here we use a grid with 512 grid points in two dimensions and 1024 in one dimension. There are also two weightings included in this calculation, the 'Weight' and 'WEIGHT_FKP' columns. The first weight refers to a correction of incompleteness in the data catalogue, while the second weight tries to optimise the signal-to-noise.

If we use this code to measure the BOSS power spectrum for the northern and southern part we get

Figure 3: Power spectrum measurements for the northern (orange)

and southern (blue) parts of the BOSS dataset.

Here we see that the power spectra of the northern and southern parts are slightly different. Whether these differences are significant is not clear though, since we don't have uncertainties on these measurements (yet).

To extract cosmological information from this power spectrum we need to get a model. We use the classy package

from classy import Class
def get_Pk(Omega_m = 0.3):

    #Start by specifying the cosmology
    Omega_b = 0.05
    Omega_cdm = Omega_m - Omega_b
    h = 0.7 #H0/100
    A_s = 2.1e-9
    n_s = 0.96

    #Create a params dictionary
    #Need to specify the max wavenumber
    k_max = 10 #UNITS: 1/Mpc

    params = {
                 'output':'mPk',
                 'non linear':'halofit',
                 'Omega_b':Omega_b,
                 'Omega_cdm':Omega_cdm,
                 'h':h,
                 'A_s':A_s,
                 'n_s':n_s,
                 'P_k_max_1/Mpc':k_max,
                 'z_max_pk':10. #Default value is 10
    }

    #Initialize the cosmology and compute everything
    cosmo = Class()
    cosmo.set(params)
    cosmo.compute()

    #Specify k and z
    k = np.logspace(-5, np.log10(k_max), num=1000) #Mpc^-1
    z = ??? Please calculate the mean redshift of your sample

    #Call these for the nonlinear and linear matter power spectra
    Pnonlin = np.array([cosmo.pk(ki, z) for ki in k])
    Plin = np.array([cosmo.pk_lin(ki, z) for ki in k])

    #NOTE: You will need to convert these to h/Mpc and (Mpc/h)^3
    #to use in the toolkit. To do this you would do:
    k /= h
    Plin *= h**3
    Pnonlin *= h**3
    return k, Plin

We can now plot the power spectrum for different values of cold dark matter density $\Omega_{cdm}$

Figure 4: Power spectrum models with different amounts of dark

matter using the class Boltzmann code.

Comparing the model with the data can allow us to determine the amount of dark matter in the Universe.

However, before we can go ahead and constrain cosmological parameters we have to get an estimate of the uncertainty on the measurements. That will be the subject of the next post.
You can find the code for this project on GitHub.
cheers
Florian

Study the Universe with Python tutorial, part 1 -- the BOSS dataset

In the next few blog posts, I will introduce a dataset of more than 1 million galaxies and I will show how easy it is to use this dataset to test cosmological models using Python. This tutorial is intended for non-experts or anybody who wants to get their hands on a very exciting dataset. If you are looking for an interesting Python project, look no further.

First, we need to download the data. We are going to use the dataset produced by the Baryon Oscillation Spectroscopic Survey (BOSS) collaboration, which represents the largest dataset of its kind to date. You can download all necessary data with the following 2 wget commands, just provide the destination path at the end. This will download a bit less than 3Gb of zipped data, which turns into 8Gb after you unzip it.

wget -N https://data.sdss.org/sas/dr12/boss/lss/galaxy_DR12v5_CMASSLOWZTOT_North.fits.gz -P path/to/folder/
wget -N https://data.sdss.org/sas/dr12/boss/lss/galaxy_DR12v5_CMASSLOWZTOT_South.fits.gz -P path/to/folder/

For this project, we need the Python package triumvirate. Please follow the install instructions on the linked page.
Now we are ready to go.

We start by reading in the data catalog using the fitsio package

from astropy.io import fits
from astropy.table import Table
def read_data(filename):
    """
    Read a FITS catalogue and:
      - convert all columns to native-endian (NumPy 2.x safe)
      - make RA, DEC, Z, NZ, and WEIGHT_* columns float64
      - ensure columns are contiguous
      - apply the redshift mask 0.5 < Z < 0.75
    Returns:
      filtered_table (Astropy Table), col_names (list)
    """
    with fits.open(filename, memmap=False) as hdul:
        data = hdul[1].data  # FITS_rec
        t = Table(data)      # work in an Astropy Table for easy column ops

    # 1) Make all columns native-endian and contiguous
    for name in t.colnames:
        col = np.asanyarray(t[name])
        dt = col.dtype
        if dt.byteorder not in ('=', '|'):
            # Convert big-endian to native-endian (NumPy 2.x compatible)
            native_dt = dt.newbyteorder('=')
            col = col.byteswap().view(native_dt)
        col = np.ascontiguousarray(col)
        t[name] = col

    # 2) Promote Triumvirate-relevant columns to float64 ("double")
    to_float64 = ['RA', 'DEC', 'Z', 'NZ', 'WEIGHT_SYSTOT', 'WEIGHT_NOZ', 'WEIGHT_CP', 'WEIGHT_FKP']
    for name in to_float64:
        if name in t.colnames:
            t[name] = np.ascontiguousarray(np.asarray(t[name], dtype=np.float64))

    # 3) Apply redshift mask
    if 'Z' not in t.colnames:
        raise KeyError("Column 'Z' not found in file.")
    redshift_mask = (t['Z'] > 0.5) & (t['Z'] < 0.75)
    t = t[redshift_mask]

    print("col_names =", t.colnames)
    print(len(t['Z']))

    return t, t.colnames

This function prints out the column headers, which should look like this

col_names =  ['RA', 'DEC', 'RUN', 'RERUN', 'CAMCOL', 'FIELD', 'ID', 'ICHUNK', 'IPOLY', 'ISECT', 'FRACPSF', 'EXPFLUX', 'DEVFLUX', 'PSFFLUX', 'MODELFLUX', 'FIBER2FLUX', 'R_DEV', 'EXTINCTION', 'PSF_FWHM', 'AIRMASS', 'SKYFLUX', 'EB_MINUS_V', 'IMAGE_DEPTH', 'IMATCH', 'Z', 'WEIGHT_FKP', 'WEIGHT_CP', 'WEIGHT_NOZ', 'WEIGHT_STAR', 'WEIGHT_SEEING', 'WEIGHT_SYSTOT', 'NZ', 'COMP', 'PLATE', 'FIBERID', 'MJD', 'FINALN', 'TILE', 'SPECTILE', 'ICOLLIDED', 'INGROUP', 'MULTGROUP']

In this tutorial, we will use only a subset of the columns given here. Namely, we will only use the positions of the galaxies given by the angles on the sky in 'RA' (right ascension) and 'DEC' (declination) as well as the redshift 'Z'. The redshift describes how much the light spectrum of the galaxy has been shifted in wavelength due to the expansion of the Universe. Most of the galaxies in this dataset are several billion light-years away from us, meaning the light travelled a significant fraction of the existence of the Universe before reaching us. Galaxies which are further away have a larger redshift and we will use the redshift later, together with the right ascension and declination, to get a 3D position for all galaxies. 

Let's explore the dataset a bit before we go any further. We can read in the data

    filename = 'galaxy_DR12v5_CMASSLOWZTOT_South.fits'
    data_south, _ = read_data(filename)
    filename = 'galaxy_DR12v5_CMASSLOWZTOT_North.fits'
    data_north, _ = read_data(filename)
    angular_plot(data_south, data_north)

and plot the angular distribution (right ascension and declination) with

import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt
def angular_plot(cat1, cat2):
    ''' Plot the catalogues '''
    cat1['RA'] -= 180
    c1 = down_sample([cat1['RA'], cat1['DEC']], 2000)
    cat2['RA'] -= 180
    c2 = down_sample([cat2['RA'], cat2['DEC']], 2000)

    plt.clf()
    plt.subplot(projection="aitoff")
    plt.title("BOSS DR12 survey footprint", y=1.1)
    plt.plot(np.radians(c1[0]), np.radians(c1[1]), '.')
    plt.plot(np.radians(c2[0]), np.radians(c2[1]), '.')
    plt.grid(True)
    plt.show()
    return 

Here we read in the two data catalogues, down-sample these catalogues (scatter plots with 1 million points will be difficult to handle), and convert the angles from degrees to radians. The down_sample() function is defined as

import random
def down_sample(in_list, N):
    n = len(in_list[0])
    if N <= n:
        idx = random.sample(range(n), N)           # no replacement
    else:
        idx = random.choices(range(n), k=N)        # with replacement

    out_list = []
    for cat in in_list:
        out_list.append([cat[i] for i in idx])
    return out_list

We make a plot by projecting the data onto a sphere using Aitoff projection. 

Figure 1: The sky distribution of galaxies in the BOSS galaxy survey using Aitoff projection.

Only 1000 points for each region are plotted.

This plot shows the sky coverage of the dataset. The galaxy survey observed in two patches on the sky, which explains why we have two files. We will call them north and south from now on even though east and west might be more intuitive from this figure. However, this is just because of the Aitoff projection. In fact, the plane of the Milky Way galaxy does go right in between the two patches, so that the southern part (blue) is indeed south of the Milky Way and the northern part (yellow) is north of it.

The northern part has about twice as many galaxies as the southern part and in total this dataset covers about 1 quarter of the sky.

We can also plot the redshift distribution

def redshift_plot(cat1, cat2):
    ''' Plot the catalogues '''
    bins = np.arange(0., 0.9, 0.01)
    weight = cat1['WEIGHT_SYSTOT'] * (cat1['WEIGHT_NOZ'] + cat1['WEIGHT_CP'] - 1.0)
    plt.hist(cat1['Z'], weights=weight, bins=bins, alpha=0.4, label='BOSS South')
    weight = cat2['WEIGHT_SYSTOT'] * (cat2['WEIGHT_NOZ'] + cat2['WEIGHT_CP'] - 1.0)
    plt.hist(cat2['Z'], weights=weight, bins=bins, alpha=0.4, label='BOSS North')
    plt.xlabel("redshift")
    plt.legend(loc=0)
    plt.show()
    return 

which looks like this

Figure 2: Redshift distribution of the BOSS galaxy sample. Redshift

0.7 corresponds to about 6 000 000 000 (6 billion) light-years. So when the light from

these galaxies was sent out, the Universe was only about 7 billion years old, 
while now it is 13 billion years old.

The next step is to turn the RA, DEC, z system into Cartesian coordinates. However, this transformation depends on the cosmological model. Using the current standard model of cosmology (also called $\Lambda$CDM) we can set up a cosmological framework by fixing a small set of parameters
(1) The amount of baryonic matter ($\Omega_b$)
(2) The amount of cold dark matter ($\Omega0_{cdm}$)

(3) The Hubble parameter $h$ (or more precisely the Hubble parameter divided by 100, $h = H/100$)

In principle, there are more free parameters, but these are the ones that matter to this analysis. Instead of fixing those parameters individually, we will use the Planck cosmology published in 2015

from mpl_toolkits import mplot3d
from astropy import units as u
from astropy.coordinates import SkyCoord, Distance
from astropy.cosmology import Planck15
def threeD_plot(cat1, cat2):
    # Get cartesian coordinates using the Planck 2018 cosmology
    c1 = down_sample([cat1['Z'], cat1['RA'], cat1['DEC']], 2000)
    dist1 = Planck15.comoving_distance(c1[0]).to_value(u.Mpc)
    dist1_h = dist1 * Planck15.h  # now in Mpc/h
    catalog1 = {}
    catalog1['x'] = dist1_h * np.cos(np.radians(c1[2])) * np.cos(np.radians(c1[1]))
    catalog1['y'] = dist1_h * np.cos(np.radians(c1[2])) * np.sin(np.radians(c1[1]))
    catalog1['z'] = dist1_h * np.sin(np.radians(c1[2])) 

    c2 = down_sample([cat2['Z'], cat2['RA'], cat2['DEC']], 2000)
    dist2 = Planck15.comoving_distance(c2[0]).to_value(u.Mpc)
    dist2_h = dist2 * Planck15.h  # now in Mpc/h
    catalog2 = {}
    catalog2['x'] = dist2_h * np.cos(np.radians(c2[2])) * np.cos(np.radians(c2[1]))
    catalog2['y'] = dist2_h * np.cos(np.radians(c2[2])) * np.sin(np.radians(c2[1]))
    catalog2['z'] = dist2_h * np.sin(np.radians(c2[2]))

    fig = plt.figure()
    ax = plt.axes(projection='3d')
    ax.scatter(catalog1['x'], catalog1['y'], catalog1['z'])
    ax.scatter(catalog2['x'], catalog2['y'], catalog2['z'])
    plt.show()
    return

The $\Lambda$CDM cosmology here assumes a flat universe with no curvature, which means that the dark energy parameter is automatically fixed through

\[
\Omega_{DarkEnergy} = 1 - \Omega_{cdm} + \Omega_b.
\]These parameters constrain the expansion history of the Universe, including its size and age. The dataset we just read with our Python script is powerful enough to constrain these parameters and, therefore, teach us something about dark matter and dark energy. We will get to that in a later post.

The code above converts the angles (RA, DEC) together with the redshift to a 3D position for the galaxies producing a plot which should look like this

Figure 3: 3D distribution of BOSS galaxies. Only 1000 randomly selected

galaxies out of more than 1 million in this dataset are included in this plot. 

Milky Way is located in the middle and the galaxies in the northern (yellow)

and southern (blue) part are observed in two opposite directions. 

And this concludes this first look at the BOSS dataset. This dataset was released in 2016 and represents cutting edge cosmological research. How to use this dataset to constrain dark matter and dark energy will be the subject of the next blog posts.

The entire Python code for this project is available on GitHub.

cheers
Florian

Elasticsearch and PubMed: Step 3, Diving in the Data

In the last two blog posts, I explained how to get a local copy of all 30 million PubMed publications and how to parse this data set into an elasticsearch database. In this blog post, I will describe a small project which leverages this data set.

We will go through all publications (title and abstract) and count how many publications contain a certain term. We will do that for different terms and plot them as a function of time. This will tell us which topics are currently trending in medical science. Here is a possible output of the program described in this post:

Figure 1: Number of PubMed publications with certain topics as a function of time.

This plot shows that there was a very steep increase in the publications about Ebola research in 2014. We can also see that the research on blood cancer or leukemia has decreased in the last 25 years, while the research on prostate cancer has increased in those years.

Elasticsearch uses a JSON based syntax, which indeed can be a bit cumbersome, but the documentation is very good.

First, we write the elasticsearch query which counts the papers. There are a couple of different queries in elasticsearch which can do that e.g. the match, term and match_phrase queries.

The match query analyzes the search term and then compares it to the stored indices.

{
    "query": {
        "match" : {
            "title" : "Cancer"
        }
    }
}

To analyze the term elasticsearch runs the same analyzer we used to index the field. You might remember from the last blog post that we used the standard analyzer, which tokenizes the text (splits it into words) and puts everything in lower case. So the match query above will return all papers which mention the term 'Cancer' or 'cancer'. 

Note that the elasticsearch approach can be limiting. For example, there is no easy way to distinguish lower case and upper case anymore because of the pre-processing of the standard analyzer. You can use the term query to prevent the search term to be analyzed

{
    "query": {
        "term" : { 
            "title" : "Cancer" 
        }
    }
}

but this search will return zero results since 'Cancer' with a capital 'C' does not exist anywhere in the index. For our case, this does not matter, but if you want to keep the ability to distinguish upper and lower case, you need to use a different analyzer when indexing your data.

In some cases, we might be interested to find multiple word terms or phrases, like 'drug resistance'. The match query is going to run the standard analyzer on this phrase and besides lower casing all words, the standard analyzer also tokenizes the phrase. So it turns 'drug resistance' into ['drug', 'resistance'] and looks for papers which contain either one of the words. We can additionally provide the and operator to enforce that documents need to contain both terms,

{
    "query": {
        "match" : {
            "title" : {
                "query": "drug resistance",
                "operator": "and"
            }
        }
    }
}

but this is still not exactly what we want since this will return documents which have the words 'resistance' and 'drug' anywhere. We want to enforce that these words need to be together. This can be achieved with the match_phrase query.

{
    "query": {
        "match_phrase":{
            "title": "drug resistance"
        }
    }
}

If you want to allow some leniency on how close the two words can appear together you can provide the slope parameter. A slope of 2 would still count documents which have not more than 2 words in between 'drug' and 'resistance'. For infinite slope one would recover the match query.

Besides making sure that the papers have the specific term, we also want to limit the publication date to a certain range. We, therefore, create a bool-must query like this

{
    "query": {
        "bool": {
            "must": [{
                "range": {
                    "arxiv_date":{
                        "gte" : low_date.strftime('%Y-%m-%d'), 
                        "lte" : up_date.strftime('%Y-%m-%d'), 
                        "format": "yyyy-MM-dd"
                    }
                }
            }]
        }
    }
}

All conditions in the bool-must query must be true. The must query in this example is a list, as shown by the [] brackets, and we can just append the match_phrase query from above.

However, in our case, we just want to find papers which have the search term in the title OR abstract. So we do not want to put this directly in a bool-must query since this would require the term to be in the title AND abstract. Instead, we create a bool-should query within the bool-must query. The bool-should query requires only one of the conditions within it to be true. Here is the final code.

def get_doc(low_date, up_date, list_of_terms=[]):
    term_doc = []
    for sub_string in list_of_terms:
        term_doc.append({
            "match_phrase":{
                "title": sub_string
            }
        })
        term_doc.append({
            "match_phrase":{
                "abstract": sub_string
            }
        })
    doc = {
        "query": {
            "bool": {
                "must": [{
                    "range": {
                        "created_date":{
                            "gte" : low_date.strftime('%Y-%m-%d'), 
                            "lte" : up_date.strftime('%Y-%m-%d'), 
                            "format": "yyyy-MM-dd"
                        }
                    }
                },
                {
                    'bool': {
                        "should": term_doc
                    }
                }]
            }
        }
    }
    return doc

This function expects a list of terms so that we can search for synonyms like ['blood cancer', 'leukemia']. For each term, we create a match_phrase query for the title and abstract fields and include this list to the bool-should query.

Now, all we have to do is to call this query for different time steps and different terms and plot it. Here is the code which loops over the last 25 years of PubMed publications and counts the number of papers. The number of papers is always normalized to the total number of papers in the same time frame.

def get_paper_count(list_of_terms, timestep):
    # We start 25 years in the past
    start_date = datetime.datetime.utcnow() - datetime.timedelta(days=365*25)

    list_of_counts = []
    list_of_dates = []
    low_date = start_date
    # loop through the data year by year
    while low_date < datetime.datetime.utcnow() - datetime.timedelta(days=10):
        up_date = low_date + datetime.timedelta(timestep)

        doc = get_doc(low_date, up_date)
        # we are only interested in the count -> size=0
        res = es.search(index=index_name, size=0, body=doc)
        norm = res['hits']['total']

        doc = get_doc(low_date, up_date, list_of_terms)
        # we are only interested in the count -> size=0
        res = es.search(index=index_name, size=0, body=doc)

        # norm should always be >0 but just in case   
        if norm > 0:
            list_of_counts.append(100.*float(res['hits']['total'])/float(norm))
            list_of_dates.append(low_date + datetime.timedelta(days=timestep/2))
        else:
            list_of_counts.append(0.)
            list_of_dates.append(low_date + datetime.timedelta(days=timestep/2))

        low_date = low_date + datetime.timedelta(timestep)
    return list_of_counts, list_of_dates 

Using this function we can create the plot I showed at the beginning of the page like this

def create_trending_plot():
    timestep = 365 # average over 365 days
    
    # Get a generic list of colors
    colors = dict(mcolors.BASE_COLORS, **mcolors.CSS4_COLORS)
    colors = [color[0] for color in colors.items()]
    # get all possible line styles
    linestyles = ['-', '--', '-.', ':']

    list_of_queries = [['prostate cancer'], 

                       ['blood cancer', 'leukemia'], 

                       ['Ebola'],
                       ['alzheimer', 'dementia']]
    timestamp = datetime.datetime.utcnow()

    plt.clf()
    # The maximum number of terms is given by the available colors
    for i, list_of_terms in enumerate(list_of_queries[:len(colors)]):
        print "i = ", i, "term = ", list_of_terms
        list_of_counts, list_of_dates = get_paper_count(list_of_terms, timestep)
        plt.plot(list_of_dates, 

                 list_of_counts, 

                 color=colors[i], 

                 label=', '.join(list_of_terms),
                 linestyle=linestyles[i%len(linestyles)])
    plt.xlabel('Date [in steps of %d days]' % timestep)
    plt.title('Relative number of papers for topic vs. time')
    plt.ylabel('Relative number of papers [%]')
    plt.legend(loc='upper left', prop={'size': 7})
    plt.savefig(OUTPUT_FOLDER + "trending_pubmed_%s.png" % 

                timestamp.strftime('%m-%d-%Y-%H-%M-%f'))
    return 

I also created a small web application using this code www.benty-fields.com/trending. The entire exercise can also be downloaded on GitHub. Let me know if you have any comments/questions below.
cheers
Florian

Elasticsearch and pubmed: Step 2 parsing the data into Elasticsearch

This is the second blog post discussing my project with the PubMed dataset. In the last post, I explained how to get a local copy of all the PubMed data. Here I will describe how I parse the data into my elasticsearch database.

The installation guidelines for an elasticsearch server are given here:
1. First download the latest version, in my case 5.5.2

    curl -L -O https://artifacts.elastic.co/downloads/elasticsearch/elasticsearch-5.5.2.tar.gz

2. Untar the file you just downloaded

    tar -xvf elasticsearch-5.5.2.tar.gz

3. Enter into the folder which was created by the previous command

    cd elasticsearch-5.5.2/bin

4. Now you can start the elasticsearch server with ./elasticsearch. In general, it is a good idea to set the swap and heap size using the corresponding environment variable. The recommendation is to use half of your available memory, which in my case is 4Gb, but I run other stuff on my machine and therefore use only 1Gb

    ES_JAVA_OPTS="-Xms1g -Xmx1g" ./elasticsearch

5. Now we need a python interface for elasticsearch. Of course, it is possible to directly interact with the elasticsearch server using curl (e.g. pycurl), but there are some interfaces which allow you to get away from the rather messy elasticsearch syntax. To be honest none of the available python packages totally convinced me, but the current standard seems to be the one which is also called elasticsearch (see e.g. here for more details)

    pip install elasticsearch

With this we can contact the elasticsearch server within our python program by including the following lines

from elasticsearch import Elasticsearch
es = Elasticsearch(hosts=['localhost:9200'])

Now we have to create an elasticsearch index (in other database systems this is called a table). This index describes the structure in which our data will be stored. I decided to store the data in 5 shards, which means that the data will be divided into 5 index files, each of which can be stored on a different machine. For now, I am not using any replicas (copies of shards). Replicas increase redundancy which allows retrieving data even if one or multiple servers are down. The decision how many shards to use is quite important since it cannot easily be changed in the future, while replicas can be added at any time.

Below I posted the function which creates the index. The settings option is used to set the number of shards and replicas together with some customised filters which I will describe in a second. The mapping describes the data vector and how it will be stored. Here I will keep things simple and only store the title, abstract and creation_date of the papers, even though the .gz files we downloaded contain much more information.

def create_pubmed_paper_index():    
    settings = {
        # changing the number of shards after the fact is not 

        # possible max Gb per shard should be 30Gb, replicas can 

        # be produced anytime
        # https://qbox.io/blog/optimizing-elasticsearch-how-many-shards-per-index
        "number_of_shards" : 5,
        "number_of_replicas": 0
    }
    mappings = {
        "pubmed-paper": {
            "properties" : {
                "title": { "type": "string", "analyzer": "standard"},
                "abstract": { "type": "string", "analyzer": "standard"},
                "created_date": {
                    "type":   "date",
                    "format": "yyyy-MM-dd"
                }
            }
        }
    }
    es.indices.delete(index=index_name, ignore=[400, 404])
    es.indices.create(index=index_name, 

                      body={ 'settings': settings,

                             'mappings': mappings }, 

                      request_timeout=30)
    return 

Elasticsearch is a non SQL database, meaning it does not follow the SQL syntax or functionality. If you are used to SQL databases, this will require some re-thinking. For example there is no way to join indices easily, which is a fundamental principle of SQL. This limits what you can do with elasticsearch significantly.

However, joins are very slow, and elasticsearch is all about speed. If you think about it, there is almost always a way around joins, all what you have to do is to store the data in the way you want to retrieve it. This can sometimes be ugly and requires a lot of memory to store redundant data, but without having to perform joins at runtime, it can be very fast.

Another way how elasticsearch saves time, is by processing the document directly when indexed. I used the standard analyzer, which lower cases and tokenizes all words. So for example the sentence

"The 2 QUICK Brown-Foxes jumped over the lazy dog's bone."

would be stored as

[ the, 2, quick, brown, foxes, jumped, over, the, lazy, dog's, bone ]

Again, this will speed up the required processing steps at runtime.


The function above creates the elasticsearch index. Now we have to read the data into this index. For that we write a small function, which unzips the files and builds an Element tree using xml.etree.cElementTree. Note that building such a tree can quickly lead to memory issues, since our .gz files are several Gb in size. So you should stay away from the often used parse() function which would eventually load the entire file into memory. Below I use the iterparse() function, which allows us to discard the elements after we have written them to the database.

import xml.etree.cElementTree as ET # C implementation of ElementTree
def fill_pubmed_papers_table(list_of_files):
    # Loop over all files, extract the information and index in bulk
    for i, f in enumerate(list_of_files):
        print("Read file %d filename = %s" % (i, f))

        time0 = time.time()
        time1 = time.time()
        inF = gzip.open(f, 'rb')
        # we have to iterate through the subtrees, ET.parse() would result
        # in memory issues
        context = ET.iterparse(inF, events=("start", "end"))
        # turn it into an iterator
        context = iter(context)

        # get the root element
        event, root = context.next()
        print("Preparing the file: %0.4fsec" % ((time.time() - time1)))

        time1 = time.time()

        documents = []
        time1 = time.time()
        for event, elem in context:
            if event == "end" and elem.tag == "PubmedArticle":
                doc, source = extract_data(elem)
                documents.append(doc)
                documents.append(source)
                elem.clear()
        root.clear()
        print("Extracting the file information: %0.4fsec" % 

              ((time.time() - time1)))
        time1 = time.time()

        res = es.bulk(index=index_name, body=documents, request_timeout=300)
        print("Indexing data: %0.4fsec" % ((time.time() - time1)))

        print("Total time spend on this file: %0.4fsec\n" % 

             ((time.time() - time0)))
        os.remove(f) # we directly remove all processed files
    return 

This function is looking for elements with the tag PubmedArticle, which we pass on to the extract_data() function. In that function, we extract the information we need. To write such a function we need to know the internal structure of the pubmed xml files. To get an idea how that structure might look like, you could print out one element using

def prettify(elem):
    from bs4 import BeautifulSoup # just for prettify
    '''Return a pretty-printed XML string for the Element.'''
    return BeautifulSoup(ET.tostring(elem, 'utf-8'), "xml").prettify()

I am using BeautifulSoup to produce a readable output since the equivalent functionality in cElementTree doesn't look as nice.

Without going into any more detail, here is the function which can extract the relevant information and store it in a class element

def extract_data(citation):
    new_pubmed_paper = Pubmed_paper()

    citation = citation.find('MedlineCitation')

    new_pubmed_paper.pm_id = citation.find('PMID').text
    new_pubmed_paper.title = citation.find('Article/ArticleTitle').text

    Abstract = citation.find('Article/Abstract')
    if Abstract is not None:
        # Here we discart information about objectives, design, 

        # results and conclusion etc.
        for text in Abstract.findall('AbstractText'):
            if text.text:
                if text.get('Label'):
                    new_pubmed_paper.abstract += '<b>' + text.get('Label') + '</b>: '
                new_pubmed_paper.abstract += text.text + '<br>'

    DateCreated = citation.find('DateCreated')
    new_pubmed_paper.created_datetime = datetime.datetime(

        int(DateCreated.find('Year').text),
        int(DateCreated.find('Month').text),
        int(DateCreated.find('Day').text)

    )
    doc, source = get_es_docs(new_pubmed_paper)
    del new_pubmed_paper
    return doc, source

where the class Pubmed_paper() is defined as

class Pubmed_paper():
    ''' Used to temporarily store a pubmed paper outside es '''
    def __init__(self):
        self.pm_id = 0

        # every paper has a created_date
        self.created_datetime = datetime.datetime.today()
        self.title = ""
        self.abstract = ""

    def __repr__(self):
        return '<Pubmed_paper %r>' % (self.pm_id)

and the function which writes the doc and source dictionaries is

def get_es_docs(paper):
    source = {
        'title': paper.title,
        'created_date': paper.created_datetime.date(),
        'abstract': paper.abstract
    }
    doc = {
        "index": {
            "_index": index_name,
            "_type": type_name,
            "_id": paper.pm_id
        }
    }
    return doc, source

To read all the files into the database will take a few hours.

This post was a bit code heavy, but now that we have written the entire dataset into elasticsearch, we can easily access it. In the next post we will start a small project using this large dataset and making use of the fast elasticsearch database. You can access the code used in this post at GitHub. Let me know if you have any comments/questions below.
cheers
Florian