#!/usr/bin/env python
# encoding: utf-8
"""
catalogs.py

Scripts to generate a coherent earthquake catalog for Oklahoma.  This catalog is derived from multiple
sources.  Each source is tracked with the data.  The data is stored in a postgis database.

Created by Austin Holland on 2012-07-12.
Copyright (c) 2012 __Oklahoma Geological Survey__. All rights reserved.
"""

import sys
sys.path.append("/Users/aholland/lib/")
import os
import shapely.geometry as sg
from osgeo import ogr
from sqlalchemy import *
from sqlalchemy.orm import *
from sqlalchemy.ext.declarative import declarative_base 
from geoalchemy import *
from datetime import *
from seisutils.eqspatial import inOK
import psycopg2
import mplstereonet
import seisutils.fmech as fmec
import numpy as np
from scipy.special import erfinv
import matplotlib.pyplot as plt
#sys.path.append("/Users/aholland/proj/hazard.d/mtoolkit/")
from mtoolkit.scientific.declustering import *
from mtoolkit.scientific.completeness import *
from mtoolkit.scientific.recurrence import *

# Global definitions
engine=create_engine('postgresql://analyst@wichita.ogs.ou.edu/analyst', echo=False)
Session = sessionmaker(bind=engine)
session = Session()
Base=declarative_base()

# File locations
SHPcounties="/Users/aholland/GIS/Geography/OKCNTYD/OKCNTYD.shp"


class EQ(Base):
  """SQLAlchemy/PostGIS earthquake dataclass and table definitions
  source_code identifies different data sources and should be a varchar(4) specifying the agency 
    OGS - Oklahoma Geological Survey
    OGSD - OGS Derived
    USGS - US Geological Survey
    ANSS - Advanced National Seismic System
    
  database analyst
  tablename proc_catalog
    id  integer primary key
    origintime datetime (fully qualified datetime if available)
    latitude float
    longitude float
    depth float
    err_lon float
    err_lat float
    err_depth float
    err_origintime float
    county string
    origin_src source_code
    m  float
    m_uncert float
    m_src source code
    mw float
    mw_src source_code
    mblg float
    mblg_src source_code
    ml float
    ml souce_code
    m3hz float
    m3hz_src source_code
    md float
    md_src source_code
    max_mmi integer
    mmi_src source_code
    geom well known binary (GeoSpatial)
    flag integer (-1 foreshock, 0 mainshock, 1 aftershock)
    """
  __tablename__='proc_catalog'
  id=Column(Integer, primary_key=True)
  origintime=Column(DateTime, nullable=False)
  latitude=Column(Float)
  longitude=Column(Float)
  depth=Column(Float)
  err_lon=Column(Float)
  err_lat=Column(Float)
  err_depth=Column(Float)
  err_origintime=Column(Float)
  county=Column(String)
  origin_src=Column(String(4))
  m=Column(Float)
  m_uncert=Column(Float)
  m_src=Column(String(4))
  mw=Column(Float)
  mw_src=Column(String(4))
  mblg_ogs=Column(Float)
  mblg_usgs=Column(Float)
  ml_ogs=Column(Float)
  m3hz_ogs=Column(Float)
  md_ogs=Column(Float)
  mb=Column(Float)
  ms=Column(Float)
  mfa=Column(Float)
  max_mmi=Column(Integer)
  mmi_src=Column(String(4))
  geom=GeometryColumn(Point(2))
  flag=Column(Integer)
    
  def __init__(self,ot,lat,lon,depth,src,epicenter):
    """obj=EQ(origintime,latitude,longitude,depth,origin_src)
    Other attributes of the object must then be specified manually"""
    self.origintime=ot
    self.latitude=lat
    self.longitude=lon
    self.depth=depth
    self.origin_src=src
    self.geom=WKTSpatialElement(epicenter.wkt)
#    self.max_mmi=None

  def usgs_magproc(self,line,includeTUL=False):
    """Process magnitudes form usgs expanded format if more than one mw take the weighted moment"""
    mag=[]
    mtype=[]
    msrc=[]
    mw=[]
    mw_uncert=[]
    if line[80:83]!='   ':
      if includeTUL:
        mag.append(float(line[75:78]))
        mtype.append(line[78:80])
        msrc.append(line[80:84])
      elif line[80:83]!='TUL':
        mag.append(float(line[75:78]))
        mtype.append(line[78:80])
        msrc.append(line[80:84])        
    if line[91:94]!='   ':
      if includeTUL:
        mag.append(float(line[86:89]))
        mtype.append(line[89:91])
        msrc.append(line[91:95])
      elif line[91:94]!='TUL':
        mag.append(float(line[86:89]))
        mtype.append(line[89:91])
        msrc.append(line[91:95])      
    if line[66:69]!='   ':
      self.mb=float(line[66:69])
    if line[70:73]!='   ':
      self.ms=float(line[70:73])
    for i in range(0,len(mag)):
      if mtype[i]=='Lg':
        if msrc!='TUL':
          self.mblg_usgs=mag[i]
        else:
          self.mblg_ogs=mag[i]
      if mtype[i]=='Md':
        if msrc=='TUL':
          self.md_ogs=mag[i]
      if mtype[i]=='MS' or mtype[i]=='Ms':
        self.ms=mag[i]
      if mtype[i]=='FA':
        self.mfa=mag[i]
      if mtype[i]=='Mw':
        mw.append(mag[i])
        mw_uncert.append(0.1)
        mw_src=msrc[i]
    if line[105]!='.':
      mmi=int(line[105])
      if hasattr(self,'max_mmi'):
        if mmi>self.max_mmi:
          self.max_mmi=mmi
          self.mmi_src='USGS'
      else:
        self.max_mmi=mmi    
    if len(mw)>1:
      self.mw,self.m_uncert=weighted_moment(np.array(mw),np.array(mw_uncert))
      self.mw_src='OGSD'
    if len(mw)==1:
      self.mw=mw[0]
      self.mw_src=mw_src

  def __repr__(self):
    return "<EQ(%s %s)" % (self.origintime,self.origin_src)


