eigenvector to tde slip matrix construction now loops over meshes and…

… is a function. Need to replace all instances of eigenvetors_two_component.
brendanjmeade · Jul 13, 2024 · 9568ffe · 9568ffe
1 parent 1eb1f52
commit 9568ffe
Showing 1 changed file with 143 additions and 72 deletions.
diff --git 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": [
+      "<class 'addict.addict.Dict'>\n",
+      "<class 'addict.addict.Dict'>\n",
+      "<class 'addict.addict.Dict'>\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": [