allfracs_final.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os \n",
    "import matplotlib\n",
    "%matplotlib widget\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import matplotlib.pyplot as plt \n",
    "import matplotlib.cm as cm\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import random\n",
    "import cutde.fullspace as FS\n",
    "import plotly\n",
    "import plotly.graph_objects as go\n",
    "from scipy.spatial.transform import Rotation as R\n",
    "# from scipy import sparse\n",
    "from scipy.linalg import svd, lstsq\n",
    "from scipy.integrate import cumulative_trapezoid\n",
    "from scipy.signal import resample, detrend,medfilt\n",
    "from scipy.stats import bootstrap\n",
    "from operator import itemgetter\n",
    "import matplotlib.dates as mdates\n",
    "from datetime import datetime\n",
    "from more_itertools import sliced\n",
    "from numpy.linalg import norm\n",
    "from numpy.linalg import solve\n",
    "# import scipy.interpolate as sp\n",
    "import cmasher as cmr\n",
    "import pickle\n",
    "from dts_funcs import * "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in DAS data\n",
    "def stack_DAS(das_dat):\n",
    "    n = das_dat.shape[1]\n",
    "    m = n // 2\n",
    "    a1 = das_dat[:,1:m+1]\n",
    "    a2 = np.fliplr(das_dat[:,m+1:])\n",
    "    return  np.divide(np.add(a1,a2),2)\n",
    "\n",
    "df_fullPDT = pd.read_pickle('/home/spri902/EGS_Collab/4850/results/maystim/processed_DAS/lpFilter/wellPDT/maystim22_26_combined_full')\n",
    "df_fullPDB = pd.read_pickle('/home/spri902/EGS_Collab/4850/results/maystim/processed_DAS/lpFilter/wellPDB/maystim22_26_combined_full')\n",
    "df_fullPST = pd.read_pickle('/home/spri902/EGS_Collab/4850/results/maystim/processed_DAS/lpFilter/wellPST/maystim22_26_combined_full')\n",
    "df_fullPSB = pd.read_pickle('/home/spri902/EGS_Collab/4850/results/maystim/processed_DAS/lpFilter/wellPSB/maystim22_26_combined_full')\n",
    "df_fullOT = pd.read_pickle('/home/spri902/EGS_Collab/4850/results/maystim/processed_DAS/lpFilter/wellOT/maystim22_26_combined_full')\n",
    "df_fullOB = pd.read_pickle('/home/spri902/EGS_Collab/4850/results/maystim/processed_DAS/lpFilter/wellOB/maystim22_26_combined_full')\n",
    "\n",
    "df_fullPDT = detrend(df_fullPDT,axis=0,type='linear')\n",
    "df_fullPDB = detrend(df_fullPDB,axis=0,type='linear')\n",
    "df_fullPST = detrend(df_fullPST,axis=0,type='linear')\n",
    "df_fullPSB = detrend(df_fullPSB,axis=0,type='linear')\n",
    "df_fullOT = detrend(df_fullOT,axis=0,type='linear')\n",
    "df_fullOB = detrend(df_fullOB,axis=0,type='linear')\n",
    "\n",
    "# df_fullPDT = stack_DAS(df_fullPDT)\n",
    "# df_fullPDB = stack_DAS(df_fullPDB)\n",
    "# df_fullPST = stack_DAS(df_fullPST)\n",
    "# df_fullPSB = stack_DAS(df_fullPSB)\n",
    "# df_fullOT = stack_DAS(df_fullOT)\n",
    "# df_fullOB = stack_DAS(df_fullOB)\n",
    "\n",
    "\n",
    "wn     = ['PDT','PDB','PST','PSB','OT','OB']\n",
    "nfile_list = sorted(os.walk('/data1/parker/EGS_iDAS'))\n",
    "nfile_list = nfile_list[1:]\n",
    "#file_list = file_list[1:]\n",
    "nfile_list = [group[2] for group in nfile_list]\n",
    "nfile_list = [item for sublist in nfile_list for item in sublist]\n",
    "# [file_list.append(f) for f in nfile_list]\n",
    "fd = [name.split(\"_\") for name in nfile_list]\n",
    "fl = [fd[file][2].split(\".\") for file in range(len(fd))]\n",
    "fl = [el[0] for el in fl]\n",
    "DASdates = [datetime.strptime(d,'%y%m%d%H%M%S') for d in sorted(fl)]\n",
    "# these are files that get skipped during the low pass filtering process and so the dates need to be removed \n",
    "ind2rem = [0, 90, 91, 257, 258, 1571, 1572, 3082, 3083, 5085, 5086, 5599, 5600, 5961, 5962, 7623, 7624, 8841, 8842, 9562]\n",
    "# remove in reverse so that the indices remain in the correct order for removal\n",
    "for index in sorted(ind2rem,reverse=True):\n",
    "    del DASdates[index]\n",
    "\n",
    "directory = sorted(os.walk('/home/spri902/mayCASSM'))\n",
    "CASSMdates = [datetime.strptime(d,'%Y%m%d%H%M%S') for d in sorted(directory[0][1])]\n",
    "chansT=np.linspace(0,df_fullPDT.shape[1] - 1,df_fullPDT.shape[1]).astype(int)\n",
    "chansB=np.linspace(0,df_fullPDB.shape[1] - 1,df_fullPDB.shape[1]).astype(int)\n",
    "# stimbeg = [96, 224, 352, 472, 507]\n",
    "# stimfin = [106, 232, 368, 475, 512]\n",
    "# stimfin = [106, 232, 368, 475, 511]\n",
    "# stimbegLines = itemgetter(*stimbeg)(CASSMdates)\n",
    "# stimfinLines = itemgetter(*stimfin)(CASSMdates)\n",
    "stimbeg = [0, 2448, 5737, 7749, 8405]\n",
    "# stimfin = [106, 232, 368, 475, 512]\n",
    "stimfin = [29, 2578, 5813, 7807, 8469]\n",
    "stimbegLines = itemgetter(*stimbeg)(DASdates)\n",
    "stimfinLines = itemgetter(*stimfin)(DASdates)\n",
    "\n",
    "dasdnums = mdates.date2num(DASdates)\n",
    "cassmdnums = mdates.date2num(CASSMdates)\n",
    "dasd = pd.date_range(start=DASdates[0],end=DASdates[-1],periods=len(cassmdnums[:-2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# scaling equation provided by Silixa (tdmsvalues/2048)*116*(fs/gaugeLength)\n",
    "\n",
    "# set up matrices to scale the data by their different Fs values\n",
    "fsvecT = np.ones_like(df_fullPDT)\n",
    "fsvecB = np.ones_like(df_fullPDB)\n",
    "fsvecOT = np.ones_like(df_fullOT)\n",
    "fsvecOB = np.ones_like(df_fullOB)\n",
    "fsvecPST = np.ones_like(df_fullPST)\n",
    "fsvecPSB = np.ones_like(df_fullPSB)\n",
    "# These are the indices where the sampling rate changes between 1kHz and 10kHz\n",
    "lst = [0,88,10000,89,1565,1000,1566,5075,10000,5076,5587,1000,5588,7607,10000,7608,9543,1000]\n",
    "lst = np.array(list(sliced((lst),3))).reshape(6,3)\n",
    "# fill the matrices with the correct sampling freqs\n",
    "for i in range(len(lst)):\n",
    "    fsvecT[lst[i][0]:lst[i][1]+1,:] = fsvecT[lst[i][0]:lst[i][1]+1,:]*lst[i][2]\n",
    "    fsvecB[lst[i][0]:lst[i][1]+1,:] = fsvecB[lst[i][0]:lst[i][1]+1,:]*lst[i][2]\n",
    "    fsvecOT[lst[i][0]:lst[i][1]+1,:] = fsvecOT[lst[i][0]:lst[i][1]+1,:]*lst[i][2]\n",
    "    fsvecOB[lst[i][0]:lst[i][1]+1,:] = fsvecOB[lst[i][0]:lst[i][1]+1,:]*lst[i][2]\n",
    "    fsvecPST[lst[i][0]:lst[i][1]+1,:] = fsvecPST[lst[i][0]:lst[i][1]+1,:]*lst[i][2]\n",
    "    fsvecPSB[lst[i][0]:lst[i][1]+1,:] = fsvecPSB[lst[i][0]:lst[i][1]+1,:]*lst[i][2]\n",
    "# divide by gauge length\n",
    "fsvecT /= 10.\n",
    "fsvecB /= 10.\n",
    "fsvecOT /= 10.\n",
    "fsvecOB /= 10.\n",
    "fsvecPST /= 10.\n",
    "fsvecPSB /= 10.\n",
    "# multiply the data by the scaling matrices and then scale by the Silixa scaling value\n",
    "scaledT = np.multiply(fsvecT,df_fullPDT)*(116/2048)\n",
    "scaledB = np.multiply(fsvecB,df_fullPDB)*(116/2048)\n",
    "scaledOT = np.multiply(fsvecOT,df_fullOT)*(116/2048)\n",
    "scaledOB = np.multiply(fsvecOB,df_fullOB)*(116/2048)\n",
    "scaledPST = np.multiply(fsvecPST,df_fullPST)*(116/2048)\n",
    "scaledPSB = np.multiply(fsvecPSB,df_fullPSB)*(116/2048)\n",
    "\n",
    "# Integrate strain rates to strain\n",
    "df_fullPDT = scaledT.copy()\n",
    "df_fullPDB = scaledB.copy()\n",
    "df_fullPST = scaledPST.copy()\n",
    "df_fullPSB = scaledPSB.copy()\n",
    "df_fullOT = scaledOT.copy()\n",
    "df_fullOB = scaledOB.copy()\n",
    "# Integrate strain rates to strain\n",
    "df_strainPDT = cumulative_trapezoid(df_fullPDT,axis=0,initial=0)\n",
    "df_strainPDB = cumulative_trapezoid(df_fullPDB,axis=0,initial=0)\n",
    "df_strainPST = cumulative_trapezoid(df_fullPST,axis=0,initial=0)\n",
    "df_strainPSB = cumulative_trapezoid(df_fullPSB,axis=0,initial=0)\n",
    "df_strainOT = cumulative_trapezoid(df_fullOT,axis=0,initial=0)\n",
    "df_strainOB = cumulative_trapezoid(df_fullOB,axis=0,initial=0)\n",
    "# create channel spacing vectors\n",
    "chansOB=np.linspace(0,df_strainOB.shape[1] - 1,df_strainOB.shape[1]).astype(int)\n",
    "chansOT=np.linspace(0,df_strainOT.shape[1] - 1,df_strainOT.shape[1]).astype(int)\n",
    "chansPSB=np.linspace(0,df_strainPSB.shape[1] - 1,df_strainPSB.shape[1]).astype(int)\n",
    "chansPST=np.linspace(0,df_strainPST.shape[1] - 1,df_strainPST.shape[1]).astype(int)\n",
    "chansPDB=np.linspace(0,df_strainPDB.shape[1] - 1,df_strainPDB.shape[1]).astype(int)\n",
    "chansPDT=np.linspace(0,df_strainPDT.shape[1] - 1,df_strainPDT.shape[1]).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in well geometry\n",
    "wellcols=np.arange(12)\n",
    "WellGeom=pd.read_excel('/home/spri902/Collab_metadata/Well_Points.xlsx',header=0,usecols=wellcols)\n",
    "WellGeom.columns = WellGeom.columns.str.replace(' ','')\n",
    "# convert from feet to meters\n",
    "WellGeom.x = WellGeom.x/3.28084\n",
    "WellGeom.y = WellGeom.y/3.28084\n",
    "WellGeom.z = WellGeom.z/3.28084\n",
    "# pull out the monitoring wells\n",
    "mwells=[]\n",
    "wellList = ['E1-OB','E1-OT','E1-PDT','E1-PDB','E1-PST','E1-PSB']\n",
    "for i in wellList:\n",
    "    tmpwell = WellGeom[WellGeom[\"HoleID\"]== i]\n",
    "    tmpwell = tmpwell.iloc[:,0:6]\n",
    "    mwells.append(tmpwell)\n",
    "mwells = pd.concat(mwells)\n",
    "# pull out the injector and production wells\n",
    "swells=[]\n",
    "wellList = ['E1-I','E1-P']\n",
    "for i in wellList:\n",
    "    tmpwell = WellGeom[WellGeom[\"HoleID\"]== i]\n",
    "    tmpwell = tmpwell.iloc[:,0:6]\n",
    "    swells.append(tmpwell)\n",
    "swells = pd.concat(swells)\n",
    "\n",
    "# take every 10th point (1 meter spacing along well)\n",
    "PDBpnts = mwells[mwells['HoleID'] == 'E1-PDB' ].iloc[0::10,:]\n",
    "PDBpnts = PDBpnts[1:df_fullPDB.shape[1]]\n",
    "PDTpnts = mwells[mwells['HoleID'] == 'E1-PDT' ].iloc[0::10,:]\n",
    "PDTpnts = PDTpnts[1:df_fullPDT.shape[1]]\n",
    "\n",
    "PSBpnts = mwells[mwells['HoleID'] == 'E1-PSB' ].iloc[0::10,:]\n",
    "PSBpnts = PSBpnts[1:df_fullPSB.shape[1]]\n",
    "PSTpnts = mwells[mwells['HoleID'] == 'E1-PST' ].iloc[0::10,:]\n",
    "PSTpnts = PSTpnts[1:df_fullPST.shape[1]]\n",
    "\n",
    "OBpnts = mwells[mwells['HoleID'] == 'E1-OB' ].iloc[0::10,:]\n",
    "OBpnts = OBpnts[1:df_fullOB.shape[1]]\n",
    "OTpnts = mwells[mwells['HoleID'] == 'E1-OT' ].iloc[0::10,:]\n",
    "OTpnts = OTpnts[1:df_fullOT.shape[1]]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# well coordinates\n",
    "\n",
    "pdtobsx = np.array(PDTpnts.x) # well PDT\n",
    "pdbobsx = np.array(PDBpnts.x) # well PDB\n",
    "pstobsx = np.array(PSTpnts.x) # well PST\n",
    "psbobsx = np.array(PSBpnts.x) # well PSB\n",
    "otobsx = np.array(OTpnts.x) # well OT\n",
    "obobsx = np.array(OBpnts.x) # well OB\n",
    "\n",
    "pdtobsy = np.array(PDTpnts.y)\n",
    "pdbobsy = np.array(PDBpnts.y)\n",
    "pstobsy = np.array(PSTpnts.y)\n",
    "psbobsy = np.array(PSBpnts.y)\n",
    "otobsy = np.array(OTpnts.y)\n",
    "obobsy = np.array(OBpnts.y)\n",
    "\n",
    "pdtobsz = np.array(PDTpnts.z)\n",
    "pdbobsz = np.array(PDBpnts.z)\n",
    "pstobsz = np.array(PSTpnts.z)\n",
    "psbobsz = np.array(PSBpnts.z)\n",
    "otobsz = np.array(OTpnts.z)\n",
    "obobsz = np.array(OBpnts.z)\n",
    "\n",
    "pdtpts = np.array([pdtobsx, pdtobsy,  pdtobsz]).reshape((3, -1)).T.copy()\n",
    "pdbpts = np.array([pdbobsx, pdbobsy,  pdbobsz]).reshape((3, -1)).T.copy()\n",
    "pstpts = np.array([pstobsx, pstobsy,  pstobsz]).reshape((3, -1)).T.copy()\n",
    "psbpts = np.array([psbobsx, psbobsy,  psbobsz]).reshape((3, -1)).T.copy()\n",
    "otpts = np.array([otobsx, otobsy,  otobsz]).reshape((3, -1)).T.copy()\n",
    "obpts = np.array([obobsx, obobsy,  obobsz]).reshape((3, -1)).T.copy()\n",
    "\n",
    "allpts = np.concatenate((pdtpts,pdbpts,pstpts,psbpts,otpts,obpts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bring in MEQ catalog\n",
    "meq = pd.read_csv('/home/spri902/EGS_Collab/4850/MEQ/MEQcat.csv',header=0,parse_dates=[11],infer_datetime_format=True)\n",
    "meq.columns = meq.columns.str.replace(' ','')\n",
    "dates = mdates.date2num(meq['time'])\n",
    "meq['datenums'] = dates\n",
    "meq = meq[(meq['time'] > pd.to_datetime('2018-05-24T21:00')) & (meq['time'] < pd.to_datetime('2018-05-25T23'))]\n",
    "# meq = meq[(meq['time'] > pd.to_datetime('2018-05-25T20:30')) & (meq['time'] < pd.to_datetime('2018-05-25T21:35'))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bring in Injection data\n",
    "os.chdir('/home/spri902/EGS_Collab/4850/stimflow/')\n",
    "injFiles = sorted(os.listdir('/home/spri902/EGS_Collab/4850/stimflow/'))\n",
    "injDat = pd.concat((pd.read_csv(f,header=1,usecols=[0,2,4,24],\\\n",
    "                           parse_dates = [0],infer_datetime_format=True) \\\n",
    "                    for f in injFiles if f.endswith('.csv')),axis=0)\n",
    "injDat.rename(columns={'hh:mm:ss':'date','LPM':'QLPM','LPM.1':'TLPM','psig.9':'psig'},inplace=True)\n",
    "injDat.reset_index(drop = True,inplace = True)\n",
    "injDat.set_index('date',inplace=True)\n",
    "injDat = injDat.iloc[14650:]\n",
    "injDat = injDat[~injDat.index.duplicated(keep='first')]\n",
    "\n",
    "cumQvol = np.cumsum(injDat.QLPM)\n",
    "cumTvol = np.cumsum(injDat.TLPM)\n",
    "flwRate = injDat.QLPM + injDat.TLPM"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Read in DTS data and calculate temperature imposed strain\n",
    "def calc_dts_strain(dts_list,well_list):\n",
    "\n",
    "\n",
    "    dts_dat = {}\n",
    "    for well in well_list:\n",
    "        time = []\n",
    "        deps = dts_list[0][well]['depth']\n",
    "        dat =  []\n",
    "        ref =  dts_list[0][well]['data'][:,0]\n",
    "        for i in range(len(dts_list)):\n",
    "            time.append(dts_list[i][well]['times'])\n",
    "            # dat.append(dts_list[i][well]['data'])\n",
    "            dat = np.concatenate([dts_list[x][well]['data']-ref[:,None] for x in range(len(dts_list))],axis=1)\n",
    "            # dat = np.concatenate([dts_list[x][well]['data'] for x in range(len(dts_list))],axis=1)\n",
    "        \n",
    "        time = np.concatenate([t.flatten() for t in time])\n",
    "    \n",
    "        dts_dat[well] = {'data': dat, 'depth': deps,\n",
    "                            'times': time}\n",
    "    return dts_dat\n",
    "\n",
    "# Create a list of dict DTS data from each pickle file\n",
    "# Each list is len(dts_list) = number of pkl files\n",
    "#  and Plot processed DTS images  \n",
    "#\n",
    "dts_list=[]\n",
    "for path,subdir,files in sorted(os.walk('/home/spri902/EGS_Collab/4850/DTS/processed/data/Exp1')):\n",
    "    for file in sorted(files):\n",
    "        with open(os.path.join(path,f'{file}'),'rb') as f:\n",
    "            well_dict = pickle.load(f)\n",
    "            dts_list.append(well_dict)\n",
    "for i in range(len(wells_4850)):    \n",
    "    plot_DTS(dts_list,well=wells_4850[i])\n",
    "\n",
    "\n",
    "# \n",
    "dts_strain = calc_dts_strain(dts_list,wn)\n",
    "dts_times = dts_strain['OT']['times']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pick a DTS record time as close as possible to a DAS record time\n",
    "dts_t0 = 423\n",
    "# dts_t0 = 523 # index for date time in DTS data\n",
    "# dts_t0 = 557\n",
    "print(dts_times[dts_t0])\n",
    "# sum every 4 channels in DTS and take average to match DAS channel spacing\n",
    "otdts = np.add.reduceat(dts_strain['OT']['data'][:,dts_t0],np.arange(0,len(dts_strain['OT']['data'][:,dts_t0]),4)) / 4\n",
    "obdts = np.add.reduceat(dts_strain['OB']['data'][:,dts_t0],np.arange(0,len(dts_strain['OB']['data'][:,dts_t0]),4)) / 4\n",
    "pstdts = np.add.reduceat(dts_strain['PST']['data'][:,dts_t0],np.arange(0,len(dts_strain['PST']['data'][:,dts_t0]),4)) / 4\n",
    "psbdts = np.add.reduceat(dts_strain['PSB']['data'][:,dts_t0],np.arange(0,len(dts_strain['PSB']['data'][:,dts_t0]),4)) / 4\n",
    "pdtdts = np.add.reduceat(dts_strain['PDT']['data'][:,dts_t0],np.arange(0,len(dts_strain['PDT']['data'][:,dts_t0]),4)) / 4\n",
    "pdbdts = np.add.reduceat(dts_strain['PDB']['data'][:,dts_t0],np.arange(0,len(dts_strain['PDB']['data'][:,dts_t0]),4)) / 4\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_patch_spines_invisible(ax):\n",
    "    ax.set_frame_on(True)\n",
    "    ax.patch.set_visible(False)\n",
    "    for sp in ax.spines.values():\n",
    "        sp.set_visible(False)\n",
    "\n",
    "def four_yaxes(ax,dastimes,chans,dasdat,flwdat,injdat,meqdat,colmap,flwcolor,prescolor,meqcol,wellnm,flag):\n",
    "    ax1 = ax.twinx()\n",
    "    ax2 = [ax,ax.twinx()]\n",
    "    ax3 = ax.twinx()\n",
    "    ax1.spines[\"right\"].set_position((\"axes\",1.1))\n",
    "    ax3.spines[\"right\"].set_position((\"axes\",1.25))\n",
    "    make_patch_spines_invisible(ax1)\n",
    "    ax1.spines[\"right\"].set_visible(True)\n",
    "    make_patch_spines_invisible(ax3)\n",
    "    ax3.spines[\"right\"].set_visible(True)\n",
    "    vm=np.nanpercentile(dasdat[7700:,:],99)\n",
    "    img = ax.pcolormesh(dastimes,chans,dasdat.T,cmap=colmap,vmin=-vm,vmax=vm)\n",
    "    ax.xaxis_date()\n",
    "    ax.xaxis.set_tick_params(which='both',labelbottom=True)\n",
    "\n",
    "    ax2[0].set_ylabel('Channel Number')\n",
    "    im2 = ax2[1].plot(flwdat,'.',color=flwcolor,markersize=2,label='Flow Rate')\n",
    "    ax2[1].set_ylabel('Flow Rate (LPM)')\n",
    "    im3 = ax3.plot(injdat,'.',color=prescolor,markersize=2,label='Inj. Pressure')\n",
    "    ax3.set_ylabel('Inj. Pressure (psig)')\n",
    "    # im2[0].yaxis.label.set_color(ax2[1].get_color())\n",
    "    ax3.yaxis.label.set_color(im3[0].get_color())\n",
    "    ax2[1].yaxis.label.set_color(im2[0].get_color())\n",
    "    im1 = ax1.plot(meqdat['datenums'],meqdat['D_OT_P'],'.',color=meqcol,label='MEQs')\n",
    "    ax1.set_ylabel('MEQs distance from OTP (m)')\n",
    "    ax1.yaxis.label.set_color(im1[0].get_color())\n",
    "    ax1.plot(meqdat['datenums'],np.zeros_like(meqdat['datenums']),'-m')\n",
    "    if flag:\n",
    "        plt.colorbar(img,orientation=\"horizontal\",label='Nanostrain Rate')\n",
    "        plt.title(f'{wellnm} Strain Rate')\n",
    "    else:\n",
    "        plt.colorbar(img,orientation=\"horizontal\",label='Nanostrain')\n",
    "        plt.title(f'{wellnm} Strain')\n",
    "    return img, im2, im3, im1\n",
    "        \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1,ax2,ax3) = plt.subplots(3,2,figsize=(12,8),sharex=True,sharey=True)\n",
    "fig.subplots_adjust(right=0.75)\n",
    "a1,a1a,a1b,a1c = four_yaxes(ax1[0],DASdates,chansOT,df_fullOT,flwRate,injDat.loc[:,(\"psig\")],meq,'seismic','cyan','green','orange',wn[4],1)\n",
    "a2,a2a,a2b,a2c = four_yaxes(ax1[1],DASdates,chansOB,df_fullOB,flwRate,injDat.loc[:,(\"psig\")],meq,'seismic','cyan','green','orange',wn[5],1)\n",
    "a3,a3a,a3b,a3c = four_yaxes(ax2[0],DASdates,chansPDT,df_fullPDT,flwRate,injDat.loc[:,(\"psig\")],meq,'seismic','cyan','green','orange',wn[0],1)\n",
    "a4,a4a,a4b,a4c = four_yaxes(ax2[1],DASdates,chansPDB,df_fullPDB,flwRate,injDat.loc[:,(\"psig\")],meq,'seismic','cyan','green','orange',wn[1],1)\n",
    "a5,a5a,a5b,a5c = four_yaxes(ax3[0],DASdates,chansPST,df_fullPST,flwRate,injDat.loc[:,(\"psig\")],meq,'seismic','cyan','green','orange',wn[2],1)\n",
    "a6,a6a,a6b,a6c = four_yaxes(ax3[1],DASdates,chansPSB,df_fullPSB,flwRate,injDat.loc[:,(\"psig\")],meq,'seismic','cyan','green','orange',wn[3],1)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1,ax2,ax3) = plt.subplots(3,2,figsize=(12,8),sharex=True,sharey=True)\n",
    "fig.subplots_adjust(right=0.75)\n",
    "a1,a1a,a1b,a1c = four_yaxes(ax1[0],DASdates,chansOT,df_strainOT,flwRate,injDat.loc[:,(\"psig\")],meq,'seismic','black','green','orange',wn[4],0)\n",
    "a2,a2a,a2b,a2c = four_yaxes(ax1[1],DASdates,chansOB,df_strainOB,flwRate,injDat.loc[:,(\"psig\")],meq,'seismic','black','green','orange',wn[5],0)\n",
    "a3,a3a,a3b,a3c = four_yaxes(ax2[0],DASdates,chansPDT,df_strainPDT,flwRate,injDat.loc[:,(\"psig\")],meq,'seismic','black','green','orange',wn[0],0)\n",
    "a4,a4a,a4b,a4c = four_yaxes(ax2[1],DASdates,chansPDB,df_strainPDB,flwRate,injDat.loc[:,(\"psig\")],meq,'seismic','black','green','orange',wn[1],0)\n",
    "a5,a5a,a5b,a5c = four_yaxes(ax3[0],DASdates,chansPST,df_strainPST,flwRate,injDat.loc[:,(\"psig\")],meq,'seismic','black','green','orange',wn[2],0)\n",
    "a6,a6a,a6b,a6c = four_yaxes(ax3[1],DASdates,chansPSB,df_strainPSB,flwRate,injDat.loc[:,(\"psig\")],meq,'seismic','black','green','orange',wn[3],0)\n",
    "\n",
    "plt.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# switch dir to where the fracture points are\n",
    "os.chdir('/home/spri902/EGS_Collab/4850/fractures/')\n",
    "# Read in xyz points for OTP-connector and hydrofracs\n",
    "with open('frac_plane_points.npy','rb') as f:\n",
    "    x = np.load(f)\n",
    "\n",
    "# Pull out OTP-connector points\n",
    "otp = x[0:80,:]\n",
    "westN = x[80:160,:]\n",
    "# westN = interparc(80,westN[:,0],westN[:,1],westN[:,2]) \n",
    "westS = x[160:240,:]\n",
    "westD = x[240:320,:]\n",
    "eastN = x[320:400,:]\n",
    "eastS = x[400:,:]\n",
    "# otp = np.array(itemgetter(19,39,59)(x))\n",
    "otpctr = otp.mean(axis=0)\n",
    "wnctr = westN.mean(axis=0)\n",
    "wsctr = westS.mean(axis=0)\n",
    "wdctr = westD.mean(axis=0)\n",
    "enctr = eastN.mean(axis=0)\n",
    "esctr = eastS.mean(axis=0)\n",
    "\n",
    "# fit the OT-P fracture plane with 4 of the 80 points\n",
    "otpseedPnts = np.array([x[0,0],x[0,1],x[0,2],\n",
    "x[19,0],x[19,1],x[19,2],\n",
    "x[39,0],x[39,1],x[39,2],\n",
    "x[59,0],x[59,1],x[59,2]]).reshape(4,3)\n",
    "wnseedPnts = np.array([x[80,0],x[80,1],x[80,2],\n",
    "x[99,0],x[99,1],x[99,2],\n",
    "x[106,0],x[106,1],x[106,2],\n",
    "x[131,0],x[131,1],x[131,2],x[141,0],x[141,1],x[141,2]]).reshape(5,3)\n",
    "wsseedPnts = np.array([x[167,0],x[167,1],x[167,2],\n",
    "x[187,0],x[187,1],x[187,2],\n",
    "x[197,0],x[197,1],x[197,2],\n",
    "x[215,0],x[215,1],x[215,2],x[225,0],x[225,1],x[225,2]]).reshape(5,3)\n",
    "wdseedPnts = np.array([x[240,0],x[240,1],x[240,2],\n",
    "x[260,0],x[260,1],x[260,2],\n",
    "x[279,0],x[279,1],x[279,2],\n",
    "x[298,0],x[298,1],x[298,2]]).reshape(4,3)\n",
    "enseedPnts = np.array([x[397,0],x[397,1],x[397,2],\n",
    "x[377,0],x[377,1],x[377,2],\n",
    "x[365,0],x[365,1],x[365,2],\n",
    "x[347,0],x[347,1],x[347,2],x[336,0],x[336,1],x[336,2]]).reshape(5,3)\n",
    "esseedPnts = np.array([x[472,0],x[472,1],x[472,2],\n",
    "x[462,0],x[462,1],x[462,2],\n",
    "x[446,0],x[446,1],x[446,2],\n",
    "x[427,0],x[427,1],x[427,2],x[417,0],x[417,1],x[417,2]]).reshape(5,3)\n",
    "\n",
    "fracList = [f'OTP','West North','West South','West Deep','East North','East South']\n",
    "nplanes = len(fracList)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write a function to fit a plane to the boundary points of the fracture\n",
    "# this actually finds the normal vector to the plane\n",
    "def planeFit(points):\n",
    "    \"\"\"\n",
    "    p, n = planeFit(points)\n",
    "\n",
    "    Given an array, points, of shape (d,n)\n",
    "    points in d-dimensional space,\n",
    "    fit a d-dimensional plane to the points.\n",
    "    Return a point, p, on the plane (the point-cloud centroid),\n",
    "    and the normal vector to the plane, n.\n",
    "    \"\"\"\n",
    "    points = np.reshape(points, (np.shape(points)[0], -1)) # Collapse trailing dimensions\n",
    "    assert points.shape[0] <= points.shape[1], \"There are only {} points in {} dimensions.\".format(points.shape[1], points.shape[0])\n",
    "    ctr = points.mean(axis=1)\n",
    "    x = points - ctr[:,np.newaxis]\n",
    "    # M = np.dot(x, x.T) # Could also use np.cov(x) here.\n",
    "    M = np.cov(x)\n",
    "    normvec = svd(M)[0][:,-1]\n",
    "    # do the dot product of the center point and the normal vector to find the offset\n",
    "    fittedPlane = -ctr.dot(normvec)\n",
    "    return ctr, normvec, fittedPlane\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get the center and the normal vector to the planes\n",
    "otpctr,otpnormvec,otpfittedPlane = planeFit(otp.T)\n",
    "wnctr,wnnormvec,wnfittedPlane = planeFit(westN.T)\n",
    "wsctr,wsnormvec,wsfittedPlane= planeFit(westS.T)\n",
    "wdctr,wdnormvec,wdfittedPlane = planeFit(westD.T)\n",
    "enctr,ennormvec,enfittedPlane  = planeFit(eastN.T)\n",
    "esctr,esnormvec,esfittedPlane = planeFit(eastS.T)\n",
    "\n",
    "# create points in the plane using the plane EQ\n",
    "# Ax + By + Cz + D = 0 ... where A, B, C and D are the normal vector components and offset\n",
    "# use x and z points in the plane to find the y points\n",
    "# this works much better than using x and y to find z \n",
    "# y = -(Ax + Cz + D) / B\n",
    "def makeplanepoints(seedPnts,normvec,fittedPlane):\n",
    "    xb = np.sort(seedPnts[:,0].copy())\n",
    "    zb = np.sort(seedPnts[:,2].copy())\n",
    "    \n",
    "    xx,zz = np.meshgrid(xb,zb)\n",
    "\n",
    "    yy = -(normvec[0] * xx + normvec[2] * zz + fittedPlane) * 1. /normvec[1]\n",
    "\n",
    "    return xx,yy,zz\n",
    "\n",
    "otpxx,otpyy,otpzz = makeplanepoints(otp,otpnormvec,otpfittedPlane)\n",
    "wnxx,wnyy,wnzz = makeplanepoints(westN,wnnormvec,wnfittedPlane)\n",
    "wsxx,wsyy,wszz = makeplanepoints(westS,wsnormvec,wsfittedPlane)\n",
    "wdxx,wdyy,wdzz = makeplanepoints(westD,wdnormvec,wdfittedPlane)\n",
    "enxx,enyy,enzz = makeplanepoints(eastN,ennormvec,enfittedPlane)\n",
    "esxx,esyy,eszz = makeplanepoints(eastS,esnormvec,esfittedPlane)\n",
    "\n",
    "# code to make OTP bigger if we want to\n",
    "xf = np.array([otpxx[i] for i in list(sliced((0,0,-1,0,-1,-1,0,-1),2))])\n",
    "yf = np.array([otpyy[i] for i in list(sliced((0,0,-1,0,-1,-1,0,-1),2))])\n",
    "zf = np.array([otpzz[i] for i in list(sliced((0,0,-1,0,-1,-1,0,-1),2))])\n",
    "corners = np.concatenate((xf[:,np.newaxis],yf[:,np.newaxis],zf[:,np.newaxis]),axis=1)\n",
    "l1vec = corners[2,:] - corners[0,:]# type: ignore\n",
    "l2vec = corners[1,:] - corners[3,:]# type: ignore\n",
    "c1 =corners[0,:] - 0.5*(l1vec)\n",
    "c2 =corners[1,:] + 0.5*(l2vec)\n",
    "c3 =corners[0,:] + 1.5*(l1vec)\n",
    "c4 =corners[1,:] - 1.5*(l2vec)\n",
    "ncorners = np.vstack((c1,c2,c3,c4))\n",
    "cnt=0\n",
    "for i in list(sliced((0,0,-1,0,-1,-1,0,-1),2)):\n",
    "    otpxx[i] = ncorners[cnt,0]\n",
    "    otpyy[i] = ncorners[cnt,1]\n",
    "    otpzz[i] = ncorners[cnt,2]\n",
    "    cnt+=1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fillPlane(x,y,z,numels):\n",
    "    # find the corners of the plane using corner indices\n",
    "    xf = np.array([x[i] for i in list(sliced((0,0,-1,0,-1,-1,0,-1),2))])\n",
    "    yf = np.array([y[i] for i in list(sliced((0,0,-1,0,-1,-1,0,-1),2))])\n",
    "    zf = np.array([z[i] for i in list(sliced((0,0,-1,0,-1,-1,0,-1),2))])\n",
    "\n",
    "    corners = np.concatenate((xf[:,np.newaxis],yf[:,np.newaxis],zf[:,np.newaxis]),axis=1)\n",
    "    # make points in between left,top,right, and bottom axes\n",
    "    left = np.linspace(corners[0,:],corners[1,:],numels)\n",
    "    top = np.linspace(corners[1,:],corners[2,:],numels)[1:-1]\n",
    "    right = np.linspace(corners[3,:],corners[2,:],numels)\n",
    "    bot = np.linspace(corners[0,:],corners[3,:],numels)[1:-1]\n",
    "    # fill the interior points in the middle\n",
    "    fp = []\n",
    "    for i in range(numels-2):\n",
    "        fp.append(np.linspace(bot[i,:],top[i,:],numels)[1:-1])\n",
    "    # fp = np.array(fp).reshape(((numels-2)**2,corners.shape[1]))\n",
    "    fp = np.array(fp)\n",
    "\n",
    "    # build array correctly\n",
    "    midarr = []\n",
    "    for i in range(numels-2):\n",
    "        midarr.append(np.vstack((bot[i,:],fp[i,:,:],top[i,:])))\n",
    "    midarr = np.array(midarr).reshape((numels*(numels-2),3))\n",
    "\n",
    "    return np.vstack((left,midarr,right))\n",
    "    # return bot,fp,top\n",
    "\n",
    "numels = [12,8,8,5,8,8]\n",
    "otpsurf_pts = fillPlane(otpxx,otpyy,otpzz,numels[0])\n",
    "wnsurf_pts = fillPlane(wnxx,wnyy,wnzz,numels[1])\n",
    "wssurf_pts = fillPlane(wsxx,wsyy,wszz,numels[2])\n",
    "wdsurf_pts = fillPlane(wdxx,wdyy,wdzz,numels[3])\n",
    "ensurf_pts= fillPlane(enxx,enyy,enzz,numels[4])\n",
    "essurf_pts = fillPlane(esxx,esyy,eszz,numels[5])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a mesh for Ben Thompson's strain matrix function in cutde\n",
    "n_els_per_dim =[]\n",
    "[n_els_per_dim.append(numels[el] - 1 ) for el in range(len(numels))]\n",
    "allsurf_pts = [otpsurf_pts,wnsurf_pts,wssurf_pts,wdsurf_pts,ensurf_pts,essurf_pts]\n",
    "\n",
    "def maketris(numels):\n",
    "    surf_tris =[]\n",
    "    nx = ny = numels + 1\n",
    "    idx = lambda i, j: i * ny + j\n",
    "    for i in range(numels):\n",
    "        for j in range(numels):\n",
    "            # x1, x2 = mesh_xs[i : i + 2]\n",
    "            # y1, y2 = mesh_ys[j : j + 2]\n",
    "            surf_tris.append([idx(i, j), idx(i + 1, j), idx(i + 1, j + 1)])\n",
    "            surf_tris.append([idx(i, j), idx(i + 1, j + 1), idx(i, j + 1)])\n",
    "\n",
    "    surf_tris = np.array(surf_tris, dtype=np.int64)\n",
    "    return surf_tris\n",
    "\n",
    "otpsurf_tris = maketris(n_els_per_dim[0])\n",
    "wnsurf_tris = maketris(n_els_per_dim[1])\n",
    "wssurf_tris = maketris(n_els_per_dim[2])\n",
    "wdsurf_tris = maketris(n_els_per_dim[3])\n",
    "ensurf_tris = maketris(n_els_per_dim[4])\n",
    "essurf_tris = maketris(n_els_per_dim[5])\n",
    "\n",
    "allsurf_tris = [otpsurf_tris,wnsurf_tris,wssurf_tris,wdsurf_tris,ensurf_tris,essurf_tris]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fig = go.Figure()\n",
    "# # fig.add_scatter3d(x=np.array(isect)[:,0], y=np.array(isect)[:,1],z=np.array(isect)[:,2],mode='markers',\n",
    "# #   marker=dict(\n",
    "# #         size=8,\n",
    "# #         color='purple',\n",
    "# #         opacity=0.8))\n",
    "# fig.add_scatter3d(x=pdtpts[:,0], y=pdtpts[:,1],z=pdtpts[:,2],mode='markers',\n",
    "#   marker=dict(\n",
    "#         size=3,\n",
    "#         color=df_fullPDT[8462,:],\n",
    "#         colorscale='rdbu',\n",
    "#         cmid = 0,\n",
    "#         opacity=0.8))\n",
    "# fig.add_scatter3d(x=pdbpts[:,0], y=pdbpts[:,1],z=pdbpts[:,2],mode='markers',\n",
    "#   marker=dict(\n",
    "#         size=3,\n",
    "#         color=df_fullPDB[8462,:],\n",
    "#         colorscale='rdbu',\n",
    "#         cmid = 0,\n",
    "#         opacity=0.8))\n",
    "# fig.add_scatter3d(x=pstpts[:,0], y=pstpts[:,1],z=pstpts[:,2],mode='markers',\n",
    "#   marker=dict(\n",
    "#         size=3,\n",
    "#         color=df_fullPST[8462,:],\n",
    "#         colorscale='rdbu',\n",
    "#         cmid = 0,\n",
    "#         opacity=0.8))\n",
    "# fig.add_scatter3d(x=psbpts[:,0], y=psbpts[:,1],z=psbpts[:,2],mode='markers',\n",
    "#   marker=dict(\n",
    "#         size=3,\n",
    "#         color=df_fullPSB[8462,:],\n",
    "#         colorscale='rdbu',\n",
    "#         cmid = 0,\n",
    "#         opacity=0.8))\n",
    "# fig.add_scatter3d(x=obpts[:,0], y=obpts[:,1],z=obpts[:,2],mode='markers',\n",
    "#   marker=dict(\n",
    "#         size=3,\n",
    "#         color=df_fullOB[8462,:],\n",
    "#         colorscale='rdbu',\n",
    "#         cmid = 0,\n",
    "#         opacity=0.8))\n",
    "# fig.add_scatter3d(x=otpts[:,0], y=otpts[:,1],z=otpts[:,2],mode='markers',\n",
    "#   marker=dict(\n",
    "#         size=3,\n",
    "#         color=df_fullOT[8462,:],\n",
    "#         colorscale='rdbu',\n",
    "#         cmid = 0,\n",
    "#         colorbar=dict(thickness=20),\n",
    "#         opacity=0.8))\n",
    "# fig.add_trace(go.Mesh3d(x=otp[:,0],y=otp[:,1],z=otp[:,2],\n",
    "#         opacity=0.5,\n",
    "#         color='grey',\n",
    "#         text=[\"OT-P Connector\"],\n",
    "#         ))\n",
    "# fig.add_trace(go.Mesh3d(x=westN[:,0],y=westN[:,1],z=westN[:,2],\n",
    "#         opacity=0.5,\n",
    "#         color='grey',\n",
    "#         text=[\"WN\"],\n",
    "#         ))\n",
    "# fig.add_trace(go.Mesh3d(x=westS[:,0],y=westS[:,1],z=westS[:,2],\n",
    "#         opacity=0.5,\n",
    "#         color='grey',\n",
    "#         text=[\"WS\"],\n",
    "#         ))\n",
    "# fig.add_trace(go.Mesh3d(x=westD[:,0],y=westD[:,1],z=westD[:,2],\n",
    "#         opacity=0.5,\n",
    "#         color='grey',\n",
    "#         text=[\"WD\"],\n",
    "#         ))\n",
    "# fig.add_trace(go.Mesh3d(x=eastN[:,0],y=eastN[:,1],z=eastN[:,2],\n",
    "#         opacity=0.5,\n",
    "#         color='grey',\n",
    "#         text=[\"EN\"],\n",
    "#         ))\n",
    "# fig.add_trace(go.Mesh3d(x=eastS[:,0],y=eastS[:,1],z=eastS[:,2],\n",
    "#         opacity=0.5,\n",
    "#         color='grey',\n",
    "#         text=[\"ES\"],\n",
    "#         ))\n",
    "# # fig.add_trace(go.Mesh3d(x=westN[:,0],y=westN[:,1],z=westN[:,2],\n",
    "# #         opacity=0.5,\n",
    "# #         color='magenta',\n",
    "# #         text=[\"wncircle\"],\n",
    "# #         ))\n",
    "# # fig.add_scatter3d(x=otpsurf_pts[:,0], y=otpsurf_pts[:,1],z=otpsurf_pts[:,2],mode='markers',\n",
    "# #   marker=dict(\n",
    "# #         size=3,\n",
    "# #         color='pink',\n",
    "# #         opacity=0.7))\n",
    "# # fig.add_scatter3d(x=ensurf_pts[:,0], y=ensurf_pts[:,1],z=ensurf_pts[:,2],mode='markers',\n",
    "# #   marker=dict(\n",
    "# #         size=3,\n",
    "# #         color='green',\n",
    "# #         opacity=0.4))\n",
    "# fig.update_layout(\n",
    "#         autosize=False,\n",
    "#         width=1200,\n",
    "#         height=800,\n",
    "#         )\n",
    "\n",
    "\n",
    "# fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create figure\n",
    "fig = go.Figure(\n",
    "    data=[go.Scatter3d(x=pdtpts[:,0], y=pdtpts[:,1],z=pdtpts[:,2],mode='markers',\n",
    "        marker=dict(\n",
    "            size=3,\n",
    "            color=df_fullPDT[0,:pdtpts.shape[0]],\n",
    "            colorscale='rdbu',\n",
    "            cmid = 0,\n",
    "            opacity=0.8)),\n",
    "          go.Scatter3d(x=pdbpts[:,0], y=pdbpts[:,1],z=pdbpts[:,2],mode='markers',\n",
    "        marker=dict(\n",
    "            size=3,\n",
    "            color=df_fullPDB[0,:],\n",
    "            colorscale='rdbu',\n",
    "            cmid = 0,\n",
    "            opacity=0.8)),\n",
    "        go.Scatter3d(x=psbpts[:,0], y=psbpts[:,1],z=psbpts[:,2],mode='markers',\n",
    "        marker=dict(\n",
    "            size=3,\n",
    "            color=df_fullPSB[0,:],\n",
    "            colorscale='rdbu',\n",
    "            cmid = 0,\n",
    "            opacity=0.8)),\n",
    "        go.Scatter3d(x=pstpts[:,0], y=pstpts[:,1],z=pstpts[:,2],mode='markers',\n",
    "        marker=dict(\n",
    "            size=3,\n",
    "            color=df_fullPST[0,:],\n",
    "            colorscale='rdbu',\n",
    "            cmid = 0,\n",
    "            opacity=0.8)),\n",
    "        go.Scatter3d(x=otpts[:,0], y=otpts[:,1],z=otpts[:,2],mode='markers',\n",
    "        marker=dict(\n",
    "            size=3,\n",
    "            color=df_fullOT[0,:otpts.shape[0]],\n",
    "            colorscale='rdbu',\n",
    "            cmid = 0,\n",
    "            opacity=0.8)),\n",
    "        go.Scatter3d(x=obpts[:,0], y=obpts[:,1],z=obpts[:,2],mode='markers',\n",
    "        marker=dict(\n",
    "            size=3,\n",
    "            color=df_fullOB[0,:],\n",
    "            colorscale='balance',\n",
    "            cmid = 0,\n",
    "            opacity=0.8)),\n",
    "                     ],\n",
    "    layout=go.Layout(\n",
    "        title_text=f'Strain Rates', hovermode=\"closest\",\n",
    "        updatemenus=[dict(type=\"buttons\",\n",
    "                          buttons=[dict(label=\"Play\",\n",
    "                                        method=\"animate\",\n",
    "                                        args=[None])])]),\n",
    "    frames=[go.Frame(\n",
    "        data=[go.Scatter3d(\n",
    "            x=pdtpts[:,0],\n",
    "            y=pdtpts[:,1],\n",
    "            z=pdtpts[:,2],\n",
    "            mode='markers',\n",
    "            marker=dict(\n",
    "            size=3,\n",
    "            color=df_fullPDT[k,:otpts.shape[0]],\n",
    "            colorscale='rdbu',\n",
    "            cmid = 0,\n",
    "            colorbar=dict(thickness=20),\n",
    "            opacity=1)),\n",
    "            ],layout=go.Layout(title_text=f'Strain Rates {DASdates[k]}'))\n",
    "        for k in np.arange(5550,df_fullOT.shape[0],50)]\n",
    "\n",
    "\n",
    ")\n",
    "fig.add_trace(go.Mesh3d(x=otp[:,0],y=otp[:,1],z=otp[:,2],\n",
    "        opacity=0.25,\n",
    "        color='grey',\n",
    "        text=[\"OT-P Connector\"],\n",
    "        ))\n",
    "fig.add_trace(go.Mesh3d(x=westN[:,0],y=westN[:,1],z=westN[:,2],\n",
    "        opacity=0.25,\n",
    "        color='grey',\n",
    "        text=[\"WN\"],\n",
    "        ))\n",
    "fig.add_trace(go.Mesh3d(x=westS[:,0],y=westS[:,1],z=westS[:,2],\n",
    "        opacity=0.25,\n",
    "        color='grey',\n",
    "        text=[\"WS\"],\n",
    "        ))\n",
    "fig.add_trace(go.Mesh3d(x=westD[:,0],y=westD[:,1],z=westD[:,2],\n",
    "        opacity=0.25,\n",
    "        color='grey',\n",
    "        text=[\"WD\"],\n",
    "        ))\n",
    "fig.add_trace(go.Mesh3d(x=eastN[:,0],y=eastN[:,1],z=eastN[:,2],\n",
    "        opacity=0.25,\n",
    "        color='grey',\n",
    "        text=[\"EN\"],\n",
    "        ))\n",
    "fig.add_trace(go.Mesh3d(x=eastS[:,0],y=eastS[:,1],z=eastS[:,2],\n",
    "        opacity=0.25,\n",
    "        color='grey',\n",
    "        text=[\"ES\"],\n",
    "        ))\n",
    "fig.update_layout(\n",
    "    autosize=False,\n",
    "    width=1000,\n",
    "    height=800,\n",
    ")\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Call strain matrix function from Ben's FullSpace code \n",
    "# output dimensions are strain_mat = (num_obs, strain_tensor_component, num_triangles, slip mode)\n",
    "# strain tensor comps: 0 = xx, 1 = yy, 2 = zz, 3 = xy, 4 = xz, 5 = yz\n",
    "# slip modes: 0 = strike slip, 1 = dip slip, 2 = tensile\n",
    "\n",
    "allstrain_mat = FS.strain_matrix(obs_pts=allpts, tris=np.concatenate([allsurf_pts[i][tri] for i,tri in enumerate(allsurf_tris)]),nu=0.24)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def genSlice(sliceLen,nummodes):\n",
    "    '''Generate a list of slices\n",
    "    for grabbing the number of \n",
    "    triangles in the mesh for each slip mode'''\n",
    "\n",
    "    first = int(sliceLen / nummodes)\n",
    "    second = int(first * 2)\n",
    "    third  = int(sliceLen)\n",
    "    \n",
    "    lst = list(sliced((int(0),first,first,second,second,third),2))\n",
    "    return lst\n",
    "def fracSlice(start,numfracs):\n",
    "    '''Generate a list of slices\n",
    "    for grabbing the number of \n",
    "    triangles in the mesh for each indivial fracture'''\n",
    "    first = start + 2*( numfracs[0]**2 )\n",
    "    second = first + 2*( numfracs[1]**2 ) \n",
    "    third  = second + 2*( numfracs[2]**2 ) \n",
    "    fourth = third + 2*( numfracs[3]**2 ) \n",
    "    fifth = fourth + 2*( numfracs[4]**2 ) \n",
    "    sixth = fifth + 2*( numfracs[5]**2 )\n",
    "    \n",
    "    lst = list(sliced((start,first,first,second,second,third,third,fourth,fourth,fifth,fifth,sixth),2))\n",
    "    return lst\n",
    "\n",
    "# slices for splitting the model estimates into their respective modes (strike, dip, tensile)\n",
    "mode = genSlice(allstrain_mat.shape[2]*allstrain_mat.shape[3],3)\n",
    "# slices for splitting elements into their respective fracture plane\n",
    "ifrac = fracSlice(int(mode[0][0]),n_els_per_dim)\n",
    "jfrac = fracSlice(int(mode[1][0]),n_els_per_dim)\n",
    "kfrac = fracSlice(int(mode[2][0]),n_els_per_dim)\n",
    "mfrac = list(zip(ifrac,jfrac,kfrac))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function for finding tangent vectors along wells to use for rotating strains \n",
    "\n",
    "def getRotVecs(points):\n",
    "    # points in N x 3 \n",
    "    # convert points to array if list or dataframe\n",
    "    points = np.array(points)\n",
    "    # find the derivative of the original points (tangent)\n",
    "    d_points = np.gradient(points,axis=0)\n",
    "    # divide by length to get unit tangent\n",
    "    tangents = d_points / norm(d_points,ord=2)\n",
    "    # find the derivative of the tangents to get normal\n",
    "    # d_tans = np.gradient(tangents,axis=1)\n",
    "    # # divide by length to get unit normal\n",
    "    # normals = d_tans / norm(d_tans,ord=2)\n",
    "    # # take cross product to get bi normal\n",
    "    # bnorms = np.cross(tangents,normals)\n",
    "    return tangents\n",
    "    # return tangents,normals,bnorms\n",
    "\n",
    "pdttans = getRotVecs(pdtpts)\n",
    "pdbtans = getRotVecs(pdbpts)\n",
    "psttans = getRotVecs(pstpts)\n",
    "psbtans = getRotVecs(psbpts)\n",
    "ottans = getRotVecs(otpts)\n",
    "obtans = getRotVecs(obpts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function for doing the strain rotations given the tangent vectors and strains \n",
    "# ( See Liam Chiang / Paul Sava paper \" High resolution multi-component DAS\" )\n",
    "# Equations in Appendix\n",
    "\n",
    "def rotStrains(tanvecs,strains):\n",
    "    R_11 = tanvecs[:,0]\n",
    "    R_12 = tanvecs[:,1]\n",
    "    R_13 = tanvecs[:,2]\n",
    "    \n",
    "    rotmat = np.zeros((strains.shape[0],strains.shape[2],strains.shape[3]))\n",
    "    G = np.array([R_11**2, R_12**2, R_13**2, 2*R_11*R_12, 2*R_11*R_13, 2*R_12*R_13]).T\n",
    "    for i in range(strains.shape[3]):\n",
    "        for j in range(strains.shape[0]):\n",
    "            \n",
    "            rotmat[j,:,i] = G[j,:] @ strains[j,:,:,i]\n",
    "    return rotmat\n",
    "\n",
    "first = len(pdtpts)\n",
    "second = len(pdtpts)+len(pdbpts)\n",
    "third = second + len(pstpts)\n",
    "fourth = third + len(psbpts)\n",
    "fifth = fourth + len(otpts)\n",
    "sixth = fifth + len(obpts)\n",
    "lenvec = list(sliced((int(0),first,first,second,second,third,third,fourth,fourth,fifth,fifth,sixth),2))\n",
    "\n",
    "# Do the rotations\n",
    "\n",
    "pdtstrain_mat_rot = rotStrains(pdttans,allstrain_mat[lenvec[0][0]:lenvec[0][1],:,:,:])\n",
    "pdbstrain_mat_rot = rotStrains(pdbtans,allstrain_mat[lenvec[1][0]:lenvec[1][1],:,:,:])\n",
    "pststrain_mat_rot = rotStrains(psttans,allstrain_mat[lenvec[2][0]:lenvec[2][1],:,:,:])\n",
    "psbstrain_mat_rot = rotStrains(psbtans,allstrain_mat[lenvec[3][0]:lenvec[3][1],:,:,:])\n",
    "otstrain_mat_rot = rotStrains(ottans,allstrain_mat[lenvec[4][0]:lenvec[4][1],:,:,:])\n",
    "obstrain_mat_rot = rotStrains(obtans,allstrain_mat[lenvec[5][0]:lenvec[5][1],:,:,:])\n",
    "\n",
    "# Combine the rotated strain matrices and make the G matrix for inversion\n",
    "allstrain_mat_rot = np.concatenate((pdtstrain_mat_rot,pdbstrain_mat_rot,pststrain_mat_rot,psbstrain_mat_rot,otstrain_mat_rot,obstrain_mat_rot))\n",
    "allG = np.hstack((allstrain_mat_rot[:,:,0],allstrain_mat_rot[:,:,1],allstrain_mat_rot[:,:,2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Definitions for Ridge Regression \n",
    "\n",
    "def ridge(X, y, l2):\n",
    "    \"\"\"Ridge Regression model with intercept term.\n",
    "    L2 penalty and intercept term included via design matrix augmentation.\n",
    "    Params:\n",
    "        X - NumPy matrix, size (N, p), of numerical predictors\n",
    "        y - NumPy array, length N, of numerical response\n",
    "        l2 - L2 penalty tuning parameter (positive scalar) \n",
    "    Returns:\n",
    "        NumPy array, length p + 1, of fitted model coefficients\n",
    "    \"\"\"\n",
    "    m, n = np.shape(X)\n",
    "    upper_half = np.hstack((np.ones((m, 1)), X))\n",
    "    # upper_half = X\n",
    "    lower = np.zeros((n, n))\n",
    "    np.fill_diagonal(lower, np.sqrt(l2))\n",
    "    lower_half = np.hstack((np.zeros((n, 1)), lower))\n",
    "    # lower_half = lower\n",
    "    X = np.vstack((upper_half, lower_half))\n",
    "    y = np.append(y, np.zeros(n))\n",
    "    if m == n:\n",
    "        soln = solve(np.dot(X.T, X), np.dot(X.T, y))\n",
    "    else:\n",
    "        soln = lstsq(np.dot(X.T, X), np.dot(X.T, y))\n",
    "        U,s,Vt = svd(np.dot(X.T, X),lapack_driver='gesvd')\n",
    "        # soln = lstsq(X, y)\n",
    "        # U,s,Vt = svd(X)\n",
    "        \n",
    "    return soln, U, s, Vt, X, y\n",
    "\n",
    "def getCond(X, y, l2):\n",
    "    \"\"\" Find condition number of augmented matrix for Ridge \n",
    "    L2 penalty and intercept term included via design matrix augmentation.\n",
    "    Params:\n",
    "        X - NumPy matrix, size (N, p), of numerical predictors\n",
    "        y - NumPy array, length N, of numerical response\n",
    "        l2 - L2 penalty tuning parameter (positive scalar) \n",
    "    Returns:\n",
    "        NumPy array, length p + 1, of fitted model coefficients\n",
    "    \"\"\"\n",
    "    m, n = np.shape(X)\n",
    "    upper_half = np.hstack((np.ones((m, 1)), X))\n",
    "    lower = np.zeros((n, n))\n",
    "    np.fill_diagonal(lower, np.sqrt(l2))\n",
    "    lower_half = np.hstack((np.zeros((n, 1)), lower))\n",
    "    X = np.vstack((upper_half, lower_half))\n",
    "    \n",
    "    condition_num = np.linalg.cond(np.dot(X.T, X))\n",
    "    return condition_num"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "das_time_ind = 5785 # 2018-05-24 22:38:09\n",
    "# das_time_ind = 7782 # 2018-05-25 15:18:48\n",
    "# das_time_ind = 8462 # 2018-05-25 20:58:48\n",
    "print(DASdates[das_time_ind])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdtdat = df_strainPDT[das_time_ind,:pdtpts.shape[0]] \n",
    "pdbdat = df_strainPDB[das_time_ind,:pdbpts.shape[0]] \n",
    "pstdat = df_strainPST[das_time_ind,:pstpts.shape[0]] \n",
    "psbdat = df_strainPSB[das_time_ind,:psbpts.shape[0]] \n",
    "otdat = df_strainOT[das_time_ind,:otpts.shape[0]] \n",
    "obdat = df_strainOB[das_time_ind,:obpts.shape[0]] \n",
    "\n",
    "pdtdat1 = df_strainPDT[das_time_ind,:pdtpts.shape[0]] \n",
    "pdbdat1 = df_strainPDB[das_time_ind,:pdbpts.shape[0]] \n",
    "pstdat1 = df_strainPST[das_time_ind,:pstpts.shape[0]] \n",
    "psbdat1 = df_strainPSB[das_time_ind,:psbpts.shape[0]] \n",
    "otdat1 = df_strainOT[das_time_ind,:otpts.shape[0]] \n",
    "obdat1 = df_strainOB[das_time_ind,:obpts.shape[0]] \n",
    "\n",
    "pdtdat = pdtdat - (1.056e-5 * pdtdts[:pdtpts.shape[0]]) * 1e9\n",
    "pdbdat = pdbdat - (1.056e-5 * pdbdts[:pdbpts.shape[0]]) * 1e9 \n",
    "pstdat = pstdat - (.0561e-5 * pstdts[:pstpts.shape[0]]) * 1e9\n",
    "psbdat = psbdat - (1.056e-5 * psbdts[:psbpts.shape[0]]) * 1e9\n",
    "otdat  = otdat - (1.056e-5 * otdts[:otpts.shape[0]]) * 1e9 \n",
    "obdat  = obdat - (1.056e-5 * obdts[:obpts.shape[0]]) * 1e9 \n",
    "Datcomb = np.hstack((detrend(pdtdat,type='constant'),\n",
    "                     detrend(pdbdat,type='constant'),\n",
    "                     detrend(pstdat,type='constant'),\n",
    "                     detrend(psbdat,type='constant'),\n",
    "                     detrend(otdat,type='constant'),\n",
    "                     detrend(obdat,type='constant')))\n",
    "Datcomb1 = np.hstack((pdtdat1,pdbdat1,pstdat1,psbdat1,otdat1,obdat1))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pdttemp= (1.056e-5 * pdtdts[:pdtpts.shape[0]]) * 1e9\n",
    "pdbtemp = (1.056e-5 * pdbdts[:pdbpts.shape[0]]) * 1e9 \n",
    "psttemp = (.0561e-5 * pstdts[:pstpts.shape[0]]) * 1e9\n",
    "psbtemp= (1.056e-5 * psbdts[:psbpts.shape[0]]) * 1e9\n",
    "ottemp  = (1.056e-5 * otdts[:otpts.shape[0]]) * 1e9 \n",
    "obtemp = (1.056e-5 * obdts[:obpts.shape[0]]) * 1e9 \n",
    "plt.figure()\n",
    "plt.plot(pdttemp,label='pdt')\n",
    "plt.plot(pdbtemp,label='pdb')\n",
    "plt.plot(psttemp,label='pst')\n",
    "plt.plot(psbtemp,label='psb')\n",
    "plt.plot(ottemp,label='ot')\n",
    "plt.plot(obtemp,label='ob')\n",
    "plt.xlabel('Channel')\n",
    "plt.ylabel('Thermo-optical Strain')\n",
    "plt.legend()\n",
    "\n",
    "stdstack = np.std(np.vstack([pdttemp,pdbtemp,psbtemp,ottemp,obtemp]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tick_labels=['30','\\n\\nPDT','90','\\n\\nPDB','144','\\n\\nPST','194','\\n\\nPSB','254','\\n\\nOT','314','\\n\\nOB']\n",
    "\n",
    "tick_locations = np.array((30,90,144,194,254,314))\n",
    "ptLen = np.cumsum((pdtpts.shape[0],pdbpts.shape[0],pstpts.shape[0],psbpts.shape[0],otpts.shape[0],obpts.shape[0]))\n",
    "new_labels = [ ''.join(x) for x in zip(tick_labels[0::2], tick_labels[1::2]) ]\n",
    "fig,ax = plt.subplots(figsize=(10,8))\n",
    "ax.plot(Datcomb,label='temp corrected')\n",
    "ax.fill_between(range(len(Datcomb)),Datcomb+stdstack,Datcomb-stdstack,color='black',label='temp corrected')\n",
    "# ax.fill_between(Datcomb-stdstack,color='black',label='temp corrected')\n",
    "ax.plot(Datcomb1,label='non-temp corrected')\n",
    "ax.set_ylabel('Strain')\n",
    "ax.set_xlabel('Well Center')\n",
    "ax.set_xticks(tick_locations, new_labels)\n",
    "ax.axvspan(0,ptLen[0],facecolor='green',alpha = 0.2)\n",
    "ax.axvspan(ptLen[0],ptLen[1],facecolor='yellow',alpha = 0.2)\n",
    "ax.axvspan(ptLen[1],ptLen[2],facecolor='red',alpha = 0.2)\n",
    "ax.axvspan(ptLen[2],ptLen[3],facecolor='blue',alpha = 0.2)\n",
    "ax.axvspan(ptLen[3],ptLen[4],facecolor='orange',alpha = 0.2)\n",
    "ax.axvspan(ptLen[4],ptLen[5],facecolor='grey',alpha = 0.2)\n",
    "ax.set_title('DAS Data Temp Corrections')\n",
    "\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "plt.figure()\n",
    "lam = np.logspace(-13,-3,25)\n",
    "cn = np.zeros_like(lam)\n",
    "sn = np.zeros_like(lam)\n",
    "SSE = np.zeros_like(lam)\n",
    "var_res = np.zeros_like(lam)\n",
    "# MSE = np.zeros_like(lam)\n",
    "for i in range(len(lam)):\n",
    "    # lam = 0.000001\n",
    "    testmest,_,_,_,X,y = ridge(allG,Datcomb1,lam[i])\n",
    "    pred_dat = X @ testmest[0]\n",
    "    # print(np.mean(Datcomb-pred_dat[:349]))\n",
    "    # print(np.std(Datcomb-pred_dat[:349]))\n",
    "    # cn[i] = getCond(allG,Datcomb,lam[i])\n",
    "    var_res[i] = np.var((X @ testmest[0] - y ))\n",
    "    # SSE[i] = np.sum(((Datcomb - allG @ testmest[0]))**2)\n",
    "    SSE[i] = norm(X @ testmest[0] - y)\n",
    "    sn[i] = norm(testmest[0])\n",
    "    # MSE[i] = np.mean(((allG @ testmest[0][1:]) - Datcomb)**2)\n",
    "# plt.semilogx(lam,MSE,'.r',label='MSE')\n",
    "# plt.semilogx(lam,SSE,'.g',label='SSE')\n",
    "# plt.loglog(lam,cn)\n",
    "plt.loglog(SSE,sn,'.')\n",
    "# plt.xlabel('log(Lambda)')\n",
    "plt.xlabel('log(Residual)')\n",
    "# plt.ylabel('Matrix Condition Number')\n",
    "plt.ylabel('log(norm(solutions))')\n",
    "# plt.title('log(Lambda) vs Condition Number')\n",
    "plt.title('L-curve')\n",
    "plt.legend()\n",
    "print(lam)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lam = np.logspace(-13,-3,25)\n",
    "fig,ax = plt.subplots()\n",
    "ax2= ax.twinx()\n",
    "fig.suptitle('Bias / Variance Tradeoff')\n",
    "ax2.loglog(lam,sn,label='solution norm')\n",
    "ax2.set_ylabel('solution norm')\n",
    "ax2.legend(loc='center right')\n",
    "ax.loglog(lam,SSE,color='orange',label='residual norm ')\n",
    "ax.set_ylabel('norm of residual')\n",
    "# ax.set_ylim(4e3,1e5)\n",
    "ax.legend(loc='center left')\n",
    "ax.set_xlabel('Lambda')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# lam = 3.5e-11\n",
    "# lam = 3.8e-9\n",
    "# lam = 8.9e-8\n",
    "lam = 1e-8\n",
    "lamT = 4.2e-8\n",
    "\n",
    "# pdtdat = df_strainPDT[das_time_ind,:pdtpts.shape[0]] \n",
    "# pdbdat = df_strainPDB[das_time_ind,:pdbpts.shape[0]] \n",
    "# pstdat = df_strainPST[das_time_ind,:pstpts.shape[0]] \n",
    "# psbdat = df_strainPSB[das_time_ind,:psbpts.shape[0]] \n",
    "# otdat = df_strainOT[das_time_ind,:otpts.shape[0]] \n",
    "# obdat = df_strainOB[das_time_ind,:obpts.shape[0]] \n",
    "# Datcomb = np.hstack((pdtdat,pdbdat,pstdat,psbdat,otdat,obdat))\n",
    "\n",
    "mest,U,s,Vt, G, daty = ridge(allG, Datcomb1,lam)\n",
    "mestTemp,Utemp,stemp,Vtemp, Gtemp, datytemp = ridge(allG, Datcomb,lamT)\n",
    "\n",
    "# U = U[:-1,:-1]\n",
    "# s = s[:-1]\n",
    "# Vt = Vt[:-1,:-1]\n",
    "# res = mest[1]\n",
    "mest = mest[0][1:]\n",
    "mestTemp = mestTemp[0][1:]\n",
    "# daty = daty[:len(mest)]\n",
    "\n",
    "# G = G[:,:-1]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do the inverse problem for several time slices\n",
    "# lam = 6.9e-8\n",
    "# strainList = [df_strainPDT,df_strainPDB,df_strainPST,df_strainPSB,df_strainOT,df_strainOB]\n",
    "# pntShapeList = [pdtpts.shape[0],pdbpts.shape[0],pstpts.shape[0],psbpts.shape[0],otpts.shape[0],obpts.shape[0]]\n",
    "# das_time_ind = np.arange(7765,7816,1)\n",
    "# def doInverse(timeEpochs,strainList,Gmat,lam,pntShapeList):\n",
    "#     mest = []\n",
    "#     for i in timeEpochs:\n",
    "#         pdtdat = strainList[0][i,:pdtpts.shape[0]] \n",
    "#         pdbdat = strainList[1][i,:pdbpts.shape[0]] \n",
    "#         pstdat = strainList[2][i,:pstpts.shape[0]] \n",
    "#         psbdat = strainList[3][i,:psbpts.shape[0]] \n",
    "#         otdat = strainList[4][i,:otpts.shape[0]] \n",
    "#         obdat = strainList[5][i,:obpts.shape[0]]\n",
    "        \n",
    "#         pdtdat = pdtdat - (1.056e-5 * pdtdts[:pdtpts.shape[0]]) * 1e9\n",
    "#         pdbdat = pdbdat - (1.056e-5 * pdbdts[:pdbpts.shape[0]]) * 1e9 \n",
    "#         pstdat = pstdat - (.0561e-5 * pstdts[:pstpts.shape[0]]) * 1e9\n",
    "#         psbdat = psbdat - (1.056e-5 * psbdts[:psbpts.shape[0]]) * 1e9\n",
    "#         otdat  = otdat - (1.056e-5 * otdts[:otpts.shape[0]]) * 1e9 \n",
    "#         obdat  = obdat - (1.056e-5 * obdts[:obpts.shape[0]]) * 1e9 \n",
    "#         Datcomb = np.hstack((detrend(pdtdat,type='constant'),\n",
    "#                      detrend(pdbdat,type='constant'),\n",
    "#                      detrend(pstdat,type='constant'),\n",
    "#                      detrend(psbdat,type='constant'),\n",
    "#                      detrend(otdat,type='constant'),\n",
    "#                      detrend(obdat,type='constant')))\n",
    "#         tmp,_,_,_,_,_ = ridge(Gmat,Datcomb,lam)\n",
    "#         mest.append(tmp[0][1:])\n",
    "#     return mest\n",
    "\n",
    "# mest = doInverse(das_time_ind,strainList,allG,lam,pntShapeList)\n",
    "# plt.figure()\n",
    "# for i in range(len(mest)):\n",
    "#     plt.plot(allG@mest[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "U,s,Vt = svd(allG,full_matrices=False)\n",
    "# G1 = Vt.T @ np.linalg.inv(np.diag(s)) @ U.T\n",
    "G1 = Vt[:100,:].T @ np.linalg.inv(np.diag(s[:100])) @ U[:,:100].T\n",
    "mest1 = G1 @ Datcomb1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(1,2)\n",
    "ax[0].plot(mest,label='Ridge')\n",
    "ax[0].legend()\n",
    "ax[1].plot(mest1,label='SVD')\n",
    "ax[1].legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(1,2)\n",
    "fig.suptitle('Data fits: predicted v. observed')\n",
    "ax[0].plot(allG @ mest,label='ridge')\n",
    "ax[0].set_ylabel('model estimates')\n",
    "ax[0].set_xlabel('data index (well) ')\n",
    "ax[0].plot(Datcomb1,label='observed')\n",
    "ax[1].plot(allG @ mest1,label='svd')\n",
    "ax[1].set_xlabel('data index (well) ')\n",
    "ax[1].plot(Datcomb1,label='observed')\n",
    "ax[0].legend()\n",
    "ax[1].legend()\n",
    "fig.tight_layout()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pred_dat1 = allG @ mest1\n",
    "pred_dat = allG @ mest\n",
    "print(f'Ridge resid mean is {np.mean(Datcomb1-pred_dat)}')\n",
    "print(f'Ridge resid std is {np.std(Datcomb1-pred_dat)}')\n",
    "print(f'SVD resid mean is {np.mean(Datcomb1-pred_dat1)}')\n",
    "print(f'SVD resid std is {np.std(Datcomb1-pred_dat1)}')\n",
    "print(f'temp corrected Ridge resid mean is {np.mean(Datcomb-pred_dat)}')\n",
    "print(f'temp correctedRidge resid std is {np.std(Datcomb-pred_dat)}')\n",
    "print(f'temp correctedSVD resid mean is {np.mean(Datcomb-pred_dat1)}')\n",
    "print(f'temp correctedSVD resid std is {np.std(Datcomb-pred_dat1)}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate data and model resolution\n",
    "\n",
    "Gg = np.linalg.inv(G[:,1:].T @ G[:,1:] ) @ allG.T \n",
    "var_res = np.var((daty - G[:,1:] @ mest ))\n",
    "W = np.linalg.inv(G[:,1:].T @ G[:,1:] ) @ (allG.T @ allG)\n",
    "covm = var_res * W @ np.linalg.inv(allG.T @ allG) @ W.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots(1,2)\n",
    "im1 = ax[0].plot(np.diag(covm))\n",
    "im2 = ax[1].imshow(covm, aspect = 'auto')\n",
    "ax[1].set_aspect('equal',adjustable='box')\n",
    "cbar1 = plt.colorbar(im2, ax=ax[1],label='Model Covariance')\n",
    "ax[1].invert_yaxis()\n",
    "# im3 = ax[2].imshow(covd, aspect = 'auto',vmin=0,vmax=0.1)\n",
    "# cbar2 = plt.colorbar(im3, ax=ax[2],label='Data Covariance')\n",
    "# im4 = ax[3].imshow(covu, aspect = 'auto')\n",
    "# cbar3 = plt.colorbar(im4, ax=ax[3],label='Unit Model Covariance')\n",
    "fig.tight_layout()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "titles=[f'Strike Slip ','Dip Slip','Tensile ']\n",
    "fig,ax = plt.subplots(nplanes,3,figsize=(12,10))\n",
    "for i in range(3):\n",
    "    for k in range(nplanes):\n",
    "\n",
    "        # vm = np.percentile(mest[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],90)\n",
    "        vm = np.percentile(mest[mode[i][0]:mode[i][1]],99)\n",
    "        # covvm = np.percentile(modcov[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],99)\n",
    "        # vm = np.percentile(mest[mode[i][0]:mode[i][1]],99)\n",
    "\n",
    "\n",
    "        # im  = ax[k,i].tripcolor(allsurf_pts[pfrac[k][0]:pfrac[k][1], 0], allsurf_pts[pfrac[k][0]:pfrac[k][1], 2],\n",
    "        #                         allsurf_tris[ifrac[k][0]:ifrac[k][1]],\n",
    "        #                         facecolors=mest[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],\n",
    "        #                         cmap='seismic',vmin=-50000,vmax=50000) # type: ignore\n",
    "        \n",
    "        im  = ax[k,i].tripcolor(np.array(allsurf_pts[k])[:,0], np.array(allsurf_pts[k])[:,2],\n",
    "                        np.array(allsurf_tris[k]),\n",
    "                        facecolors=mest[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],\n",
    "                        cmap='seismic',\n",
    "                        vmin = -vm,\n",
    "                        vmax = vm)\n",
    "                        # vmin=min(mest[mode[i][0]:mode[i][1]]),\n",
    "                        # vmax=max(mest[mode[i][0]:mode[i][1]])) # type: ignore\n",
    "        # im  = ax[k,i].tripcolor(np.array(allsurf_pts[k])[:,0], np.array(allsurf_pts[k])[:,2],\n",
    "        #         np.array(allsurf_tris[k]),\n",
    "        #         facecolors=mest[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],\n",
    "        #         cmap='cmr.lilac_r') # type: ignore\n",
    "        # im  = ax[k,i].tripcolor(np.array(allsurf_pts[k])[:,0], np.array(allsurf_pts[k])[:,2],\n",
    "        #         np.array(allsurf_tris[k]),\n",
    "        #         facecolors=modcov[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],\n",
    "        #         cmap='cmr.lilac_r',vmin=94000,vmax=95000) # type: ignore\n",
    "        ax[k,i].set_xlabel('X direction (m)')\n",
    "        ax[k,i].set_ylabel('Z direction (m)')\n",
    "        # ax[i].axis('equal')\n",
    "        # ax[i].set_aspect('equal','box')\n",
    "        plt.colorbar(im,ax=ax[k,i])\n",
    "        ax[k,i].set_title(f'{fracList[k]} {titles[i]}')\n",
    "        fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine model estimates over each mode of slip to get total slip magnitude\n",
    "\n",
    "ss = []\n",
    "ds = []\n",
    "ts = []\n",
    "mode = genSlice(len(mest),3)\n",
    "ifrac = fracSlice(int(mode[0][0]),n_els_per_dim)\n",
    "for k in range(nplanes):\n",
    "\n",
    "    ss.append(mest[mode[0][0]:mode[0][1]][ifrac[k][0]:ifrac[k][1]])\n",
    "    ds.append(mest[mode[1][0]:mode[1][1]][ifrac[k][0]:ifrac[k][1]])\n",
    "    ts.append(mest[mode[2][0]:mode[2][1]][ifrac[k][0]:ifrac[k][1]])\n",
    "\n",
    "zipped_list = list(zip(ss,ds,ts))\n",
    "allslip = []\n",
    "i=0\n",
    "for frac in range(len(zipped_list)):\n",
    "    tmp2 = []\n",
    "    for elem in range(len(zipped_list[frac][i])):\n",
    "        tmp = []\n",
    "        for mode in range(len(zipped_list[0])):\n",
    "            tmp.append(zipped_list[frac][mode][elem]**2)\n",
    "        tmp2.append(np.sqrt(sum(tmp)))\n",
    "        i =+ 1\n",
    "    allslip.append(tmp2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look at some confidence intervals\n",
    "# calc standard devs\n",
    "modcov = np.diag(covm)\n",
    "frac_var = []\n",
    "frac_stds = []\n",
    "maxslip = []\n",
    "maxslipInd = []\n",
    "covfracs = []\n",
    "mode = genSlice(len(mest),3)\n",
    "ifrac = fracSlice(int(mode[0][0]),n_els_per_dim)\n",
    "for i in range(len(mode)):\n",
    "    for k in range(nplanes):\n",
    "        covfracs.append(modcov[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]])\n",
    "        frac_stds.append(np.std(mest[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]]))\n",
    "        frac_var.append(np.var(mest[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]]))\n",
    "        maxslip.append(np.max(np.abs(mest[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]])))\n",
    "        maxslipInd.append(np.argmax(np.abs(mest[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]])))\n",
    "        print(f'maxslip is {maxslip[k]}')\n",
    "        print(f'maxslipInd is {maxslipInd[k]}')\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CI for highest slip on each fracture for each mode\n",
    "\n",
    "maxslip_CI = []\n",
    "maxslip_CIstd = []\n",
    "# mode = genSlice(len(mest),3)\n",
    "# ifrac = fracSlice(int(mode[0][0]),mode[0][1],nplanes)\n",
    "for i in range(len(maxslip)):\n",
    "    maxslip_CI.append((maxslip[i] - covfracs[i][maxslipInd[i]], maxslip[i] + covfracs[i][maxslipInd[i]]))\n",
    "    maxslip_CIstd.append((maxslip[i] - frac_stds[i], maxslip[i] + frac_stds[i]))\n",
    "print(maxslip_CI)\n",
    "print(maxslip_CIstd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# CIs for integrated slip on each fracture plane\n",
    "\n",
    "# mode = genSlice(len(mest),3)\n",
    "# ifrac = fracSlice(int(mode[0][0]),mode[0][1],nplanes)\n",
    "# avg_open = []\n",
    "# for i in range(len(mode)):\n",
    "#     for k in range(nplanes):\n",
    "#         sumvar = np.sum(modcov[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]])\n",
    "#         avg_open.append(sumvar/num_tri_frac)\n",
    "fracind = list(sliced((0,1,2,3,4,5,0,1,2,3,4,5,0,1,2,3,4,5),1))\n",
    "num_tri_frac = [2*(els**2) for els in n_els_per_dim]\n",
    "avg_open = []\n",
    "for i in range(len(maxslip_CI)):\n",
    "    sumvar = maxslip_CI[i][0] + maxslip_CI[i][1]\n",
    "    avg_open.append(sumvar/num_tri_frac[fracind[i][0]])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Average opening on each fault plane for each mode in microns\n",
    "[print(avg_open[i]/1e3) for i in range(len(avg_open))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_dat = mest\n",
    "mode = genSlice(len(mest),3)\n",
    "ifrac = fracSlice(int(mode[0][0]),n_els_per_dim)\n",
    "for i in range(len(mode)):\n",
    "        fig = go.Figure(data=[go.Scatter3d(x=meq.x,y=meq.y,z=meq.z,mode='markers'\n",
    "                                        ,marker=dict(size=3,color=meq.datenums,colorscale='RdBu',\n",
    "                                                ))])\n",
    "        fig.add_scatter3d(x=swells[swells[\"HoleID\"]=='E1-I'].x, # type: ignore\n",
    "                           y=swells[swells[\"HoleID\"]=='E1-I'].y, # type: ignore\n",
    "                           z=swells[swells[\"HoleID\"]=='E1-I'].z, # type: ignore\n",
    "                           mode='lines',\n",
    "        marker=dict(\n",
    "                size=3,\n",
    "                color='black',\n",
    "                opacity=0.8))\n",
    "        fig.add_scatter3d(x=swells[swells[\"HoleID\"]=='E1-P'].x, # type: ignore\n",
    "                          y=swells[swells[\"HoleID\"]=='E1-P'].y, # type: ignore\n",
    "                          z=swells[swells[\"HoleID\"]=='E1-P'].z, # type: ignore\n",
    "                          mode='lines',\n",
    "        marker=dict(\n",
    "                size=3,\n",
    "                color='black',\n",
    "                opacity=0.8))\n",
    "        fig.add_scatter3d(x=pdtpts[:,0], y=pdtpts[:,1],z=pdtpts[:,2],mode='lines',\n",
    "        marker=dict(\n",
    "                size=3,\n",
    "                color='red',\n",
    "                opacity=0.8))\n",
    "        fig.add_scatter3d(x=pdbpts[:,0], y=pdbpts[:,1],z=pdbpts[:,2],mode='lines',\n",
    "        marker=dict(\n",
    "                size=3,\n",
    "                color='red',\n",
    "                opacity=0.8))\n",
    "        fig.add_scatter3d(x=pstpts[:,0], y=pstpts[:,1],z=pstpts[:,2],mode='lines',\n",
    "        marker=dict(\n",
    "                size=3,\n",
    "                color='blue',\n",
    "                opacity=0.8))\n",
    "        fig.add_scatter3d(x=psbpts[:,0], y=psbpts[:,1],z=psbpts[:,2],mode='lines',\n",
    "        marker=dict(\n",
    "                size=3,\n",
    "                color='blue',\n",
    "                opacity=0.8))\n",
    "        fig.add_scatter3d(x=obpts[:,0], y=obpts[:,1],z=obpts[:,2],mode='lines',\n",
    "        marker=dict(\n",
    "                size=3,\n",
    "                color='green',\n",
    "                opacity=0.8))\n",
    "        fig.add_scatter3d(x=otpts[:,0], y=otpts[:,1],z=otpts[:,2],mode='lines',\n",
    "        marker=dict(\n",
    "                size=3,\n",
    "                color='green',\n",
    "                opacity=0.8))\n",
    "        for k in range(nplanes):\n",
    "                # normV = Normalize(vmin=0,vmax=10000)\n",
    "                fig.add_trace(go.Mesh3d(x=allsurf_pts[k][:,0], y=allsurf_pts[k][:,1], z=allsurf_pts[k][:,2],\n",
    "                i=np.array(allsurf_tris[k])[:,0],j=np.array(allsurf_tris[k])[:,1],k=np.array(allsurf_tris[k])[:,2],\n",
    "                colorscale='thermal',facecolor=plot_dat[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],\n",
    "                intensity=plot_dat[mode[i][0]:mode[i][1]],\n",
    "                intensitymode='cell'))\n",
    "        # for k in range(nplanes):\n",
    "        #         # normV = Normalize(vmin=0,vmax=10000)\n",
    "        #         fig.add_trace(go.Mesh3d(x=allsurf_pts[k][:,0], y=allsurf_pts[k][:,1], z=allsurf_pts[k][:,2],\n",
    "        #         i=np.array(allsurf_tris[k])[:,0],j=np.array(allsurf_tris[k])[:,1],k=np.array(allsurf_tris[k])[:,2],\n",
    "        #         colorscale='thermal',facecolor=plot_dat[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],\n",
    "        #         intensity=plot_dat[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],\n",
    "        #         intensitymode='cell',\n",
    "        #         cmin=np.min(np.array([plot_dat[mfrac[k][0][0]:mfrac[k][0][1]]] + [plot_dat[mfrac[k][1][0]:mfrac[k][1][1]]] + [plot_dat[mfrac[k][2][0]:mfrac[k][2][1]]])),\n",
    "        #         cmax=np.max(np.array([plot_dat[mfrac[k][0][0]:mfrac[k][0][1]]] + [plot_dat[mfrac[k][1][0]:mfrac[k][1][1]]] + [plot_dat[mfrac[k][2][0]:mfrac[k][2][1]]]))))\n",
    "        # for k in range(nplanes):\n",
    "        #         # normV = Normalize(vmin=0,vmax=10000)\n",
    "        #         fig.add_trace(go.Mesh3d(x=allsurf_pts[k][:,0], y=allsurf_pts[k][:,1], z=allsurf_pts[k][:,2],\n",
    "        #         i=np.array(allsurf_tris[k])[:,0],j=np.array(allsurf_tris[k])[:,1],k=np.array(allsurf_tris[k])[:,2],\n",
    "        #         colorscale='thermal',facecolor=plot_dat[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],\n",
    "        #         intensity=plot_dat[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],\n",
    "        #         intensitymode='cell',\n",
    "        #         cmin=-np.percentile(plot_dat[mode[i][0]:mode[i][1]],99),\n",
    "        #         cmax=np.percentile(plot_dat[mode[i][0]:mode[i][1]],99)))\n",
    "        fig.update_layout(\n",
    "                autosize=False,\n",
    "                width=1200,\n",
    "                height=800,)\n",
    "        fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "minslip = []\n",
    "maxslip = []\n",
    "[minslip.append(min(allslip[i])) for i in range(len(allslip))]\n",
    "[maxslip.append(max(allslip[i])) for i in range(len(allslip))]\n",
    "\n",
    "\n",
    "fig = go.Figure(data=[go.Scatter3d(x=meq.x,y=meq.y,z=meq.z,mode='markers'\n",
    "                                ,marker=dict(size=3,color=meq.datenums,colorscale='RdBu',\n",
    "                                        ))])\n",
    "fig.add_scatter3d(x=swells[swells[\"HoleID\"]=='E1-I'].x, # type: ignore\n",
    "                    y=swells[swells[\"HoleID\"]=='E1-I'].y, # type: ignore\n",
    "                    z=swells[swells[\"HoleID\"]=='E1-I'].z, # type: ignore\n",
    "                    mode='lines',\n",
    "marker=dict(\n",
    "        size=3,\n",
    "        color='black',\n",
    "        opacity=0.8))\n",
    "fig.add_scatter3d(x=swells[swells[\"HoleID\"]=='E1-P'].x, # type: ignore\n",
    "                    y=swells[swells[\"HoleID\"]=='E1-P'].y, # type: ignore\n",
    "                    z=swells[swells[\"HoleID\"]=='E1-P'].z, # type: ignore\n",
    "                    mode='lines',\n",
    "marker=dict(\n",
    "        size=3,\n",
    "        color='black',\n",
    "        opacity=0.8))\n",
    "fig.add_scatter3d(x=pdtpts[:,0], y=pdtpts[:,1],z=pdtpts[:,2],mode='lines',\n",
    "marker=dict(\n",
    "        size=3,\n",
    "        color='red',\n",
    "        opacity=0.8))\n",
    "fig.add_scatter3d(x=pdbpts[:,0], y=pdbpts[:,1],z=pdbpts[:,2],mode='lines',\n",
    "marker=dict(\n",
    "        size=3,\n",
    "        color='red',\n",
    "        opacity=0.8))\n",
    "fig.add_scatter3d(x=pstpts[:,0], y=pstpts[:,1],z=pstpts[:,2],mode='lines',\n",
    "marker=dict(\n",
    "        size=3,\n",
    "        color='blue',\n",
    "        opacity=0.8))\n",
    "fig.add_scatter3d(x=psbpts[:,0], y=psbpts[:,1],z=psbpts[:,2],mode='lines',\n",
    "marker=dict(\n",
    "        size=3,\n",
    "        color='blue',\n",
    "        opacity=0.8))\n",
    "fig.add_scatter3d(x=obpts[:,0], y=obpts[:,1],z=obpts[:,2],mode='lines',\n",
    "marker=dict(\n",
    "        size=3,\n",
    "        color='green',\n",
    "        opacity=0.8))\n",
    "fig.add_scatter3d(x=otpts[:,0], y=otpts[:,1],z=otpts[:,2],mode='lines',\n",
    "marker=dict(\n",
    "        size=3,\n",
    "        color='green',\n",
    "        opacity=0.8))\n",
    "\n",
    "# for k in range(nplanes):\n",
    "#         # normV = Normalize(vmin=0,vmax=10000)\n",
    "#         fig.add_trace(go.Mesh3d(x=allsurf_pts[k][:,0], y=allsurf_pts[k][:,1], z=allsurf_pts[k][:,2],\n",
    "#         i=allsurf_tri_arr[k,:,0],j=allsurf_tri_arr[k,:,1],k=allsurf_tri_arr[k,:,2],\n",
    "#         colorscale='thermal',facecolor=mest[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],\n",
    "#         intensity=mest[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],\n",
    "#         intensitymode='cell',\n",
    "#         cmin=np.min(np.array([mest[mfrac[k][0][0]:mfrac[k][0][1]]] + [mest[mfrac[k][1][0]:mfrac[k][1][1]]] + [mest[mfrac[k][2][0]:mfrac[k][2][1]]])),\n",
    "#         cmax=np.max(np.array([mest[mfrac[k][0][0]:mfrac[k][0][1]]] + [mest[mfrac[k][1][0]:mfrac[k][1][1]]] + [mest[mfrac[k][2][0]:mfrac[k][2][1]]]))))\n",
    "# for k in range(nplanes):\n",
    "#         # normV = Normalize(vmin=0,vmax=10000)\n",
    "#         fig.add_trace(go.Mesh3d(x=allsurf_pts[k][:,0], y=allsurf_pts[k][:,1], z=allsurf_pts[k][:,2],\n",
    "#         i=allsurf_tri_arr[k,:,0],j=allsurf_tri_arr[k,:,1],k=allsurf_tri_arr[k,:,2],\n",
    "#         colorscale='thermal',facecolor=mest[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],\n",
    "#         intensity=mest[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],\n",
    "#         intensitymode='cell',\n",
    "#         cmin=min(mest[mode[i][0]:mode[i][1]]),\n",
    "#         cmax=max(mest[mode[i][0]:mode[i][1]])))\n",
    "        # cmin=-np.percentile(mest[mode[i][0]:mode[i][1]],99),\n",
    "        # cmax=np.percentile(mest[mode[i][0]:mode[i][1]],99)))\n",
    "for k in range(nplanes):\n",
    "    # normV = Normalize(vmin=0,vmax=10000)\n",
    "    fig.add_trace(go.Mesh3d(x=allsurf_pts[k][:,0], y=allsurf_pts[k][:,1], z=allsurf_pts[k][:,2],\n",
    "                i=np.array(allsurf_tris[k])[:,0],j=np.array(allsurf_tris[k])[:,1],k=np.array(allsurf_tris[k])[:,2],\n",
    "                colorscale='thermal',facecolor=allslip[k],intensity=allslip[k],intensitymode='cell',\n",
    "                cmin=-np.min(np.percentile(np.array(allslip),99)),\n",
    "                cmax=np.max(np.percentile(np.array(allslip),99))))\n",
    "# for k in range(nplanes):\n",
    "#     # normV = Normalize(vmin=0,vmax=10000)\n",
    "#     fig.add_trace(go.Mesh3d(x=allsurf_pts[k][:,0], y=allsurf_pts[k][:,1], z=allsurf_pts[k][:,2],\n",
    "#                 i=allsurf_tri_arr[k,:,0],j=allsurf_tri_arr[k,:,1],k=allsurf_tri_arr[k,:,2],\n",
    "#                 colorscale='YlOrRd',facecolor=covslip[k],intensity=covslip[k],intensitymode='cell',cmin=1000,cmax=100000))\n",
    "# for k in range(nplanes):\n",
    "#     # normV = Normalize(vmin=0,vmax=10000)\n",
    "#     fig.add_trace(go.Mesh3d(x=allsurf_pts[k][:,0], y=allsurf_pts[k][:,1], z=allsurf_pts[k][:,2],\n",
    "#                 i=allsurf_tri_arr[k,:,0],j=allsurf_tri_arr[k,:,1],k=allsurf_tri_arr[k,:,2],\n",
    "#                 colorscale='YlOrRd',facecolor=np.array(modcov[k]),intensity=np.array(modcov[k]),intensitymode='cell',cauto=True))\n",
    "fig.update_layout(\n",
    "        autosize=False,\n",
    "        width=1200,\n",
    "        height=800,)\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for j in range(len(mest)):\n",
    "    nplanes = 6\n",
    "    fracList = [f'OTP','West North','West South','West Deep','East North','East South']\n",
    "    titles=[f'Strike Slip ','Dip Slip','Tensile ']\n",
    "    fig,ax = plt.subplots(nplanes,3,figsize=(12,10))\n",
    "    mode = genSlice(len(mest[j]),3)\n",
    "    ifrac = fracSlice(mode[0][1],nplanes)\n",
    "    pfrac = fracSlice(len(allsurf_pts),nplanes)\n",
    "    for i in range(3):\n",
    "        for k in range(nplanes):\n",
    "        # mode = list(sliced((0,8,8,16,16,24),2))\n",
    "\n",
    "            # vm = np.percentile(mest[j][mode[i][0]:mode[i][1]],99)\n",
    "            vm = np.percentile(mest[j][mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],99)\n",
    "\n",
    "\n",
    "            # im  = ax[k,i].tripcolor(allsurf_pts[pfrac[k][0]:pfrac[k][1], 0], allsurf_pts[pfrac[k][0]:pfrac[k][1], 2],\n",
    "            #                         allsurf_tris[ifrac[k][0]:ifrac[k][1]],\n",
    "            #                         facecolors=mest[mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],\n",
    "            #                         cmap='seismic',vmin=-50000,vmax=50000) # type: ignore\n",
    "            \n",
    "            im  = ax[k,i].tripcolor(np.array(allsurf_pts[k])[:,0], np.array(allsurf_pts[k])[:,2],\n",
    "                            np.array(allsurf_tris[k]),\n",
    "                            facecolors=mest[j][mode[i][0]:mode[i][1]][ifrac[k][0]:ifrac[k][1]],\n",
    "                            cmap='seismic',vmin=-vm,vmax=vm) # type: ignore\n",
    "            ax[k,i].set_xlabel('X direction (m)')\n",
    "            ax[k,i].set_ylabel('Z direction (m)')\n",
    "            # ax[i].axis('equal')\n",
    "            # ax[i].set_aspect('equal','box')\n",
    "            plt.colorbar(im,ax=ax[k,i])\n",
    "            ax[k,i].set_title(f'{fracList[k]} {titles[i]}')\n",
    "            plt.suptitle(f'{DASdates[timeEpochs[j]]}')\n",
    "            fig.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def lineSearch(x, g, p, func, t=0.5):\n",
    "\t\"\"\"\n",
    "\tLine Search Function\n",
    "\tsimple descent line search\n",
    "\n",
    "\tinput\n",
    "\t-----\n",
    "\tx : array_like\n",
    "\t\tBase point.\n",
    "\tg : array_like\n",
    "\t\tGradient on the base point.\n",
    "\tp : array_like\n",
    "\t\tGiven descent direction.\n",
    "\tfunc : function\n",
    "\t\tInput x and return funciton value.\n",
    "\tt : float, optional\n",
    "\t\tstep size shrink ratio\n",
    "\n",
    "\toutput\n",
    "\t------\n",
    "\tstep_size : float or None\n",
    "\t\tWhen sucess return the step_size, otherwise return None.\n",
    "\t\"\"\"\n",
    "\t# initial step size\n",
    "\tstep_size = 1.0\n",
    "\t# \n",
    "\tm = g.dot(p)\n",
    "\tif m < 0:\n",
    "\t\tprint('Line search: not a descent direction.')\n",
    "\t\treturn None\n",
    "\t#\n",
    "\tf = func(x)\n",
    "\ty = x - step_size*p\n",
    "\t#\n",
    "\twhile func(y) > f:\n",
    "\t\tstep_size *= t\n",
    "\t\tif step_size < 1e-15:\n",
    "\t\t\tprint('Line search: step size too small.')\n",
    "\t\t\treturn None\n",
    "\t\ty = x - step_size*p\n",
    "\t#\n",
    "\treturn step_size\n",
    "\n",
    "# line search function - armijo line search\n",
    "# -----------------------------------------------------------------------------\n",
    "def lineSearch_armijo(x, g, p, func, c=0.01, t=0.5):\n",
    "\t\"\"\"\n",
    "\tLine Search Function\n",
    "\tarmijo line search\n",
    "\n",
    "\tinput\n",
    "\t-----\n",
    "\tx : array_like\n",
    "\t\tBase point.\n",
    "\tg : array_like\n",
    "\t\tGradient on the base point.\n",
    "\tp : array_like\n",
    "\t\tGiven descent direction.\n",
    "\tfunc : function\n",
    "\t\tInput x and return funciton value.\n",
    "\tc : float, optional\n",
    "\t\thas to strictly be between 0 and 1\n",
    "\tt : float, optional\n",
    "\t\tstep size shrink ratio\n",
    "\n",
    "\toutput\n",
    "\t------\n",
    "\tstep_size : float or None\n",
    "\t\tWhen sucess return the step_size, otherwise return None.\n",
    "\t\"\"\"\n",
    "\t# save guard for c\n",
    "\tassert 0 < c < 1, 'c needs to strictly be in 0 and 1'\n",
    "\t# initial step size\n",
    "\tstep_size = 1.0\n",
    "\t# \n",
    "\tm = c*g.dot(p)\n",
    "\tif m < 0:\n",
    "\t\tprint('Line search: not a descent direction.')\n",
    "\t\treturn None\n",
    "\t#\n",
    "\tf = func(x)\n",
    "\ty = x - step_size*p\n",
    "\t#\n",
    "\twhile func(y) > f - step_size*m:\n",
    "\t\tstep_size *= t\n",
    "\t\tif step_size < 1e-15:\n",
    "\t\t\tprint('Line search: step size too small.')\n",
    "\t\t\treturn None\n",
    "\t\ty = x - step_size*p\n",
    "\t#\n",
    "\treturn step_size\n",
    "\n",
    "\n",
    "def optimizeWithGD(x0, func, grad, step_size, tol=1e-7, max_iter=10000, use_line_search=False):\n",
    "    \"\"\"\n",
    "    Optimize with Gradient Descent\n",
    "    \n",
    "    input\n",
    "    -----\n",
    "    x0 : array_like\n",
    "        Starting point for the solver.\n",
    "    func : function\n",
    "        Takes x and return the obj function value.\n",
    "    grad : function\n",
    "        Takes x and returns the gradient of \"func\".\n",
    "    step_size : float or None\n",
    "        If it is a float number and `use_line_search=False`, it will be used as the step size.\n",
    "        Otherwise, line search will be used\n",
    "    tol : float, optional\n",
    "        Gradient tolerance for terminating the solver.\n",
    "    max_iter : int, optional\n",
    "        Maximum number of iterations for terminating the solver.\n",
    "    use_line_search : bool, optional\n",
    "        When it is true a line search will be used, other wise `step_size` has to be provided.\n",
    "        \n",
    "    output\n",
    "    ------\n",
    "    x : array_like\n",
    "        Final solution\n",
    "    obj_his : array_like\n",
    "        Objective function's values convergence history\n",
    "    err_his : array_like\n",
    "        Norm of gradient convergence history\n",
    "    exit_flag : int\n",
    "        0, norm of gradient below `tol`\n",
    "        1, exceed maximum number of iteration\n",
    "        2, line search fail\n",
    "        3, other\n",
    "    \"\"\"\n",
    "    # safeguard\n",
    "    if not use_line_search and step_size is None:\n",
    "        print('Please specify the step_size or use the line search.')\n",
    "        return x0, np.array([]), np.array([]), 3\n",
    "    \n",
    "    # initial step\n",
    "    x = np.copy(x0)\n",
    "    g = grad(x)\n",
    "    #\n",
    "    obj = func(x)\n",
    "    err = norm(g)\n",
    "    #\n",
    "    obj_his = np.zeros(max_iter + 1)\n",
    "    err_his = np.zeros(max_iter + 1)\n",
    "    #\n",
    "    obj_his[0] = obj\n",
    "    err_his[0] = err\n",
    "    \n",
    "    # start iterations\n",
    "    iter_count = 0\n",
    "    while err >= tol:\n",
    "        if use_line_search:\n",
    "            step_size = lineSearch(x, g, g, func)\n",
    "        #\n",
    "        # if line search fail step_size will be None\n",
    "        if step_size is None:\n",
    "            print('Gradient descent line search fail.')\n",
    "            return x, obj_his[:iter_count+1], err_his[:iter_count+1], 2\n",
    "        #\n",
    "        # gradient descent step\n",
    "        x -= step_size * g\n",
    "        # update obj function and gradient and error\n",
    "        g = grad(x)\n",
    "        #\n",
    "        obj = func(x)\n",
    "        err = norm(g)\n",
    "        #\n",
    "        iter_count += 1\n",
    "        obj_his[iter_count] = obj\n",
    "        err_his[iter_count] = err\n",
    "        #\n",
    "        # check if exceed maximum number of iteration\n",
    "        if iter_count >= max_iter:\n",
    "            print('Gradient descent reached maximum number of iteration.')\n",
    "            return x, obj_his[:iter_count+1], err_his[:iter_count+1], 1\n",
    "    #\n",
    "    return x, obj_his[:iter_count+1], err_his[:iter_count+1], 0\n",
    "\n",
    "def optimizeWithPGD(x0, func_f, func_g, grad_f, prox_g, beta_f, tol=1e-5, max_iter=50000):\n",
    "    \"\"\"\n",
    "    Optimize with Proximal Gradient Descent Method\n",
    "        min_x f(x) + g(x)\n",
    "    where f is beta smooth and g is proxiable.\n",
    "    \n",
    "    input\n",
    "    -----\n",
    "    x0 : array_like\n",
    "        Starting point for the solver\n",
    "    func_f : function\n",
    "        Input x and return the function value of f\n",
    "    func_g : function\n",
    "        Input x and return the function value of g\n",
    "    grad_f : function\n",
    "        Input x and return the gradient of f\n",
    "    prox_g : function\n",
    "        Input x and a constant float number and return the prox solution\n",
    "    beta_f : float\n",
    "        beta smoothness constant for f\n",
    "    tol : float, optional\n",
    "        Gradient tolerance for terminating the solver.\n",
    "    max_iter : int, optional\n",
    "        Maximum number of iteration for terminating the solver.\n",
    "        \n",
    "    output\n",
    "    ------\n",
    "    x : array_like\n",
    "        Final solution\n",
    "    obj_his : array_like\n",
    "        Objective function value convergence history\n",
    "    err_his : array_like\n",
    "        Norm of gradient convergence history\n",
    "    exit_flag : int\n",
    "        0, norm of gradient below `tol`\n",
    "        1, exceed maximum number of iteration\n",
    "        2, others\n",
    "    \"\"\"\n",
    "    # initial information\n",
    "    x = x0.copy()\n",
    "    g = grad_f(x)\n",
    "    #\n",
    "    step_size = 1.0/beta_f\n",
    "    # not recording the initial point since we do not have measure of the optimality\n",
    "    obj_his = np.zeros(max_iter)\n",
    "    err_his = np.zeros(max_iter)\n",
    "    \n",
    "    # start iteration\n",
    "    iter_count = 0\n",
    "    err = tol + 1.0\n",
    "    while err >= tol:\n",
    "        #####\n",
    "        # the proximal gradient step\n",
    "        x_new = prox_g(x - step_size*g,step_size)\n",
    "        #####\n",
    "        # update information\n",
    "        obj = func_f(x_new) + func_g(x_new)\n",
    "        err = norm(x - x_new)/step_size\n",
    "        #\n",
    "        np.copyto(x, x_new)\n",
    "        g = grad_f(x)\n",
    "        #\n",
    "        obj_his[iter_count] = obj\n",
    "        err_his[iter_count] = err\n",
    "        #\n",
    "        # check if exceed maximum number of iteration\n",
    "        iter_count += 1\n",
    "        if iter_count >= max_iter:\n",
    "            print('Proximal gradient descent reach maximum of iteration')\n",
    "            return x, obj_his[:iter_count], err_his[:iter_count], 1\n",
    "    #\n",
    "    return x, obj_his[:iter_count], err_his[:iter_count], 0\n",
    "# Newton's Method\n",
    "# -----------------------------------------------------------------------------\n",
    "def optimizeWithNT(x0, func, grad, hess, tol=1e-6, max_iter=100):\n",
    "    \"\"\"\n",
    "    Optimize with Newton's Method\n",
    "    \n",
    "    input\n",
    "    -----\n",
    "    x0 : array_like\n",
    "        Starting point for the solver.\n",
    "    func : function\n",
    "        Input x and return the function value.\n",
    "    grad : function\n",
    "        Input x and return the gradient.\n",
    "    hess : function\n",
    "        Input x and return the Hessian matrix.\n",
    "    tol : float, optional\n",
    "        Gradient tolerance for terminating the solver.\n",
    "    max_iter : int, optional\n",
    "        Maximum number of iteration for terminating the solver.\n",
    "        \n",
    "    output\n",
    "    ------\n",
    "    x : array_like\n",
    "        Final solution\n",
    "    obj_his : array_like\n",
    "        Objective function value convergence history\n",
    "    err_his : array_like\n",
    "        Norm of gradient convergence history\n",
    "    exit_flag : int\n",
    "        0, norm of gradient below `tol`\n",
    "        1, exceed maximum number of iteration\n",
    "        2, others\n",
    "    \"\"\"\n",
    "    # initial step\n",
    "    x = np.copy(x0)\n",
    "    g = grad(x)\n",
    "    H = hess(x)\n",
    "    #\n",
    "    obj = func(x)\n",
    "    err = norm(g)\n",
    "    #\n",
    "    obj_his = np.zeros(max_iter + 1)\n",
    "    err_his = np.zeros(max_iter + 1)\n",
    "    #\n",
    "    obj_his[0] = obj\n",
    "    err_his[0] = err\n",
    "    \n",
    "    # start iteration\n",
    "    iter_count = 0\n",
    "    while err >= tol:\n",
    "        # Newton's step\n",
    "        x -= solve(H, g)\n",
    "        #\n",
    "        # update function, gradient and Hessian\n",
    "        g = grad(x)\n",
    "        H = hess(x)\n",
    "        #\n",
    "        obj = func(x)\n",
    "        err = norm(g)\n",
    "        #\n",
    "        iter_count += 1\n",
    "        obj_his[iter_count] = obj\n",
    "        err_his[iter_count] = err\n",
    "        #\n",
    "        # check if exceed maximum number of iteration\n",
    "        if iter_count >= max_iter:\n",
    "            print('Gradient descent reach maximum number of iteration.')\n",
    "            return x, obj_his[:iter_count+1], err_his[:iter_count+1], 1\n",
    "    #\n",
    "    return x, obj_his[:iter_count+1], err_his[:iter_count+1], 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the function, prox and the beta constant\n",
    "def func_f_cs(x):\n",
    "    f_cs = (0.5)*norm(A@x - b,2)**2\n",
    "    return f_cs\n",
    "def func_g_cs(x):\n",
    "    g_cs = lam_cs*norm(x,2)\n",
    "    return g_cs\n",
    "def grad_f_cs(x):\n",
    "    Ax = A@x\n",
    "    grad_f = A.T@(Ax - b) \n",
    "    return grad_f\n",
    "def prox_g_cs(x, t):\n",
    "    # prox of 1 norm\n",
    "    prox_gx = np.sign(x)*np.maximum(np.abs(x) - lam_cs*t, 0)\n",
    "    # prox of 2 norm\n",
    "    # prox_gx = np.maximum(0,(1 - (t / norm(x,2)))) * x\n",
    "    # if norm(x,2) <= t:\n",
    "    #     prox_gx = 0\n",
    "    # elif norm(x,2) > t:\n",
    "    #     prox_gx = ((norm(x,2) - t) / norm(x,2)) * x \n",
    "        \n",
    "    return prox_gx\n",
    "def prox_g_ps(x,t):\n",
    "    # prox of 2 norm\n",
    "    prox_gx = np.maximum(0,(1 - (t / norm(x,2)))) * x\n",
    "    return prox_gx\n",
    "beta_f_cs = norm(A.T, 2)**2\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "obspy",
   "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.9.13"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}