# Create the database tables so they can be populated and such
GeometryDDL(EQ.__table__)
Base.metadata.create_all(engine)     

def addcat_usgs(dts=10.):
  usgsfiles=['eqcatalogs.d/usgspde_expanded.txt','eqcatalogs.d/usgshis_expanded.txt']
  dt=timedelta(seconds=dts)
  for filename in usgsfiles:
    fh=open(filename,'r')

  
    for line in fh:
      oterr=0.
      if line[1:4]=='PDE' or line[1:4]=='USH': # Earthquake line process
#        try:
          lat=float(line[33:39])
          lon=float(line[40:48])
          if line[51:53]!='  ':
            z=float(line[51:53])
          county=inOK(lon,lat,SHPcounties,'DESCRIP')
          if county: # Skip event if not in Oklahoma
            year=int(line[9:13])
            month=int(line[14:16])
            day=int(line[17:19])
            hour=int(line[20:22])
            minute=int(line[22:24])
            if line[24:26]!='  ':
              sec=int(line[24:26].rstrip())
            else:
              sec=30
              oterr=30.
            if line[27:29]!='  ':
              dsec=float(line[26:29].rstrip())
              msec=int(dsec*1.e6)
            else:
              msec=0
            ot=datetime(year,month,day,hour,minute,sec,msec)
            mag=[]
            mtype=[]
            msrc=[]
            mmi=0
            res=session.query(EQ).filter(and_(EQ.origintime>ot-dt,EQ.origintime<ot+dt)).all()
            if len(res)==1:  #this earthquake is in the OGS database add magnitudes as appropriate
              #print "Add USGS magstuff \n\t",line
              res=res[0]
              res.usgs_magproc(line)
              session.add(res)    
            elif len(res)==0:  #add this earthquake it has no record in the OGS database
              print "Adding event from USGS PDE: \n\t",line

              res=EQ(ot,lat,lon,z,'USGS',sg.Point(lon,lat))
              res.err_origintime=oterr
              res.usgs_magproc(line,includeTUL=True)
              session.add(res)
            else:
              print "USGS ",ot," event matches multiple records"
#        except:
#            print "Trouble processing, line ignored\n",line
    fh.close()
  session.commit()


def EQfrom_ogscat(line):
  monthstr={'JAN':1,'FEB':2,'MAR':3,'APR':4,'MAY':5,'JUN':6,'JUL':7,'AUG':8,'SEP':9,'OCT':10,
    'NOV':11,'DEC':12}

  if line[12:14]=='  ' or len(line)<50:  #Ignore earthquakes with no hour how well do we really know them
    return None
  else:
    yr=int(line[0:4])
    mo=monthstr[line[5:8]]
    dy=int(line[9:11])
    hr=int(line[12:14])
    mn=int(line[15:17])
    sec=int(line[18:20])
    
    if line[21:23]!='  ':
      dsec=float(line[20:23].rstrip())
      msec=int(dsec*1.e6)
    else:
      msec=0
    ot=datetime(yr,mo,dy,hr,mn,sec,msec)
    
    lat=float(line[53:59])
    lon=float(line[61:69])
    z=float(line[72:76])
    county=inOK(lon,lat,SHPcounties,'DESCRIP')
   
    if county:
      eq=EQ(ot,lat,lon,z,'OGS',sg.Point(lon,lat))
      eq.county=county
      #print line[39]
      if line[39]=='F':
        eq.max_mmi=2
        eq.mmi_src="OGS"
      elif line[39]!=' ':
        eq.max_mmi=int(line[39])
        eq.mmi_src="OGS"

      if line[49:52]!='   ':
       eq.md_ogs=float(line[49:52])
      if line[45:48]!='   ':
        eq.mblg_ogs=float(line[45:48])   
      if line[41:44]!='   ':
        eq.m3hz_ogs=float(line[41:44])


      return eq
    else:
      print ot, lon, lat,
      return None


def populate_ogsearthquakes():  
  # Populate our earthquakes
  eqs=[]
