From 6e9b13c0372e1d7eb0b8fcee7ecc62b86ca6b8ca Mon Sep 17 00:00:00 2001 From: Landung 'Don' Setiawan Date: Tue, 4 Apr 2023 12:12:09 -0700 Subject: [PATCH 01/11] Add function for computing harmonic mean --- src/seagap/configs/solver.py | 3 ++ src/seagap/harmonic_mean.py | 80 ++++++++++++++++++++++++++++++++++++ 2 files changed, 83 insertions(+) create mode 100644 src/seagap/harmonic_mean.py diff --git a/src/seagap/configs/solver.py b/src/seagap/configs/solver.py index a3b9a06..de029e9 100644 --- a/src/seagap/configs/solver.py +++ b/src/seagap/configs/solver.py @@ -148,6 +148,9 @@ class Solver(BaseModel): "Options: ship/SV3 = 0.0s, WG = 0.1s" ), ) + harmonic_mean_start_depth: float = Field( + 0.0, description="Start depth in meters for harmonic mean computation." + ) input_files: SolverInputs = Field( ..., description="Input files data path specifications." ) diff --git a/src/seagap/harmonic_mean.py b/src/seagap/harmonic_mean.py new file mode 100644 index 0000000..aee0657 --- /dev/null +++ b/src/seagap/harmonic_mean.py @@ -0,0 +1,80 @@ +import numba +import math + +@numba.njit +def _compute_hm(dd, sv, start_depth, end_depth, start_index): + """ + Computes harmonic mean. + It's a direct translation from the original Fortran code found in + src/cal_sv_harmonic_mean/get_sv_harmonic_mean.F called + subroutine `sv_harmon_mean` + """ + # TODO: Find a way to vectorize this computation + zs = start_depth + ze = end_depth + + z1=dd[start_index] + z2=dd[start_index+1] + + c_z1 = sv[start_index] + c_z2 = sv[start_index+1] + + zi = zs + if(z2 >= ze): + zf = ze + else: + zf = z2 + + cumsum = 0.0 + + for i in range(start_index+1, len(dd)): + b = ( c_z2 - c_z1) / ( z2 - z1 ) + wi = zi - z1 + c_z1/b + wf = zf - z1 + c_z1/b + + wi = math.log( (zi - z1)*b + c_z1)/b + wf = math.log( (zf - z1)*b + c_z1)/b + + delta = (wf-wi) + cumsum = cumsum + delta + z1=zf + z2=dd[i+1] + c_z1 = c_z2 + c_z2 = sv[i+1] + zi=zf + + if ze > zi and ze < z2: + zf = ze + else: + zf = z2 + + if z1 >= ze: break + + return (ze-zs)/cumsum + +def sv_harmonic_mean(svdf, start_depth, end_depth): + """ + Computes harmonic mean from a sound profile + containing depth (dd) and sound speed (sv) + + Parameters + ---------- + svdf : pd.DataFrame + Sound speed profile data as dataframe + start_depth : int or float + The start depth for harmonic mean to be computed + end_depth : int or float + The end depth for harmonic mean to be computed + + Returns + ------- + float + The sound speed harmonic mean value + """ + abs_start = abs(start_depth) + abs_end = abs(end_depth) + abs_sv = abs(svdf) + # Get the index for the start of depth closest to specified start depth + start_index = abs_sv[(abs_sv['dd'].round() >= start_depth)].index[0] + + return _compute_hm(abs_sv['dd'].values, abs_sv['sv'].values, abs_start, abs_end, start_index) From 6a34aafbf6330d11522f0a8711792498940aceb4 Mon Sep 17 00:00:00 2001 From: Landung 'Don' Setiawan Date: Thu, 6 Apr 2023 12:50:28 -0700 Subject: [PATCH 02/11] Add function to compute harmonic mean from sound speed profile --- src/seagap/configs/solver.py | 12 +++++------ src/seagap/harmonic_mean.py | 42 ++++++++++++++++++++++++++++++++---- 2 files changed, 44 insertions(+), 10 deletions(-) diff --git a/src/seagap/configs/solver.py b/src/seagap/configs/solver.py index de029e9..376dfda 100644 --- a/src/seagap/configs/solver.py +++ b/src/seagap/configs/solver.py @@ -141,10 +141,10 @@ class Solver(BaseModel): travel_times_variance: float = Field( 1e-10, description="VARIANCE (s**2) PXP two-way travel time measurement" ) - transponder_wait_time: float = Field( + transducer_delay_time: float = Field( 0.0, description=( - "Transponder Wait Time - delay at surface transducer (secs.). " + "Transducer Delay Time - delay at surface transducer (secs.). " "Options: ship/SV3 = 0.0s, WG = 0.1s" ), ) @@ -171,9 +171,9 @@ class Solver(BaseModel): ), ) - @validator("transponder_wait_time") - def check_transponder_wait_time(cls, v): - """Validate transponder wait time value""" + @validator("transducer_delay_time") + def check_transducer_delay_time(cls, v): + """Validate transducer wait time value""" if v not in [0.0, 0.1]: - raise ValueError("Transponder wait time must either be 0.0s or 0.1s") + raise ValueError("Transducer wait time must either be 0.0s or 0.1s") return v diff --git a/src/seagap/harmonic_mean.py b/src/seagap/harmonic_mean.py index aee0657..f9ec5e6 100644 --- a/src/seagap/harmonic_mean.py +++ b/src/seagap/harmonic_mean.py @@ -1,8 +1,17 @@ +import pandas as pd +import numpy as np +from typing import Union import numba import math @numba.njit -def _compute_hm(dd, sv, start_depth, end_depth, start_index): +def _compute_hm( + dd: np.ndarray, + sv: np.ndarray, + start_depth: Union[int, float], + end_depth: Union[int, float], + start_index: int +): """ Computes harmonic mean. It's a direct translation from the original Fortran code found in @@ -10,24 +19,42 @@ def _compute_hm(dd, sv, start_depth, end_depth, start_index): subroutine `sv_harmon_mean` """ # TODO: Find a way to vectorize this computation + + # Assign start depth and end depth to zs and ze zs = start_depth ze = end_depth + # Extract the first and second depth values z1=dd[start_index] z2=dd[start_index+1] + # Extract the first and second sound speed values c_z1 = sv[start_index] c_z2 = sv[start_index+1] + # Set start depth to initial depth to keep track of it zi = zs - if(z2 >= ze): + + # If the second depth value (z2) is greater than and equal to the + # end depth value (ze) then the final depth (zf) should be assigned + # to the end depth value (ze), otherwise, the final depth + # should be assigned to the second depth values + if z2 >= ze: zf = ze else: zf = z2 - + + # Start cumulative sum as 0.0 cumsum = 0.0 + # Loop over the whole array for depth and sound speed, + # starting from the index of the second value to the whole + # depth data, which is assumed to be the same exact number + # as the sound speed data for i in range(start_index+1, len(dd)): + # calculate the slope of the two points + # for sound speed and depth + # slope = (sv2 - sv1) / (d2 - d1) b = ( c_z2 - c_z1) / ( z2 - z1 ) wi = zi - z1 + c_z1/b wf = zf - z1 + c_z1/b @@ -52,7 +79,11 @@ def _compute_hm(dd, sv, start_depth, end_depth, start_index): return (ze-zs)/cumsum -def sv_harmonic_mean(svdf, start_depth, end_depth): +def sv_harmonic_mean( + svdf: pd.DataFrame, + start_depth: Union[int, float], + end_depth: Union[int, float] +): """ Computes harmonic mean from a sound profile containing depth (dd) and sound speed (sv) @@ -71,6 +102,9 @@ def sv_harmonic_mean(svdf, start_depth, end_depth): float The sound speed harmonic mean value """ + # Clean up the sound speed value, ensuring that there's no negative value + svdf = svdf[svdf['sv'] < 0] + # Make all of the values absolute values, so we're only dealing with positives abs_start = abs(start_depth) abs_end = abs(end_depth) abs_sv = abs(svdf) From b56c1577b5f7ce013ba0e21d3822a85caf42daff Mon Sep 17 00:00:00 2001 From: Landung 'Don' Setiawan Date: Thu, 6 Apr 2023 14:02:50 -0700 Subject: [PATCH 03/11] Add harmonic mean tests --- tests/data/sound_profile_hm.csv | 985 ++++++++++++++++++++++++++++++++ tests/test_harmonic_mean.py | 50 ++ 2 files changed, 1035 insertions(+) create mode 100644 tests/data/sound_profile_hm.csv create mode 100644 tests/test_harmonic_mean.py diff --git a/tests/data/sound_profile_hm.csv b/tests/data/sound_profile_hm.csv new file mode 100644 index 0000000..8d24dd7 --- /dev/null +++ b/tests/data/sound_profile_hm.csv @@ -0,0 +1,985 @@ +dd,sv +-1.984,1502.5 +-2.976,1502.54 +-3.967,1502.57 +-4.959,1502.51 +-5.951,1502.49 +-6.943,1502.36 +-7.934999999999999,1502.3 +-8.926999999999998,1502.22 +-9.918,1501.84 +-10.91,1501.55 +-11.902,1501.48 +-12.894,1501.07 +-13.886,1500.89 +-14.877,1500.73 +-15.869,1500.65 +-16.861,1500.7 +-17.853,1500.38 +-18.844,1500.09 +-19.836,1499.68 +-20.828,1499.01 +-21.82,1498.3 +-22.811,1497.3 +-23.803,1496.32 +-24.795,1495.76 +-25.787,1495.74 +-26.778,1495.12 +-27.77,1494.3 +-28.762,1494.18 +-29.753,1493.57 +-30.745,1492.59 +-31.737,1491.65 +-32.729,1491.14 +-33.72,1491.03 +-34.712,1490.46 +-35.704,1490.02 +-36.695,1489.85 +-37.687,1489.87 +-38.679,1489.93 +-39.67,1489.85 +-40.662,1489.44 +-41.654,1489.38 +-42.645,1489.22 +-43.637,1489.13 +-44.629,1489.1 +-45.62,1489.02 +-46.612,1488.98 +-47.604,1488.96 +-48.595,1488.92 +-49.587,1488.74 +-50.578,1488.71 +-51.57,1488.68 +-52.562,1488.58 +-53.553,1488.48 +-54.545,1488.42 +-56.528,1488.43 +-57.52,1488.3 +-58.511,1488.02 +-59.503,1487.89 +-60.494,1487.77 +-61.486,1487.74 +-62.477,1487.53 +-63.469,1487.45 +-64.46,1487.43 +-65.452,1487.29 +-66.444,1487.12 +-67.435,1486.91 +-68.427,1486.8 +-69.418,1486.73 +-70.41,1486.64 +-71.401,1486.62 +-72.393,1486.55 +-75.367,1486.6 +-76.35899999999998,1486.65 +-77.34999999999998,1486.73 +-78.342,1486.79 +-79.333,1486.78 +-80.325,1486.79 +-81.316,1486.81 +-82.307,1486.85 +-83.299,1486.91 +-84.29,1486.94 +-85.28199999999998,1486.92 +-86.27299999999998,1486.94 +-87.265,1486.96 +-89.248,1486.95 +-90.239,1486.94 +-91.23,1486.95 +-92.222,1486.94 +-93.213,1486.83 +-94.205,1486.77 +-95.196,1486.74 +-96.187,1486.75 +-98.17,1486.76 +-99.161,1486.73 +-100.153,1486.68 +-101.144,1486.64 +-102.136,1486.61 +-103.127,1486.47 +-104.118,1486.31 +-105.11,1486.15 +-106.101,1486.11 +-107.092,1486.09 +-108.084,1485.98 +-109.075,1486.07 +-110.066,1485.97 +-111.058,1485.91 +-112.049,1485.92 +-114.032,1485.93 +-116.014,1485.9 +-117.005,1485.91 +-118.988,1485.82 +-119.979,1485.8 +-121.962,1485.79 +-125.927,1485.8 +-126.918,1485.66 +-127.909,1485.63 +-128.901,1485.64 +-129.892,1485.57 +-130.883,1485.52 +-131.874,1485.51 +-132.865,1485.45 +-133.857,1485.39 +-134.848,1485.28 +-135.839,1485.24 +-137.821,1485.2 +-138.813,1485.18 +-139.804,1485.17 +-140.795,1485.16 +-141.786,1485.11 +-142.777,1485.09 +-143.769,1485.1 +-144.76,1485.11 +-145.751,1485.12 +-147.733,1485.07 +-148.724,1485.01 +-149.715,1484.97 +-150.707,1484.91 +-151.698,1484.89 +-153.68,1484.82 +-154.671,1484.8 +-155.662,1484.75 +-156.653,1484.64 +-157.644,1484.62 +-158.635,1484.59 +-159.626,1484.46 +-160.618,1484.43 +-162.6,1484.41 +-163.591,1484.34 +-164.582,1484.29 +-165.573,1484.26 +-166.564,1484.24 +-167.555,1484.14 +-168.546,1484.05 +-171.519,1483.95 +-172.51,1483.89 +-173.501,1483.88 +-175.483,1483.78 +-176.474,1483.74 +-179.447,1483.72 +-180.438,1483.7 +-181.429,1483.61 +-182.42,1483.5 +-183.411,1483.39 +-184.402,1483.37 +-187.375,1483.38 +-190.348,1483.37 +-191.339,1483.26 +-192.33,1483.21 +-193.32,1483.18 +-194.311,1483.16 +-195.302,1483.13 +-196.293,1483.11 +-197.284,1483.0 +-198.275,1482.95 +-199.266,1482.86 +-200.257,1482.84 +-201.248,1482.8 +-202.239,1482.79 +-203.229,1482.76 +-204.22,1482.73 +-205.211,1482.79 +-206.202,1482.74 +-207.193,1482.72 +-208.184,1482.69 +-209.175,1482.68 +-211.156,1482.67 +-212.147,1482.68 +-213.138,1482.7 +-214.129,1482.69 +-215.12,1482.62 +-216.11,1482.6 +-217.101,1482.59 +-218.092,1482.57 +-219.083,1482.56 +-220.074,1482.57 +-221.064,1482.55 +-222.055,1482.47 +-223.046,1482.43 +-224.037,1482.4 +-225.027,1482.34 +-226.018,1482.28 +-227.009,1482.24 +-228.0,1482.21 +-228.99,1482.15 +-229.981,1482.14 +-230.972,1482.05 +-232.953,1482.07 +-233.944,1482.05 +-234.935,1481.96 +-235.926,1481.9 +-236.916,1481.91 +-237.907,1481.89 +-239.888,1481.87 +-240.879,1481.86 +-242.86,1481.85 +-243.851,1481.82 +-244.842,1481.78 +-245.832,1481.74 +-246.823,1481.68 +-247.814,1481.53 +-248.804,1481.47 +-249.795,1481.49 +-250.786,1481.55 +-252.767,1481.54 +-253.758,1481.43 +-254.748,1481.3 +-256.73,1481.32 +-258.711,1481.28 +-259.701,1481.27 +-261.683,1481.28 +-263.664,1481.3 +-264.654,1481.17 +-265.645,1481.08 +-266.635,1481.07 +-267.626,1481.04 +-268.617,1481.03 +-270.598,1481.01 +-271.588,1480.91 +-272.579,1480.93 +-273.569,1480.96 +-274.56,1480.99 +-275.55,1481.01 +-276.541,1480.99 +-277.531,1480.98 +-278.522,1480.96 +-279.512,1480.85 +-280.503,1480.71 +-281.493,1480.58 +-282.484,1480.54 +-283.474,1480.55 +-284.465,1480.54 +-285.455,1480.55 +-286.446,1480.57 +-287.436,1480.59 +-288.427,1480.52 +-289.417,1480.49 +-290.408,1480.46 +-292.389,1480.42 +-293.379,1480.34 +-294.369,1480.31 +-295.36,1480.32 +-297.341,1480.28 +-298.331,1480.24 +-299.322,1480.22 +-301.302,1480.24 +-302.293,1480.25 +-304.273,1480.24 +-306.254,1480.26 +-308.235,1480.22 +-310.216,1480.2 +-312.196,1480.19 +-313.187,1480.18 +-314.177,1480.12 +-315.167,1480.07 +-316.158,1480.05 +-317.148,1480.03 +-318.138,1480.01 +-319.129,1479.99 +-320.119,1479.97 +-321.109,1479.98 +-322.1,1479.99 +-323.09,1480.01 +-326.061,1480.0 +-327.051,1479.99 +-328.041,1479.95 +-329.032,1479.93 +-331.012,1479.94 +-332.002,1479.95 +-332.993,1479.96 +-333.983,1479.95 +-335.963,1479.94 +-337.944,1479.95 +-340.914,1479.96 +-342.895,1479.95 +-343.885,1479.96 +-344.875,1479.95 +-345.865,1479.93 +-347.846,1479.92 +-348.836,1479.89 +-349.826,1479.87 +-350.816,1479.85 +-351.806,1479.82 +-354.777,1479.78 +-355.767,1479.77 +-356.757,1479.72 +-357.747,1479.73 +-358.737,1479.71 +-359.728,1479.66 +-360.718,1479.65 +-361.708,1479.64 +-362.698,1479.62 +-363.688,1479.61 +-364.678,1479.62 +-365.668,1479.63 +-366.658,1479.64 +-367.648,1479.65 +-368.638,1479.64 +-369.628,1479.68 +-374.579,1479.57 +-375.569,1479.44 +-376.559,1479.41 +-377.549,1479.39 +-378.539,1479.42 +-379.529,1479.35 +-380.519,1479.32 +-381.509,1479.31 +-382.499,1479.33 +-383.489,1479.34 +-384.479,1479.3 +-385.469,1479.29 +-386.459,1479.25 +-388.439,1479.28 +-389.429,1479.31 +-390.419,1479.35 +-391.409,1479.38 +-392.399,1479.36 +-393.389,1479.34 +-394.379,1479.35 +-395.369,1479.36 +-396.359,1479.37 +-397.349,1479.39 +-400.318,1479.37 +-401.308,1479.39 +-403.288,1479.3 +-404.278,1479.26 +-406.258,1479.28 +-407.248,1479.3 +-408.238,1479.32 +-409.227,1479.3 +-410.217,1479.28 +-411.207,1479.26 +-412.197,1479.17 +-413.187,1479.13 +-414.177,1479.09 +-415.167,1479.05 +-417.146,1479.03 +-418.136,1479.0 +-419.126,1478.93 +-420.116,1478.95 +-421.106,1478.96 +-422.096,1478.99 +-423.085,1479.11 +-424.075,1479.07 +-425.065,1479.04 +-426.055,1479.08 +-427.045,1479.11 +-429.024,1479.12 +-430.014,1479.09 +-431.004,1479.05 +-431.993,1478.99 +-432.983,1478.91 +-433.973,1478.94 +-434.963,1478.89 +-435.953,1478.86 +-436.942,1478.88 +-437.932,1478.9 +-438.922,1478.93 +-440.901,1478.91 +-441.891,1478.92 +-442.881,1478.93 +-443.87,1478.92 +-445.85,1478.93 +-446.839,1478.94 +-447.829,1478.95 +-448.819,1478.97 +-452.778,1478.95 +-454.757,1478.96 +-455.747,1478.97 +-457.726,1478.95 +-458.716,1478.96 +-460.695,1478.97 +-461.684,1478.98 +-462.674,1478.95 +-463.664,1478.94 +-465.643,1478.93 +-468.612,1478.91 +-470.591,1478.93 +-471.58,1478.92 +-473.56,1478.93 +-474.549,1478.92 +-475.539,1478.88 +-477.518,1478.86 +-478.507,1478.85 +-479.497,1478.81 +-480.487,1478.8 +-482.466,1478.81 +-483.455,1478.82 +-484.445,1478.83 +-485.434,1478.84 +-486.424,1478.83 +-488.403,1478.82 +-489.392,1478.8 +-491.371,1478.78 +-492.361,1478.75 +-493.35,1478.73 +-494.34,1478.72 +-496.319,1478.73 +-497.308,1478.74 +-498.297,1478.75 +-499.287,1478.76 +-501.266,1478.77 +-502.255,1478.76 +-503.245,1478.78 +-505.223,1478.76 +-507.202,1478.78 +-509.181,1478.8 +-511.16,1478.82 +-512.149,1478.83 +-513.139,1478.84 +-515.117,1478.83 +-516.107,1478.8 +-519.075,1478.78 +-521.054,1478.77 +-522.043,1478.8 +-523.032,1478.75 +-526.99,1478.74 +-527.979,1478.73 +-528.968,1478.72 +-529.957,1478.7 +-531.936,1478.72 +-532.925,1478.73 +-533.915,1478.74 +-536.882,1478.75 +-538.861,1478.76 +-543.807,1478.75 +-544.796,1478.76 +-546.775,1478.77 +-548.753,1478.76 +-550.732,1478.74 +-551.721,1478.73 +-552.71,1478.74 +-553.699,1478.75 +-554.689,1478.76 +-555.678,1478.78 +-556.667,1478.79 +-557.656,1478.81 +-558.645,1478.82 +-559.634,1478.83 +-560.624,1478.8 +-561.613,1478.78 +-564.58,1478.8 +-567.548,1478.81 +-568.537,1478.82 +-569.526,1478.83 +-570.515,1478.82 +-572.493,1478.83 +-573.482,1478.8 +-574.471,1478.73 +-575.46,1478.72 +-577.438,1478.69 +-578.428,1478.66 +-580.406,1478.67 +-581.395,1478.68 +-584.362,1478.7 +-585.351,1478.71 +-586.34,1478.72 +-587.329,1478.73 +-588.318,1478.74 +-589.307,1478.76 +-590.296,1478.75 +-591.285,1478.76 +-592.274,1478.77 +-593.263,1478.78 +-594.252,1478.79 +-595.241,1478.8 +-596.23,1478.82 +-598.208,1478.8 +-599.197,1478.77 +-600.186,1478.76 +-601.175,1478.78 +-602.164,1478.8 +-603.153,1478.81 +-605.131,1478.82 +-606.12,1478.83 +-607.109,1478.85 +-608.097,1478.86 +-610.075,1478.87 +-613.042,1478.88 +-614.031,1478.89 +-615.02,1478.91 +-616.998,1478.89 +-617.987,1478.86 +-619.964,1478.88 +-621.942,1478.89 +-622.931,1478.91 +-623.92,1478.92 +-624.909,1478.91 +-625.897,1478.92 +-626.886,1478.93 +-628.864,1478.92 +-629.853,1478.9 +-630.842,1478.86 +-632.819,1478.87 +-633.808,1478.89 +-634.797,1478.9 +-635.786,1478.92 +-636.774,1478.93 +-637.763,1478.95 +-638.752,1478.97 +-639.741,1478.98 +-640.73,1478.99 +-641.718,1479.0 +-642.707,1478.99 +-643.696,1478.98 +-644.685,1478.97 +-647.651,1478.99 +-648.64,1479.0 +-649.628,1479.01 +-650.617,1479.02 +-651.606,1479.01 +-653.583,1479.0 +-654.572,1479.01 +-658.527,1479.03 +-660.504,1479.04 +-661.493,1479.05 +-662.481,1479.06 +-664.459,1479.05 +-666.436,1479.06 +-667.425,1479.05 +-668.413,1479.06 +-669.402,1479.07 +-670.391,1479.09 +-672.368,1479.11 +-673.356,1479.09 +-674.345,1479.1 +-675.334,1479.11 +-676.322,1479.12 +-677.311,1479.14 +-679.288,1479.15 +-680.277,1479.16 +-681.265,1479.17 +-682.254,1479.18 +-683.242,1479.19 +-684.231,1479.2 +-685.219,1479.22 +-686.208,1479.23 +-687.197,1479.24 +-691.151,1479.23 +-692.139,1479.24 +-693.128,1479.26 +-694.116,1479.27 +-695.105,1479.28 +-696.093,1479.3 +-697.082,1479.31 +-698.07,1479.32 +-699.059,1479.34 +-700.047,1479.35 +-701.036,1479.36 +-702.024,1479.38 +-704.001,1479.39 +-705.978,1479.41 +-706.967,1479.42 +-707.955,1479.4 +-708.943,1479.39 +-709.932,1479.41 +-711.909,1479.42 +-712.897,1479.44 +-713.886,1479.45 +-715.863,1479.44 +-716.851,1479.43 +-717.839,1479.45 +-718.828,1479.46 +-719.816,1479.47 +-720.8049999999998,1479.48 +-721.793,1479.49 +-722.7809999999998,1479.5 +-723.77,1479.51 +-724.758,1479.49 +-726.735,1479.5 +-727.7229999999998,1479.51 +-728.711,1479.52 +-729.7,1479.53 +-730.688,1479.54 +-731.677,1479.55 +-732.6649999999998,1479.57 +-734.6409999999998,1479.58 +-735.63,1479.59 +-736.618,1479.61 +-738.595,1479.63 +-739.5829999999999,1479.64 +-741.5599999999998,1479.66 +-742.548,1479.67 +-744.524,1479.66 +-746.501,1479.64 +-747.489,1479.61 +-748.4779999999998,1479.62 +-753.419,1479.63 +-754.407,1479.64 +-755.395,1479.66 +-756.383,1479.65 +-757.3719999999998,1479.66 +-759.3479999999998,1479.68 +-762.313,1479.69 +-763.301,1479.7 +-764.289,1479.71 +-765.277,1479.72 +-766.265,1479.73 +-767.254,1479.74 +-769.23,1479.75 +-771.206,1479.76 +-774.171,1479.74 +-778.123,1479.75 +-779.111,1479.76 +-780.099,1479.78 +-781.087,1479.79 +-782.076,1479.8 +-783.0639999999999,1479.81 +-785.0399999999998,1479.82 +-786.028,1479.83 +-787.0159999999998,1479.84 +-788.004,1479.85 +-788.9919999999998,1479.86 +-789.98,1479.87 +-790.9679999999998,1479.88 +-791.956,1479.87 +-792.9439999999998,1479.88 +-793.932,1479.89 +-794.9199999999998,1479.9 +-795.908,1479.91 +-796.897,1479.93 +-797.885,1479.94 +-798.873,1479.95 +-800.849,1479.96 +-801.837,1479.97 +-802.825,1479.98 +-803.813,1479.99 +-804.801,1480.0 +-806.777,1480.01 +-807.765,1480.03 +-809.74,1480.04 +-810.7279999999998,1480.02 +-812.7039999999998,1480.03 +-813.692,1480.04 +-814.6799999999998,1480.05 +-815.668,1480.04 +-816.6559999999998,1480.05 +-817.644,1480.07 +-818.6319999999998,1480.08 +-819.62,1480.09 +-820.6079999999998,1480.1 +-821.596,1480.11 +-823.572,1480.12 +-824.5589999999999,1480.13 +-825.547,1480.14 +-826.5349999999999,1480.16 +-827.523,1480.17 +-828.5109999999999,1480.18 +-831.475,1480.19 +-832.462,1480.2 +-833.45,1480.22 +-835.426,1480.21 +-836.414,1480.22 +-837.402,1480.23 +-838.389,1480.24 +-839.3769999999998,1480.23 +-841.3529999999998,1480.24 +-842.341,1480.25 +-843.3289999999998,1480.26 +-844.316,1480.28 +-845.3039999999999,1480.29 +-846.292,1480.3 +-847.2799999999999,1480.32 +-849.255,1480.33 +-850.243,1480.34 +-851.231,1480.35 +-852.219,1480.36 +-853.206,1480.37 +-854.1939999999998,1480.39 +-855.182,1480.41 +-856.1699999999998,1480.42 +-857.157,1480.43 +-859.133,1480.44 +-860.121,1480.45 +-861.1079999999998,1480.46 +-862.096,1480.45 +-863.0839999999998,1480.46 +-864.071,1480.47 +-868.022,1480.48 +-869.01,1480.49 +-871.9729999999998,1480.5 +-872.961,1480.51 +-873.948,1480.52 +-875.924,1480.53 +-879.874,1480.54 +-880.8619999999999,1480.55 +-881.849,1480.56 +-882.837,1480.57 +-883.825,1480.59 +-884.812,1480.6 +-885.7999999999998,1480.61 +-887.775,1480.62 +-888.763,1480.63 +-889.75,1480.64 +-894.688,1480.63 +-896.663,1480.64 +-897.6509999999998,1480.65 +-898.638,1480.66 +-899.626,1480.65 +-900.613,1480.62 +-901.601,1480.59 +-902.589,1480.6 +-903.576,1480.62 +-905.551,1480.64 +-906.5390000000002,1480.66 +-908.514,1480.67 +-911.476,1480.69 +-912.464,1480.7 +-913.451,1480.72 +-914.4380000000002,1480.73 +-915.426,1480.74 +-916.413,1480.75 +-917.401,1480.76 +-919.376,1480.75 +-921.351,1480.74 +-922.338,1480.75 +-923.325,1480.76 +-924.3130000000002,1480.77 +-925.3,1480.78 +-926.288,1480.79 +-928.263,1480.8 +-929.25,1480.81 +-930.237,1480.82 +-932.2120000000002,1480.83 +-933.2,1480.84 +-935.174,1480.86 +-936.162,1480.88 +-937.149,1480.89 +-938.136,1480.9 +-939.124,1480.91 +-940.1110000000001,1480.93 +-942.086,1480.91 +-944.06,1480.92 +-946.035,1480.93 +-947.022,1480.94 +-948.0100000000001,1480.95 +-948.997,1480.97 +-949.984,1480.98 +-950.972,1481.0 +-952.946,1480.99 +-958.87,1481.01 +-959.857,1481.03 +-960.844,1481.0 +-963.806,1481.01 +-964.793,1481.02 +-965.781,1481.03 +-966.768,1481.05 +-967.755,1481.06 +-968.742,1481.01 +-969.729,1481.02 +-970.717,1481.04 +-971.704,1481.06 +-972.691,1481.07 +-973.678,1481.09 +-974.665,1481.11 +-975.653,1481.12 +-976.6400000000002,1481.14 +-977.627,1481.15 +-978.614,1481.17 +-979.601,1481.18 +-980.588,1481.15 +-981.576,1481.16 +-983.55,1481.18 +-984.537,1481.19 +-985.524,1481.2 +-986.511,1481.22 +-987.498,1481.24 +-989.473,1481.26 +-990.46,1481.27 +-991.447,1481.29 +-992.434,1481.3 +-993.421,1481.31 +-994.408,1481.32 +-995.395,1481.34 +-997.369,1481.35 +-998.3560000000001,1481.36 +-999.343,1481.37 +-1000.33,1481.38 +-1002.3,1481.39 +-1003.29,1481.4 +-1004.28,1481.41 +-1005.27,1481.42 +-1006.25,1481.43 +-1007.24,1481.44 +-1008.23,1481.45 +-1009.21,1481.46 +-1012.17,1481.47 +-1013.16,1481.49 +-1018.1,1481.5 +-1019.08,1481.51 +-1020.07,1481.53 +-1021.06,1481.55 +-1022.04,1481.56 +-1023.03,1481.57 +-1024.02,1481.59 +-1025.01,1481.61 +-1025.99,1481.62 +-1026.98,1481.64 +-1027.97,1481.65 +-1028.95,1481.66 +-1029.94,1481.67 +-1030.93,1481.69 +-1031.91,1481.7 +-1032.9,1481.72 +-1033.89,1481.73 +-1034.88,1481.75 +-1035.86,1481.76 +-1036.85,1481.78 +-1037.84,1481.8 +-1038.82,1481.81 +-1039.81,1481.82 +-1040.8,1481.84 +-1042.77,1481.85 +-1043.76,1481.86 +-1044.74,1481.87 +-1045.73,1481.89 +-1047.7,1481.91 +-1048.69,1481.92 +-1049.68,1481.93 +-1050.66,1481.95 +-1051.65,1481.96 +-1052.64,1481.98 +-1053.62,1481.99 +-1054.61,1482.01 +-1055.6,1482.03 +-1056.59,1482.04 +-1057.57,1482.06 +-1058.56,1482.07 +-1059.55,1482.05 +-1060.53,1482.06 +-1061.52,1482.07 +-1062.51,1482.08 +-1063.49,1482.09 +-1065.47,1482.1 +-1066.45,1482.11 +-1067.44,1482.1 +-1069.41,1482.11 +-1070.4,1482.12 +-1072.37,1482.1 +-1073.36,1482.07 +-1075.33,1482.08 +-1076.32,1482.1 +-1077.31,1482.12 +-1078.29,1482.14 +-1079.28,1482.15 +-1080.27,1482.17 +-1081.25,1482.18 +-1082.24,1482.19 +-1083.23,1482.2 +-1084.21,1482.21 +-1085.2,1482.23 +-1086.19,1482.24 +-1087.17,1482.25 +-1088.16,1482.26 +-1089.15,1482.28 +-1090.13,1482.3 +-1091.12,1482.31 +-1092.11,1482.33 +-1093.09,1482.34 +-1094.08,1482.36 +-1095.07,1482.37 +-1096.05,1482.39 +-1100.0,1482.38 +-1101.97,1482.39 +-1102.96,1482.41 +-1103.95,1482.42 +-1104.93,1482.43 +-1105.92,1482.44 +-1106.91,1482.45 +-1107.89,1482.46 +-1108.88,1482.47 +-1109.87,1482.48 +-1110.85,1482.49 +-1111.84,1482.5 +-1113.81,1482.51 +-1114.8,1482.49 +-1116.77,1482.47 +-1117.76,1482.41 +-1119.73,1482.39 +-1120.72,1482.43 +-1121.7,1482.41 +-1123.68,1482.42 +-1124.66,1482.44 +-1125.65,1482.42 +-1126.64,1482.44 +-1127.62,1482.45 +-1128.61,1482.46 +-1129.6,1482.47 +-1130.58,1482.46 +-1131.57,1482.47 +-1132.56,1482.48 +-1133.54,1482.5 +-1134.53,1482.51 +-1137.49,1482.53 +-1138.47,1482.55 +-1140.45,1482.56 +-1141.43,1482.58 +-1142.42,1482.59 +-1143.4,1482.61 +-1145.38,1482.59 +-1146.37,1482.6 +-1147.35,1482.58 +-1148.34,1482.59 +-1149.32,1482.6 +-1150.31,1482.61 +-1151.3,1482.55 +-1152.28,1482.5 +-1153.27,1482.51 +-1154.26,1482.52 +-1155.24,1482.54 +-1156.23,1482.56 +-1157.21,1482.58 +-1158.2,1482.6 +-1159.19,1482.61 +-1160.17,1482.63 +-1161.16,1482.64 +-1162.15,1482.66 +-1163.13,1482.68 +-1164.12,1482.7 +-1165.11,1482.71 +-1166.09,1482.72 +-1168.06,1482.74 +-1169.05,1482.76 +-1170.04,1482.78 +-1171.02,1482.8 +-1172.99,1482.82 +-1175.95,1482.84 +-1177.93,1482.86 +-1178.91,1482.88 +-1179.9,1482.89 +-1180.88,1482.91 +-1182.86,1482.92 +-1183.84,1482.94 +-1184.83,1482.97 +-1186.8,1482.99 +-1187.79,1483.01 +-1188.77,1483.03 +-1189.76,1483.04 +-1190.75,1483.06 +-1191.73,1483.08 +-1193.7,1483.09 +-1194.69,1483.11 +-1197.65,1483.13 +-1198.64,1483.15 +-1199.62,1483.16 +-1200.61,1483.18 +-1202.58,1483.21 +-1203.57,1483.23 +-1204.55,1483.24 +-1205.54,1483.26 +-1206.53,1483.28 +-1207.51,1483.29 +-1208.5,1483.31 +-1209.48,1483.33 +-1210.47,1483.35 +-1211.45,1483.36 +-1212.44,1483.39 +-1213.43,1483.4 +-1214.41,1483.42 +-1215.4,1483.44 +-1216.38,1483.46 +-1217.37,1483.47 +-1218.36,1483.49 +-1219.34,1483.5 +-1220.33,1483.51 diff --git a/tests/test_harmonic_mean.py b/tests/test_harmonic_mean.py new file mode 100644 index 0000000..a18dc23 --- /dev/null +++ b/tests/test_harmonic_mean.py @@ -0,0 +1,50 @@ +import pandas as pd +import pytest +import numpy as np + +from . import TEST_DATA_FOLDER +from seagap.harmonic_mean import sv_harmonic_mean, _compute_hm + +@pytest.fixture() +def sound_profile_data() -> pd.DataFrame: + sv_file = TEST_DATA_FOLDER / 'sound_profile_hm.csv' + return pd.read_csv(sv_file) + + +@pytest.mark.parametrize( + "end_depth,expected_hm", + [ + (-1176.5866, 1481.551), + (-1146.5881, 1481.521), + (-1133.7305, 1481.509) + ] +) +def test_sv_harmonic_mean(end_depth, expected_hm, sound_profile_data): + svdf = sound_profile_data + start_depth = -4 + harmonic_mean = round(sv_harmonic_mean(svdf, start_depth, end_depth), 3) + + assert harmonic_mean == expected_hm + +@pytest.mark.parametrize( + "test_idx,expected_hm", + [ + (3, 1501.69), + (6, 1501.225) + ] +) +def test__compute_hm(test_idx, expected_hm): + dd = np.arange(7) * 10 + sv = np.arange(1502, 1500, step=-0.31) + + # Get partial of the data from test index + part_dd = dd[:test_idx] + part_sv = sv[:test_idx] + start_depth, end_depth = part_dd[0], part_dd[-1] + result_hm = round(_compute_hm(dd, sv, start_depth, end_depth, 0), 3) + + # Check for result to match expected + assert result_hm == expected_hm + + # Result should be same as regular mean as the data is linear here + assert result_hm == round(sum(part_sv) / test_idx, 3) From fe67dea8bc721360a4c6f7812aa311094c53046a Mon Sep 17 00:00:00 2001 From: Landung 'Don' Setiawan Date: Thu, 6 Apr 2023 14:07:22 -0700 Subject: [PATCH 04/11] Fix lint and add test folder path --- src/seagap/harmonic_mean.py | 76 +++++++++++++++++++++---------------- tests/__init__.py | 5 +++ 2 files changed, 48 insertions(+), 33 deletions(-) diff --git a/src/seagap/harmonic_mean.py b/src/seagap/harmonic_mean.py index f9ec5e6..b887618 100644 --- a/src/seagap/harmonic_mean.py +++ b/src/seagap/harmonic_mean.py @@ -1,8 +1,10 @@ -import pandas as pd -import numpy as np +import math from typing import Union + import numba -import math +import numpy as np +import pandas as pd + @numba.njit def _compute_hm( @@ -10,7 +12,7 @@ def _compute_hm( sv: np.ndarray, start_depth: Union[int, float], end_depth: Union[int, float], - start_index: int + start_index: int, ): """ Computes harmonic mean. @@ -19,23 +21,23 @@ def _compute_hm( subroutine `sv_harmon_mean` """ # TODO: Find a way to vectorize this computation - + # Assign start depth and end depth to zs and ze zs = start_depth ze = end_depth - + # Extract the first and second depth values - z1=dd[start_index] - z2=dd[start_index+1] + z1 = dd[start_index] + z2 = dd[start_index + 1] # Extract the first and second sound speed values c_z1 = sv[start_index] - c_z2 = sv[start_index+1] + c_z2 = sv[start_index + 1] # Set start depth to initial depth to keep track of it zi = zs - # If the second depth value (z2) is greater than and equal to the + # If the second depth value (z2) is greater than and equal to the # end depth value (ze) then the final depth (zf) should be assigned # to the end depth value (ze), otherwise, the final depth # should be assigned to the second depth values @@ -43,51 +45,54 @@ def _compute_hm( zf = ze else: zf = z2 - + # Start cumulative sum as 0.0 cumsum = 0.0 - + # Loop over the whole array for depth and sound speed, # starting from the index of the second value to the whole # depth data, which is assumed to be the same exact number # as the sound speed data - for i in range(start_index+1, len(dd)): + for i in range(start_index + 1, len(dd)): # calculate the slope of the two points # for sound speed and depth # slope = (sv2 - sv1) / (d2 - d1) - b = ( c_z2 - c_z1) / ( z2 - z1 ) - wi = zi - z1 + c_z1/b - wf = zf - z1 + c_z1/b + b = (c_z2 - c_z1) / (z2 - z1) + wi = zi - z1 + c_z1 / b + wf = zf - z1 + c_z1 / b - wi = math.log( (zi - z1)*b + c_z1)/b - wf = math.log( (zf - z1)*b + c_z1)/b + wi = math.log((zi - z1) * b + c_z1) / b + wf = math.log((zf - z1) * b + c_z1) / b - delta = (wf-wi) + delta = wf - wi cumsum = cumsum + delta - z1=zf - z2=dd[i+1] + z1 = zf + z2 = dd[i + 1] c_z1 = c_z2 - c_z2 = sv[i+1] - zi=zf + c_z2 = sv[i + 1] + zi = zf if ze > zi and ze < z2: zf = ze else: zf = z2 - if z1 >= ze: break + if z1 >= ze: + break - return (ze-zs)/cumsum + if cumsum == 0: + # If cumulative sum is 0, most likely only one value + return sv[start_index] + return (ze - zs) / cumsum + def sv_harmonic_mean( - svdf: pd.DataFrame, - start_depth: Union[int, float], - end_depth: Union[int, float] + svdf: pd.DataFrame, start_depth: Union[int, float], end_depth: Union[int, float] ): """ Computes harmonic mean from a sound profile containing depth (dd) and sound speed (sv) - + Parameters ---------- svdf : pd.DataFrame @@ -96,19 +101,24 @@ def sv_harmonic_mean( The start depth for harmonic mean to be computed end_depth : int or float The end depth for harmonic mean to be computed - + Returns ------- float The sound speed harmonic mean value """ # Clean up the sound speed value, ensuring that there's no negative value - svdf = svdf[svdf['sv'] < 0] + svdf = svdf[svdf["sv"] > 0].reset_index(drop=True) # Make all of the values absolute values, so we're only dealing with positives abs_start = abs(start_depth) abs_end = abs(end_depth) abs_sv = abs(svdf) # Get the index for the start of depth closest to specified start depth - start_index = abs_sv[(abs_sv['dd'].round() >= start_depth)].index[0] + try: + start_index = abs_sv[(abs_sv["dd"].round() >= abs_start)].index[0] + except IndexError: + raise ValueError("Dataframe is empty! Please check your data inputs.") - return _compute_hm(abs_sv['dd'].values, abs_sv['sv'].values, abs_start, abs_end, start_index) + return _compute_hm( + abs_sv["dd"].values, abs_sv["sv"].values, abs_start, abs_end, start_index + ) diff --git a/tests/__init__.py b/tests/__init__.py index e69de29..c154d79 100644 --- a/tests/__init__.py +++ b/tests/__init__.py @@ -0,0 +1,5 @@ +from pathlib import Path + +HERE = Path(__file__).parent + +TEST_DATA_FOLDER = HERE / "data" From 016030f0a8a60c819f3e456555c2f8a930c2207e Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 6 Apr 2023 21:08:11 +0000 Subject: [PATCH 05/11] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/seagap/harmonic_mean.py | 2 +- tests/test_harmonic_mean.py | 25 +++++++++---------------- 2 files changed, 10 insertions(+), 17 deletions(-) diff --git a/src/seagap/harmonic_mean.py b/src/seagap/harmonic_mean.py index b887618..b4d16aa 100644 --- a/src/seagap/harmonic_mean.py +++ b/src/seagap/harmonic_mean.py @@ -79,7 +79,7 @@ def _compute_hm( if z1 >= ze: break - + if cumsum == 0: # If cumulative sum is 0, most likely only one value return sv[start_index] diff --git a/tests/test_harmonic_mean.py b/tests/test_harmonic_mean.py index a18dc23..8fffe0c 100644 --- a/tests/test_harmonic_mean.py +++ b/tests/test_harmonic_mean.py @@ -1,23 +1,21 @@ +import numpy as np import pandas as pd import pytest -import numpy as np + +from seagap.harmonic_mean import _compute_hm, sv_harmonic_mean from . import TEST_DATA_FOLDER -from seagap.harmonic_mean import sv_harmonic_mean, _compute_hm + @pytest.fixture() def sound_profile_data() -> pd.DataFrame: - sv_file = TEST_DATA_FOLDER / 'sound_profile_hm.csv' + sv_file = TEST_DATA_FOLDER / "sound_profile_hm.csv" return pd.read_csv(sv_file) @pytest.mark.parametrize( "end_depth,expected_hm", - [ - (-1176.5866, 1481.551), - (-1146.5881, 1481.521), - (-1133.7305, 1481.509) - ] + [(-1176.5866, 1481.551), (-1146.5881, 1481.521), (-1133.7305, 1481.509)], ) def test_sv_harmonic_mean(end_depth, expected_hm, sound_profile_data): svdf = sound_profile_data @@ -26,13 +24,8 @@ def test_sv_harmonic_mean(end_depth, expected_hm, sound_profile_data): assert harmonic_mean == expected_hm -@pytest.mark.parametrize( - "test_idx,expected_hm", - [ - (3, 1501.69), - (6, 1501.225) - ] -) + +@pytest.mark.parametrize("test_idx,expected_hm", [(3, 1501.69), (6, 1501.225)]) def test__compute_hm(test_idx, expected_hm): dd = np.arange(7) * 10 sv = np.arange(1502, 1500, step=-0.31) @@ -42,7 +35,7 @@ def test__compute_hm(test_idx, expected_hm): part_sv = sv[:test_idx] start_depth, end_depth = part_dd[0], part_dd[-1] result_hm = round(_compute_hm(dd, sv, start_depth, end_depth, 0), 3) - + # Check for result to match expected assert result_hm == expected_hm From 766c131a77cf6b6d918fa7c5d67dec72d01f180e Mon Sep 17 00:00:00 2001 From: Landung 'Don' Setiawan Date: Thu, 6 Apr 2023 15:21:46 -0700 Subject: [PATCH 06/11] Remove catching single value --- src/seagap/harmonic_mean.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/seagap/harmonic_mean.py b/src/seagap/harmonic_mean.py index b4d16aa..6e74a8a 100644 --- a/src/seagap/harmonic_mean.py +++ b/src/seagap/harmonic_mean.py @@ -80,9 +80,6 @@ def _compute_hm( if z1 >= ze: break - if cumsum == 0: - # If cumulative sum is 0, most likely only one value - return sv[start_index] return (ze - zs) / cumsum From 418eb5de02826a7bffbfce571c926435edcde4d3 Mon Sep 17 00:00:00 2001 From: Landung 'Don' Setiawan Date: Thu, 6 Apr 2023 16:42:01 -0700 Subject: [PATCH 07/11] Add docs and return type annotations as suggested --- src/seagap/harmonic_mean.py | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/src/seagap/harmonic_mean.py b/src/seagap/harmonic_mean.py index 6e74a8a..15bf622 100644 --- a/src/seagap/harmonic_mean.py +++ b/src/seagap/harmonic_mean.py @@ -1,5 +1,4 @@ import math -from typing import Union import numba import numpy as np @@ -10,15 +9,30 @@ def _compute_hm( dd: np.ndarray, sv: np.ndarray, - start_depth: Union[int, float], - end_depth: Union[int, float], + start_depth: float, + end_depth: float, start_index: int, -): +) -> float: """ Computes harmonic mean. It's a direct translation from the original Fortran code found in src/cal_sv_harmonic_mean/get_sv_harmonic_mean.F called subroutine `sv_harmon_mean` + + Parameters + ---------- + dd : np.ndarray + Depth values array dd[n] + sv : np.ndarray + Sound speed values array sv[n] + start_depth : float + The start depth for calculation + end_depth : float + The end depth for calculation + start_index : int + The start index for the first value + computed by start depth + """ # TODO: Find a way to vectorize this computation @@ -84,8 +98,8 @@ def _compute_hm( def sv_harmonic_mean( - svdf: pd.DataFrame, start_depth: Union[int, float], end_depth: Union[int, float] -): + svdf: pd.DataFrame, start_depth: float, end_depth: float +) -> float: """ Computes harmonic mean from a sound profile containing depth (dd) and sound speed (sv) From 03cda3adead05d47a4d1f26d8d558e77b0504b2a Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 6 Apr 2023 23:42:12 +0000 Subject: [PATCH 08/11] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/seagap/harmonic_mean.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/seagap/harmonic_mean.py b/src/seagap/harmonic_mean.py index 15bf622..cf8e331 100644 --- a/src/seagap/harmonic_mean.py +++ b/src/seagap/harmonic_mean.py @@ -97,9 +97,7 @@ def _compute_hm( return (ze - zs) / cumsum -def sv_harmonic_mean( - svdf: pd.DataFrame, start_depth: float, end_depth: float -) -> float: +def sv_harmonic_mean(svdf: pd.DataFrame, start_depth: float, end_depth: float) -> float: """ Computes harmonic mean from a sound profile containing depth (dd) and sound speed (sv) From 3c084cf12941f9afdcc84874c7c4851963b5ed26 Mon Sep 17 00:00:00 2001 From: Landung 'Don' Setiawan Date: Fri, 7 Apr 2023 16:49:47 -0700 Subject: [PATCH 09/11] Remove validation for transducer_delay_time --- src/seagap/configs/solver.py | 12 +----------- 1 file changed, 1 insertion(+), 11 deletions(-) diff --git a/src/seagap/configs/solver.py b/src/seagap/configs/solver.py index 376dfda..0f3ff96 100644 --- a/src/seagap/configs/solver.py +++ b/src/seagap/configs/solver.py @@ -143,10 +143,7 @@ class Solver(BaseModel): ) transducer_delay_time: float = Field( 0.0, - description=( - "Transducer Delay Time - delay at surface transducer (secs.). " - "Options: ship/SV3 = 0.0s, WG = 0.1s" - ), + description="Transducer Delay Time - delay at surface transducer (secs). ", ) harmonic_mean_start_depth: float = Field( 0.0, description="Start depth in meters for harmonic mean computation." @@ -170,10 +167,3 @@ class Solver(BaseModel): "which data points will be excluded from solution" ), ) - - @validator("transducer_delay_time") - def check_transducer_delay_time(cls, v): - """Validate transducer wait time value""" - if v not in [0.0, 0.1]: - raise ValueError("Transducer wait time must either be 0.0s or 0.1s") - return v From cd6c6d5c2a41a385cede262b555dd3b5a1f4ee91 Mon Sep 17 00:00:00 2001 From: Landung 'Don' Setiawan Date: Fri, 7 Apr 2023 16:52:37 -0700 Subject: [PATCH 10/11] Remove unused module --- src/seagap/configs/solver.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/seagap/configs/solver.py b/src/seagap/configs/solver.py index 0f3ff96..3eb36f9 100644 --- a/src/seagap/configs/solver.py +++ b/src/seagap/configs/solver.py @@ -6,7 +6,7 @@ from typing import Any, List, Literal, Optional from uuid import uuid4 -from pydantic import BaseModel, Field, PrivateAttr, validator +from pydantic import BaseModel, Field, PrivateAttr from .io import InputData From 4b488e9c9fff552f02d51d00a975dce7bcf10cdd Mon Sep 17 00:00:00 2001 From: Landung 'Don' Setiawan Date: Mon, 10 Apr 2023 11:47:37 -0700 Subject: [PATCH 11/11] Update hmean computation to use scipy's hmean --- pyproject.toml | 1 + src/seagap/harmonic_mean.py | 123 ++++++++++-------------------------- tests/test_harmonic_mean.py | 32 +++++++--- 3 files changed, 56 insertions(+), 100 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index 0f9c5b8..1c4d5da 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -28,6 +28,7 @@ dependencies = [ "pandas>=1.5.3", "pydantic>=1.10.6", "pluggy>=1.0.0", + "scipy>=1.10.1", "typer>=0.7.0" ] diff --git a/src/seagap/harmonic_mean.py b/src/seagap/harmonic_mean.py index cf8e331..eeda430 100644 --- a/src/seagap/harmonic_mean.py +++ b/src/seagap/harmonic_mean.py @@ -1,100 +1,45 @@ -import math - -import numba -import numpy as np import pandas as pd +import scipy.stats -@numba.njit -def _compute_hm( - dd: np.ndarray, - sv: np.ndarray, - start_depth: float, - end_depth: float, - start_index: int, -) -> float: +def _compute_hm(svdf: pd.DataFrame, start_depth: float, end_depth: float) -> float: """ - Computes harmonic mean. - It's a direct translation from the original Fortran code found in - src/cal_sv_harmonic_mean/get_sv_harmonic_mean.F called - subroutine `sv_harmon_mean` + Computes harmonic mean using `scipy's hmean `_ method. + It takes the sound speed 'sv' as the input array and the depth 'dd' differences as the weights. + + Note that this function assumes absolute values for depth array, + sound speed array, start depth, and end depth. + + The underlying formula is + + H = (w1+...+wn) / ((w1/x1)+...+(wn/xn)) + + H is the resulting harmonic mean + w is the weight value, in this case, the depth differences + x is the input value, in this case, the sound speed Parameters ---------- - dd : np.ndarray - Depth values array dd[n] - sv : np.ndarray - Sound speed values array sv[n] + svdf : pd.DataFrame + Sound speed profile data as dataframe with columns 'dd' and 'sv' start_depth : float The start depth for calculation end_depth : float The end depth for calculation - start_index : int - The start index for the first value - computed by start depth - """ - # TODO: Find a way to vectorize this computation - - # Assign start depth and end depth to zs and ze - zs = start_depth - ze = end_depth - - # Extract the first and second depth values - z1 = dd[start_index] - z2 = dd[start_index + 1] - - # Extract the first and second sound speed values - c_z1 = sv[start_index] - c_z2 = sv[start_index + 1] - - # Set start depth to initial depth to keep track of it - zi = zs - - # If the second depth value (z2) is greater than and equal to the - # end depth value (ze) then the final depth (zf) should be assigned - # to the end depth value (ze), otherwise, the final depth - # should be assigned to the second depth values - if z2 >= ze: - zf = ze - else: - zf = z2 - - # Start cumulative sum as 0.0 - cumsum = 0.0 - - # Loop over the whole array for depth and sound speed, - # starting from the index of the second value to the whole - # depth data, which is assumed to be the same exact number - # as the sound speed data - for i in range(start_index + 1, len(dd)): - # calculate the slope of the two points - # for sound speed and depth - # slope = (sv2 - sv1) / (d2 - d1) - b = (c_z2 - c_z1) / (z2 - z1) - wi = zi - z1 + c_z1 / b - wf = zf - z1 + c_z1 / b - - wi = math.log((zi - z1) * b + c_z1) / b - wf = math.log((zf - z1) * b + c_z1) / b - - delta = wf - wi - cumsum = cumsum + delta - z1 = zf - z2 = dd[i + 1] - c_z1 = c_z2 - c_z2 = sv[i + 1] - zi = zf - - if ze > zi and ze < z2: - zf = ze - else: - zf = z2 - - if z1 >= ze: - break - - return (ze - zs) / cumsum + """ # noqa + for col in ["dd", "sv"]: + if col not in svdf.columns: + raise ValueError(f"{col} column must exist in the input dataframe!") + + filtdf = svdf[ + (svdf["dd"].round() >= start_depth) & (svdf["dd"].round() <= end_depth) + ] + + # Get weights + weights = filtdf["dd"].diff() + + return scipy.stats.hmean(filtdf["sv"], weights=weights, nan_policy="omit") def sv_harmonic_mean(svdf: pd.DataFrame, start_depth: float, end_depth: float) -> float: @@ -123,11 +68,7 @@ def sv_harmonic_mean(svdf: pd.DataFrame, start_depth: float, end_depth: float) - abs_end = abs(end_depth) abs_sv = abs(svdf) # Get the index for the start of depth closest to specified start depth - try: - start_index = abs_sv[(abs_sv["dd"].round() >= abs_start)].index[0] - except IndexError: + if len(abs_sv) == 0: raise ValueError("Dataframe is empty! Please check your data inputs.") - return _compute_hm( - abs_sv["dd"].values, abs_sv["sv"].values, abs_start, abs_end, start_index - ) + return _compute_hm(abs_sv, abs_start, abs_end) diff --git a/tests/test_harmonic_mean.py b/tests/test_harmonic_mean.py index 8fffe0c..b0a4e02 100644 --- a/tests/test_harmonic_mean.py +++ b/tests/test_harmonic_mean.py @@ -15,7 +15,7 @@ def sound_profile_data() -> pd.DataFrame: @pytest.mark.parametrize( "end_depth,expected_hm", - [(-1176.5866, 1481.551), (-1146.5881, 1481.521), (-1133.7305, 1481.509)], + [(-1176.5866, 1481.542), (-1146.5881, 1481.513), (-1133.7305, 1481.5)], ) def test_sv_harmonic_mean(end_depth, expected_hm, sound_profile_data): svdf = sound_profile_data @@ -25,19 +25,33 @@ def test_sv_harmonic_mean(end_depth, expected_hm, sound_profile_data): assert harmonic_mean == expected_hm -@pytest.mark.parametrize("test_idx,expected_hm", [(3, 1501.69), (6, 1501.225)]) -def test__compute_hm(test_idx, expected_hm): +@pytest.mark.parametrize( + "start_idx,end_idx,expected_hm", + [(0, 3, 1501.535), (0, 6, 1501.07), (2, 5, 1500.915), (4, 7, 1500.295)], +) +def test__compute_hm(start_idx, end_idx, expected_hm): dd = np.arange(7) * 10 sv = np.arange(1502, 1500, step=-0.31) + svdf = pd.DataFrame(dict(dd=dd, sv=sv)) + # Get partial of the data from test index - part_dd = dd[:test_idx] - part_sv = sv[:test_idx] - start_depth, end_depth = part_dd[0], part_dd[-1] - result_hm = round(_compute_hm(dd, sv, start_depth, end_depth, 0), 3) + partdf = svdf[start_idx:end_idx].copy() + start_depth, end_depth = partdf.iloc[0]["dd"], partdf.iloc[-1]["dd"] + result_hm = round(_compute_hm(svdf, start_depth, end_depth), 3) # Check for result to match expected assert result_hm == expected_hm - # Result should be same as regular mean as the data is linear here - assert result_hm == round(sum(part_sv) / test_idx, 3) + # Result should be same as manual computation of weighted harmonic mean + # https://en.wikipedia.org/wiki/Harmonic_mean#Weighted_harmonic_mean + # + # H = (w1+...+wn) / ((w1/x1)+...+(wn/xn)) + # + # H is the resulting harmonic mean + # w is the weight value, in this case, the depth differences + # x is the input value, in this case, the sound speed + w = partdf["dd"].diff() + x = partdf["sv"] + H = w.dropna().sum() / (w / x).dropna().sum() + assert result_hm == round(H, 3)