From 8d2764360da4ae7ec80f1683bb72e6e76239c75b Mon Sep 17 00:00:00 2001 From: Tobias Hienzsch Date: Tue, 3 Sep 2024 05:11:29 +0200 Subject: [PATCH] [python] Remove dead code --- .../common/{myasserts.py => asserts.py} | 14 +- src/python/pffdtd/common/misc.py | 32 +++ src/python/pffdtd/common/wavfile.py | 18 ++ src/python/pffdtd/geometry/box.py | 157 ++++++++------ src/python/pffdtd/geometry/math.py | 194 ------------------ .../pffdtd/geometry/tri_ray_intersection.py | 2 +- src/python/pffdtd/sim3d/engine.py | 3 +- src/python/pffdtd/sim3d/process_outputs.py | 3 +- src/python/pffdtd/voxelizer/vox_grid.py | 3 +- src/python/pffdtd/voxelizer/vox_grid_base.py | 2 +- src/python/pffdtd/voxelizer/vox_scene.py | 8 +- src/python/test/test_sim3d_infinite_baffle.py | 2 +- .../test/test_sim3d_locate_sound_source.py | 4 +- 13 files changed, 164 insertions(+), 278 deletions(-) rename src/python/pffdtd/common/{myasserts.py => asserts.py} (56%) create mode 100644 src/python/pffdtd/common/misc.py create mode 100644 src/python/pffdtd/common/wavfile.py diff --git a/src/python/pffdtd/common/myasserts.py b/src/python/pffdtd/common/asserts.py similarity index 56% rename from src/python/pffdtd/common/myasserts.py rename to src/python/pffdtd/common/asserts.py index 6ccae8d..4ac98df 100644 --- a/src/python/pffdtd/common/myasserts.py +++ b/src/python/pffdtd/common/asserts.py @@ -13,15 +13,13 @@ ############################################################################## import numpy as np + + def assert_np_array_float(x): - assert isinstance(x,np.ndarray) - assert x.dtype in [np.dtype('float32'),np.dtype('float64')] + assert isinstance(x, np.ndarray) + assert x.dtype in [np.dtype('float32'), np.dtype('float64')] -def assert_np_array_complex(x): - assert isinstance(x,np.ndarray) - #allow float or complex - assert x.dtype in [np.dtype('complex64'),np.dtype('complex128'),np.dtype('float32'),np.dtype('float64')] -#python 'int' is arbitrary size +# python 'int' is arbitrary size def assert_is_int(x): - assert isinstance(x,(int,np.integer)) + assert isinstance(x, (int, np.integer)) diff --git a/src/python/pffdtd/common/misc.py b/src/python/pffdtd/common/misc.py new file mode 100644 index 0000000..9025801 --- /dev/null +++ b/src/python/pffdtd/common/misc.py @@ -0,0 +1,32 @@ +import multiprocessing as mp +from pathlib import Path + + +def get_default_nprocs(): + return max(1, int(0.8*mp.cpu_count())) + + +def clear_dat_folder(dat_folder_str=None): + # clear dat folder + dat_path = Path(dat_folder_str) + if not dat_path.exists(): + dat_path.mkdir(parents=True) + else: + assert dat_path.is_dir() + + assert dat_path.exists() + assert dat_path.is_dir() + for f in dat_path.glob('*'): + try: + f.unlink() + except OSError as e: + print("Error: %s : %s" % (f, e.strerror)) + + +def yes_or_no(question): + while "the answer is invalid": + reply = str(input(question+' (y/n): ')).lower().strip() + if reply[:1] == 'y': + return True + if reply[:1] == 'n': + return False diff --git a/src/python/pffdtd/common/wavfile.py b/src/python/pffdtd/common/wavfile.py new file mode 100644 index 0000000..77356b3 --- /dev/null +++ b/src/python/pffdtd/common/wavfile.py @@ -0,0 +1,18 @@ +import numpy as np +import scipy.io.wavfile + + +def wavread(fname): + fs, data = scipy.io.wavfile.read(fname) # reads in (Nsamples,Nchannels) + if data.dtype == np.int16: + data = data/32768.0 + fs = np.float64(fs) + return fs, np.float64(data.T) + + +def wavwrite(fname, fs, data): + # expects (Nchannels,Nsamples), this will also assert that + data = np.atleast_2d(data) + # reads in (Nsamples,Nchannels) + scipy.io.wavfile.write(fname, int(fs), np.float32(data.T)) + print(f'wrote {fname} at fs={fs/1000:.2f} kHz') diff --git a/src/python/pffdtd/geometry/box.py b/src/python/pffdtd/geometry/box.py index 4d55660..1121403 100644 --- a/src/python/pffdtd/geometry/box.py +++ b/src/python/pffdtd/geometry/box.py @@ -1,61 +1,75 @@ # SPDX-License-Identifier: MIT +from typing import Any import numpy as np -from numpy import array as npa import numpy.random as npr -from pffdtd.geometry.math import mydefault,rotmatrix_ax_ang +from pffdtd.geometry.math import rotmatrix_ax_ang + + +def mydefault(a: Any, default: Any): + # only used in box + if a is None: + a = default + return a + class Box: """This is a class for a box. Used in a few places (vox_grid, tri_box). """ - def __init__(self,Lx=None,Ly=None,Lz=None,Rax=None,Rang=None,shift=None,centered=True): - #defaults - Lx = mydefault(Lx,1.0) - Ly = mydefault(Ly,1.0) - Lz = mydefault(Lz,1.0) - Rax = mydefault(Rax,npa([1.,1.,1.])) - Rang = mydefault(Rang,0.) - shift = mydefault(shift,npa([0.,0.,0.])) + + def __init__(self, Lx=None, Ly=None, Lz=None, Rax=None, Rang=None, shift=None, centered=True): + # defaults + Lx = mydefault(Lx, 1.0) + Ly = mydefault(Ly, 1.0) + Lz = mydefault(Lz, 1.0) + Rax = mydefault(Rax, np.array([1., 1., 1.])) + Rang = mydefault(Rang, 0.) + shift = mydefault(shift, np.array([0., 0., 0.])) self.centered = centered - self.init(Lx,Ly,Lz,Rax,Rang,shift) + self.init(Lx, Ly, Lz, Rax, Rang, shift) - def init(self,Lx,Ly,Lz,Rax,Rang,shift): + def init(self, Lx, Ly, Lz, Rax, Rang, shift): - verts = npa([[0.,0.,0.],[0.,0.,Lz],[0.,Ly,0.],[0.,Ly,Lz],[Lx,0.,0.],[Lx,0.,Lz],[Lx,Ly,0.],[Lx,Ly,Lz]]) + verts = np.array([[0., 0., 0.], [0., 0., Lz], [0., Ly, 0.], [0., Ly, Lz], [ + Lx, 0., 0.], [Lx, 0., Lz], [Lx, Ly, 0.], [Lx, Ly, Lz]]) if self.centered: - verts -= 0.5*npa([Lx,Ly,Lz]) + verts -= 0.5*np.array([Lx, Ly, Lz]) - A = npa([[-1.,0.,0.],[0.,-1.,0.],[0.,0.,-1.],[1.,0.,0.],[0.,1.,0.],[0.,0.,1.]]) + A = np.array([[-1., 0., 0.], [0., -1., 0.], [0., 0., -1.], + [1., 0., 0.], [0., 1., 0.], [0., 0., 1.]]) if self.centered: - b = 0.5*npa([Lx,Ly,Lz,Lx,Ly,Lz]) + b = 0.5*np.array([Lx, Ly, Lz, Lx, Ly, Lz]) else: - b = npa([0,0,0,Lx,Ly,Lz]) + b = np.array([0, 0, 0, Lx, Ly, Lz]) - #rotate verts - R = rotmatrix_ax_ang(Rax,Rang) + # rotate verts + R = rotmatrix_ax_ang(Rax, Rang) verts = verts.dot(R.T) + shift - #get AABB min/max - bmin = np.amin(verts,0) - bmax = np.amax(verts,0) + # get AABB min/max + bmin = np.amin(verts, 0) + bmax = np.amax(verts, 0) - #update A,b + # update A,b A = A.dot(R.T) b += shift.dot(A.T) - #set member variables + # set member variables self.A = A self.b = b self.verts = verts self.bmin = bmin self.bmax = bmax - self.edges = npa([[0,1],[0,2],[0,4],[1,3],[1,5],[2,3],[2,6],[4,5],[4,6],[3,7],[5,7],[6,7]]) - self.tris = npa([[0,1,3],[0,3,2],[1,7,3],[1,5,7],[0,2,6],[0,6,4],[4,7,5],[4,6,7],[2,3,7],[2,7,6],[0,5,1],[0,4,5]]) - self.quads = npa([[0,1,3,2],[0,4,5,1],[4,6,7,5],[1,5,7,3],[2,3,7,6],[0,2,6,4]]) + self.edges = np.array([[0, 1], [0, 2], [0, 4], [1, 3], [1, 5], [2, 3], [ + 2, 6], [4, 5], [4, 6], [3, 7], [5, 7], [6, 7]]) + self.tris = np.array([[0, 1, 3], [0, 3, 2], [1, 7, 3], [1, 5, 7], [0, 2, 6], [ + 0, 6, 4], [4, 7, 5], [4, 6, 7], [2, 3, 7], [2, 7, 6], [0, 5, 1], [0, 4, 5]]) + self.quads = np.array([[0, 1, 3, 2], [0, 4, 5, 1], [4, 6, 7, 5], [ + 1, 5, 7, 3], [2, 3, 7, 6], [0, 2, 6, 4]]) def randomise(self): Lx = 10*npr.random() @@ -64,23 +78,24 @@ def randomise(self): Rax = npr.rand(3) Rang = (-1.0 + 2.0*npr.random())*90 shift = npr.randn(3)*2 - 1 - self.init(Lx,Ly,Lz,Rax,Rang,shift) + self.init(Lx, Ly, Lz, Rax, Rang, shift) - #just draw (no show) - def _draw(self,backend='mayavi',r=None,color=(0,1,0)): + # just draw (no show) + def _draw(self, backend='mayavi', r=None, color=(0, 1, 0)): box = self if backend == 'mayavi': from mayavi import mlab for edge in box.edges: - mlab.plot3d(box.verts[edge,0],box.verts[edge,1],box.verts[edge,2],color=color,tube_radius=r) + mlab.plot3d(box.verts[edge, 0], box.verts[edge, 1], + box.verts[edge, 2], color=color, tube_radius=r) mlab.draw() print('drawing box..') else: raise - #draw and show - def draw(self,backend='mayavi',r=None,color=(0,1,0)): + # draw and show + def draw(self, backend='mayavi', r=None, color=(0, 1, 0)): box = self if backend == 'mayavi': from mayavi import mlab @@ -88,15 +103,17 @@ def draw(self,backend='mayavi',r=None,color=(0,1,0)): fig = mlab.gcf() box._draw() box._draw_faces() - fake_verts = np.c_[box.bmin,box.bmax] #stacks as column vectors 2x3 - mlab.plot3d(fake_verts[0],fake_verts[1],fake_verts[2],transparent=True,opacity=0) + # stacks as column vectors 2x3 + fake_verts = np.c_[box.bmin, box.bmax] + mlab.plot3d(fake_verts[0], fake_verts[1], + fake_verts[2], transparent=True, opacity=0) mlab.axes(xlabel='x', ylabel='y', zlabel='z', color=(0., 0., 0.)) fig.scene.interactor.interactor_style = tvtk.InteractorStyleTerrain() mlab.show() else: raise - def show(self,backend='mayavi'): + def show(self, backend='mayavi'): if backend == 'mayavi': from mayavi import mlab mlab.show() @@ -104,71 +121,85 @@ def show(self,backend='mayavi'): import matplotlib.pyplot as plt plt.show() - #draw faces without show() - def _draw_faces(self,backend='mayavi'): + # draw faces without show() + def _draw_faces(self, backend='mayavi'): box = self if backend == 'mayavi': from mayavi import mlab - x = box.verts[:,0]; y = box.verts[:,1]; z = box.verts[:,2] - mlab.triangular_mesh(x,y,z,box.tris,color=(0,0,0),opacity=0.6) + x = box.verts[:, 0] + y = box.verts[:, 1] + z = box.verts[:, 2] + mlab.triangular_mesh( + x, y, z, box.tris, color=(0, 0, 0), opacity=0.6) mlab.draw() print('drawing box faces..') else: raise + def main(): - #no shift, no rotation + # no shift, no rotation Lx = npr.random() Ly = npr.random() Lz = npr.random() print(f'Lx = {Lx:.2f}, Ly = {Ly:.2f}, Lz = {Lz:.2f}') - box = Box(Lx,Ly,Lz) - print('box 0 ... bmin = %.2f,%.2f,%.2f' % (box.bmin[0],box.bmin[1],box.bmin[2])) - print('box 0 ... bmax = %.2f,%.2f,%.2f' % (box.bmax[0],box.bmax[1],box.bmax[2])) - assert np.all(np.mean(box.verts,0).dot(box.A.T) < box.b) + box = Box(Lx, Ly, Lz) + print('box 0 ... bmin = %.2f,%.2f,%.2f' % + (box.bmin[0], box.bmin[1], box.bmin[2])) + print('box 0 ... bmax = %.2f,%.2f,%.2f' % + (box.bmax[0], box.bmax[1], box.bmax[2])) + assert np.all(np.mean(box.verts, 0).dot(box.A.T) < box.b) assert np.all((box.verts*(1-1e-4)).dot(box.A.T) < box.b) assert np.any((box.verts*(1+1e-4)).dot(box.A.T) > box.b) - #rotated + # rotated Lx = npr.random() Ly = npr.random() Lz = npr.random() - print('Lx = %.2f, Ly = %.2f, Lz = %.2f' % (Lx,Ly,Lz)) + print('Lx = %.2f, Ly = %.2f, Lz = %.2f' % (Lx, Ly, Lz)) Rax = npr.rand(3) Rang = (-1.0 + 2.0*npr.random())*90 - print('Rax = %.2f,%.2f,%.2f, Rang = %.2f degrees' % (Rax[0],Rax[1],Rax[2],Rang)) + print('Rax = %.2f,%.2f,%.2f, Rang = %.2f degrees' % + (Rax[0], Rax[1], Rax[2], Rang)) print() - box = Box(Lx,Ly,Lz,Rax,Rang) - print('Lx = %.2f, Ly = %.2f, Lz = %.2f' % (Lx,Ly,Lz)) - print('Rax = %.2f,%.2f,%.2f, Rang = %.2f degrees' % (Rax[0],Rax[1],Rax[2],Rang)) - print('box 1 ... bmin = %.2f,%.2f,%.2f' % (box.bmin[0],box.bmin[1],box.bmin[2])) - print('box 1 ... bmax = %.2f,%.2f,%.2f' % (box.bmax[0],box.bmax[1],box.bmax[2])) - assert np.all(np.mean(box.verts,0).dot(box.A.T) < box.b) + box = Box(Lx, Ly, Lz, Rax, Rang) + print('Lx = %.2f, Ly = %.2f, Lz = %.2f' % (Lx, Ly, Lz)) + print('Rax = %.2f,%.2f,%.2f, Rang = %.2f degrees' % + (Rax[0], Rax[1], Rax[2], Rang)) + print('box 1 ... bmin = %.2f,%.2f,%.2f' % + (box.bmin[0], box.bmin[1], box.bmin[2])) + print('box 1 ... bmax = %.2f,%.2f,%.2f' % + (box.bmax[0], box.bmax[1], box.bmax[2])) + assert np.all(np.mean(box.verts, 0).dot(box.A.T) < box.b) assert np.all((box.verts*(1-1e-4)).dot(box.A.T) < box.b) assert np.any((box.verts*(1+1e-4)).dot(box.A.T) > box.b) print() - #shifted and rotated + # shifted and rotated Lx = npr.random() Ly = npr.random() Lz = npr.random() - print('Lx = %.2f, Ly = %.2f, Lz = %.2f' % (Lx,Ly,Lz)) + print('Lx = %.2f, Ly = %.2f, Lz = %.2f' % (Lx, Ly, Lz)) Rax = npr.rand(3) Rang = (-1.0 + 2.0*npr.random())*90 - print('Rax = %.2f,%.2f,%.2f, Rang = %.2f degrees' % (Rax[0],Rax[1],Rax[2],Rang)) + print('Rax = %.2f,%.2f,%.2f, Rang = %.2f degrees' % + (Rax[0], Rax[1], Rax[2], Rang)) shift = npr.rand(3)*100 - print('shift = %.2f,%.2f,%.2f' % (shift[0],shift[1],shift[2])) - - box = Box(Lx,Ly,Lz,Rax,Rang,shift) - print('box 2 ... bmin = %.2f,%.2f,%.2f' % (box.bmin[0],box.bmin[1],box.bmin[2])) - print('box 2 ... bmax = %.2f,%.2f,%.2f' % (box.bmax[0],box.bmax[1],box.bmax[2])) - assert np.all(np.mean(box.verts,0).dot(box.A.T) < box.b) + print('shift = %.2f,%.2f,%.2f' % (shift[0], shift[1], shift[2])) + + box = Box(Lx, Ly, Lz, Rax, Rang, shift) + print('box 2 ... bmin = %.2f,%.2f,%.2f' % + (box.bmin[0], box.bmin[1], box.bmin[2])) + print('box 2 ... bmax = %.2f,%.2f,%.2f' % + (box.bmax[0], box.bmax[1], box.bmax[2])) + assert np.all(np.mean(box.verts, 0).dot(box.A.T) < box.b) assert np.all((((box.verts-shift)*(1-1e-4))+shift).dot(box.A.T) < box.b) assert np.any((((box.verts-shift)*(1+1e-4))+shift).dot(box.A.T) > box.b) box.draw() + if __name__ == '__main__': main() diff --git a/src/python/pffdtd/geometry/math.py b/src/python/pffdtd/geometry/math.py index 4697600..7f48441 100644 --- a/src/python/pffdtd/geometry/math.py +++ b/src/python/pffdtd/geometry/math.py @@ -1,21 +1,13 @@ # SPDX-License-Identifier: MIT """Miscellaneous python/numpy functions. Not all used or useful. """ -import argparse -import hashlib -import multiprocessing as mp -import os -from pathlib import Path from typing import Any -import struct import numpy as np import numpy.linalg as npl from numpy import cos, sin, pi -import scipy.io.wavfile -from pffdtd.common.myasserts import assert_np_array_complex EPS = np.finfo(np.float64).eps @@ -76,38 +68,6 @@ def rotate_az_el_deg(az_d, el_d): return R, Raz, Rel -def mydefault(a: Any, default: Any): - # only used in box - if a is None: - a = default - return a - - -def s2dhms(time: int): - # seconds to day/hour/min/s - # give time in seconds - day = time // (24 * 3600) - time = time % (24 * 3600) - hour = time // 3600 - time %= 3600 - minutes = time // 60 - time %= 60 - seconds = time - return (day, hour, minutes, seconds) - - -def clamp(x: Any, xmin: Any, xmax: Any): - assert isinstance(x, (np.ndarray, int, np.integer, float)) - assert isinstance(xmin, (np.ndarray, int, np.integer, float)) - assert isinstance(xmax, (np.ndarray, int, np.integer, float)) - return np.where(x < xmin, xmin, np.where(x > xmax, xmax, x)) - - -def vclamp(x, xmin, xmax): - # only this works for vectorized xmin,xmax - return np.where(x < xmin, xmin, np.where(x > xmax, xmax, x)) - - def dot2(v): return dotv(v, v) @@ -120,51 +80,10 @@ def vecnorm(v1): return np.sqrt(dot2(v1)) -def vecnorm2(v1): - return np.sum(v1**2, axis=-1) - - def normalise(v1, eps=EPS): return (v1.T/(vecnorm(v1)+eps)).T -def ceilint(x: float): - assert isinstance(x, float) - return np.int_(np.ceil(x)) - - -def roundint(x: float): - assert isinstance(x, float) - return np.int_(np.round(x)) - - -def floorint(x: float): - assert isinstance(x, float) - return np.int_(np.floor(x)) - - -def maxabs(x: Any): - assert_np_array_complex(x) - if x.ndim == 2: - return np.amax(np.abs(x), axis=-1)[:, None] - else: - return np.amax(np.abs(x)) - - -def to_col_2d(x: Any): - assert x.ndim < 3 - x = x.reshape(x.shape[-1], -1) - assert x.shape[0] > x.shape[1] - return x - - -def to_row_2d(x: Any): - assert x.ndim < 3 - x = x.reshape(-1, x.shape[-1]) - assert x.shape[1] > x.shape[0] - return x - - def ind2sub3d(ii, Nx, Ny, Nz): iz = ii % Nz iy = (ii - iz)//Nz % Ny @@ -177,103 +96,6 @@ def rel_diff(x0, x1): return (x0-x1)/(2.0**np.floor(np.log2(x0))) -def get_default_nprocs(): - return max(1, int(0.8*mp.cpu_count())) - - -def str2bool(v: str): - # taken from SO, for argparse bool problem - assert isinstance(v, str) - if v.lower() in ('true'): - return True - elif v.lower() in ('false'): - return False - else: - raise argparse.ArgumentTypeError('Boolean value expected.') - - -def clear_dat_folder(dat_folder_str=None): - # clear dat folder - dat_path = Path(dat_folder_str) - if not dat_path.exists(): - dat_path.mkdir(parents=True) - else: - assert dat_path.is_dir() - - assert dat_path.exists() - assert dat_path.is_dir() - for f in dat_path.glob('*'): - try: - f.unlink() - except OSError as e: - print("Error: %s : %s" % (f, e.strerror)) - - -def clear_folder(folder_str): - print(f'clearing {folder_str}...') - os.system('rm -rf ' + str(Path(folder_str)) + '/*') - print(f'{folder_str} cleared') - - -def double_to_hex(d): - return hex(struct.unpack('