#  fh=open("earthquakes.d/1897-2007-OGS_KVL.csv",'r')
#  fh.readline()
#  for line in fh:
#    res=EQfrom_kvlcsv(line)
#    if isinstance(res,EQE):
#      eqs.append(res)
#  fh.close()
  ogscats=['1897-1969','1970-1976','1977','1978','1979','1980','1981','1982','1983','1984','1985',
  '1986','1987','1988','1989','1990','1991','1992','1993','1994','1995','1996','1997','1998','1999',
  '2000','2001','2002','2003','2004','2005','2006','2007','2008']
  mmistr_dict={'I':1,'II':2,'III':3,'IV':4,'V':5,'VI':6,'VII':7,'VIII':8,'IX':9,'X':10}
  for name in ogscats:
    print name
    fh=open("eqcatalogs.d/ogs.d/"+name+".txt",'r')
    fh.readline()
    fh.readline()
    fh.readline()
    fh.readline()
    for line in fh:
      res=EQfrom_ogscat(line)
      if isinstance(res,EQ):
        eqs.append(res)
    fh.close()

  session.add_all(eqs)
  session.commit()
  # Add Data in yearly postgis catalogs
  eqs=[]
  try:
     conn = psycopg2.connect("dbname='analyst' user='analyst' host='wichita.ogs.ou.edu' password='OKPunt*db'");
     cur=conn.cursor()
  except:
      print "Unable to connect to database"
      sys.exit()
  years=[2010,2011,2012]
  for year in years:
    sql="SELECT origintime,latitude,longitude,depth,error_lat,error_lon,error_depth,error_ot,magnitude,magnitude_type,magnitude_agency,max_mmi,county FROM eqcat_%s order by origintime;" % (year)
    cur.execute(sql)
    rows=cur.fetchall()
    for row in rows:
      if row[12]: # If county is defined ie the eq occurred in Oklahoma
        eq=EQ(row[0],row[1],row[2],row[3],'OGS',sg.Point(row[2],row[1]))
        eq.county=row[12]
        eq.err_lat=row[4]
        eq.err_lon=row[5]
        eq.err_depth=row[6]
        eq.err_origintime=row[7]
        if row[11]:
          eq.max_mmi=mmistr_dict[row[11]]
        for i in range(0,len(row[8])):
          if row[10][i]=='OGS':
            if row[9][i]=='ML':
              eq.ml_ogs=row[8][i]
            if row[9][i]=="mbLg":
              eq.mblg_ogs=row[8][i]
            if row[9][i]=="Md":
              eq.md_ogs=row[8][i]
        eqs.append(eq)
    
  session.add_all(eqs)
  session.commit()
  return len(eqs)
  
def ml_mw_relationship():
  eqs=session.query(EQ).filter(EQ.ml_ogs!=None).filter(EQ.mw!=None).all()
  ml=[]
  mw=[]
  mw_avg=[]
  ml_xavg=[]
  for eq in eqs:   
    mw.append(eq.mw)
    ml.append(eq.ml_ogs)
    if eq.mw_src=='OGSD':
      ml_xavg.append(eq.ml_ogs)
      mw_avg.append(eq.mw)
  x=np.arange(0.,6.0,0.2)
  mw_miao=.817+.7603*x
  fig=plt.figure()
  ax = fig.add_subplot(111, aspect='equal')
  #print mw,ml
  plt.plot(ml,mw,'or',label='RMT')
  plt.plot(ml_xavg,mw_avg,'ok',label='RMT & CMT')
  plt.plot(x,x,'-b',label='Hanks & Kanamori')
  plt.plot(x,mw_miao,':b',label='Miao & Langston')
  plt.xlabel('$M_L$')
  plt.ylabel('$M_W$')
  plt.legend(loc='upper left')
  plt.savefig('figures.d/mw_ml_rel_published.pdf',bbox='tight')
  # Determine our own linear relationship
  import scipy.stats
  slope,intercept,r_value, p_value, slope_std_error=scipy.stats.linregress(ml,mw)
  print "R: ",r_value**2
  print "N=",len(ml)
  # Calculate some additional outputs
  y = intercept + slope * x
  predict_y = intercept + slope * ml
  pred_error = mw - predict_y
  degrees_of_freedom = len(ml) - 2
  residual_std_error = np.sqrt(np.sum(pred_error**2) / degrees_of_freedom)
  sigma=np.sqrt(np.var(pred_error))
  fig=plt.figure()
  ax = fig.add_subplot(111, aspect='equal')
  plt.plot(ml,mw,'or',label='RMT')
  plt.plot(ml_xavg,mw_avg,'ok',label='RMT & CMT')
  plt.plot(x,y,'-b')
  plt.text(0.5,5.,"$M=%.3f+%.3fM_L$" %(intercept,slope))
  plt.text(0.5,4.7,"$\sigma=%.3f$" % (sigma))
  plt.xlabel('$M_L$')
  plt.ylabel('$M_W$')
  plt.savefig('figures.d/mw_ml_rel.pdf',bbox='tight')
  return True
  
