From 9568ffe7779db4d06ec0e0f7b5aebfde9024350d Mon Sep 17 00:00:00 2001 From: Brendan Meade Date: Fri, 12 Jul 2024 23:09:45 -0400 Subject: [PATCH] eigenvector to tde slip matrix construction now loops over meshes and is a function. Need to replace all instances of eigenvetors_two_component. --- .../celeri_western_north_america_qp_04.ipynb | 215 ++++++++++++------ 1 file changed, 143 insertions(+), 72 deletions(-) diff --git a/notebooks/celeri_western_north_america_qp_04.ipynb b/notebooks/celeri_western_north_america_qp_04.ipynb index 6e41129..446e59f 100644 --- a/notebooks/celeri_western_north_america_qp_04.ipynb +++ b/notebooks/celeri_western_north_america_qp_04.ipynb @@ -40,18 +40,18 @@ "name": "stdout", "output_type": "stream", "text": [ - "\u001b[32m2024-07-12 17:48:51.947\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-12 17:48:51.947\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mget_logger\u001b[0m:\u001b[36m6344\u001b[0m - \u001b[1mRUN_NAME: 0000000022\u001b[0m\n", - "\u001b[32m2024-07-12 17:48:51.948\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/0000000022/0000000022.log\u001b[0m\n", - "\u001b[32m2024-07-12 17:48:51.948\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-12 17:48:51.952\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-12 17:48:51.954\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-12 17:48:51.954\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", + "\u001b[32m2024-07-12 18:23:51.544\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-12 18:23:51.545\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mget_logger\u001b[0m:\u001b[36m6344\u001b[0m - \u001b[1mRUN_NAME: 0000000029\u001b[0m\n", + "\u001b[32m2024-07-12 18:23:51.545\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/0000000029/0000000029.log\u001b[0m\n", + "\u001b[32m2024-07-12 18:23:51.546\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-12 18:23:51.550\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-12 18:23:51.552\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-12 18:23:51.553\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-12 17:48:52.066\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-12 17:48:52.069\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-12 17:48:52.070\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-12 17:48:52.071\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" + "\u001b[32m2024-07-12 18:23:51.665\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-12 18:23:51.669\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-12 18:23:51.670\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-12 18:23:51.671\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" ] } ], @@ -141,9 +141,9 @@ "name": "stdout", "output_type": "stream", "text": [ - "\u001b[32m2024-07-12 17:48:52.822\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-12 17:48:53.492\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-12 17:48:53.496\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" + "\u001b[32m2024-07-12 18:23:52.422\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-12 18:23:52.990\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-12 18:23:52.994\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" ] } ], @@ -219,7 +219,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# KL machinery (panel)" + "# KL functions" ] }, { @@ -441,61 +441,7 @@ " # 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" - ] - }, - { - "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "n_eigen=30.0\n", - "eigenvectors_two_component.shape=(3682, 60)\n", - "operator.shape=(3519, 153)\n", - "10.0\n" - ] - } - ], - "source": [ - "# Calculate eigenvales and eigenvectors\n", - "eigenvalues, 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", - "\n", - "# TODO: Push to function\n", - "# Create eigenmodes to TDE slip matrix (panel style)\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", - "\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", - "data_vector = celeri.get_data_vector(assembly, index, meshes)\n", - "data_vector_eigen = get_data_vector_eigen(assembly, index)\n", - "weighting_vector = celeri.get_weighting_vector(command, station, meshes, index)\n", - "weighting_vector_eigen = get_weighting_vector_eigen(command, station, meshes, index)\n", - "operator = celeri.get_full_dense_operator(operators, meshes, index)\n", - "operator_eigen = get_full_dense_operator_eigen(operators, meshes, index)\n", + " return operator\n", "\n", "\n", "def post_process_estimation_eigen(estimation_eigen, operators, station, index):\n", @@ -611,6 +557,97 @@ " estimation_eigen.north_vel_tde = estimation_eigen.vel_tde[1::2]" ] }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\u001b[32m2024-07-12 23:04:37.967\u001b[0m | \u001b[1mINFO \u001b[0m | \u001b[36m__main__\u001b[0m:\u001b[36mget_eigenvectors_to_tde_slip\u001b[0m:\u001b[36m43\u001b[0m - \u001b[1mStart: Modes to TDE slip for mesh: ../data/mesh/cascadia.msh\u001b[0m\n", + "\u001b[32m2024-07-12 23:04:38.694\u001b[0m | \u001b[32m\u001b[1mSUCCESS \u001b[0m | \u001b[36m__main__\u001b[0m:\u001b[36mget_eigenvectors_to_tde_slip\u001b[0m:\u001b[36m72\u001b[0m - \u001b[32m\u001b[1mFinish: Modes to TDE slip for mesh: ../data/mesh/cascadia.msh\u001b[0m\n", + "n_eigen=30.0\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", + "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": {}, @@ -620,7 +657,7 @@ }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 23, "metadata": {}, "outputs": [ { @@ -693,6 +730,40 @@ "test_agreement(estimation_eigen_cvxopt_bounded)" ] }, + { + "cell_type": "code", + "execution_count": 25, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "\n", + "\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 25, + "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": {}, @@ -702,7 +773,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 18, "metadata": {}, "outputs": [], "source": [