add a gaussian_filter and find_peaks for scalarfield and vectorfield #1120
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,168 @@ | ||
| import numpy as np | ||
|
|
||
| from pyphare.pharesee.hierarchy import ScalarField, VectorField | ||
Check noticeCode scanning / CodeQL Unused import Note |
||
| from pyphare.pharesee.hierarchy.hierarchy_utils import compute_hier_from | ||
Check noticeCode scanning / CodeQL Unused import Note |
||
|
|
||
| from scipy.ndimage import gaussian_filter | ||
|
|
||
|
|
||
|
|
||
| def _compute_gaussian_filter_on_scalarfield(patch_datas, **kwargs): | ||
| from scipy.ndimage import gaussian_filter | ||
|
|
||
| ndim = patch_datas["value"].box.ndim | ||
| nb_ghosts = kwargs["nb_ghosts"] | ||
| sigma = kwargs["sigma"] | ||
| ds = np.asarray(patch_datas["value"][:]) | ||
|
|
||
| ds_ = np.full(list(ds.shape), np.nan) | ||
|
|
||
| gf_ = gaussian_filter(ds, sigma=sigma) | ||
| select = tuple([slice(nb_ghosts, -nb_ghosts) for _ in range(ndim)]) | ||
| ds_[select] = np.asarray(gf_[select]) | ||
|
|
||
| return ( | ||
| {"name": "value", "data": ds_, "centering": patch_datas["value"].centerings}, | ||
| ) | ||
|
|
||
|
|
||
| def _compute_gaussian_filter_on_vectorfield(patch_datas, **kwargs): | ||
| from scipy.ndimage import gaussian_filter | ||
|
|
||
| ref_name = next(iter(patch_datas.keys())) | ||
|
|
||
| ndim = patch_datas[ref_name].box.ndim | ||
| nb_ghosts = kwargs["nb_ghosts"] | ||
| sigma = kwargs["sigma"] | ||
| ds_x = np.asarray(patch_datas["x"][:]) | ||
| ds_y = np.asarray(patch_datas["y"][:]) | ||
| ds_z = np.asarray(patch_datas["z"][:]) | ||
|
|
||
| dsx_ = np.full(list(ds_x.shape), np.nan) | ||
| dsy_ = np.full(list(ds_y.shape), np.nan) | ||
| dsz_ = np.full(list(ds_z.shape), np.nan) | ||
|
|
||
| gfx_ = gaussian_filter(ds_x, sigma=sigma) | ||
| gfy_ = gaussian_filter(ds_y, sigma=sigma) | ||
| gfz_ = gaussian_filter(ds_z, sigma=sigma) | ||
|
|
||
| select = tuple([slice(nb_ghosts, -nb_ghosts) for _ in range(ndim)]) | ||
|
|
||
| dsx_[select] = np.asarray(gfx_[select]) | ||
| dsy_[select] = np.asarray(gfy_[select]) | ||
| dsz_[select] = np.asarray(gfz_[select]) | ||
|
|
||
| return ( | ||
| {"name": "x", "data": dsx_, "centering": patch_datas["x"].centerings}, | ||
| {"name": "y", "data": dsy_, "centering": patch_datas["y"].centerings}, | ||
| {"name": "z", "data": dsz_, "centering": patch_datas["z"].centerings}, | ||
| ) | ||
|
|
||
|
|
||
| def gFilt(hier, **kwargs): | ||
|
github-advanced-security[bot] marked this conversation as resolved.
Fixed
|
||
| sigma = kwargs.get("sigma", 2) | ||
|
|
||
| # time0 = list(hier.times())[0] | ||
| # level0 = 0 | ||
| # p0 = 0 | ||
| # pd0 = hier.levels(time0)[level0].patches[p0].patch_datas | ||
| # key0 = list(pd0.keys())[0] | ||
| # nb_ghosts = np.max(pd0[key0].ghosts_nbr) | ||
|
|
||
|
|
||
| nb_ghosts = np.max(list(hier.level(0).patches[0].patch_datas.values())[0].ghosts_nbr) | ||
|
|
||
| if nb_ghosts < sigma : | ||
| print("nb_ghosts ({0}) < sigma ({1}) : your gaussian filter might be dirty".format(nb_ghosts, sigma)) | ||
|
|
||
| if hier.ndim == 1: | ||
| if isinstance(hier, ScalarField) : | ||
| h = compute_hier_from(_compute_gaussian_filter_on_scalarfield, hier, nb_ghosts=nb_ghosts, sigma=sigma) | ||
| return ScalarField(h) | ||
| elif isinstance(hier, VectorField) : | ||
| h = compute_hier_from(_compute_gaussian_filter_on_vectorfield, hier, nb_ghosts=nb_ghosts, sigma=sigma) | ||
| return VectorField(h) | ||
| else: | ||
| return NotImplemented | ||
|
coderabbitai[bot] marked this conversation as resolved.
Outdated
|
||
| else: | ||
| return NotImplemented | ||
|
|
||
|
|
||
| def gF(hier, **kwargs): | ||
| sigma = kwargs.get("sigma", 4) | ||
| if sigma == 1: | ||
| raise ValueError("sigma value has to be at least 2") | ||
| h_ = hier.__deepcopy__(memo={}) | ||
| ndim = hier.ndim | ||
| n_pad = 4*sigma+1 | ||
|
Comment on lines
+90
to
+96
There was a problem hiding this comment. Sigma validation is incomplete and
Proposed fix def gF(hier, **kwargs):
sigma = kwargs.get("sigma", 4)
- if sigma == 1:
- raise ValueError("sigma value has to be at least 2")
+ if not np.isscalar(sigma):
+ raise TypeError("sigma must be a scalar value")
+ if sigma < 2:
+ raise ValueError("sigma value has to be at least 2")
h_ = hier.__deepcopy__(memo={})
ndim = hier.ndim
- n_pad = 4*sigma+1
+ n_pad = int(4*sigma+1)🧰 Tools🪛 Ruff (0.15.0)[warning] 93-93: Avoid specifying long messages outside the exception class (TRY003) 🤖 Prompt for AI Agents |
||
| # The gaussian filter is calculated on the box extended by | ||
| # n_pad. Hence, the number of points is large enough so that the value | ||
| # at the last point of the real box is as equal as possible to the one | ||
| # at the first point of the next box... | ||
|
|
||
| for time in h_.times(): | ||
| interp_ = hier.interpol(time) | ||
| for lvl in h_.levels(time).values(): | ||
| for patch in lvl.patches: | ||
| names = list(patch.patch_datas.keys()) | ||
| box = patch.box | ||
|
|
||
| for name in names: | ||
| pdata = patch.patch_datas[name] | ||
| nb_ghosts = pdata.ghosts_nbr | ||
| if not n_pad > nb_ghosts: | ||
| raise ValueError('sigma value is too small') | ||
|
Comment on lines
+112
to
+113
There was a problem hiding this comment. Inverted guard condition: The check Consider clarifying the intent and using - if not n_pad > nb_ghosts:
- raise ValueError('sigma value is too small')
+ if np.any(n_pad <= nb_ghosts):
+ raise ValueError('sigma value is too small: n_pad must exceed nb_ghosts')🧰 Tools🪛 Ruff (0.15.0)[warning] 113-113: Avoid specifying long messages outside the exception class (TRY003) 🤖 Prompt for AI Agents |
||
|
|
||
| r_ = [] | ||
| for i in range(ndim): | ||
| s_ = np.arange(box.lower[i]-n_pad, box.upper[i]+2+n_pad)*pdata.dl[i] | ||
| r_.append(s_) | ||
|
|
||
| func, _ = interp_[name] | ||
|
|
||
| if ndim == 1: | ||
| data = func(r_[0]) | ||
| elif ndim == 2: | ||
| data = func(r_[0], r_[1]) | ||
| elif ndim == 3: | ||
| data = func(r_[0], r_[1], r_[2]) | ||
| else: | ||
| raise ValueError('unvalid dimension') | ||
|
|
||
|
coderabbitai[bot] marked this conversation as resolved.
|
||
| gf_ = gaussian_filter(data, sigma=sigma) | ||
| select = tuple([slice(n_pad-nb_ghosts[i], -n_pad+nb_ghosts[i]) for i in range(ndim)]) | ||
| pdata.dataset = np.asarray(gf_[select]) | ||
|
|
||
| return h_ | ||
|
|
||
|
Comment on lines
+90
to
+136
There was a problem hiding this comment. 🧩 Analysis chain🌐 Web query:
💡 Result: In
This is exactly how SciPy’s docs define it: “scalar or sequence of scalars … as a sequence [per axis], or as a single number, in which case it is equal for all axes.” [1] Example: for a 2D array, Source: SciPy API docs for [1] https://docs.scipy.org/doc/scipy-1.11.4/reference/generated/scipy.ndimage.gaussian_filter.html 🏁 Script executed: # Search for sigma validation patterns in the codebase
rg "sigma" pyphare/pyphare/core/ufuncs.py -B 2 -A 2
# Look for tests of the gF function
find . -name "*.py" -type f | xargs grep -l "gF" | grep -i testRepository: PHAREHUB/PHARE Length of output: 2272 🏁 Script executed: # Check for any existing sigma type validation elsewhere in the codebase
rg "np.isscalar.*sigma" --type py
rg "isinstance.*sigma" --type pyRepository: PHAREHUB/PHARE Length of output: 40 🏁 Script executed: # Understand how numpy handles float indices in slicing
python3 << 'EOF'
import numpy as np
# Test 1: Does numpy allow float indices in slicing?
try:
arr = np.array([1, 2, 3, 4, 5])
result = arr[2.5:]
print("Float index in slice: works")
except TypeError as e:
print(f"Float index in slice: {e}")
# Test 2: What about float arithmetic in slice construction?
n_pad = 2 * 4.5 # Simulating float sigma
sigma = 4.5
try:
print(f"n_pad (float): {n_pad}, type: {type(n_pad)}")
ndim = 1
select = tuple([slice(n_pad, -n_pad)])
print(f"Slice construction with float: {select}")
arr = np.zeros((10,))
result = arr[select]
print(f"Slicing with float: works, result shape: {result.shape}")
except (TypeError, ValueError) as e:
print(f"Slicing with float: {e}")
# Test 3: Check scipy.ndimage.gaussian_filter sigma behavior
from scipy.ndimage import gaussian_filter
test_data = np.random.rand(10, 10)
try:
result = gaussian_filter(test_data, sigma=4.5)
print(f"gaussian_filter with float sigma: works")
except Exception as e:
print(f"gaussian_filter with float sigma: {e}")
try:
result = gaussian_filter(test_data, sigma=[2.0, 3.0])
print(f"gaussian_filter with sequence sigma: works")
except Exception as e:
print(f"gaussian_filter with sequence sigma: {e}")
EOFRepository: PHAREHUB/PHARE Length of output: 144 🏁 Script executed: # Verify what happens when float values are used in slice objects with numpy arrays
python3 << 'EOF'
import numpy as np
# Test float indices in slice with numpy array
arr = np.zeros((10,))
n_pad = 4.5 # float value
try:
select = tuple([slice(n_pad, -n_pad)])
result = arr[select]
print(f"Float in slice works. Result shape: {result.shape}")
print(f"Slice object: {select}")
except TypeError as e:
print(f"Float in slice fails: {e}")
# Test with actual usage from code
ndim = 1
nb_ghosts = np.array([1])
n_pad_float = 4.5
n_pad_int = int(np.ceil(4.5))
select_float = tuple([slice(n_pad_float-nb_ghosts[i], -n_pad_float+nb_ghosts[i]) for i in range(ndim)])
select_int = tuple([slice(n_pad_int-nb_ghosts[i], -n_pad_int+nb_ghosts[i]) for i in range(ndim)])
arr = np.zeros((15,))
try:
result_float = arr[select_float]
result_int = arr[select_int]
print(f"Float slice result: {result_float.shape}")
print(f"Int slice result: {result_int.shape}")
print(f"Are results identical? {np.array_equal(result_float, result_int)}")
except Exception as e:
print(f"Error: {e}")
# Test with sequence sigma
from scipy.ndimage import gaussian_filter
data = np.random.rand(10, 10)
try:
# This should work according to scipy docs
result = gaussian_filter(data, sigma=(2.0, 3.0))
print(f"gaussian_filter with sequence sigma works")
# But what happens if we treat it like the gF code does?
sigma = (2.0, 3.0)
n_pad = 2 * sigma # This will create a tuple
print(f"n_pad with sequence sigma: {n_pad}, type: {type(n_pad)}")
# Try to use it in slicing
ndim = 2
nb_ghosts = np.array([1, 1])
select = tuple([slice(n_pad[i]-nb_ghosts[i], -n_pad[i]+nb_ghosts[i]) for i in range(ndim)])
result = data[select]
print(f"Slicing with sequence-derived n_pad works: {result.shape}")
except Exception as e:
print(f"Error with sequence sigma: {type(e).__name__}: {e}")
EOFRepository: PHAREHUB/PHARE Length of output: 144 🏁 Script executed: # Check if there are any tests for the gF function to see expected behavior
find . -name "*.py" -type f -exec grep -l "def test.*gF\|gF.*test" {} \;Repository: PHAREHUB/PHARE Length of output: 40 🏁 Script executed: # Search for actual usage of gF function in the codebase
rg "gF\(" --type py -B 2 -A 2Repository: PHAREHUB/PHARE Length of output: 292 🏁 Script executed: # Check if gF is ever called with sequence sigma in tests or examples
rg "gF\s*\(" --type py -A 5 | head -50
# Look for test files that might call gF
fd "test.*\.py" --type f | head -20Repository: PHAREHUB/PHARE Length of output: 1424 🏁 Script executed: # Simulate what happens with sequence sigma in the actual code path
python3 << 'EOF'
import numpy as np
# Simulate the actual code behavior with sequence sigma
sigma_sequence = (2.0, 3.0)
n_pad = 2 * sigma_sequence
print(f"n_pad with sequence sigma: {n_pad}")
print(f"Type of n_pad: {type(n_pad)}")
# Try the slice arithmetic
ndim = 2
nb_ghosts = np.array([1, 1])
try:
select = tuple([slice(n_pad[i]-nb_ghosts[i], -n_pad[i]+nb_ghosts[i]) for i in range(ndim)])
print(f"Slice construction works: {select}")
except (TypeError, IndexError) as e:
print(f"Slice construction fails: {type(e).__name__}: {e}")
# Now test with float scalar sigma
sigma_float = 2.5
n_pad_float = 2 * sigma_float
print(f"\nn_pad with float scalar sigma: {n_pad_float}")
print(f"Type of n_pad: {type(n_pad_float)}")
try:
select = tuple([slice(n_pad_float-nb_ghosts[i], -n_pad_float+nb_ghosts[i]) for i in range(ndim)])
print(f"Slice construction works with float: {select}")
# Test actual slicing
data = np.zeros((10, 10))
result = data[select]
print(f"Actual slicing works with float indices: {result.shape}")
except (TypeError, IndexError) as e:
print(f"Error: {type(e).__name__}: {e}")
EOFRepository: PHAREHUB/PHARE Length of output: 175 🏁 Script executed: # Simple test without numpy of the sequence sigma issue
python3 << 'EOF'
# Test 1: What happens when we multiply a tuple by 2?
sigma_sequence = (2.0, 3.0)
n_pad = 2 * sigma_sequence
print(f"n_pad with sequence sigma: {n_pad}")
print(f"Type: {type(n_pad)}")
# Test 2: Can we use tuple elements in arithmetic?
try:
result = n_pad[0] - 1
print(f"Can index tuple n_pad: {result}")
except TypeError as e:
print(f"Error indexing: {e}")
# Test 3: What about scalar float multiplication?
sigma_scalar = 2.5
n_pad_scalar = 2 * sigma_scalar
print(f"\nn_pad with float scalar sigma: {n_pad_scalar}")
print(f"Type: {type(n_pad_scalar)}")
# Test 4: Using float in range operations (simulating slice indices)
try:
n_pad_float = 4.5
idx = int(n_pad_float)
print(f"Can convert float to int: {idx}")
# But what about direct float use?
print(f"Float value {n_pad_float} - int {1} = {n_pad_float - 1}, type: {type(n_pad_float - 1)}")
except Exception as e:
print(f"Error: {e}")
EOFRepository: PHAREHUB/PHARE Length of output: 293 🏁 Script executed: # Check the actual scipy.ndimage.gaussian_filter source or behavior with different sigma types
# by looking at how it's used in other parts of the codebase
rg "gaussian_filter" --type py -B 2 -A 2 | grep -E "(sigma|gaussian_filter)" | head -30Repository: PHAREHUB/PHARE Length of output: 2006 🏁 Script executed: # Look at how sigma is validated in the gFilt function mentioned earlier
rg "def gFilt" --type py -A 20Repository: PHAREHUB/PHARE Length of output: 1585 Add validation for scalar 🐛 Proposed fix def gF(hier, **kwargs):
- sigma = kwargs.get("sigma", 4)
- if sigma == 1:
- raise ValueError("sigma value has to be at least 2")
+ sigma = kwargs.get("sigma", 4)
+ if not np.isscalar(sigma):
+ raise ValueError("sigma must be a scalar")
+ if sigma < 2:
+ raise ValueError("sigma value has to be at least 2")
h_ = hier.__deepcopy__(memo={})
ndim = hier.ndim
- n_pad = 2*sigma
+ n_pad = int(np.ceil(2 * sigma))🧰 Tools🪛 Ruff (0.14.13)93-93: Avoid specifying long messages outside the exception class (TRY003) 127-127: Avoid specifying long messages outside the exception class (TRY003) 🤖 Prompt for AI Agents |
||
|
|
||
| def peakIds(hier, **kwargs): | ||
| from scipy.signal import find_peaks | ||
|
|
||
| times = list(hier.times()) | ||
| if len(times) == 1: | ||
| time = times[0] | ||
| else: | ||
| raise ValueError('multiple time is not yet implemented') | ||
|
|
||
| pks_ = np.array([]) | ||
| hgs_ = np.array([]) | ||
|
|
||
| names_ = kwargs.pop("names", None) | ||
| if names_ is None: | ||
| names_ = list(hier.levels(time)[0].patches[0].patch_datas.keys()) | ||
|
|
||
| ph_ = kwargs.get('peak_heights', None) | ||
| if ph_ is None: | ||
| raise ValueError("the kwarg 'peak_heights' is mandatory for now...") | ||
|
|
||
| for lvl in hier.levels(time).values(): | ||
| for patch in lvl.patches: | ||
| for name in names_: | ||
| pdata = patch.patch_datas[name] | ||
| ds = np.asarray(pdata.dataset) | ||
| pks = find_peaks(ds, **kwargs) | ||
| for pk, hg in zip(pks[0], pks[1]['peak_heights']): | ||
| pks_ = np.append(pks_, np.add(np.multiply(pk, patch.dl), patch.origin)) | ||
| hgs_ = np.append(hgs_, hg) | ||
|
coderabbitai[bot] marked this conversation as resolved.
Outdated
|
||
|
|
||
| return pks_, hgs_ | ||
|
There was a problem hiding this comment. 🧩 Analysis chain🌐 Web query:
💡 Result:
About the
Sources: [1] 🏁 Script executed: cat -n pyphare/pyphare/core/ufuncs.py | sed -n '136,166p'Repository: PHAREHUB/PHARE Length of output: 1306 Pass 🐛 Proposed fix names_ = kwargs.pop("names", None)
if names_ is None:
names_ = list(hier.levels(time)[0].patches[0].patch_datas.keys())
- ph_ = kwargs.get('peak_heights', None)
- if ph_ is None:
- raise ValueError("the kwarg 'peak_heights' is mandatory for now...")
+ find_peaks_kwargs = {k: v for k, v in kwargs.items() if k != "peak_heights"}
+ if "height" not in find_peaks_kwargs:
+ raise ValueError("the kwarg 'height' is mandatory for now...")
for lvl in hier.levels(time).values():
for patch in lvl.patches:
for name in names_:
pdata = patch.patch_datas[name]
ds = np.asarray(pdata.dataset)
- pks = find_peaks(ds, **kwargs)
+ pks = find_peaks(ds, **find_peaks_kwargs)
for pk, hg in zip(pks[0], pks[1]['peak_heights']):
pks_ = np.append(pks_, np.add(np.multiply(pk, patch.dl), patch.origin))
hgs_ = np.append(hgs_, hg)🧰 Tools🪛 Ruff (0.14.13)143-143: Avoid specifying long messages outside the exception class (TRY003) 154-154: Avoid specifying long messages outside the exception class (TRY003) 162-162: Add explicit value for parameter (B905) 🤖 Prompt for AI Agents |
||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
is there something to do here ?