def plot_mblg():
  eqs=session.query(EQ).filter(EQ.mblg_ogs!=None).filter(EQ.mw!=None).all()
  mblg_ogs=[]
  mw_ogs=[]
  for eq in eqs:   
    mw_ogs.append(eq.mw)
    mblg_ogs.append(eq.mblg_ogs)
  eqs=session.query(EQ).filter(EQ.mblg_usgs!=None).filter(EQ.mw!=None).all()
  mblg_usgs=[]
  mw_usgs=[]
  for eq in eqs:   
    mw_usgs.append(eq.mw)
    mblg_usgs.append(eq.mblg_usgs) 
  x=np.arange(0.,6.0,0.2)
  mw_johnson=np.zeros(len(x))
  mw_ceusssc=np.zeros(len(x))
  for i in np.arange(0,len(x)):
    mw_johnson[i]=1.14+0.24*x[i]+0.093*(x[i]**2)
    mw_ceusssc[i]=x[i]-0.316
  print "Total #mblg ",len(mw_usgs)+len(mw_ogs)
  fig=plt.figure()
  ax = fig.add_subplot(111, aspect='equal')
  plt.plot(mblg_ogs,mw_ogs,'ok',label='$M_{bLg}$ OGS')
  plt.plot(mblg_usgs,mw_usgs,'or',label='$M_{bLg}$ USGS')
  plt.plot(x,mw_johnson,'-b',label="Johnson (1996a)")
  plt.plot(x,mw_ceusssc,':b',label="CEUS-SSC (2012)")
  plt.xlabel('$M_{bLg}$')
  plt.ylabel('$M_W$')
  plt.legend(loc='upper left')
  plt.savefig('figures.d/mblg_rel.pdf',bbox='tight')
  
  eqs=session.query(EQ).filter(EQ.mblg_ogs!=None).filter(EQ.ml_ogs!=None).all()
  mblg_ogs=[]
  mw_ogs=[]
  for eq in eqs:   
    mw_ogs.append(eq.ml_ogs*1.096-0.0583)
    mblg_ogs.append(eq.mblg_ogs)
  eqs=session.query(EQ).filter(EQ.mblg_usgs!=None).filter(EQ.ml_ogs!=None).all()
  mblg_usgs=[]
  mw_usgs=[]
  for eq in eqs:   
    mw_usgs.append(eq.ml_ogs*1.096-0.0583)
    mblg_usgs.append(eq.mblg_usgs) 
  print "Total #mblg for MLM",len(mw_usgs)+len(mw_ogs)
  fig=plt.figure()
  ax = fig.add_subplot(111, aspect='equal')
  plt.plot(mblg_ogs,mw_ogs,'ok',label='$M_{bLg}$ OGS')
  plt.plot(mblg_usgs,mw_usgs,'or',label='$M_{bLg}$ USGS')
  plt.plot(x,mw_johnson,'-b',label="Johnson (1996a)")
  plt.plot(x,mw_ceusssc,':b',label="CEUS-SSC (2012)")
  plt.xlabel('$M_{bLg}$')
  plt.ylabel('$M (M_L)$')
  plt.legend(loc='upper left')
  plt.savefig('figures.d/mblg_MLM_rel.pdf',bbox='tight')
  
def weighted_moment(m,m_sigma):
  N=len(m)
  numsum=0.
  denomsum=0.
  m_var=(m_sigma**2)
  for i in np.arange(0,N):
    numsum+=m[i]/m_var[i]
    denomsum+=1/m_var[i]

  me=numsum/denomsum
  m_uncert=np.sqrt(1./denomsum)

  return me,m_uncert

def m_mblg(mblg):
  return 1.14+0.24*mblg+0.093*(mblg**2)
def m_ms(ms):
  return 5.74-0.772*ms+0.128*(ms**2)
def m_mb(mb):
  return 1.487+.4527*mb+.0513*(mb**2)
def m_io(mmi):
  m={2:2.9,3:3.3,4:3.8,5:4.2,6:4.7,7:5.3,8:5.8,9:6.4}
  m_sigma={2:0.59,3:0.55,4:.53,5:.52,6:.52,7:.52,8:.52,9:.53}
  return m[mmi],m_sigma[mmi]

def estimate_moment():
  # Estimate moment magnitude based on CEUS-SSC conversions
  for eq in session.query(EQ).order_by(EQ.id).all():
#    if eq.origintime.year<1962:
#      # Assume the magnitude is determined from the intensity 
#      # if no intensity use magnitude
#      if eq.max_mmi:
#        if eq.max_mmi<=6:
#          eq.m=0.017+0.666*eq.max_mmi
#        if eq.max_mmi>6:
#          eq.m=4.008+3.411*np.sqrt(2)*erfinv((eq.max_mmi-6.)/6.5)
#        eq.m_uncert=0.5
#        eq.m_src='OGSD'
#      else:
#        if eq.mblg_ogs:
#          eq.m=eq.mblg_ogs-0.316
#          eq.m_uncert=0.24
#          eq.m_src='OGSD'
#    else:
     if eq.mw:
        eq.m=eq.mw
        if not eq.m_uncert:
          eq.m_uncert=0.1
          eq.m_src=eq.mw_src
     else:
        m=[]
        m_sigma=[]
        if eq.m3hz_ogs:
          m.append(m_mblg(eq.m3hz_ogs))
          m_sigma.append(0.332)
        if eq.mb:
          m.append(m_mb(eq.mb))
          m_sigma.append(0.394)
        if eq.ms:
          m.append(m_ms(eq.ms))
          m_sigma.append(0.271)
        if eq.mfa:
          m.append(eq.mfa)
          m_sigma.append(0.3)
        if eq.mblg_ogs:
          m.append(m_mblg(eq.mblg_ogs))
          m_sigma.append(0.332)
        if eq.mblg_usgs:
          m.append(m_mblg(eq.mblg_usgs))
          m_sigma.append(0.332)          
        if eq.ml_ogs:
          m.append(-.0583+1.096*eq.ml_ogs)
          m_sigma.append(0.21)
        if eq.md_ogs:
          m.append(0.869+0.762*eq.md_ogs)
          m_sigma.append(0.25)
        if eq.max_mmi!=None:
          if eq.max_mmi>1:
            me,me_uncert=m_io(eq.max_mmi)
            m.append(me)
            m_sigma.append(me_uncert)
      
        if len(m)>0:
          m=np.array(m)
          m_sigma=np.array(m_sigma)
          print m,m_sigma
          eq.m,eq.m_uncert=weighted_moment(m,m_sigma)
          print eq.m,eq.m_uncert
          eq.m_src="OGSD"
          
        session.add(eq)
        session.commit()
def decluster():
  """Decluster the catalog using gardner_knopoff_decluster from OpenQuake's declustering"""
  eqs=session.query(EQ).order_by(EQ.id).all()
  cat=np.zeros(shape=(len(eqs),6))
  for i in np.arange(0,len(eqs)):
    cat[i][0]=eqs[i].origintime.year
    cat[i][1]=eqs[i].origintime.month
    cat[i][2]=eqs[i].origintime.day
    cat[i][3]=eqs[i].longitude
    cat[i][4]=eqs[i].latitude
    cat[i][5]=eqs[i].m
    
  vcl,vmain,flag=gardner_knopoff_decluster(cat,fs_time_prop=0.2)
  print len(vcl),len(vmain),len(flag)
  
  # Update database for earthquake dependence
  for i in np.arange(0,len(eqs)):  
    eqs[i].flag=flag[i]
    
  session.add_all(eqs)
  session.commit()
  
  return flag
  
def completeness(ttol=0.2,dm=.4):
  """Generate stepp plots for the different catalog cuts.  Couldn't get the completeness determinations
  in OpenQuake to give reasonable results, but it is easy enough to identify them from stepp plots."""
  #1882-2008.filter(EQ.flag==0)
  eqs=session.query(EQ).filter(EQ.m>=2.2).filter(EQ.flag==0).filter(EQ.origintime<datetime(2009,1,1,0,0,0)).all()
  year=np.zeros(len(eqs))
  m=np.zeros(len(eqs))
  for i in np.arange(0,len(eqs)):  
    year[i]=eqs[i].origintime.year
    m[i]=eqs[i].m
  
  stepp_plot(year,m,dt=1,xlabel='Years before Jan. 1, 2009',figname='figures.d/stepp_2008declustered.pdf')

  #1882-2008 Complete Catalog
  eqs=session.query(EQ).filter(EQ.m>=2.2).filter(EQ.origintime<datetime(2009,1,1,0,0,0)).all()
  year=np.zeros(len(eqs))
  m=np.zeros(len(eqs))
  for i in np.arange(0,len(eqs)):  
    year[i]=eqs[i].origintime.year
    m[i]=eqs[i].m
  
  stepp_plot(year,m,dt=1,xlabel='Years before Jan. 1, 2009',figname='figures.d/stepp_2008complete.pdf')

  #1882-present declustered
  eqs=session.query(EQ).filter(EQ.m>=2.2).filter(EQ.flag==0).filter(EQ.origintime<datetime(2012,1,1,0,0,0)).all()
  year=np.zeros(len(eqs))
  m=np.zeros(len(eqs))
  for i in np.arange(0,len(eqs)):  
    year[i]=eqs[i].origintime.year
    m[i]=eqs[i].m
  
  stepp_plot(year,m,dt=1,xlabel='Years before Jan. 1, 2012',figname='figures.d/stepp_2011declustered.pdf')
  
  #1882-present declustered
  eqs=session.query(EQ).filter(EQ.m>=2.2).filter(EQ.origintime<datetime(2012,1,1,0,0,0)).all()
  year=np.zeros(len(eqs))
  m=np.zeros(len(eqs))
  for i in np.arange(0,len(eqs)):  
    year[i]=eqs[i].origintime.year
    m[i]=eqs[i].m
  
  stepp_plot(year,m,dt=1,xlabel='Years before Jan. 1, 2012',figname='figures.d/stepp_2011complete.pdf')
  
  return True

