From 7ea1c25910dd684490d8d3191fcb1836713ca314 Mon Sep 17 00:00:00 2001 From: Brendan Meade Date: Sat, 13 Jul 2024 16:14:19 -0400 Subject: [PATCH] Starting to remove eigenvectors_two_components and build index_eigen --- .../celeri_western_north_america_qp_05.ipynb | 951 ++++++++++++++++++ 1 file changed, 951 insertions(+) create mode 100644 notebooks/celeri_western_north_america_qp_05.ipynb diff --git a/notebooks/celeri_western_north_america_qp_05.ipynb b/notebooks/celeri_western_north_america_qp_05.ipynb new file mode 100644 index 0000000..a737473 --- /dev/null +++ b/notebooks/celeri_western_north_america_qp_05.ipynb @@ -0,0 +1,951 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "%config InlineBackend.figure_format = \"retina\"" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import addict\n", + "import scipy\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "import celeri\n", + "\n", + "plt.rcParams[\"text.usetex\"] = False" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001b[32m2024-07-13 16:12:21.717\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mget_logger\u001b[0m:\u001b[36m6343\u001b[0m - \u001b[1mRead: ../data/command/western_north_america_command.json\u001b[0m\n", + "\u001b[32m2024-07-13 16:12:21.718\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mget_logger\u001b[0m:\u001b[36m6344\u001b[0m - \u001b[1mRUN_NAME: 0000000041\u001b[0m\n", + "\u001b[32m2024-07-13 16:12:21.718\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mget_logger\u001b[0m:\u001b[36m6345\u001b[0m - \u001b[1mWrite log file: ../runs/0000000041/0000000041.log\u001b[0m\n", + "\u001b[32m2024-07-13 16:12:21.718\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mread_data\u001b[0m:\u001b[36m302\u001b[0m - \u001b[1mReading data files\u001b[0m\n", + "\u001b[32m2024-07-13 16:12:21.723\u001b[0m | \u001b[32m\u001b[1mSUCCESS \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mread_data\u001b[0m:\u001b[36m306\u001b[0m - \u001b[32m\u001b[1mRead: ../data/segment/western_north_america_segment.csv\u001b[0m\n", + "\u001b[32m2024-07-13 16:12:21.725\u001b[0m | \u001b[32m\u001b[1mSUCCESS \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mread_data\u001b[0m:\u001b[36m311\u001b[0m - \u001b[32m\u001b[1mRead: ../data/block/western_north_america_block.csv\u001b[0m\n", + "\u001b[32m2024-07-13 16:12:21.726\u001b[0m | \u001b[32m\u001b[1mSUCCESS \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mread_data\u001b[0m:\u001b[36m318\u001b[0m - \u001b[32m\u001b[1mRead: ../data/mesh/western_north_america_mesh_parameters.json\u001b[0m\n", + "\n", + "\u001b[32m2024-07-13 16:12:21.836\u001b[0m | \u001b[32m\u001b[1mSUCCESS \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mread_data\u001b[0m:\u001b[36m464\u001b[0m - \u001b[32m\u001b[1mRead: ../data/mesh/cascadia.msh\u001b[0m\n", + "\u001b[32m2024-07-13 16:12:21.839\u001b[0m | \u001b[32m\u001b[1mSUCCESS \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mread_data\u001b[0m:\u001b[36m500\u001b[0m - \u001b[32m\u001b[1mRead: ../data/station/western_north_america_station.csv\u001b[0m\n", + "\u001b[32m2024-07-13 16:12:21.840\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mread_data\u001b[0m:\u001b[36m515\u001b[0m - \u001b[1mNo mogi_file_name\u001b[0m\n", + "\u001b[32m2024-07-13 16:12:21.841\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mread_data\u001b[0m:\u001b[36m537\u001b[0m - \u001b[1mNo sar_file_name\u001b[0m\n" + ] + } + ], + "source": [ + "# TODO: CHANGE THIS to 1/ 1e3\n", + "\n", + "DEG_PER_MYR_TO_RAD_PER_YR = 1 / 1e6\n", + "COMMAND_FILE_NAME = \"../data/command/western_north_america_command.json\"\n", + "command = celeri.get_command(COMMAND_FILE_NAME)\n", + "celeri.create_output_folder(command)\n", + "logger = celeri.get_logger(command)\n", + "segment, block, meshes, station, mogi, sar = celeri.read_data(command)\n", + "\n", + "# NOTE: This modification of the number of eigenvalues is for experiments\n", + "meshes[0].n_modes_strike_slip = 10\n", + "meshes[0].n_modes_dip_slip = 50\n", + "meshes[0].n_modes_total = meshes[0].n_modes_strike_slip + meshes[0].n_modes_dip_slip\n", + "\n", + "# Set n_modes to the great of strike-slip or dip slip modes\n", + "# This overrides meshes[i].n_modes in mesh files\n", + "meshes[0].n_modes = np.max([meshes[0].n_modes_strike_slip, meshes[0].n_modes_dip_slip])\n", + "\n", + "station = celeri.process_station(station, command)\n", + "segment = celeri.process_segment(segment, command, meshes)\n", + "sar = celeri.process_sar(sar, command)\n", + "closure, block = celeri.assign_block_labels(segment, station, block, mogi, sar)\n", + "assembly = addict.Dict()\n", + "operators = addict.Dict()\n", + "operators.meshes = [addict.Dict()] * len(meshes)\n", + "assembly = celeri.merge_geodetic_data(assembly, station, sar)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Utilities" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "def test_agreement(estimation):\n", + " # np.savez(\n", + " # \"test_eigen_arrays.npz\",\n", + " # estimation_eigen_cvxopt_bounded_slip_rates=estimation_eigen_cvxopt_bounded.slip_rates,\n", + " # estimation_eigen_cvxopt_bounded_tde_rates=estimation_eigen_cvxopt_bounded.tde_rates,\n", + " # estimation_eigen_cvxopt_bounded_east_vel_residual=estimation_eigen_cvxopt_bounded.east_vel_residual,\n", + " # estimation_eigen_cvxopt_bounded_north_vel_residual=estimation_eigen_cvxopt_bounded.north_vel_residual,\n", + " # )\n", + "\n", + " testing_eigen_arrays = np.load(\"test_eigen_arrays.npz\")\n", + " print(\n", + " np.allclose(\n", + " estimation.slip_rates,\n", + " testing_eigen_arrays[\"estimation_eigen_cvxopt_bounded_slip_rates\"],\n", + " )\n", + " )\n", + "\n", + " print(\n", + " np.allclose(\n", + " estimation.tde_rates,\n", + " testing_eigen_arrays[\"estimation_eigen_cvxopt_bounded_tde_rates\"],\n", + " )\n", + " )\n", + " print(\n", + " np.allclose(\n", + " estimation.east_vel_residual,\n", + " testing_eigen_arrays[\"estimation_eigen_cvxopt_bounded_east_vel_residual\"],\n", + " )\n", + " )\n", + " print(\n", + " np.allclose(\n", + " estimation.north_vel_residual,\n", + " testing_eigen_arrays[\"estimation_eigen_cvxopt_bounded_north_vel_residual\"],\n", + " )\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001b[32m2024-07-13 16:12:22.574\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mget_elastic_operators\u001b[0m:\u001b[36m1687\u001b[0m - \u001b[1mUsing precomputed elastic operators\u001b[0m\n", + "\u001b[32m2024-07-13 16:12:23.148\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mget_slip_rate_constraints\u001b[0m:\u001b[36m2918\u001b[0m - \u001b[1mFound 1 slip rate constraints\u001b[0m\n", + "\u001b[32m2024-07-13 16:12:23.151\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mget_slip_rate_constraints\u001b[0m:\u001b[36m2921\u001b[0m - \u001b[1mStrike-slip rate constraint on cfm_san_andreas_mojave_extruded_trace_part1_sa: rate = -50.00 (mm/yr), 1-sigma uncertainty = +/-1.00 (mm/yr)\u001b[0m\n" + ] + } + ], + "source": [ + "# Get all elastic operators for segments and TDEs\n", + "celeri.get_elastic_operators(operators, meshes, segment, station, command)\n", + "\n", + "# Get TDE smoothing operators\n", + "celeri.get_all_mesh_smoothing_matrices(meshes, operators)\n", + "\n", + "n_blocks = len(block)\n", + "operators.rotation_to_velocities = celeri.get_rotation_to_velocities_partials(\n", + " station, n_blocks\n", + ")\n", + "operators.global_float_block_rotation = celeri.get_global_float_block_rotation_partials(\n", + " station\n", + ")\n", + "assembly, operators.block_motion_constraints = celeri.get_block_motion_constraints(\n", + " assembly, block, command\n", + ")\n", + "assembly, operators.slip_rate_constraints = celeri.get_slip_rate_constraints(\n", + " assembly, segment, block, command\n", + ")\n", + "operators.rotation_to_slip_rate = celeri.get_rotation_to_slip_rate_partials(\n", + " segment, block\n", + ")\n", + "(\n", + " operators.block_strain_rate_to_velocities,\n", + " strain_rate_block_index,\n", + ") = celeri.get_block_strain_rate_to_velocities_partials(block, station, segment)\n", + "operators.mogi_to_velocities = celeri.get_mogi_to_velocities_partials(\n", + " mogi, station, command\n", + ")\n", + "celeri.get_tde_slip_rate_constraints(meshes, operators)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Full direct dense block model solve\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "index = celeri.get_index(assembly, station, block, meshes, mogi)\n", + "estimation = addict.Dict()\n", + "estimation.data_vector = celeri.get_data_vector(assembly, index, meshes)\n", + "estimation.weighting_vector = celeri.get_weighting_vector(\n", + " command, station, meshes, index\n", + ")\n", + "estimation.operator = celeri.get_full_dense_operator(operators, meshes, index)\n", + "\n", + "# Solve the overdetermined linear system\n", + "estimation.state_covariance_matrix = np.linalg.inv(\n", + " estimation.operator.T * estimation.weighting_vector @ estimation.operator\n", + ")\n", + "estimation.state_vector = (\n", + " estimation.state_covariance_matrix\n", + " @ estimation.operator.T\n", + " * estimation.weighting_vector\n", + " @ estimation.data_vector\n", + ")\n", + "\n", + "celeri.post_process_estimation(estimation, operators, station, index)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# KL functions" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "index_eigen.n_modes_total=60\n" + ] + } + ], + "source": [ + "index_eigen = addict.Dict()\n", + "\n", + "# Get total number of modes\n", + "index_eigen.n_modes_total = 0\n", + "for i in range(len(meshes)):\n", + " index_eigen.n_modes_total += meshes[i].n_modes_strike_slip\n", + " index_eigen.n_modes_total += meshes[i].n_modes_dip_slip\n", + "\n", + "print(f\"{index_eigen.n_modes_total=}\")\n", + "\n", + "# TODO: Add strain elements\n", + "n_strain_blocks = np.sum(np.where(block.strain_rate_flag == 1)[0])\n", + "n_mogi = len(mogi)\n", + "n_state_vector_eigen = (\n", + " 3 * n_blocks + index_eigen.n_modes_total + n_mogi + n_strain_blocks\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "def get_eigenvalues_and_eigenvectors(n_eigenvalues, x, y, z, distance_exponent):\n", + " n_tde = x.size\n", + "\n", + " # Calculate Cartesian distances between triangle centroids\n", + " centroid_coordinates = np.array([x, y, z]).T\n", + " distance_matrix = scipy.spatial.distance.cdist(\n", + " centroid_coordinates, centroid_coordinates, \"euclidean\"\n", + " )\n", + "\n", + " # Rescale distance matrix to the range 0-1\n", + " distance_matrix = (distance_matrix - np.min(distance_matrix)) / np.ptp(\n", + " distance_matrix\n", + " )\n", + "\n", + " # Calculate correlation matrix\n", + " correlation_matrix = np.exp(-(distance_matrix**distance_exponent))\n", + " # correlation_matrix = np.exp(-distance_matrix)\n", + "\n", + " # https://stackoverflow.com/questions/12167654/fastest-way-to-compute-k-largest-eigenvalues-and-corresponding-eigenvectors-with\n", + " eigenvalues, eigenvectors = scipy.linalg.eigh(\n", + " correlation_matrix,\n", + " subset_by_index=[n_tde - n_eigenvalues, n_tde - 1],\n", + " )\n", + " eigenvalues = np.real(eigenvalues)\n", + " eigenvectors = np.real(eigenvectors)\n", + " ordered_index = np.flip(np.argsort(eigenvalues))\n", + " eigenvalues = eigenvalues[ordered_index]\n", + " eigenvectors = eigenvectors[:, ordered_index]\n", + " return eigenvalues, eigenvectors\n", + "\n", + "\n", + "def get_data_vector_eigen(assembly, index):\n", + " data_vector = np.zeros(\n", + " 2 * index.n_stations\n", + " + 3 * index.n_block_constraints\n", + " + index.n_slip_rate_constraints\n", + " + index.n_tde_constraints_total\n", + " )\n", + "\n", + " # Add GPS stations to data vector\n", + " data_vector[index.start_station_row : index.end_station_row] = celeri.interleave2(\n", + " assembly.data.east_vel, assembly.data.north_vel\n", + " )\n", + "\n", + " # Add block motion constraints to data vector\n", + " data_vector[index.start_block_constraints_row : index.end_block_constraints_row] = (\n", + " DEG_PER_MYR_TO_RAD_PER_YR * assembly.data.block_constraints\n", + " )\n", + "\n", + " # Add slip rate constraints to data vector\n", + " data_vector[\n", + " index.start_slip_rate_constraints_row : index.end_slip_rate_constraints_row\n", + " ] = assembly.data.slip_rate_constraints\n", + " return data_vector\n", + "\n", + "\n", + "def get_weighting_vector_eigen(command, station, meshes, index):\n", + " # Initialize and build weighting matrix\n", + " weighting_vector = np.ones(\n", + " 2 * index.n_stations\n", + " + 3 * index.n_block_constraints\n", + " + index.n_slip_rate_constraints\n", + " + index.n_tde_constraints_total\n", + " )\n", + "\n", + " weighting_vector[index.start_station_row : index.end_station_row] = (\n", + " celeri.interleave2(1 / (station.east_sig**2), 1 / (station.north_sig**2))\n", + " )\n", + "\n", + " weighting_vector[\n", + " index.start_block_constraints_row : index.end_block_constraints_row\n", + " ] = command.block_constraint_weight\n", + "\n", + " weighting_vector[\n", + " index.start_slip_rate_constraints_row : index.end_slip_rate_constraints_row\n", + " ] = command.slip_constraint_weight * np.ones(index.n_slip_rate_constraints)\n", + "\n", + " for i in range(len(meshes)):\n", + " # TODO: This is too hacky to keep. Revise with no reference to smoothing\n", + " weighting_vector[\n", + " index.start_tde_smoothing_row[i] : index.start_tde_smoothing_row[i]\n", + " + index.n_tde_constraints[i]\n", + " ] = command.tri_con_weight * np.ones(index.n_tde_constraints[i])\n", + "\n", + " return weighting_vector\n", + "\n", + "\n", + "def get_full_dense_operator_eigen(operators, meshes, index):\n", + " # Initialize linear operator\n", + " operator = np.zeros(\n", + " (\n", + " 2 * index.n_stations\n", + " + 3 * index.n_block_constraints\n", + " + index.n_slip_rate_constraints\n", + " + index.n_tde_constraints_total,\n", + " 3 * index.n_blocks + index_eigen.n_modes_total,\n", + " )\n", + " )\n", + "\n", + " # Insert block rotations and elastic velocities from fully locked segments\n", + " operators.rotation_to_slip_rate_to_okada_to_velocities = (\n", + " operators.slip_rate_to_okada_to_velocities @ operators.rotation_to_slip_rate\n", + " )\n", + " operator[\n", + " index.start_station_row : index.end_station_row,\n", + " index.start_block_col : index.end_block_col,\n", + " ] = (\n", + " operators.rotation_to_velocities[index.station_row_keep_index, :]\n", + " - operators.rotation_to_slip_rate_to_okada_to_velocities[\n", + " index.station_row_keep_index, :\n", + " ]\n", + " )\n", + "\n", + " # Insert block motion constraints\n", + " operator[\n", + " index.start_block_constraints_row : index.end_block_constraints_row,\n", + " index.start_block_col : index.end_block_col,\n", + " ] = operators.block_motion_constraints\n", + "\n", + " # Insert slip rate constraints\n", + " operator[\n", + " index.start_slip_rate_constraints_row : index.end_slip_rate_constraints_row,\n", + " index.start_block_col : index.end_block_col,\n", + " ] = operators.slip_rate_constraints\n", + "\n", + " # Insert eigenvector to velocity matrix\n", + " # TODO: Generalize over multiple meshes\n", + " # Eliminate vertical elastic velocities\n", + " tde_keep_row_index = celeri.get_keep_index_12(\n", + " operators.tde_to_velocities[0].shape[0]\n", + " )\n", + " tde_keep_col_index = celeri.get_keep_index_12(\n", + " operators.tde_to_velocities[0].shape[1]\n", + " )\n", + " eigen_to_velocities = (\n", + " -operators.tde_to_velocities[0][tde_keep_row_index, :][:, tde_keep_col_index]\n", + " @ eigenvectors_two_component\n", + " )\n", + "\n", + " # TODO: Update index here with eigen specific terms\n", + " # TODO: Generalize over multiple meshes\n", + " operator[\n", + " index.start_station_row : index.end_station_row,\n", + " index.start_tde_col[0] : index.start_tde_col[0] + 2 * meshes[0].n_modes,\n", + " ] = eigen_to_velocities\n", + "\n", + " # Insert eigenvector to TDE constraints matrix\n", + " # TODO: Generalize over multiple meshes\n", + " print(meshes[0].mesh_tde_modes_bc_weight)\n", + " eigen_to_tde_slip_rate_constraints = (\n", + " meshes[0].mesh_tde_modes_bc_weight\n", + " * operators.tde_slip_rate_constraints[0]\n", + " @ eigenvectors_two_component\n", + " )\n", + "\n", + " operator[\n", + " index.start_slip_rate_constraints_row : index.start_slip_rate_constraints_row\n", + " + index.n_tde_constraints_total,\n", + " index.start_tde_col[0] : index.start_tde_col[0] + 2 * meshes[0].n_modes,\n", + " ] = eigen_to_tde_slip_rate_constraints\n", + "\n", + " # # Insert TDE to velocity matrix\n", + " # for i in range(len(meshes)):\n", + " # # Insert TDE to velocity matrix\n", + " # tde_keep_row_index = celeri.get_keep_index_12(operators.tde_to_velocities[i].shape[0])\n", + " # tde_keep_col_index = celeri.get_keep_index_12(operators.tde_to_velocities[i].shape[1])\n", + " # operator[\n", + " # index.start_station_row : index.end_station_row,\n", + " # index.start_tde_col[i] : index.end_tde_col[i],\n", + " # ] = -operators.tde_to_velocities[i][tde_keep_row_index, :][\n", + " # :, tde_keep_col_index\n", + " # ]\n", + "\n", + " # # # Insert TDE smoothing matrix\n", + " # # smoothing_keep_index = celeri.get_keep_index_12(\n", + " # # operators.tde_to_velocities[i].shape[1]\n", + " # # )\n", + " # # operator[\n", + " # # index.start_tde_smoothing_row[i] : index.end_tde_smoothing_row[i],\n", + " # # index.start_tde_col[i] : index.end_tde_col[i],\n", + " # # ] = operators.smoothing_matrix[i].toarray()[smoothing_keep_index, :][\n", + " # # :, smoothing_keep_index\n", + " # # ]\n", + "\n", + " # # Insert TDE slip rate constraints into estimation operator\n", + " # operator[\n", + " # index.start_tde_constraint_row[i] : index.end_tde_constraint_row[i],\n", + " # index.start_tde_col[i] : index.end_tde_col[i],\n", + " # ] = operators.tde_slip_rate_constraints[i]\n", + " return operator\n", + "\n", + "\n", + "def post_process_estimation_eigen(estimation_eigen, operators, station, index):\n", + " \"\"\"Calculate derived values derived from the block model linear estimate (e.g., velocities, undertainties)\n", + "\n", + " Args:\n", + " estimation (Dict): Estimated state vector and model covariance\n", + " operators (Dict): All linear operators\n", + " station (pd.DataFrame): GPS station data\n", + " idx (Dict): Indices and counts of data and array sizes\n", + " \"\"\"\n", + "\n", + " estimation_eigen.eigenvalues = estimation_eigen.state_vector[3 * n_blocks : :]\n", + " estimation_eigen.predictions = (\n", + " estimation_eigen.operator @ estimation_eigen.state_vector\n", + " )\n", + " estimation_eigen.vel = estimation_eigen.predictions[0 : 2 * index.n_stations]\n", + " estimation_eigen.east_vel = estimation_eigen.vel[0::2]\n", + " estimation_eigen.north_vel = estimation_eigen.vel[1::2]\n", + "\n", + " # Estimate slip rate uncertainties\n", + " estimation_eigen.slip_rate_sigma = np.sqrt(\n", + " np.diag(\n", + " operators.rotation_to_slip_rate\n", + " @ estimation_eigen.state_covariance_matrix[\n", + " 0 : 3 * index.n_blocks, 0 : 3 * index.n_blocks\n", + " ]\n", + " @ operators.rotation_to_slip_rate.T\n", + " )\n", + " ) # I don't think this is correct because for the case when there is a rotation vector a priori\n", + " estimation_eigen.strike_slip_rate_sigma = estimation_eigen.slip_rate_sigma[0::3]\n", + " estimation_eigen.dip_slip_rate_sigma = estimation_eigen.slip_rate_sigma[1::3]\n", + " estimation_eigen.tensile_slip_rate_sigma = estimation_eigen.slip_rate_sigma[2::3]\n", + "\n", + " # Calculate mean squared residual velocity\n", + " estimation_eigen.east_vel_residual = estimation_eigen.east_vel - station.east_vel\n", + " estimation_eigen.north_vel_residual = estimation_eigen.north_vel - station.north_vel\n", + "\n", + " # Extract TDE slip rates from state vector\n", + " estimation_eigen.tde_rates = (\n", + " eigenvectors_two_component @ estimation_eigen.state_vector[3 * n_blocks : :]\n", + " )\n", + "\n", + " # Isolate strike- and dip-slip rates\n", + " estimation_eigen.tde_strike_slip_rates = estimation_eigen.tde_rates[0::2]\n", + " estimation_eigen.tde_dip_slip_rates = estimation_eigen.tde_rates[1::2]\n", + "\n", + " # Create a pseudo state vector that is the length of a TDE state vector\n", + " estimation_eigen.pseudo_tde_state_vector = np.zeros(\n", + " 3 * n_blocks + 2 * index.n_tde_total\n", + " )\n", + " estimation_eigen.pseudo_tde_state_vector[0 : 3 * index.n_blocks] = (\n", + " estimation_eigen.state_vector[0 : 3 * index.n_blocks]\n", + " )\n", + " estimation_eigen.pseudo_tde_state_vector[3 * index.n_blocks : :] = (\n", + " estimation_eigen.tde_rates\n", + " )\n", + "\n", + " # Extract segment slip rates from state vector\n", + " estimation_eigen.slip_rates = (\n", + " operators.rotation_to_slip_rate\n", + " @ estimation_eigen.state_vector[0 : 3 * index.n_blocks]\n", + " )\n", + " estimation_eigen.strike_slip_rates = estimation_eigen.slip_rates[0::3]\n", + " estimation_eigen.dip_slip_rates = estimation_eigen.slip_rates[1::3]\n", + " estimation_eigen.tensile_slip_rates = estimation_eigen.slip_rates[2::3]\n", + "\n", + " # Calculate rotation only velocities\n", + " estimation_eigen.vel_rotation = (\n", + " operators.rotation_to_velocities[index.station_row_keep_index, :]\n", + " @ estimation_eigen.state_vector[0 : 3 * index.n_blocks]\n", + " )\n", + " estimation_eigen.east_vel_rotation = estimation_eigen.vel_rotation[0::2]\n", + " estimation_eigen.north_vel_rotation = estimation_eigen.vel_rotation[1::2]\n", + "\n", + " # Calculate fully locked segment velocities\n", + " estimation_eigen.vel_elastic_segment = (\n", + " operators.rotation_to_slip_rate_to_okada_to_velocities[\n", + " index.station_row_keep_index, :\n", + " ]\n", + " @ estimation_eigen.state_vector[0 : 3 * index.n_blocks]\n", + " )\n", + " estimation_eigen.east_vel_elastic_segment = estimation_eigen.vel_elastic_segment[\n", + " 0::2\n", + " ]\n", + " estimation_eigen.north_vel_elastic_segment = estimation_eigen.vel_elastic_segment[\n", + " 1::2\n", + " ]\n", + "\n", + " # TODO: Calculate block strain rate velocities\n", + " estimation_eigen.east_vel_block_strain_rate = np.zeros(len(station))\n", + " estimation_eigen.north_vel_block_strain_rate = np.zeros(len(station))\n", + "\n", + " # Calculate TDE velocities\n", + " estimation_eigen.vel_tde = np.zeros(2 * index.n_stations)\n", + " for i in range(len(operators.tde_to_velocities)):\n", + "\n", + " tde_keep_row_index = celeri.get_keep_index_12(\n", + " operators.tde_to_velocities[i].shape[0]\n", + " )\n", + " tde_keep_col_index = celeri.get_keep_index_12(\n", + " operators.tde_to_velocities[i].shape[1]\n", + " )\n", + " estimation_eigen.vel_tde += (\n", + " operators.tde_to_velocities[i][tde_keep_row_index, :][:, tde_keep_col_index]\n", + " @ eigenvectors_two_component\n", + " @ estimation_eigen.state_vector[\n", + " index.start_tde_col[i] : index.end_tde_col[i]\n", + " ]\n", + " # @ estimation_eigen.eigenvalues\n", + " )\n", + " estimation_eigen.east_vel_tde = estimation_eigen.vel_tde[0::2]\n", + " estimation_eigen.north_vel_tde = estimation_eigen.vel_tde[1::2]" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001b[32m2024-07-13 16:12:26.547\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36m__main__\u001b[0m:\u001b[36mget_eigenvectors_to_tde_slip\u001b[0m:\u001b[36m32\u001b[0m - \u001b[1mStart: Eigenvectors to TDE slip for mesh: ../data/mesh/cascadia.msh\u001b[0m\n", + "\u001b[32m2024-07-13 16:12:27.284\u001b[0m | \u001b[32m\u001b[1mSUCCESS \u001b[0m | \u001b[36m__main__\u001b[0m:\u001b[36mget_eigenvectors_to_tde_slip\u001b[0m:\u001b[36m61\u001b[0m - \u001b[32m\u001b[1mFinish: Eigenvectors to TDE slip for mesh: ../data/mesh/cascadia.msh\u001b[0m\n", + "eigenvectors_two_component.shape=(3682, 60)\n", + "operator.shape=(3519, 153)\n", + "10.0\n" + ] + } + ], + "source": [ + "# Get eigenvectors\n", + "_, eigenvectors = get_eigenvalues_and_eigenvectors(\n", + " meshes[0].n_modes,\n", + " meshes[0].x_centroid,\n", + " meshes[0].y_centroid,\n", + " meshes[0].z_centroid,\n", + " distance_exponent=1.0,\n", + ")\n", + "\n", + "# TODO: Push to function\n", + "# Create eigevectors to TDE slip matrix\n", + "# This shouls be in operators ???\n", + "eigenvectors_two_component = np.zeros(\n", + " (\n", + " 2 * eigenvectors.shape[0],\n", + " meshes[0].n_modes_strike_slip + meshes[0].n_modes_dip_slip,\n", + " )\n", + ")\n", + "eigenvectors_two_component[0::2, 0 : meshes[0].n_modes_strike_slip] = eigenvectors[\n", + " :, 0 : meshes[0].n_modes_strike_slip\n", + "]\n", + "eigenvectors_two_component[\n", + " 1::2,\n", + " meshes[0].n_modes_strike_slip : meshes[0].n_modes_strike_slip\n", + " + meshes[0].n_modes_dip_slip,\n", + "] = eigenvectors[:, 0 : meshes[0].n_modes_dip_slip]\n", + "\n", + "\n", + "# TODO: Replace all mentions of eigenvectors_two_component with operators.eigenvectors_to_tde_slip\n", + "def get_eigenvectors_to_tde_slip(operators, meshes):\n", + " for i in range(len(meshes)):\n", + " logger.info(f\"Start: Eigenvectors to TDE slip for mesh: {meshes[i].file_name}\")\n", + " # Get eigenvectors for curren mesh\n", + " _, eigenvectors = get_eigenvalues_and_eigenvectors(\n", + " meshes[i].n_modes,\n", + " meshes[i].x_centroid,\n", + " meshes[i].y_centroid,\n", + " meshes[i].z_centroid,\n", + " distance_exponent=1.0, # Make this something set in mesh_parameters.json\n", + " )\n", + "\n", + " # Create eigenvectors to TDE slip matrix\n", + " operators.eigenvectors_to_tde_slip[i] = np.zeros(\n", + " (\n", + " 2 * eigenvectors.shape[0],\n", + " meshes[i].n_modes_strike_slip + meshes[i].n_modes_dip_slip,\n", + " )\n", + " )\n", + "\n", + " # Place strike-slip panel\n", + " operators.eigenvectors_to_tde_slip[i][\n", + " 0::2, 0 : meshes[i].n_modes_strike_slip\n", + " ] = eigenvectors[:, 0 : meshes[i].n_modes_strike_slip]\n", + "\n", + " # Place dip-slip panel\n", + " operators.eigenvectors_to_tde_slip[i][\n", + " 1::2,\n", + " meshes[i].n_modes_strike_slip : meshes[i].n_modes_strike_slip\n", + " + meshes[i].n_modes_dip_slip,\n", + " ] = eigenvectors[:, 0 : meshes[i].n_modes_dip_slip]\n", + " logger.success(\n", + " f\"Finish: Eigenvectors to TDE slip for mesh: {meshes[i].file_name}\"\n", + " )\n", + "\n", + "\n", + "get_eigenvectors_to_tde_slip(operators, meshes)\n", + "\n", + "\n", + "data_vector_eigen = get_data_vector_eigen(assembly, index)\n", + "weighting_vector_eigen = get_weighting_vector_eigen(command, station, meshes, index)\n", + "operator_eigen = get_full_dense_operator_eigen(operators, meshes, index)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# KL + QP solve" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "True\n", + "True\n", + "True\n", + "True\n" + ] + } + ], + "source": [ + "# Options dictionary for lsqlin_qp\n", + "opts = {\"show_progress\": False}\n", + "\n", + "lower_bound = -np.ones_like(estimation.state_vector)\n", + "upper_bound = np.ones_like(estimation.state_vector)\n", + "\n", + "# TDE strike-slip\n", + "lower_bound[index.start_tde_col[0] : index.end_tde_col[0] : 2] = -5\n", + "upper_bound[index.start_tde_col[0] : index.end_tde_col[0] : 2] = 5\n", + "\n", + "# TDE dip-slip\n", + "lower_bound[index.start_tde_col[0] + 1 : index.end_tde_col[0] : 2] = 0\n", + "upper_bound[index.start_tde_col[0] + 1 : index.end_tde_col[0] : 2] = 30\n", + "\n", + "lower_bound_tde_only = lower_bound[3 * n_blocks : :]\n", + "upper_bound_tde_only = upper_bound[3 * n_blocks : :]\n", + "\n", + "A = np.zeros((2 * eigenvectors_two_component.shape[0], n_state_vector_eigen))\n", + "A[0 : eigenvectors_two_component.shape[0], 3 * n_blocks :] = -eigenvectors_two_component\n", + "A[eigenvectors_two_component.shape[0] :, 3 * n_blocks :] = eigenvectors_two_component\n", + "\n", + "# TODO: Need better use of index for shape here\n", + "b = np.zeros(2 * eigenvectors_two_component.shape[0])\n", + "b[0 : eigenvectors_two_component.shape[0]] = -lower_bound_tde_only\n", + "b[eigenvectors_two_component.shape[0] :] = upper_bound_tde_only\n", + "\n", + "\n", + "# TODO: Better name for ret\n", + "ret = celeri.lsqlin_qp(\n", + " operator_eigen * np.sqrt(weighting_vector_eigen[:, None]),\n", + " data_vector_eigen * np.sqrt(weighting_vector_eigen),\n", + " 0,\n", + " A, # A=None, inequality operator\n", + " b, # b=None, inequality vector\n", + " None,\n", + " None,\n", + " None,\n", + " None,\n", + " None,\n", + " opts,\n", + ")\n", + "\n", + "estimation_eigen_cvxopt_bounded = addict.Dict()\n", + "estimation_eigen_cvxopt_bounded.state_covariance_matrix = np.linalg.inv(\n", + " operator_eigen.T * weighting_vector_eigen @ operator_eigen\n", + ")\n", + "estimation_eigen_cvxopt_bounded.state_vector = np.array(ret[\"x\"]).flatten()\n", + "forward_data_vector_eigen = (\n", + " operator_eigen @ estimation_eigen_cvxopt_bounded.state_vector\n", + ")\n", + "estimation_eigen_cvxopt_bounded.operator = operator_eigen\n", + "post_process_estimation_eigen(\n", + " estimation_eigen_cvxopt_bounded, operators, station, index\n", + ")\n", + "\n", + "test_agreement(estimation_eigen_cvxopt_bounded)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "\n", + "\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "operators.tde_to_velocities[0]\n", + "operators.eigenvectors_two_component[0] = eigenvectors_two_component\n", + "print(type(operators.tde_to_velocities))\n", + "print(type(operators.eigenvectors_two_component))\n", + "print(type(operators.eigenvectors_to_tde_slip))\n", + "np.allclose(operators.eigenvectors_to_tde_slip[0], eigenvectors_two_component)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Eigen specific values `index`" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "# def update_index_eigen(index, meshes):\n", + "# index.eigen_start_tde_constraint_row = np.zeros(len(meshes), dtype=int)\n", + "# index.eigen_end_tde_constraint_row = np.zeros(len(meshes), dtype=int)\n", + "# index.eigen_start_tde_constraint_row[0] = index.end_slip_rate_constraints_row\n", + "# index.eigen_end_tde_constraint_row[0] = (\n", + "# index.eigen_start_tde_constraint_row[0] + index.n_tde_constraints[0]\n", + "# )\n", + "# # Starting at 1 because the 0 case is already done (above)\n", + "# for i in range(1, len(meshes)):\n", + "# index.eigen_start_tde_constraint_row[i] = index.eigen_end_tde_constraint_row[\n", + "# i - 1\n", + "# ]\n", + "# index.eigen_end_tde_constraint_row[i] = (\n", + "# index.eigen_start_tde_constraint_row[i - 1] + index.n_tde_constraints[i]\n", + "# )\n", + "\n", + "\n", + "# def get_weighting_vector_eigen(command, station, meshes, index):\n", + "# # Initialize and build weighting matrix\n", + "# weighting_vector_eigen = np.ones(\n", + "# 2 * index.n_stations\n", + "# + 3 * index.n_block_constraints\n", + "# + index.n_slip_rate_constraints\n", + "# + index.n_tde_constraints_total\n", + "# )\n", + "\n", + "# # Velocity uncertainties\n", + "# weighting_vector_eigen[index.start_station_row : index.end_station_row] = (\n", + "# celeri.interleave2(1 / (station.east_sig**2), 1 / (station.north_sig**2))\n", + "# )\n", + "\n", + "# # Soft prior block rotation weights\n", + "# weighting_vector_eigen[\n", + "# index.start_block_constraints_row : index.end_block_constraints_row\n", + "# ] = command.block_constraint_weight\n", + "\n", + "# # Soft prior slip rate weights\n", + "# weighting_vector_eigen[\n", + "# index.start_slip_rate_constraints_row : index.end_slip_rate_constraints_row\n", + "# ] = command.slip_constraint_weight * np.ones(index.n_slip_rate_constraints)\n", + "\n", + "# # Soft prior TDE weights\n", + "# for i in range(len(meshes)):\n", + "# weighting_vector_eigen[\n", + "# index.eigen_start_tde_constraint_row[i] : index.end_tde_constraint_row[i]\n", + "# ] = command.tri_con_weight * np.ones(index.n_tde_constraints[i])\n", + "# return weighting_vector_eigen\n", + "\n", + "\n", + "# update_index_eigen(index, meshes)\n", + "# wve = get_weighting_vector_eigen(command, station, meshes, index)\n", + "# np.allclose(weighting_vector_eigen, wve)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "celeri", + "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.10.9" + }, + "vscode": { + "interpreter": { + "hash": "4e02521ee8166e85f0cdf9248b501b87197c4fbf1c25b3c3121662d555f974cc" + } + }, + "widgets": { + "application/vnd.jupyter.widget-state+json": { + "state": { + "4d250c5d35aa493295ca814fb3eaa1ee": { + "model_module": "@jupyter-widgets/base", + "model_module_version": "1.2.0", + "model_name": "LayoutModel", + "state": {} + }, + "6faf75ca5f3b41388f284e98ec2cf803": { + "model_module": "jupyter-matplotlib", + "model_module_version": "^0.9.0", + "model_name": "ToolbarModel", + "state": { + "layout": "IPY_MODEL_9b061db2dc65459ca586b9b9f73c2362", + "toolitems": [ + [ + "Home", + "Reset original view", + "home", + "home" + ], + [ + "Back", + "Back to previous view", + "arrow-left", + "back" + ], + [ + "Forward", + "Forward to next view", + "arrow-right", + "forward" + ], + [ + "Pan", + "Left button pans, Right button zooms\nx/y fixes axis, CTRL fixes aspect", + "arrows", + "pan" + ], + [ + "Zoom", + "Zoom to rectangle\nx/y fixes axis, CTRL fixes aspect", + "square-o", + "zoom" + ], + [ + "Download", + "Download plot", + "floppy-o", + "save_figure" + ] + ] + } + }, + "9b061db2dc65459ca586b9b9f73c2362": { + "model_module": "@jupyter-widgets/base", + "model_module_version": "1.2.0", + "model_name": "LayoutModel", + "state": {} + }, + "c25a38234e8f4e818670d9767f95a430": { + "model_module": "jupyter-matplotlib", + "model_module_version": "^0.9.0", + "model_name": "MPLCanvasModel", + "state": { + "_cursor": "default", + "_figure_label": "Figure 1", + "_height": 708, + "_width": 1746, + "layout": "IPY_MODEL_4d250c5d35aa493295ca814fb3eaa1ee", + "toolbar": "IPY_MODEL_6faf75ca5f3b41388f284e98ec2cf803", + "toolbar_position": "left" + } + } + }, + "version_major": 2, + "version_minor": 0 + } + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}