eigenvectors_two_component eliminated!

brendanjmeade · Jul 13, 2024 · 7bf5c43 · 7bf5c43
1 parent 7ea1c25
commit 7bf5c43
Showing 1 changed file with 39 additions and 95 deletions.
diff --git a/notebooks/celeri_western_north_america_qp_05.ipynb b/notebooks/celeri_western_north_america_qp_05.ipynb
@@ -40,18 +40,18 @@
      "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",
+      "\u001b[32m2024-07-13 16:36:50.175\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:36:50.176\u001b[0m | \u001b[1mINFO    \u001b[0m | \u001b[36mceleri.celeri\u001b[0m:\u001b[36mget_logger\u001b[0m:\u001b[36m6344\u001b[0m - \u001b[1mRUN_NAME: 0000000047\u001b[0m\n",
+      "\u001b[32m2024-07-13 16:36:50.177\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/0000000047/0000000047.log\u001b[0m\n",
+      "\u001b[32m2024-07-13 16:36:50.177\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:36:50.182\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:36:50.184\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:36:50.185\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"
+      "\u001b[32m2024-07-13 16:36:50.296\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:36:50.299\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:36:50.300\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:36:50.301\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"
      ]
     }
    ],
@@ -143,9 +143,9 @@
      "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"
+      "\u001b[32m2024-07-13 16:36:51.108\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:36:51.680\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:36:51.683\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"
      ]
     }
    ],
@@ -398,7 +398,7 @@
     "    )\n",
     "    eigen_to_velocities = (\n",
     "        -operators.tde_to_velocities[0][tde_keep_row_index, :][:, tde_keep_col_index]\n",
-    "        @ eigenvectors_two_component\n",
+    "        @ operators.eigenvectors_to_tde_slip[0]\n",
     "    )\n",
     "\n",
     "    # TODO: Update index here with eigen specific terms\n",
@@ -414,7 +414,7 @@
     "    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",
+    "        @ operators.eigenvectors_to_tde_slip[0]\n",
     "    )\n",
     "\n",
     "    operator[\n",
@@ -492,7 +492,8 @@
     "\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",
+    "        operators.eigenvectors_to_tde_slip[0]\n",
+    "        @ estimation_eigen.state_vector[3 * n_blocks : :]\n",
     "    )\n",
     "\n",
     "    # Isolate strike- and dip-slip rates\n",
@@ -557,7 +558,7 @@
     "        )\n",
     "        estimation_eigen.vel_tde += (\n",
     "            operators.tde_to_velocities[i][tde_keep_row_index, :][:, tde_keep_col_index]\n",
-    "            @ eigenvectors_two_component\n",
+    "            @ operators.eigenvectors_to_tde_slip[0]\n",
     "            @ estimation_eigen.state_vector[\n",
     "                index.start_tde_col[i] : index.end_tde_col[i]\n",
     "            ]\n",
@@ -576,44 +577,19 @@
      "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",
+      "\u001b[32m2024-07-13 16:36:54.348\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"
+     ]
+    },
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "\u001b[32m2024-07-13 16:36:55.109\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",
       "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",
@@ -651,8 +627,6 @@
     "\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)"
@@ -699,15 +673,19 @@
     "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",
+    "# Write linear inequality constraints for bounds on TDE slip\n",
+    "A = np.zeros((2 * operators.eigenvectors_to_tde_slip[0].shape[0], n_state_vector_eigen))\n",
+    "A[0 : operators.eigenvectors_to_tde_slip[0].shape[0], 3 * n_blocks :] = (\n",
+    "    -operators.eigenvectors_to_tde_slip[0]\n",
+    ")\n",
+    "A[operators.eigenvectors_to_tde_slip[0].shape[0] :, 3 * n_blocks :] = (\n",
+    "    operators.eigenvectors_to_tde_slip[0]\n",
+    ")\n",
     "\n",
+    "b = np.zeros(2 * operators.eigenvectors_to_tde_slip[0].shape[0])\n",
+    "b[0 : operators.eigenvectors_to_tde_slip[0].shape[0]] = -lower_bound_tde_only\n",
+    "b[operators.eigenvectors_to_tde_slip[0].shape[0] :] = upper_bound_tde_only\n",
     "\n",
     "# TODO: Better name for ret\n",
     "ret = celeri.lsqlin_qp(\n",
@@ -740,40 +718,6 @@
     "test_agreement(estimation_eigen_cvxopt_bounded)"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 11,
-   "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "<class 'addict.addict.Dict'>\n",
-      "<class 'addict.addict.Dict'>\n",
-      "<class 'addict.addict.Dict'>\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": {},
@@ -783,7 +727,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [],
    "source": [