def stepp_plot(year,mw,dt=5,xlabel='',figname='temp.pdf'):
  tbins=np.arange(np.max(year)+1,np.min(year),-dt)
  mbins=[2.2,2.9,3.6,4.3,5.0,5.8]
#  stepp=np.zeros(shape=(len(mbins)-1,len(tbins)-1))
  n=np.zeros(shape=(len(mbins)-1,len(tbins)-1))  
  q=np.zeros(len(tbins)-1)
  for i in np.arange(0,len(year)):
    for j in np.arange(1,len(tbins)):
      for k in np.arange(1,len(mbins)):
        if year[i]<tbins[j-1] and year[i]>=tbins[j]:
          if mw[i]>=mbins[k-1] and mw[i]<mbins[k]:
            #print len(tbins)-1,j
            n[k-1][j-1]+=1
    
  t=np.abs(tbins[1:]-(np.max(tbins)))
  #print n
  #rec=np.cumsum(n)/t
  #print t       
  #print rec  
  plt.figure()
  colors=['red','blue','green','cyan','magenta']
  for k in np.arange(0,len(mbins)-1):
    #qt=[]
    st=[]
    for j in np.arange(0,len(t)):
      st.append(np.sum(n[k][j:])/(t[-1]))
      #qt.append(np.sum(q[j:])/(t[-1]))
    print st
    plt.loglog(t,st,'-',color=colors[k],linewidth=1.5)
  plt.grid(True,which="both",ls="-")
  plt.legend(['2.2$\leq$M<2.9','2.9$\leq$M<3.6','3.6$\leq$M<4.3','4.3$\leq$M<5.0','5.0$\leq$M<5.8'])
  plt.ylabel('Annual Rate')
  plt.xlabel(xlabel)
  plt.savefig(figname,bbox='tight')
  plt.show()
  
