Skip to content

linear_tf() frequency grid stops one step short of fmax for common inputs, including the defaults (fmax=30, freq_resolution=0.05) #98

Open

Description

@jsh9

Found while reviewing #74. Pre-existing on main (b764057); #74 does not change the behavior.

Problem

linear_tf() sets the number of frequencies with int(np.floor_divide(fmax, freq_resolution)). np.floor_divide() rounds down the exact quotient of the two binary floats. 0.05 and 0.1 are stored as slightly more than 1/20 and 1/10, so np.floor_divide(30.0, 0.05) is 599.0, not 600.0. The grid then ends at 29.95 Hz, not at fmax. The docstring calls fmax the "Maximum frequency of interest".

# length of transfer function
TF_size = int(np.floor_divide(fmax, freq_resolution))
freq_array = np.linspace(
freq_resolution, freq_resolution * TF_size, num=TF_size
)

freq_resolution : float, default=0.05
Frequency resolution of the frequency spectrum.
fmax : float, default=30.0
Maximum frequency of interest.

Vs_Profile.get_ampl_function() and Vs_Profile.get_transfer_function() use this grid, with the same defaults. Vs_Profile.get_f0_RO() and get_f0_BH() call get_ampl_function() with the defaults, so they never look at 30 Hz:

def get_ampl_function(
self,
*,
show_fig: bool = False,
freq_resolution: float = 0.05,
fmax: float = 30.0,
) -> tuple[Frequency_Spectrum, Frequency_Spectrum, Frequency_Spectrum]:

def get_f0_RO(self) -> float:
"""
Return the rock-outcrop fundamental frequency.
Returns
-------
f0_RO : float
Rock-outcrop fundamental frequency.
"""
return self.get_ampl_function(show_fig=False)[0].get_f0()
def get_f0_BH(self) -> float:
"""
Return the borehole fundamental frequency.
Returns
-------
f0_BH : float
Borehole fundamental frequency.
"""
return self.get_ampl_function(show_fig=False)[1].get_f0()

import numpy as np
from PySeismoSoil.class_Vs_profile import Vs_Profile

print(np.floor_divide(30.0, 0.05))  # 599.0  (expected 600.0)

prof = Vs_Profile(np.array([[10.0, 0.0], [200.0, 1000.0]]).T)
for fmax, freq_resolution in [(30.0, 0.05), (30.0, 0.1), (10.0, 0.1), (20.0, 0.2), (15.0, 0.5)]:
    af_RO = prof.get_ampl_function(fmax=fmax, freq_resolution=freq_resolution)[0]
    freq = af_RO.amplitude_2col[:, 0]
    print(fmax, freq_resolution, len(freq), round(freq[-1], 6))
# 30.0 0.05 599 29.95
# 30.0 0.1 299 29.9
# 10.0 0.1 99 9.9
# 20.0 0.2 99 19.8
# 15.0 0.5 30 15.0
fmax freq_resolution points last frequency expected
30 0.05 (default) 599 29.95 600 points, 30
30 0.1 299 29.9 300 points, 30
10 0.1 99 9.9 100 points, 10
20 0.2 99 19.8 100 points, 20
15 0.5 30 15 correct

test_get_amplif_function() uses fmax=15 and freq_resolution=0.5, which divide exactly, so it does not catch this.

Fix

  1. Use true division with a small tolerance, then round down:

    TF_size = int(np.floor(fmax / freq_resolution + 1e-6))

    A plain int(fmax / freq_resolution) is not enough, because the quotient can also land just below a whole number (0.3 / 0.1 is 2.9999999999999996). With only this change on main, the repro gives 600, 300, 100, 100 and 30 points, each ending at fmax, and all 240 tests pass. linear_site_resp() also calls linear_tf(); the tests that run it (tests/test_helper_simulations.py, and Linear_Simulation in tests/test_class_simulation.py) still pass.

  2. Add a test that Vs_Profile.get_ampl_function() with the defaults returns 600 frequencies ending at 30.0, and that fmax=10, freq_resolution=0.1 returns 100 frequencies ending at 10.0.

Metadata

Metadata

Assignees

No one assigned

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions