Merge pull request #28 from IAA-CSIC/dev

Dev: add Flare analysis to MUTIS

docs/index.rst

@@ -12,6 +12,7 @@ A Python package for the analysis of correlations of light curves and their stat
 
    recipes/mutis_workflow.ipynb
    recipes/Synthethic_generation_of_light_curves.ipynb
+   recipes/Flare_analysis.ipynb
    recipes/PSD_tests.ipynb
    recipes/sin_tests.ipynb
    recipes/stochastics_tests.ipynb
@@ -23,4 +24,5 @@ A Python package for the analysis of correlations of light curves and their stat
 
    reference/correlation
    reference/signal
+   reference/flares
    reference/lib

docs/recipes/Flare_analysis.ipynb

@@ -0,0 +1,127 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Flare analysis"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%load_ext autoreload\n",
+    "%autoreload 2\n",
+    "\n",
+    "import numpy as np\n",
+    "import scipy as sp\n",
+    "import pandas as pd\n",
+    "\n",
+    "import matplotlib as mplt\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "import astropy\n",
+    "import astropy.units as u\n",
+    "\n",
+    "import mutis\n",
+    "from mutis.flares import BayesianBlocks"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Load data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "data = pd.read_csv('data/gamma-I.dat', comment='!')            \n",
+    "data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Build signal and plot its Bayesian Block representation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sig = mutis.Signal(data['jyear'], 1e6*data['CFlux'], 1e6*data['CFluxErr'])\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(10,6))\n",
+    "bayblocks = BayesianBlocks(sig, p=1e-8)\n",
+    "bayblocks.get_flares()\n",
+    "bayblocks.signal.plot()\n",
+    "bayblocks.plot()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Get flares from BayesianBlocks, and plot them"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig, ax = plt.subplots(figsize=(10,6))\n",
+    "bayblocks.signal.plot()\n",
+    "\n",
+    "bayblocks.get_flares()\n",
+    "\n",
+    "bayblocks.plot('color.strict')\n",
+    "\n",
+    "plt.axhline(y=np.median(sig.values), color='r', linestyle='--', alpha=1, label='1 MAD')\n",
+    "plt.axhline(y=np.median(sig.values)+np.sqrt(np.median((sig.values-np.median(sig.values))**2)), color='m', linestyle='--', alpha=1, label='Mad')\n",
+    "\n",
+    "for flare in bayblocks.get_flare_list():\n",
+    "    flare.plot()\n",
+    "\n",
+    "plt.legend()\n",
+    "    \n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.1"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

docs/reference/flares.rst

@@ -0,0 +1,3 @@
+.. automodapi:: mutis.flares
+    :no-inheritance-diagram:
+    :include-all-objects:

mutis/flares/__init__.py

@@ -0,0 +1,5 @@
+# Licensed under a 3-clause BSD style license - see LICENSE
+""" Flare analysis """
+
+from .bayblocks import BayesianBlocks
+from .flare import Flare

mutis/flares/bayblocks.py