def recurrence(dmag=.5,norm_mag=2.9,interval=5.0):
  """Calculate recurrence statistics using completeness values from stepp plots"""
  completeness=np.array([[1980.,2.2],
                [1980.,2.9],
                [1960.,3.6],
                [1960.,4.3],
                [1897.,5.0],
                [1897.,5.7]])
  #.filter(EQ.flag==0)
  eqs=session.query(EQ).filter(EQ.m>=2.9).filter(EQ.flag==0).filter(EQ.origintime<datetime(2009,1,1,0,0,0)).all()
  year=np.zeros(len(eqs))
  m=np.zeros(len(eqs))
  for i in np.arange(0,len(eqs)):  
    year[i]=eqs[i].origintime.year
    m[i]=eqs[i].m
  
  # 2.64 at 1 yr interval and dmag=.5

  b,sigma_b,a,sigma_a=recurrence_analysis(year,m,completeness,dmag,'Weichert',norm_mag,interval)
  #print "Nnorm=",10**(a-b*norm_mag),"N0.=",b*norm_mag+10**(a-b*norm_mag)
  a=b*norm_mag+10**(a-b*norm_mag)
  print "2008 declustered b=%f +/-%f a=%f +/- %f, weichert norm mag=%f" % (b,sigma_b,a,sigma_a,norm_mag)
  
  rec_table=recurrence_table(m,dmag,year)
  n_pred=[]
  neg_sigma=[]
  pos_sigma=[]
  for m_est in rec_table[:,0]:
    n_pred.append(10**(a-b*m_est))
    pos_sigma.append(10**(a+sigma_a-(b-sigma_b)*m_est))
    neg_sigma.append(10**(a-sigma_a-(b+sigma_b)*m_est))
  n_pred=np.array(n_pred)
  neg_sigma=np.array(neg_sigma)
  pos_sigma=np.array(pos_sigma)
  plt.figure()
  plt.semilogy(rec_table[:,0],rec_table[:,4],'ok')
  plt.semilogy(rec_table[:,0],n_pred,'-r')
  plt.semilogy(rec_table[:,0],neg_sigma,':r')
  plt.semilogy(rec_table[:,0],pos_sigma,':r')
  plt.xlabel("Magnitude (M)")
  plt.ylabel("# Earthquakes/year M>")
  plt.grid(True,which="both",ls=":")
  plt.savefig('figures.d/b-value_2008declustered.pdf',bbox='tight')
  plt.show()
  
  eqs=session.query(EQ).filter(EQ.m>=2.9).filter(EQ.flag==0).filter(EQ.origintime<datetime(2012,1,1,0,0,0)).all()
  year=np.zeros(len(eqs))
  m=np.zeros(len(eqs))
  for i in np.arange(0,len(eqs)):  
    year[i]=eqs[i].origintime.year
    m[i]=eqs[i].m
  # 2.64 at 1 yr interval and dmag=.5

  b,sigma_b,a,sigma_a=recurrence_analysis(year,m,completeness,dmag,'Weichert',norm_mag,interval)
  a=b*norm_mag+10**(a-b*norm_mag)
  print "2012 declustered b=%f +/-%f a=%f +/- %f, weichert norm mag=%f" % (b,sigma_b,a,sigma_a,norm_mag)
  
  rec_table=recurrence_table(m,dmag,year)
  n_pred=[]
  neg_sigma=[]
  pos_sigma=[]
  for m_est in rec_table[:,0]:
    n_pred.append(10**(a-b*m_est))
    pos_sigma.append(10**(a+sigma_a-(b-sigma_b)*m_est))
    neg_sigma.append(10**(a-sigma_a-(b+sigma_b)*m_est))
  n_pred=np.array(n_pred)
  neg_sigma=np.array(neg_sigma)
  pos_sigma=np.array(pos_sigma)
  plt.figure()
  plt.semilogy(rec_table[:,0],rec_table[:,4],'ok')
  plt.semilogy(rec_table[:,0],n_pred,'-r')
  plt.semilogy(rec_table[:,0],neg_sigma,':r')
  plt.semilogy(rec_table[:,0],pos_sigma,':r')
  plt.xlabel("Magnitude (M)")
  plt.ylabel("# Earthquakes/year M>")
  plt.grid(True,which="both",ls=":")
  plt.savefig('figures.d/b-value_2012declustered.pdf',bbox='tight')
  plt.show()
  
  
  eqs=session.query(EQ).filter(EQ.m>=2.2).filter(EQ.flag==0).\
  filter(EQ.origintime<datetime(2012,1,1,0,0,0)).filter(EQ.origintime>datetime(2008,12,31,23,59,59)).all()
  year=np.zeros(len(eqs))
  m=np.zeros(len(eqs))
  for i in np.arange(0,len(eqs)):  
    year[i]=eqs[i].origintime.year
    m[i]=eqs[i].m
  # 2.64 at 1 yr interval and dmag=.5
  norm_mag=2.74
  rec_table=recurrence_table(m,dmag,year)
  b,sigma_b=b_max_likelihood(rec_table[:,0],rec_table[:,1],dmag=dmag,m_c=2.2)
  a,sigma_a=a_value(b,np.array(rec_table[:,0]),np.array(rec_table[:,4]))
  print "2009-2011 declustered b=%f +/-%f a=%f +/- %f, mle norm mag=%f" % (b,sigma_b,a,sigma_a,norm_mag)
  
  
  n_pred=[]
  neg_sigma=[]
  pos_sigma=[]
  for m_est in rec_table[:,0]:
    n_pred.append(10**(a-b*m_est))
    pos_sigma.append(10**(a+sigma_a-(b-sigma_b)*m_est))
    neg_sigma.append(10**(a-sigma_a-(b+sigma_b)*m_est))
  n_pred=np.array(n_pred)
  neg_sigma=np.array(neg_sigma)
  pos_sigma=np.array(pos_sigma)
  plt.figure()
  plt.semilogy(rec_table[:,0],rec_table[:,4],'ok')
  plt.semilogy(rec_table[:,0],n_pred,'-r')
  plt.semilogy(rec_table[:,0],neg_sigma,':r')
  plt.semilogy(rec_table[:,0],pos_sigma,':r')
  plt.xlabel("Magnitude (M)")
  plt.ylabel("# Earthquakes/year M>")
  plt.grid(True,which="both",ls=":")
  plt.savefig('figures.d/b-value_2009-2012declustered.pdf',bbox='tight')
  plt.show()
  
  eqs=session.query(EQ).filter(EQ.m>=2.2).\
  filter(EQ.origintime<datetime(2012,1,1,0,0,0)).filter(EQ.origintime>datetime(2008,12,31,23,59,59)).all()
  year=np.zeros(len(eqs))
  m=np.zeros(len(eqs))
  for i in np.arange(0,len(eqs)):  
    year[i]=eqs[i].origintime.year
    m[i]=eqs[i].m
  # 2.64 at 1 yr interval and dmag=.5
  rec_table=recurrence_table(m,dmag,year)
  b,sigma_b=b_max_likelihood(rec_table[:,0],rec_table[:,1],dmag=dmag,m_c=2.2)
  a,sigma_a=a_value(b,np.array(rec_table[:,0]),np.array(rec_table[:,4]))
  print "2009-2011 clustered b=%f +/-%f a=%f +/- %f, mle norm mag=%f" % (b,sigma_b,a,sigma_a,norm_mag)
  
  
  n_pred=[]
  neg_sigma=[]
  pos_sigma=[]
  for m_est in rec_table[:,0]:
    n_pred.append(10**(a-b*m_est))
    pos_sigma.append(10**(a+sigma_a-(b-sigma_b)*m_est))
    neg_sigma.append(10**(a-sigma_a-(b+sigma_b)*m_est))
  n_pred=np.array(n_pred)
  neg_sigma=np.array(neg_sigma)
  pos_sigma=np.array(pos_sigma)
  plt.figure()
  plt.semilogy(rec_table[:,0],rec_table[:,4],'ok')
  plt.semilogy(rec_table[:,0],n_pred,'-r')
  plt.semilogy(rec_table[:,0],neg_sigma,':r')
  plt.semilogy(rec_table[:,0],pos_sigma,':r')
  plt.xlabel("Magnitude (M)")
  plt.ylabel("# Earthquakes/year M>")
  plt.grid(True,which="both",ls=":")
  plt.savefig('figures.d/b-value_2009-2012clustered.pdf',bbox='tight')
  plt.show()

