364 lines
14 KiB
Plaintext
364 lines
14 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"id": "initial_id",
|
|
"metadata": {
|
|
"collapsed": true
|
|
},
|
|
"source": [
|
|
"%load_ext autoreload\n",
|
|
"%autoreload 2"
|
|
],
|
|
"outputs": [],
|
|
"execution_count": null
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"source": [
|
|
"%reload_ext autoreload\n",
|
|
"%autoreload 2\n",
|
|
"import time\n",
|
|
"import pickle\n",
|
|
"import numpy as np\n",
|
|
"from numpy import nan\n",
|
|
"import winkelumrechnungen as wu\n",
|
|
"import os\n",
|
|
"from contextlib import contextmanager, redirect_stdout, redirect_stderr\n",
|
|
"import pandas as pd\n",
|
|
"import plotly.graph_objects as go\n",
|
|
"import math\n",
|
|
"\n",
|
|
"from ellipsoide import EllipsoidTriaxial\n",
|
|
"from GHA_triaxial.panou import alpha_ell2para, alpha_para2ell\n",
|
|
"\n",
|
|
"from GHA_triaxial.panou import gha1_num, gha1_ana\n",
|
|
"from GHA_triaxial.approx_gha1 import gha1_approx\n",
|
|
"\n",
|
|
"from GHA_triaxial.panou_2013_2GHA_num import gha2_num\n",
|
|
"from GHA_triaxial.ES_gha2 import gha2_ES\n",
|
|
"from GHA_triaxial.approx_gha2 import gha2_approx\n",
|
|
"\n",
|
|
"from GHA_triaxial.numeric_examples_panou import get_tables as get_tables_panou\n",
|
|
"from GHA_triaxial.numeric_examples_karney import get_random_examples as get_examples_karney"
|
|
],
|
|
"id": "961cb22764c5bcb9",
|
|
"outputs": [],
|
|
"execution_count": null
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"source": [
|
|
"@contextmanager\n",
|
|
"def suppress_print():\n",
|
|
" with open(os.devnull, 'w') as fnull:\n",
|
|
" with redirect_stdout(fnull), redirect_stderr(fnull):\n",
|
|
" yield"
|
|
],
|
|
"id": "e4740625a366ac13",
|
|
"outputs": [],
|
|
"execution_count": null
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"source": [
|
|
"# steps_gha1_num = [2000, 5000, 10000, 20000]\n",
|
|
"# maxM_gha1_ana = [20, 50]\n",
|
|
"# parts_gha1_ana = [4, 8]\n",
|
|
"# dsPart_gha1_approx = [600, 1250, 6000, 60000] # entspricht bei der Erde ca. 10000, 5000, 1000, 100\n",
|
|
"#\n",
|
|
"# steps_gha2_num = [2000, 5000, 10000, 20000]\n",
|
|
"# dsPart_gha2_ES = [600, 1250, 6000] # entspricht bei der Erde ca. 10000, 5000, 1000\n",
|
|
"# dsPart_gha2_approx = [600, 1250, 6000, 60000] # entspricht bei der Erde ca. 10000, 5000, 1000, 100\n",
|
|
"\n",
|
|
"steps_gha1_num = [2000]\n",
|
|
"maxM_gha1_ana = [10, 20, 40, 60, 80]\n",
|
|
"parts_gha1_ana = [4, 8, 12, 16]\n",
|
|
"dsPart_gha1_approx = [600]\n",
|
|
"\n",
|
|
"steps_gha2_num = [16000]\n",
|
|
"dsPart_gha2_ES = [20]\n",
|
|
"dsPart_gha2_approx = [600]"
|
|
],
|
|
"id": "96093cdde03f8d57",
|
|
"outputs": [],
|
|
"execution_count": null
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"source": [
|
|
"# test = \"Karney\"\n",
|
|
"test = \"Panou\"\n",
|
|
"if test == \"Karney\":\n",
|
|
" ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(\"KarneyTest2024\")\n",
|
|
" examples = get_examples_karney(10, 42)\n",
|
|
"elif test == \"Panou\":\n",
|
|
" ell: EllipsoidTriaxial = EllipsoidTriaxial.init_name(\"BursaSima1980round\")\n",
|
|
" tables = get_tables_panou()\n",
|
|
" table_indices = []\n",
|
|
" examples = []\n",
|
|
" for i, table in enumerate(tables):\n",
|
|
" for example in table:\n",
|
|
" table_indices.append(i+1)\n",
|
|
" examples.append(example)"
|
|
],
|
|
"id": "6e384cc01c2dbe",
|
|
"outputs": [],
|
|
"execution_count": null
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"source": [
|
|
"results = {}\n",
|
|
"for example in examples:\n",
|
|
" example_results = {}\n",
|
|
"\n",
|
|
" beta0, lamb0, alpha0_ell, beta1, lamb1, alpha1_ell, s = example\n",
|
|
" P0 = ell.ell2cart(beta0, lamb0)\n",
|
|
" P1 = ell.ell2cart(beta1, lamb1)\n",
|
|
" _, _, alpha0_para = alpha_ell2para(ell, beta0, lamb0, alpha0_ell)\n",
|
|
"\n",
|
|
" # for steps in steps_gha1_num:\n",
|
|
" # start = time.perf_counter()\n",
|
|
" # try:\n",
|
|
" # P1_num, alpha1_num_1 = gha1_num(ell, P0, alpha0_ell, s, num=steps)\n",
|
|
" # end = time.perf_counter()\n",
|
|
" # beta1_num, lamb1_num = ell.cart2ell(P1_num)\n",
|
|
" # d_beta1 = abs(wu.rad2deg(beta1_num - beta1)) / 3600\n",
|
|
" # d_lamb1 = abs(wu.rad2deg(lamb1_num - lamb1)) / 3600\n",
|
|
" # d_alpha1 = abs(wu.rad2deg(alpha1_num_1 - alpha1_ell)) / 3600\n",
|
|
" # d_time = end - start\n",
|
|
" # example_results[f\"GHA1_num_{steps}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n",
|
|
" # except Exception as e:\n",
|
|
" # print(e)\n",
|
|
" # example_results[f\"GHA1_num_{steps}\"] = (nan, nan, nan, nan)\n",
|
|
" #\n",
|
|
" for maxM in maxM_gha1_ana:\n",
|
|
" for parts in parts_gha1_ana:\n",
|
|
" start = time.perf_counter()\n",
|
|
" try:\n",
|
|
" P1_ana, alpha1_ana_para = gha1_ana(ell, P0, alpha0_para, s, maxM=maxM, maxPartCircum=parts)\n",
|
|
" end = time.perf_counter()\n",
|
|
" beta1_ana, lamb1_ana = ell.cart2ell(P1_ana)\n",
|
|
" _, _, alpha1_ana_ell = alpha_para2ell(ell, beta1_ana, lamb1_ana, alpha1_ana_para)\n",
|
|
" d_beta1 = abs(wu.rad2deg(beta1_ana - beta1)) / 3600\n",
|
|
" d_lamb1 = abs(wu.rad2deg(lamb1_ana - lamb1)) / 3600\n",
|
|
" d_alpha1 = abs(wu.rad2deg(alpha1_ana_ell - alpha1_ell)) / 3600\n",
|
|
" d_time = end - start\n",
|
|
" example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n",
|
|
" except Exception as e:\n",
|
|
" print(e)\n",
|
|
" example_results[f\"GHA1_ana_{maxM}_{parts}\"] = (nan, nan, nan, nan)\n",
|
|
" #\n",
|
|
" # for dsPart in dsPart_gha1_approx:\n",
|
|
" # ds = ell.ax/dsPart\n",
|
|
" # start = time.perf_counter()\n",
|
|
" # try:\n",
|
|
" # P1_approx, alpha1_approx = gha1_approx(ell, P0, alpha0_ell, s, ds=ds)\n",
|
|
" # end = time.perf_counter()\n",
|
|
" # beta1_approx, lamb1_approx = ell.cart2ell(P1_approx)\n",
|
|
" # d_beta1 = abs(wu.rad2deg(beta1_approx - beta1)) / 3600\n",
|
|
" # d_lamb1 = abs(wu.rad2deg(lamb1_approx - lamb1)) / 3600\n",
|
|
" # d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600\n",
|
|
" # d_time = end - start\n",
|
|
" # example_results[f\"GHA1_approx_{ds:.3f}\"] = (d_beta1, d_lamb1, d_alpha1, d_time)\n",
|
|
" # except Exception as e:\n",
|
|
" # print(e)\n",
|
|
" # example_results[f\"GHA1_approx_{ds:.3f}\"] = (nan, nan, nan, nan)\n",
|
|
"\n",
|
|
" # for steps in steps_gha2_num:\n",
|
|
" # start = time.perf_counter()\n",
|
|
" # try:\n",
|
|
" # alpha0_num, alpha1_num_2, s_num = gha2_num(ell, beta0, lamb0, beta1, lamb1, n=steps)\n",
|
|
" # end = time.perf_counter()\n",
|
|
" # d_alpha0 = abs(wu.rad2deg(alpha0_num - alpha0_ell)) / 3600\n",
|
|
" # d_alpha1 = abs(wu.rad2deg(alpha1_num_2 - alpha1_ell)) / 3600\n",
|
|
" # d_s = abs(s_num - s) / 1000\n",
|
|
" # d_time = end - start\n",
|
|
" # example_results[f\"GHA2_num_{steps}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n",
|
|
" # print(\"konvergiert\")\n",
|
|
" # except Exception as e:\n",
|
|
" # print(e)\n",
|
|
" # print(example)\n",
|
|
" # example_results[f\"GHA2_num_{steps}\"] = (nan, nan, nan, nan)\n",
|
|
"\n",
|
|
" # for dsPart in dsPart_gha2_ES:\n",
|
|
" # ds = ell.ax/dsPart\n",
|
|
" # start = time.perf_counter()\n",
|
|
" # try:\n",
|
|
" # with suppress_print():\n",
|
|
" # alpha0_ES, alpha1_ES, s_ES = gha2_ES(ell, P0, P1, maxSegLen=ds)\n",
|
|
" # end = time.perf_counter()\n",
|
|
" # d_alpha0 = abs(wu.rad2deg(alpha0_ES - alpha0_ell)) / 3600\n",
|
|
" # d_alpha1 = abs(wu.rad2deg(alpha1_ES - alpha1_ell)) / 3600\n",
|
|
" # d_s = abs(s_ES - s) / 1000\n",
|
|
" # d_time = end - start\n",
|
|
" # example_results[f\"GHA2_ES_{ds:.3f}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n",
|
|
" # except Exception as e:\n",
|
|
" # print(e)\n",
|
|
" # example_results[f\"GHA2_ES_{ds:.3f}\"] = (nan, nan, nan, nan)\n",
|
|
" #\n",
|
|
" # for dsPart in dsPart_gha2_approx:\n",
|
|
" # ds = ell.ax/dsPart\n",
|
|
" # start = time.perf_counter()\n",
|
|
" # try:\n",
|
|
" # alpha0_approx, alpha1_approx, s_approx = gha2_approx(ell, P0, P1, ds=ds)\n",
|
|
" # end = time.perf_counter()\n",
|
|
" # d_alpha0 = abs(wu.rad2deg(alpha0_approx - alpha0_ell)) / 3600\n",
|
|
" # d_alpha1 = abs(wu.rad2deg(alpha1_approx - alpha1_ell)) / 3600\n",
|
|
" # d_s = abs(s_approx - s) / 1000\n",
|
|
" # d_time = end - start\n",
|
|
" # example_results[f\"GHA2_approx_{ds:.3f}\"] = (d_alpha0, d_alpha1, d_s, d_time)\n",
|
|
" # except Exception as e:\n",
|
|
" # print(e)\n",
|
|
" # example_results[f\"GHA2_approx_{ds:.3f}\"] = (nan, nan, nan, nan)\n",
|
|
"\n",
|
|
" results[f\"beta0: {wu.rad2deg(beta0):.3f}, lamb0: {wu.rad2deg(lamb0):.3f}, alpha0: {wu.rad2deg(alpha0_ell):.3f}, s: {s}\"] = example_results"
|
|
],
|
|
"id": "ef8849908ed3231e",
|
|
"outputs": [],
|
|
"execution_count": null
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"source": [
|
|
"# with open(f\"gha_results{test}.pkl\", \"wb\") as f:\n",
|
|
"# pickle.dump(results, f)"
|
|
],
|
|
"id": "664105ea35d50a7b",
|
|
"outputs": [],
|
|
"execution_count": null
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"source": [
|
|
"# with open(f\"gha_results{test}.pkl\", \"rb\") as f:\n",
|
|
"# results = pickle.load(f)"
|
|
],
|
|
"id": "ad8aa9dcc8af4a05",
|
|
"outputs": [],
|
|
"execution_count": null
|
|
},
|
|
{
|
|
"metadata": {},
|
|
"cell_type": "code",
|
|
"source": [
|
|
"def format_max(values, is_angle=False):\n",
|
|
" arr = np.array(values, dtype=float)\n",
|
|
" if arr.size==0:\n",
|
|
" return np.nan\n",
|
|
" maxi = np.nanmax(np.abs(arr))\n",
|
|
" if maxi is None or (isinstance(maxi,float) and (math.isnan(maxi))):\n",
|
|
" return \"nan\"\n",
|
|
" i = 2\n",
|
|
" while maxi > np.nanmean(np.abs(arr)):\n",
|
|
" maxi = np.sort(np.abs(arr[~np.isnan(arr)]))[-i]\n",
|
|
" i += 1\n",
|
|
" if is_angle:\n",
|
|
" maxi = wu.rad2deg(maxi)*3600\n",
|
|
" if f\"{maxi:.3g}\" == 0:\n",
|
|
" pass\n",
|
|
" return f\"{maxi:.3g}\"\n",
|
|
"\n",
|
|
"\n",
|
|
"def build_max_table(gha_prefix, title, group_value = None):\n",
|
|
" # metrics\n",
|
|
" if gha_prefix==\"GHA1\":\n",
|
|
" metrics = ['dBeta [\"]', 'dLambda [\"]', 'dAlpha1 [\"]', 'time [s]']\n",
|
|
" angle_mask = [True, True, True, False]\n",
|
|
" else:\n",
|
|
" metrics = ['dAlpha0 [\"]', 'dAlpha1 [\"]', 'dStrecke [m]', 'time [s]']\n",
|
|
" angle_mask = [True, True, False, False]\n",
|
|
"\n",
|
|
" # collect keys in this group\n",
|
|
" if group_value is None:\n",
|
|
" example_keys = [example_key for example_key in results.keys()]\n",
|
|
" else:\n",
|
|
" example_keys = [example_key for example_key, group_index in zip(results.keys(), table_indices) if group_index==group_value]\n",
|
|
" # all variant keys (inner dict keys) matching prefix\n",
|
|
" algorithms = sorted({algorithm for example_key in example_keys for algorithm in results[example_key].keys() if algorithm.startswith(gha_prefix)})\n",
|
|
"\n",
|
|
" header = [\"Algorithmus\", \"Parameter\"] + list(metrics)\n",
|
|
" cells = [[] for i in range(len(metrics) + 2)]\n",
|
|
" for algorithm in algorithms:\n",
|
|
" if algorithm == \"GHA1_ana_80_16\":\n",
|
|
" pass\n",
|
|
" ghaNr, variant, params = algorithm.split(\"_\", 2)\n",
|
|
" cells[0].append(variant)\n",
|
|
" cells[1].append(params)\n",
|
|
" for i, metric in enumerate(metrics):\n",
|
|
" values = []\n",
|
|
" for example_key in example_keys:\n",
|
|
" values.append(results[example_key][algorithm][i])\n",
|
|
" cells[i+2].append(format_max(values, is_angle=angle_mask[i]))\n",
|
|
"\n",
|
|
" header = dict(\n",
|
|
" values=header,\n",
|
|
" align=\"center\",\n",
|
|
" fill_color=\"lightgrey\",\n",
|
|
" font=dict(size=13)\n",
|
|
" )\n",
|
|
" cells = dict(\n",
|
|
" values=cells,\n",
|
|
" align=\"center\"\n",
|
|
" )\n",
|
|
"\n",
|
|
" fig = go.Figure(data=[go.Table(header=header, cells=cells)])\n",
|
|
" fig.update_layout(title=title,\n",
|
|
" template=\"simple_white\",\n",
|
|
" width=800,\n",
|
|
" height=280,\n",
|
|
" margin=dict(l=20, r=20, t=60, b=20))\n",
|
|
" return fig\n",
|
|
"\n",
|
|
"figs = []\n",
|
|
"if test == \"Panou\":\n",
|
|
" for table_index in sorted(set(table_indices))[:-1]:\n",
|
|
" fig1 = build_max_table(\"GHA1\", f\"{test} - Gruppe {table_index} - GHA1\", table_index)\n",
|
|
" fig2 = build_max_table(\"GHA2\", f\"{test} - Gruppe {table_index} - GHA2\", table_index)\n",
|
|
" figs.append(fig1)\n",
|
|
" figs.append(fig2)\n",
|
|
"elif test == \"Karney\":\n",
|
|
" fig1 = build_max_table(\"GHA1\", f\"{test} - GHA1\")\n",
|
|
" fig2 = build_max_table(\"GHA2\", f\"{test} - GHA2\")\n",
|
|
" figs.append(fig1)\n",
|
|
" figs.append(fig2)\n",
|
|
"\n",
|
|
"for fig in figs:\n",
|
|
" fig.show()"
|
|
],
|
|
"id": "b46d57fc0d794e28",
|
|
"outputs": [],
|
|
"execution_count": null
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python 3",
|
|
"language": "python",
|
|
"name": "python3"
|
|
},
|
|
"language_info": {
|
|
"codemirror_mode": {
|
|
"name": "ipython",
|
|
"version": 2
|
|
},
|
|
"file_extension": ".py",
|
|
"mimetype": "text/x-python",
|
|
"name": "python",
|
|
"nbconvert_exporter": "python",
|
|
"pygments_lexer": "ipython2",
|
|
"version": "2.7.6"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 5
|
|
}
|