@@ -0,0 +1,232 @@
+""" Flare analysis """
+
+import logging
+
+import numpy as np
+import scipy as sp
+import pandas as pd
+
+import matplotlib as mplt
+import matplotlib.pyplot as plt
+
+import sys
+import os
+import glob
+
+import re
+from datetime import datetime
+
+from astropy.stats import bayesian_blocks
+
+import warnings 
+
+import mutis.flares.flare as flare
+
+log = logging.getLogger(__name__)
+
+class BayesianBlocks:
+    """ Return a Bayesian Block representation of the signal
+    
+    Attributes
+    ----------
+    edges: list
+            list of edges defining the blocks.
+    values: list
+            list of values defining the height of the blocks.
+    signal: mutis.Signal
+            the mutis.Signal() instance that this BayesianBlocks() represents.
+    inflare: list
+            list of boolean representing whether a block is in flare or not.
+    """
+
+    def __init__(self, signal, p=0.1):
+        edges = bayesian_blocks(signal.times, signal.values, signal.dvalues, fitness='measures', p0=p)
+
+        values = list()
+        for i in range( len(edges)-1 ):
+            msk = (edges[i] < signal.times) & (signal.times < edges[i+1])
+            if np.sum(msk) > 0:
+                value = np.average( signal.values[msk], weights=1/signal.dvalues[msk] )
+                values.append( value )
+            else:
+                values.append( 0 )
+
+        self.edges = np.array(edges)
+        self.values = np.array(values)
+        self.signal = signal
+
+        self.inflare = None
+
+    def plot(self, style='color.strict', ax=None, **kwargs):
+        """ Plot the bayesian block representation
+        
+        E.g:
+        >>> signal = mutis.Signal(data['jyear'], data['Flux'], data['Flux_err'])
+        >>> BayesianBlocks(signal).plot(color='k')
+
+        """
+
+        ax = plt.gca() if ax is None else ax
+
+        #self.signal.plot()
+
+        x = self.edges
+        y = self.values[np.r_[0:len(self.values),-1]]
+
+
+        if self.inflare is None or style == 'none':
+
+            ax.step(x, y, 'k', where='post', **kwargs)
+
+        else:
+
+            if style == 'area':
+
+                ax.step(x, y, 'k', where='post')
+                for i in range(len(self.inflare)):
+                    if self.inflare[i]:
+                        ax.axvspan(x[i], x[i+1], facecolor='r', edgecolor=None, alpha=0.2, **kwargs)
+
+            elif style == 'color.simple':
+
+                for i in range(len(self.edges)-2):
+                    ax.step(self.edges[[i,i+1]], self.values[[i,i+1]], 'r' if self.inflare[i] else 'k', where='post', **kwargs)
+                ax.step(self.edges[[-2,-1]], self.values[[-1,-1]], 'r' if self.inflare[-1] else 'k', where='post', **kwargs)
+
+            elif style == 'color.loose':
+
+                for i in range(len(self.edges)-2):
+                    ax.step(self.edges[[i,i+1]], self.values[[i,i]], 'r' if self.inflare[i] else 'k', where='post', **kwargs)
+                    ax.step(self.edges[[i+1,i+1]], self.values[[i,i+1]], 'r' if self.inflare[i] or self.inflare[i+1] else 'k', where='post', **kwargs)
+                ax.step(self.edges[[-2,-1]], self.values[[-1,-1]], 'r' if self.inflare[-1] else 'k', where='post', **kwargs)
+
+            elif style == 'color.strict':
+
+                for i in range(len(self.edges)-2):
+                    ax.step(self.edges[[i,i+1]], self.values[[i,i]], 'r' if self.inflare[i] else 'k', where='post', **kwargs)
+                    ax.step(self.edges[[i+1,i+1]], self.values[[i,i+1]], 'r' if self.inflare[i] and self.inflare[i+1] else 'k', where='post', **kwargs)
+                ax.step(self.edges[[-2,-1]], self.values[[-1,-1]], 'r' if self.inflare[-1] else 'k', where='post', **kwargs)
+
+            else:
+                raise Exception("Unknown style. Available styles are 'none', 'area', 'color.simple', 'color.strict', 'color.loose'.")
+
+
+
+    def get_flares(self, thresh=1):
+        """ Get a list of flares following the algorithm proposed in 
+[Meyer, Scargle, Blandford (2019)]
+(https://iopscience.iop.org/article/10.3847/1538-4357/ab1651/pdf):
+        
+```There is no generally accepted consensus on the best way to
+determine which data points belong to a flaring state and which
+characterize the quiescent level. Nalewajko (2013)suggested
+the simple definition that a flare is a continuous time interval
+associated with a flux peak in which the flux is larger than half
+the peak flux value. This definition is intuitive, however, and it
+is unclear how to treat overlapping flares and identify flux
+peaks in an objective way. Here we use a simple two-step
+procedure tailored to the block representation: (1)identify a
+block that is higher than both the previous and subsequent
+blocks as a peak, and (2)proceed downward from the peak in
+both directions as long as the blocks are successively lower.```
+        """
+        di_max = 5
+
+        # get list of local maxima
+        # bad beheaviour at 0 and -1:
+        #imaxL = sp.signal.argrelextrema(self.values, np.greater)[0]
+        imaxL = list()
+        for i in range(0,len(self.values)):
+            if i == 0 and thresh*self.values[0] > self.values[1]:
+                imaxL.append(i)
+            elif i == len(self.values)-1 and self.values[i-1] < thresh*self.values[i]:
+                imaxL.append(i)
+            elif self.values[i-1] < thresh*self.values[i] > self.values[i+1]:
+                imaxL.append(i)
+        imaxL = np.array(imaxL)
+
+
+        inflareL = np.full(self.values.shape, False)
+
+        # get list of local maxima over a threshold (flare peaks)
+        for imax in imaxL:
+            inflare = True
+
+            if imax == 0:
+                pass
+            elif not (self.values[imax-1] < thresh*self.values[imax]):
+                inflare = False
+
+
+            if imax == len(self.values)-1:
+                pass
+            elif  not (thresh*self.values[imax] > self.values[imax+1]):
+                inflare = False
+
+            inflareL[imax] = inflare
+
+        # extend the flare to adyacent blocks
+        for imax in np.argwhere(inflareL):
+
+            if imax == 0:
+                pass
+            else:
+                di = 1
+                stop = False
+                while not stop and di < di_max:
+                    if imax-di-1 < 0:
+                        break
+
+                    if self.values[imax-di-1] < thresh*self.values[imax-di]:
+                        inflareL[imax-di] = True
+                    else:
+                        stop = True
+
+
+                    di = di + 1
+
+            if imax == len(self.values)-1:
+                pass
+            else:
+                di = 1
+                stop = False
+                while not stop and di < di_max:
+                    if imax+di+1 > len(self.values)-1:
+                        break
+
+                    if thresh*self.values[imax+di] > self.values[imax+di+1]:
+                        inflareL[imax+di] = True
+                    else:
+                        stop = True
+
+
+                    di = di + 1
+
+
+        self.inflare = np.array(inflareL)
+
+
+    def get_flare_list(self):
+        """ Join all flares into a list of mutis.flares.Flare() """
+
+        groups = np.split(np.arange(len(self.inflare)), 
+                          np.where((np.abs(np.diff(np.asarray(self.inflare, dtype=int))) == 1))[0]+1)
+
+        # (groups contains the groups of indices of self.inflare/self.values 
+        # that corresond to flares or no flares)
+
+        #print(self.inflare)
+        #print([list(grp) for grp in groups])
+
+        flareL = list()
+        for i, grp in enumerate(groups):
+            if np.all(self.inflare[grp]):
+                #print(grp)
+                tstart = self.edges[grp[0]]
+                tstop = self.edges[grp[-1]+1] # flare ends when no flare begins
+
+                flareL.append(flare.Flare(tstart,tstop))
+
+        #print(flareL)
+
+        return flareL