def recurrence2010(dmag=.5):
  eqs=session.query(EQ).filter(EQ.m>=2.2).filter(EQ.flag==0).\
  filter(EQ.origintime>datetime(2009,12,31,23,59,59)).all()
  year=np.zeros(len(eqs))
  m=np.zeros(len(eqs))
  for i in np.arange(0,len(eqs)):  
    year[i]=eqs[i].origintime.year
    m[i]=eqs[i].m
  # 2.64 at 1 yr interval and dmag=.5
  norm_mag=2.74
  rec_table=recurrence_table(m,dmag,year)
  b,sigma_b=b_max_likelihood(rec_table[:,0],rec_table[:,1],dmag=dmag,m_c=2.2)
  a,sigma_a=a_value(b,np.array(rec_table[:,0]),np.array(rec_table[:,4]))
  print "2010 declustered b=%f +/-%f a=%f +/- %f, mle norm mag=%f" % (b,sigma_b,a,sigma_a,norm_mag)
  
  
  n_pred=[]
  neg_sigma=[]
  pos_sigma=[]
  for m_est in rec_table[:,0]:
    n_pred.append(10**(a-b*m_est))
    pos_sigma.append(10**(a+sigma_a-(b-sigma_b)*m_est))
    neg_sigma.append(10**(a-sigma_a-(b+sigma_b)*m_est))
  n_pred=np.array(n_pred)
  neg_sigma=np.array(neg_sigma)
  pos_sigma=np.array(pos_sigma)
  plt.figure()
  plt.semilogy(rec_table[:,0],rec_table[:,4],'ok')
  plt.semilogy(rec_table[:,0],n_pred,'-r')
  plt.semilogy(rec_table[:,0],neg_sigma,':r')
  plt.semilogy(rec_table[:,0],pos_sigma,':r')
  plt.xlabel("Magnitude (M)")
  plt.ylabel("# Earthquakes/year M>")
  xlim=plt.xlim()
  ylim=plt.ylim()
  plt.text(xlim[0]+.1*(xlim[1]-xlim[0]),ylim[0]*11,"$log_{10}(N)=%.4f+%.4fM$" % (a,b))
  plt.grid(True,which="both",ls=":")
  plt.savefig('figures.d/b-value_2010declustered.pdf',bbox='tight')
  plt.show()
  
  eqs=session.query(EQ).filter(EQ.m>=2.2).\
  filter(EQ.origintime>datetime(2008,12,31,23,59,59)).all()
  year=np.zeros(len(eqs))
  m=np.zeros(len(eqs))
  for i in np.arange(0,len(eqs)):  
    year[i]=eqs[i].origintime.year
    m[i]=eqs[i].m
  # 2.64 at 1 yr interval and dmag=.5
  rec_table=recurrence_table(m,dmag,year)
  b,sigma_b=b_max_likelihood(rec_table[:,0],rec_table[:,1],dmag=dmag,m_c=2.2)
  a,sigma_a=a_value(b,np.array(rec_table[:,0]),np.array(rec_table[:,4]))
  print "2010 clustered b=%f +/-%f a=%f +/- %f, mle norm mag=%f" % (b,sigma_b,a,sigma_a,norm_mag)
  
  
  n_pred=[]
  neg_sigma=[]
  pos_sigma=[]
  for m_est in rec_table[:,0]:
    n_pred.append(10**(a-b*m_est))
    pos_sigma.append(10**(a+sigma_a-(b-sigma_b)*m_est))
    neg_sigma.append(10**(a-sigma_a-(b+sigma_b)*m_est))
  n_pred=np.array(n_pred)
  neg_sigma=np.array(neg_sigma)
  pos_sigma=np.array(pos_sigma)
  plt.figure()
  plt.semilogy(rec_table[:,0],rec_table[:,4],'ok')
  plt.semilogy(rec_table[:,0],n_pred,'-r')
  plt.semilogy(rec_table[:,0],neg_sigma,':r')
  plt.semilogy(rec_table[:,0],pos_sigma,':r')
  plt.xlabel("Magnitude (M)")
  plt.ylabel("# Earthquakes/year M>")
  xlim=plt.xlim()
  ylim=plt.ylim()
  plt.text(xlim[0]+.1*(xlim[1]-xlim[0]),ylim[0]*11,"$log_{10}(N)=%.4f+%.4fM$" % (a,b))
  plt.grid(True,which="both",ls=":")
  plt.savefig('figures.d/b-value_2010clustered.pdf',bbox='tight')
  plt.show()

def a_value(b,mval,number_obs):
  nlog=[]
  m=[]
  for i in np.arange(0,len(number_obs)):
    if number_obs[i]>0.:
      nlog.append(np.log10(number_obs[i]))
      m.append(mval[i])
  nlog=np.array(nlog)
  m=np.array(m)
  a=nlog+b*m
  #print a
  return np.mean(a),np.std(a)
  
def main():
  pass


if __name__ == '__main__':
  main()

