motifs_stats_energy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "\n",
    "\n",
    "from scipy.optimize import curve_fit\n",
    "from config.sqlcollect import getTable\n",
    "#from config.network import graphCreation\n",
    "from config.cubes import makeCubes\n",
    "from config.seismicZones import query\n",
    "from config.network import graphCreation3\n",
    "\n",
    "\n",
    "from motifs_areas_volumes_energy import*\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# For which region networks do you want to analyze motifs ?\n",
    "region = input('Input region : vrancea / romania / california / italy / japan : ')\n",
    "\n",
    "# What motifs do you want to discover ?\n",
    "motif = input('Input motif: Triangles (for areas) / Tetrahedrons (for volumes) : ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not os.path.exists(f'./results/{region}/motifStatistics'):\n",
    "    os.makedirs(f'results/{region}/motifStatistics')\n",
    "\n",
    "\n",
    "sql_query= query(region)\n",
    "plt.rcParams['text.usetex'] = True\n",
    "\n",
    "for mag in (1,2,3):\n",
    "    \n",
    "    # Magnitude windows for the condition that collects the database through mySQL\n",
    "    sql_query+=f\" AND magnitude>={mag}\"\n",
    "\n",
    "    # Collect the database and create the graph\n",
    "    # Collect the earthquakes\n",
    "    quakes = getTable(sql_query)\n",
    "\n",
    "\n",
    "    for side in (5,10):\n",
    "\n",
    "        # Make the Cubes\n",
    "        quakes = makeCubes(quakes,region,side,energyRelease=True)\n",
    "\n",
    "        # Get maximum magnitude for displaying on the plot\n",
    "        realmagMax = max(quakes['magnitude'])\n",
    "\n",
    "\n",
    "        # Make the network\n",
    "        quakesGraph = graphCreation3(quakes)\n",
    "\n",
    "\n",
    "        # Get the motif nodes and the energy in motifs    \n",
    "        motifNodes,energyMotif = totalMeanEnergyMotif(region,side,mag,motif,quakesGraph,quakes)\n",
    "        \n",
    "        # TRIANGLES \n",
    "        if motif == 'Triangles':\n",
    "            motifTotalEnergy,motifMeanEnergy = areasInTriangles(motifNodes,energyMotif,quakesGraph,quakes)\n",
    "        # TETRAHEDRONS\n",
    "        else:\n",
    "            motifTotalEnergy,motifMeanEnergy = volumesInSquares(motifNodes,energyMotif,quakesGraph,quakes)\n",
    "\n",
    "\n",
    "        # TOTAL ENERGY / MEAN ENERGY PLOTS\n",
    "        for motifstats in [motifTotalEnergy,motifMeanEnergy]:\n",
    "\n",
    "            # hist and bins for the plots\n",
    "            hist, bins = np.histogram(motifstats,bins=round(np.sqrt(len(motifstats))))\n",
    "        \n",
    "            # Create the x as hist with zeros, force into floats ! \n",
    "            x = np.zeros_like(hist.astype(float))\n",
    "            for i in range(1,len(bins)):\n",
    "                x[i-1]=((bins[i]+bins[i-1])/2)\n",
    "\n",
    "            for i in range(len(hist)):\n",
    "                if hist[i]==0:\n",
    "                    y=np.array(hist[:i])\n",
    "                    x=np.array(x[:i])\n",
    "                    break\n",
    "                else:\n",
    "                    y=hist\n",
    "\n",
    "            # normalize y axis \n",
    "            y_norm = [float(i)/sum(y) for i in y]\n",
    "\n",
    "            # the curve fit of the data\n",
    "            pars, cov = curve_fit(f=power_law,xdata=x,ydata=y_norm,maxfev=5000)\n",
    "\n",
    "            # Compute the chi_squared goodness of fit = sum( ( observed - expected )^2 / expected )\n",
    "            chi_squared = np.sum((y_norm-power_law(x,*pars))**2/power_law(x,*pars))\n",
    "            \n",
    "\n",
    "            # Plot the results\n",
    "            fig = plt.figure(figsize=(15,10))\n",
    "            ax = fig.add_subplot(111)\n",
    "\n",
    "            # The data, scattered\n",
    "            ax.scatter(x,y_norm)\n",
    "            ax.set_xscale('log')\n",
    "            ax.set_yscale('log')\n",
    "\n",
    "            # The fit\n",
    "            ax.plot(x,power_law(x,*pars),\n",
    "                        label=f'$\\gamma$ = {np.round(pars[1],4)}\\n$\\chi^2$ = {np.round(chi_squared,4)}',\n",
    "                        color='red')\n",
    "\n",
    "            # Legend : gamma coefficient of fit and chi_squared goodness of fit\n",
    "            plt.legend(loc='upper right',fontsize=26, frameon=True)\n",
    "\n",
    "            #r'$\\mathbf{Longitude}$'\n",
    "\n",
    "            \n",
    "            plt.setp(ax.get_xticklabels(), fontsize=26)\n",
    "            plt.setp(ax.get_yticklabels(), fontsize=26)\n",
    "\n",
    "\n",
    "            ax.tick_params(axis='both', which='major', pad=7)\n",
    "\n",
    "            # Title of connectivity distribution ( data + fit )\n",
    "            #ax.set_title(f'cube size = {side} km ', fontsize=26, fontweight='bold')\n",
    "\n",
    "            # Magnitude range:\n",
    "            ax.text(0.05,0.1, f'Magnitude Range:\\n {mag}'+r'$<$' + 'mag' + r'$<$' +f'{realmagMax}', fontsize=20,\n",
    "                    bbox=dict(facecolor='white', boxstyle='round', alpha=0.3),transform=ax.transAxes) \n",
    "\n",
    "            # TOTAL MAGNITUDE IN MOTIF\n",
    "            if motifstats == motifTotalEnergy:\n",
    "\n",
    "                # TRIANGLES LABELS\n",
    "                if motif == \"Triangles\":\n",
    "                    ax.set_xlabel(r'$\\mathbf{S_{TE}}$', fontsize=26,labelpad=20)\n",
    "                    ax.set_ylabel(r'$\\mathbf{P(S_{TE})}$', fontsize=26,labelpad=20)\n",
    "                # TETRAHEDRONS LABELS\n",
    "                else:\n",
    "                    ax.set_xlabel(r'$\\mathbf{V_{TE}}$', fontsize=26,labelpad=20)\n",
    "                    ax.set_ylabel(r'$\\mathbf{P(V_{TE})}$', fontsize=26,labelpad=20)\n",
    "\n",
    "\n",
    "                fig.savefig(f'./results/{region}/motifStatistics/motif{motif}_TotalEnergy_network{region}_{side}km_{mag}mag.png',\n",
    "                            facecolor='white', transparent=False, bbox_inches='tight')           \n",
    "\n",
    "            # MEAN MAGNITUDE IN MOTIF\n",
    "            else:\n",
    "\n",
    "                # TRIANGLE LABELS\n",
    "                if motif == \"Triangles\":\n",
    "                    plt.xlabel(r'$\\mathbf{S_{ME}}$', fontsize=26,labelpad=20)\n",
    "                    plt.ylabel(r'$\\mathbf{P(S_{ME})}$', fontsize=26,labelpad=20)\n",
    "                # TETRAHEDRONS LABELS\n",
    "                else:\n",
    "                    plt.xlabel(r'$\\mathbf{V_{ME}}$', fontsize=26,labelpad=20)\n",
    "                    plt.ylabel(r'$\\mathbf{P(V_{ME})}$', fontsize=26,labelpad=20)\n",
    "\n",
    "\n",
    "                fig.savefig(f'./results/{region}/motifStatistics/motif{motif}_MeanEnergy_network{region}_{side}km_{mag}mag.png',\n",
    "                                facecolor='white', transparent=False, bbox_inches='tight')\n",
    "\n",
    "            plt.close();\n",
    "\n",
    "        # Extract the magnitude restrictions from the condition \n",
    "        sql_query = sql_query.replace(f\" AND magnitude>={mag}\", '');\n",
    "\n",
    "        # Use break to get only the 5km sides\n",
    "        break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "3885b2b68169bbcc5c24138d8a594f3460f3ddac7a4132a4d1d891044f50ddd3"
  },
  "kernelspec": {
   "display_name": "Python 3.8.10 ('venv': venv)",
   "language": "python",
   "name": "python3"
  },
  "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.8.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}