From a758f159c4b3f1331f85dcdf7a7be620ee51287d Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Wed, 18 Sep 2024 13:24:52 +0100 Subject: [PATCH 1/3] upgrade Fourier methods and cosmological handling --- exptool/analysis/centering.py | 526 +++++++++++++++++---------------- exptool/observables/bulge.py | 59 ++++ exptool/observables/fourier.py | 374 ++++++++++++----------- 3 files changed, 510 insertions(+), 449 deletions(-) create mode 100644 exptool/observables/bulge.py diff --git a/exptool/analysis/centering.py b/exptool/analysis/centering.py index f4291cb..7e44d42 100644 --- a/exptool/analysis/centering.py +++ b/exptool/analysis/centering.py @@ -1,270 +1,274 @@ -########################################################################################### -# -# centering.py -# Tools to -# -# -# 09 Jul 2021: First introduction -# -# -''' - -centering.py (part of exptool.basis) - Tools for wrangling unruly n-body data - Originally used for NewHorizon bar detection - - -# read in the test data -P = np.genfromtxt('testbar.small.part',names=True,dtype=None) - -rtest,fpower,fangle,fvpower = fourier_tabulate(P['x'],P['y'],P['vx'],P['vy'],P['mass']) - -# the simplest test is the ratio of m=4 to m=2 Fourier velocities - -plt.plot(rtest,fvpower[4]/fvpower[2]) -# when this ratio is over 1, it's the end of the bar. - - - -''' - -import numpy as np - - - -import numpy as np - -def xnorm(vec): return np.sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]) +""" +centering.py +Tools for wrangling unruly n-body data. + - recenter particle data + - rotate particle data + - compute density profiles + - compute Fourier moments + - compute velocity-weighted Fourier moments + - compute Fourier moments in radial bins -def compute_fourier(x,y,w,harmonic=2): - phi = np.arctan2(y,x) - A = np.nansum(w*np.cos(harmonic*phi))/np.nansum(w) - B = np.nansum(w*np.sin(harmonic*phi))/np.nansum(w) - mod = np.sqrt(A*A + B*B) - angle = np.arctan2(B,A)/2. - return mod,angle +Originally used for NewHorizon bar detection (Reddish et al. 2022) + +09 Jul 2021: First introduction +18 Sep 2024: Class structure overhaul +""" +import numpy as np -def recentre_pos_and_vel(x,y,z,vx,vy,vz,mass,r_max): - """compute center of mass and mass weighted velocity center within r_max and recentre position and velocity - """ - mask = (x**2 + y**2 + z**2) < r_max**2 - mass_tot = np.sum(mass[mask]) - x_cm = np.sum(x[mask]*mass[mask])/mass_tot - y_cm = np.sum(y[mask]*mass[mask])/mass_tot - z_cm = np.sum(z[mask]*mass[mask])/mass_tot - vx_cm = np.sum(vx[mask]*mass[mask])/mass_tot - vy_cm = np.sum(vy[mask]*mass[mask])/mass_tot - vz_cm = np.sum(vz[mask]*mass[mask])/mass_tot - return x-x_cm,y-y_cm,z-z_cm,vx-vx_cm,vy-vy_cm,vz-vz_cm - -def rotate(x,y,z,axis,angle): - # HERE: - x_rot=x - y_rot=y - z_rot=z - axisx=axis[0]*np.ones(len(x)) - axisy=axis[1]*np.ones(len(x)) - axisz=axis[2]*np.ones(len(x)) - - dot=x*axisx+y*axisy+z*axisz - - crossx=axisy*z-axisz*y - crossy=axisz*x-axisx*z - crossz=axisx*y-axisy*x - - cosa=np.cos(angle) - sina=np.sin(angle) - - x_rot=x*cosa+crossx*sina+axis[0]*dot*(1-cosa) - y_rot=y*cosa+crossy*sina+axis[1]*dot*(1-cosa) - z_rot=z*cosa+crossz*sina+axis[2]*dot*(1-cosa) - - return x_rot,y_rot,z_rot - -def compute_rotation_to_vec(x,y,z,vx,vy,vz,mass,vec): - """compute rotation axis and anglr to vector vec +class ParticleAlignment: """ - Lxtot = np.sum(y*vz*mass - z*vy*mass) - Lytot = np.sum(z*vx*mass - x*vz*mass) - Lztot = np.sum(x*vy*mass - y*vx*mass) - - L = np.zeros(3) - L = np.hstack([Lxtot,Lytot,Lztot]) + Example usage: + + import numpy as np + + # Sample particle data (replace with your actual data) + N = 1000 # Number of particles + x = np.random.randn(N) # Random x positions + y = np.random.randn(N) # Random y positions + z = np.random.randn(N) # Random z positions + vx = np.random.randn(N) # Random velocities in x + vy = np.random.randn(N) # Random velocities in y + vz = np.random.randn(N) # Random velocities in z + mass = np.random.rand(N) # Random masses + + # Create an instance of ParticleAlignment + particle_system = ParticleAlignment(x, y, z, vx, vy, vz, mass) + + # Recenter positions and velocities within a maximum radius of 10 + x_recentered, y_recentered, z_recentered, vx_recentered, vy_recentered, vz_recentered = particle_system.recentre_pos_and_vel(r_max=10) + + # Rotate particles around a specified axis (e.g., z-axis) by 45 degrees + axis = np.array([0, 0, 1.]) # Rotate around z-axis + angle = np.pi / 4 # 45 degrees in radians + x_rot, y_rot, z_rot = particle_system.rotate(axis, angle) + + # Compute the rotation axis and angle to align the angular momentum with a given vector + vec = np.array([1., 0, 0]) # Target vector (e.g., x-axis) + axis, angle = particle_system.compute_rotation_to_vec(vec) + + # Use the shrinking sphere method to recenter and shrink + w = np.random.rand(N) # Random weighting factors for the particles + x_shrunk, y_shrunk, z_shrunk, vx_shrunk, vy_shrunk, vz_shrunk = particle_system.shrinking_sphere(w) + + # Compute the density profile for the system + R = np.sqrt(x**2 + y**2 + z**2) # Radial distances + W = mass # Weights as the particle masses + density, enclosed_mass, potential = particle_system.compute_density_profile(R, W) + + # Print results + print("Recentered positions:", x_recentered, y_recentered, z_recentered) + print("Rotated positions:", x_rot, y_rot, z_rot) + print("Rotation axis and angle:", axis, angle) + print("Density profile:", density) - L = L/np.linalg.norm(L) - vec = vec/np.linalg.norm(vec) - axis = np.zeros(3) - axis = np.cross(L,vec) - axis = axis/np.linalg.norm(axis) - - c = np.dot(L,vec) # cosine of the angle - angle = np.arccos(np.clip(c, -1, 1)) # angle between L and vec - #@@@angle = np.arccos(c) # angle between L and vec - - return axis, angle - - - - -def shrinking_sphere(x,y,z,vx,vy,vz,w,rmin=1.,stepsize=0.5,tol=0.001,verbose=0): - # - if stepsize >= 1.: - print('reverting to default step size:') - stepsize=0.5 - # - tshiftx = np.nanmedian(np.nansum(w*x)/np.nansum(w)) - tshifty = np.nanmedian(np.nansum(w*y)/np.nansum(w)) - tshiftz = np.nanmedian(np.nansum(w*z)/np.nansum(w)) - # - # first guess and normalisation - x -= tshiftx - y -= tshifty - z -= tshiftz - rval = np.sqrt(x*x + y*y + z*z) - rmax = np.nanmax(rval) - rmax0 = np.nanmax(rval) - # - if verbose: print('initial guess: {0:5.0f},{1:5.0f},{2:5.0f}'.format(tshiftx,tshifty,tshiftz)) - while rmax > rmin: - #print(rmax) - u = np.where(rvalrbinstmp[indx]) & (R rmin: + u = np.where(rval < stepsize * rmax)[0] + if float(u.size) / float(self.x.size) < tol: + print(f'Too few particles to continue at radius ratio {stepsize * rmax / rmax0}') + break + + comx = np.nanmedian(np.nansum(w[u] * self.x[u]) / np.nansum(w[u])) + comy = np.nanmedian(np.nansum(w[u] * self.y[u]) / np.nansum(w[u])) + comz = np.nanmedian(np.nansum(w[u] * self.z[u]) / np.nansum(w[u])) + + self.x -= comx + self.y -= comy + self.z -= comz + + tshiftx += comx + tshifty += comy + tshiftz += comz + + rval = np.sqrt(self.x**2 + self.y**2 + self.z**2) + rmax *= stepsize + + comvx = np.nanmedian(np.nansum(w[u] * self.vx[u]) / np.nansum(w[u])) + comvy = np.nanmedian(np.nansum(w[u] * self.vy[u]) / np.nansum(w[u])) + comvz = np.nanmedian(np.nansum(w[u] * self.vz[u]) / np.nansum(w[u])) + + if verbose: + print(f'Final shift: {tshiftx:.0f}, {tshifty:.0f}, {tshiftz:.0f}') + print(f'Final velocity shift: {comvx:.0f}, {comvy:.0f}, {comvz:.0f}') + + self.vx -= comvx + self.vy -= comvy + self.vz -= comvz + + return self.x, self.y, self.z, self.vx, self.vy, self.vz + + def compute_density_profile(self, R, W, rbins=10.**np.linspace(-3.7, 0.3, 100)): + """ + Compute the density profile, enclosed mass, and potential for the particle system. + + Parameters: + R : array + Radial positions of the particles. + W : array + Weights of the particles. + rbins : array, optional + Radial bins to compute the profile (default is logarithmic bins from 10^-3.7 to 10^0.3). + + Returns: + dens : array + Density profile. + menc : array + Enclosed mass profile. + potp : array + Potential profile. + """ + dens = np.zeros(rbins.size) + menc = np.zeros(rbins.size) + potp = np.zeros(rbins.size) + + astronomicalG = 0.0000043009125 + rbinstmp = np.concatenate([rbins, [2. * rbins[-1] - rbins[-2]]]) + + for indx, val in enumerate(rbinstmp[:-1]): + w = np.where((R > rbinstmp[indx]) & (R < rbinstmp[indx + 1]))[0] + wenc = np.where((R < rbinstmp[indx + 1]))[0] + shellsize = (4 / 3.) * np.pi * (rbinstmp[indx + 1]**3 - rbinstmp[indx]**3) + dens[indx] = np.nansum(W[w]) / shellsize + menc[indx] = np.nansum(W[wenc]) + potp[indx] = np.sqrt(astronomicalG * menc[indx] / (rbinstmp[indx + 1])) + + return dens, menc, potp diff --git a/exptool/observables/bulge.py b/exptool/observables/bulge.py new file mode 100644 index 0000000..7bde4d8 --- /dev/null +++ b/exptool/observables/bulge.py @@ -0,0 +1,59 @@ +""" +bulge.py +Tools for measuring the bulge mass of a galaxy. + - measure_bulge_mass: Compute the bulge mass of a galaxy + +18 Sep 2024: First introduction + +""" +import numpy as np + +def measure_bulge_mass(x, y, vx, vy, mass): + """ + Compute the bulge mass of a galaxy based on particle positions, velocities, and masses. + + The function assumes that the bulge is composed of stars with negative tangential velocities + (i.e., retrograde orbits), based on the assumption that prograde stars are part of the disc. + It computes the radius that contains 50% of the bulge stars and returns the total bulge mass. + + Parameters: + - x: Array of particle x-positions. + - y: Array of particle y-positions. + - vx: Array of particle velocities in the x-direction. + - vy: Array of particle velocities in the y-direction. + - mass: Array of particle masses. + + Returns: + - bulge_mass: The estimated total bulge mass (double the sum of negative velocity stars' masses). + + TODO: + - Extend to 3D: Currently, this function works in 2D (x and y plane). If ywe want a more accurate 3D mass measurement, we might want to include the z-coordinate for the radial position or compute velocities more generally in 3D. + + - Handling edge cases: we might want to include additional checks for cases where no retrograde particles are found, or where the tangential velocity is undefined due to division by zero (e.g., at the galaxy center where rpos = 0). + + """ + # Compute the radial distance from the center in the x-y plane + rpos = np.sqrt(x**2 + y**2) + + # Compute the tangential velocity + vtan = (x * vy - y * vx) / rpos + + # Determine the overall rotational direction (positive for prograde, negative for retrograde) + direction = np.sign(np.nanmean(vtan)) + + # Align all tangential velocities with the prograde direction + vtan *= direction + + # Identify particles with retrograde orbits (vtan < 0) + negvel = np.where(vtan < 0)[0] + + # Compute the radius containing 50% of the bulge stars (those with negative velocities) + bulge_radius_50 = np.nanpercentile(rpos[negvel], 50) + print(f"50% bulge radius: {bulge_radius_50:3.2f} kpc") + + # Compute the total bulge mass as double the mass of retrograde stars + bulge_mass = 2.0 * np.nansum(mass[negvel]) + + return bulge_mass + + diff --git a/exptool/observables/fourier.py b/exptool/observables/fourier.py index f121cd9..293458f 100644 --- a/exptool/observables/fourier.py +++ b/exptool/observables/fourier.py @@ -1,203 +1,201 @@ -############ -# -# CODE FOR SHADOW BAR PLOTS -# -HALO DENSITY OVER TIME, HALO DENSITY COMPARED TO NORMAL -# -HALO VELOCITY COMPARED TO NORMAL, ON/OFF BAR -# -# -''' -Fourier methods for images - -''' - - - -def excise_center(xpos,ypos,zpos,xvel,yvel,zvel,pote,mass,cutoff): - r = (xpos*xpos+ypos*ypos+zpos*zpos)**0.5 - lowr = np.where(r < cutoff)[0] - return xpos[lowr],ypos[lowr],zpos[lowr],xvel[lowr],yvel[lowr],zvel[lowr],pote[lowr],mass[lowr] - - - -def quik_transform(xpos,ypos,time,barfile): - tst = SimReadThree() - tst.X = np.zeros([1,len(xpos)]) - tst.X[0,:] = xpos - tst.Y = np.zeros([1,len(ypos)]) - tst.Y[0,:] = ypos - tst.T = time - tst.bartransform(barfile) - return tst.TX[0,:],tst.TY[0,:] - - -def add_file(infile,runtag,samplei,rcut): - filein = infile+'.'+runtag+'.%05.i' %samplei - times,massi,xposi,yposi,zposi,xveli,yveli,zveli,potei,nbodies,comptitle = psp_full_read(filein,'dark',2,10000000) - xpos,ypos,zpos,xvel,yvel,zvel,pote,mass = excise_center(xposi,yposi,zposi,xveli,yveli,zveli,potei,massi,rcut) - return xpos,ypos,zpos,xvel,yvel,zvel,pote,mass,times - - -def within_annulus(xpos,ypos,zpos,rcentral,rwidth,zheight): - # calculate whether a particle is within a given annulus (and some given height cut) - # - rval = (xpos**2.+ypos**2.)**0.5 - if ( (rval > rcentral-rwidth) & (rval < rcentral+rwidth) & (abs(zpos) < zheight) ): return 1 - else: return 0 - - - """ -# make the 3d density - -binz = 50 -tdbins = np.linspace(-0.07,0.07,binz) -dtd = tdbins[1]-tdbins[0] -tddist = np.zeros([binz,binz,binz]) - -for i in range(0,len(XPOS)): - indx = int(np.floor( (XPOS[i] - tdbins[0])/dtd)) - indy = int(np.floor( (YPOS[i] - tdbins[0])/dtd)) - indz = int(np.floor( (ZPOS[i] - tdbins[0])/dtd)) - if (indx>=0) and (indx=0) and (indy=0) and (indz=0) and (indx=0) and (indy=0) and (indz rcentral - rwidth) & (rval < rcentral + rwidth) & (np.abs(self.zpos) < zheight) + return within_annulus + + def within_radius(self, radius): + """ + Determine which particles are within a spherical radius. + + Parameters: + - radius: Radius to filter particles within. + + Returns: + - Boolean array indicating particles within the radius. + """ + r = np.sqrt(self.xpos**2 + self.ypos**2 + self.zpos**2) + return r < radius -# flatten these out -XPOS = np.concatenate(xpos) -YPOS = np.concatenate(ypos) -ZPOS = np.concatenate(zpos) -XVEL = np.concatenate(xvel) -YVEL = np.concatenate(yvel) -ZVEL = np.concatenate(zvel) -#POTE = np.concatenate(pote) -MASS = np.concatenate(mass) -#TIMES = np.concatenate(times) - -#print 'It took %3.2f seconds to build the particle list.' % (time.time()-t1) - - - - -xpos_tmp,ypos_tmp,zpos_tmp,xvel_tmp,yvel_tmp,zvel_tmp,pote_tmp,mass_tmp,times_tmp = add_file(infile,runtag,startimg+i,0.06,comp) - -times,massi,xposi,yposi,zposi,xveli,yveli,zveli,potei,nbodies,comptitle - - - -TI,MASS,XPOS,YPOS,ZPOS,XVEL,YVEL,ZVEL,POTE,NP,CT = psp_full_read('/Users/mpetersen/Research/NBody/Disk064a/OUT.run064a.01000','star',2,10000000) - - - -rbins = np.linspace(0.0,0.1,20) + def compute_fourier(self, harmonic=2): + """ + Compute Fourier moments for the particle system. + + Parameters: + - harmonic: The harmonic number for the Fourier moment (default is 2). + + Returns: + - mod: The modulus of the Fourier moment. + - angle: The phase angle of the Fourier moment. + """ + phi = np.arctan2(self.y, self.x) + A = np.nansum(self.mass * np.cos(harmonic * phi)) / np.nansum(self.mass) + B = np.nansum(self.mass * np.sin(harmonic * phi)) / np.nansum(self.mass) + mod = np.sqrt(A * A + B * B) + angle = np.arctan2(B, A) / 2. + return mod, angle + + def compute_fourier_vel(self, harmonic=2): + """ + Compute velocity-weighted Fourier moments. + + Parameters: + - harmonic: The harmonic number for the Fourier moment (default is 2). + + Returns: + - mod: The modulus of the velocity-weighted Fourier moment. + """ + phi = np.arctan2(self.y, self.x) + # Compute the tangential velocity + vel = (self.x * self.vy - self.y * self.vx) / np.sqrt(self.x**2 + self.y**2) + A = np.nansum(self.mass * vel * np.cos(harmonic * phi)) / np.nansum(self.mass) + B = np.nansum(self.mass * vel * np.sin(harmonic * phi)) / np.nansum(self.mass) + mod = np.sqrt(A * A + B * B) + return mod + + def fourier_tabulate(self, bins=75): + """ + Compute the Fourier moments and velocity-weighted Fourier moments in radial bins. + + Parameters: + - bins: Number of radial bins to use (default is 75). + + Returns: + - rtest: The radial bin centres. + - fpower: The Fourier moments for harmonics 0 to 4. + - fangle: The Fourier angles for harmonics 0 to 4. + - fvpower: The velocity-weighted Fourier moments for harmonics 0 to 4. + """ + rval = np.sqrt(self.x**2 + self.y**2) + rtest = np.linspace(0., np.nanpercentile(rval, 75), bins) + dr = rtest[1] - rtest[0] + rbin = (np.floor(rval / dr)).astype(int) + + fpower = np.zeros((5, rtest.size)) + fangle = np.zeros((5, rtest.size)) + fvpower = np.zeros((5, rtest.size)) + + for ir, rv in enumerate(rtest): + w = np.where(rbin == ir)[0] + for h in range(5): + fpower[h, ir], fangle[h, ir] = self.compute_fourier(harmonic=h) + fvpower[h, ir] = self.compute_fourier_vel(harmonic=h) + + return rtest, fpower, fangle, fvpower -RVAL = (YPOS**2.+XPOS**2.)**0.5 + def compute_fourier_in_annuli(self, rcentral, rwidth, zheight, harmonic=2): + """ + Compute Fourier moments for particles within a specified annulus. -lowx = np.where(RVAL < 0.1)[0] -xpos = XPOS[lowx] -ypos = YPOS[lowx] -mass = MASS[lowx] + Parameters: + - rcentral: Central radius of the annulus. + - rwidth: Width of the annulus. + - zheight: Height of the cylindrical region. + - harmonic: Harmonic number for Fourier analysis. + Returns: + - mod: Amplitude of the Fourier component within the annulus. + - angle: Angle of the Fourier component within the annulus. + """ + within = self.within_annulus(rcentral, rwidth, zheight) + return self.compute_fourier(self.xpos[within], self.ypos[within], self.mass[within], harmonic) -def compute_fourier(XPOS,YPOS,MASS,RBINS,mmax): - r_dig = np.digitize( (YPOS**2.+XPOS**2.)**0.5,rbins,right=True) - aval = np.zeros([mmax,len(rbins)]) - bval = np.zeros([mmax,len(rbins)]) - for indx,r in enumerate(rbins): - yes = np.where( r_dig == indx)[0] - print len(yes) - aval[0,indx] = np.sum(MASS[yes]) - for m in range(1,mmax): - aval[m,indx] = np.sum(MASS[yes] * np.cos(float(m)*np.arctan2(YPOS[yes],XPOS[yes]))) - bval[m,indx] = np.sum(MASS[yes] * np.sin(float(m)*np.arctan2(YPOS[yes],XPOS[yes]))) - return aval,bval + def compute_fourier_in_radius(self, radius, harmonic=2): + """ + Compute Fourier moments for particles within a specified spherical radius. + Parameters: + - radius: Radius to filter particles within. + - harmonic: Harmonic number for Fourier analysis. + Returns: + - mod: Amplitude of the Fourier component within the radius. + - angle: Angle of the Fourier component within the radius. + """ + within = self.within_radius(radius) + return self.compute_fourier(self.xpos[within], self.ypos[within], self.mass[within], harmonic) -aval,bval = compute_fourier(xpos,ypos,mass,rbins,20) -mvals = np.linspace(0,20,1) -imval = (aval**2.+bval**2.)**0.5 -""" From 37e78375c93802ee4893c4e5b3f29bd2cdb5669b Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Wed, 18 Sep 2024 13:32:35 +0100 Subject: [PATCH 2/3] cherry pick improvements from stale cosmological branch --- exptool/basis/bsplinebasis.py | 149 +++++++++++++++++++++++++------ exptool/io/outcoef.py | 20 ++--- exptool/observables/visualize.py | 2 +- exptool/utils/kde_3d.py | 10 +-- 4 files changed, 139 insertions(+), 42 deletions(-) diff --git a/exptool/basis/bsplinebasis.py b/exptool/basis/bsplinebasis.py index 66b78d6..e5ec70b 100644 --- a/exptool/basis/bsplinebasis.py +++ b/exptool/basis/bsplinebasis.py @@ -3,41 +3,138 @@ Set up bspline basis 31-Mar-2021: First written +04 Nov 2023: revised to improve playground vision is to follow Vasiliev's technique for deriving radial basis functions using spline representation -built on scipy's implementation, but may use the naive implementation -as well. +This version has edges clamped at 0. + TODO 1. investigate loss functions 2. investigate penalised fitting +3. explore monotone cubic spline + +import numpy as np +import matplotlib.pyplot as plt + +# Generate random control points for the B-spline curve +np.random.seed(0) +control_points = np.random.rand(10) * 10 # 10 random control points + +# Generate knot vector (assuming cubic B-spline with open uniform knots) +num_points = len(control_points) +num_knots = num_points + 4 # 3 for cubic B-spline +knots = np.linspace(0, 10, num_knots) + +# Create a B-spline curve object +bspline_curve = BSpline(control_points) + +# Evaluate the B-spline curve at various points for plotting +x_values = np.linspace(0, 10, 1000) +y_values = np.array([bspline_curve.evaluate(x, knots) for x in x_values]) + +# Plot the B-spline curve +plt.figure(figsize=(8, 6)) +plt.plot(x_values, y_values, color='blue', label='B-spline Curve') +plt.scatter(knots[self.degree:-self.degree], control_points, color='red', label='Control Points') +plt.xlabel('X') +plt.ylabel('Y') +plt.title('B-spline Curve with Control Points') +plt.legend() +plt.grid(True) +plt.show() + +# Evaluate the error between B-spline curve and observed data +error = bspline_curve.evaluate_error(observed_data, knots) +print("Error:", error) + +# Function to minimize (difference between B-spline curve and points) +def objective_function(control_points): + curve_values = np.array([bspline(x, t, control_points, k) for x in points]) + error = np.sum((curve_values - points) ** 2) + return error + +# Optimize control points to minimize the difference between B-spline curve and points +from scipy.optimize import minimize +result = minimize(objective_function, control_points, method='BFGS') + +# Fitted control points +fitted_control_points = result.x + ''' import numpy as np -from scipy.interpolate import BSpline - - -def B(x, k, i, t): - ''' - the pedagogical example from https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.BSpline.html - - ''' - if k == 0: - return 1.0 if t[i] <= x < t[i+1] else 0.0 - if t[i+k] == t[i]: - c1 = 0.0 - else: - c1 = (x - t[i])/(t[i+k] - t[i]) * B(x, k-1, i, t) - if t[i+k+1] == t[i+1]: - c2 = 0.0 - else: - c2 = (t[i+k+1] - x)/(t[i+k+1] - t[i+1]) * B(x, k-1, i+1, t) - return c1 + c2 - -def bspline(x, t, c, k): - n = len(t) - k - 1 - assert (n >= k+1) and (len(c) >= n) - return sum(c[i] * B(x, k, i, t) for i in range(n)) + +class BSpline: + def __init__(self, control_points, degree=3): + """ + Initialize a B-spline curve with given control points and degree. + + Parameters: + - control_points: List of control points defining the shape of the B-spline curve. + - degree: Degree of the B-spline curve (default is 3 for cubic B-spline). + """ + self.control_points = control_points + self.degree = degree + + def _basis_function(self, x, k, i, t): + """ + Compute the value of a B-spline basis function. + + Parameters: + - x: Evaluation point. + - k: Degree of the B-spline basis function. + - i: Index of the basis function within the B-spline curve. + - t: Knot vector defining the position and multiplicity of knots. + + Returns: + - The value of the B-spline basis function at the given point x. + """ + # Base case: B-spline basis function of degree 0 is 1 within the interval [t[i], t[i+1]) + if k == 0: + return 1.0 if t[i] <= x < t[i+1] else 0.0 + + # Recursive calculation for higher degree basis functions + c1 = 0.0 if t[i+k] == t[i] else (x - t[i]) / (t[i+k] - t[i]) * B(x, k-1, i, t) + c2 = 0.0 if t[i+k+1] == t[i+1] else (t[i+k+1] - x) / (t[i+k+1] - t[i+1]) * B(x, k-1, i+1, t) + + return c1 + c2 + + def evaluate(self, x, knots): + """ + Evaluate the B-spline curve at a specific point. + + Parameters: + - x: Evaluation point. + - knots: Knot vector. + + Returns: + - The value of the B-spline curve at the given point x. + """ + n = len(knots) - self.degree - 1 + assert (n >= self.degree + 1) and (len(self.control_points) >= n), "Invalid number of control points or knots." + + # Calculate the B-spline curve value by summing contributions from control points + curve_value = sum(self.control_points[i] * self._basis_function(x, self.degree, i, knots) for i in range(n)) + return curve_value + + def evaluate_error(self, observed_data, knots): + """ + Evaluate the error of the B-spline curve compared to observed data. + + Parameters: + - observed_data: List or array of observed data points. + - knots: Knot vector. + + Returns: + - The sum of squared differences between the B-spline curve and observed data points. + """ + # Evaluate the B-spline curve at each observed data point + bspline_values = np.array([self.evaluate(x, knots) for x in observed_data]) + + # Calculate the sum of squared differences (error) between B-spline curve and observed data + error = np.sum((bspline_values - observed_data) ** 2) + return error \ No newline at end of file diff --git a/exptool/io/outcoef.py b/exptool/io/outcoef.py index fb83ba6..73f8f82 100644 --- a/exptool/io/outcoef.py +++ b/exptool/io/outcoef.py @@ -162,7 +162,7 @@ def read_binary_eof_coefficients_old(self): # return to beginning f.seek(0) - [time0] = np.fromfile(f, dtype=np.float,count=1) + [time0] = np.fromfile(f, dtype='float',count=1) [mmax,nmax] = np.fromfile(f, dtype=np.uint32,count=2) # hard-coded to match specifications. @@ -178,17 +178,17 @@ def read_binary_eof_coefficients_old(self): for tt in range(0,n_outputs): - [time0] = np.fromfile(f, dtype=np.float,count=1) + [time0] = np.fromfile(f, dtype='float',count=1) [dummym,dummyn] = np.fromfile(f, dtype=np.uint32,count=2) times[tt] = time0 for mm in range(0,mmax+1): - coef_array[tt,0,mm,:] = np.fromfile(f, dtype=np.float,count=nmax) + coef_array[tt,0,mm,:] = np.fromfile(f, dtype='float',count=nmax) if mm > 0: - coef_array[tt,1,mm,:] = np.fromfile(f, dtype=np.float,count=nmax) + coef_array[tt,1,mm,:] = np.fromfile(f, dtype='float',count=nmax) #return times,coef_array @@ -228,7 +228,7 @@ def read_binary_sl_coefficients_old(self): [string1] = np.fromfile(f, dtype='a64',count=1) - [time0,scale] = np.fromfile(f, dtype=np.float,count=2) + [time0,scale] = np.fromfile(f, dtype='float',count=2) [nmax,lmax] = np.fromfile(f, dtype=np.uint32,count=2) # hard-coded to match specifications. @@ -247,14 +247,14 @@ def read_binary_sl_coefficients_old(self): for tt in range(0,n_outputs): [string1] = np.fromfile(f, dtype='a64',count=1) - [time0,scale] = np.fromfile(f, dtype=np.float,count=2) + [time0,scale] = np.fromfile(f, dtype='float',count=2) [nmax,lmax] = np.fromfile(f, dtype=np.uint32,count=2) times[tt] = time0 for nn in range(0,nmax): - coef_array[tt,:,nn] = np.fromfile(f, dtype=np.float,count=lmax*(lmax+2)+1) + coef_array[tt,:,nn] = np.fromfile(f, dtype='float',count=lmax*(lmax+2)+1) #return times,coef_array self.T = times @@ -327,7 +327,7 @@ def read_binary_sl_coefficients(self): for nn in range(0,D['nmax']): - coef_array[tt,:,nn] = np.fromfile(f, dtype=np.float,count=D['lmax']*(D['lmax']+2)+1) + coef_array[tt,:,nn] = np.fromfile(f, dtype='float',count=D['lmax']*(D['lmax']+2)+1) tt += 1 f.close() @@ -399,10 +399,10 @@ def read_binary_eof_coefficients(self): for mm in range(0,D['mmax']+1): - coef_array[tt,0,mm,:] = np.fromfile(f, dtype=np.float,count=D['nmax']) + coef_array[tt,0,mm,:] = np.fromfile(f, dtype='float',count=D['nmax']) if mm > 0: - coef_array[tt,1,mm,:] = np.fromfile(f, dtype=np.float,count=D['nmax']) + coef_array[tt,1,mm,:] = np.fromfile(f, dtype='float',count=D['nmax']) tt += 1 diff --git a/exptool/observables/visualize.py b/exptool/observables/visualize.py index 39282b7..27f2ec6 100644 --- a/exptool/observables/visualize.py +++ b/exptool/observables/visualize.py @@ -449,7 +449,7 @@ def show_dump(infile,comp,nout=-1,type='pos',transform=True,\ ax2.xaxis.labelpad = 18 # XZ - ax3.contourf(kdeXZx,kdeXZz,XZ,levels_edge,cmap=cwheel) + ax3.contourf(kdeXZx,-kdeXZz,XZ,levels_edge,cmap=cwheel) ax3.axis([-0.95*face_extents,0.95*face_extents,-0.95*edge_extents,0.95*edge_extents]) ax3.set_xticklabels(()) for label in ax3.get_yticklabels(): diff --git a/exptool/utils/kde_3d.py b/exptool/utils/kde_3d.py index 7c4ad8b..6dbee85 100644 --- a/exptool/utils/kde_3d.py +++ b/exptool/utils/kde_3d.py @@ -146,9 +146,9 @@ def fast_kde(x, y, z, gridsize=(200, 200, 200), extents=None, nocorrelation=Fals inv_cov = np.linalg.inv(cov * scotts_factor**3.) # x & y (pixel) coords of the kernel grid, with = <0,0> in center - xx = np.arange(kern_nx, dtype=np.float) - kern_nx / 2.0 - yy = np.arange(kern_ny, dtype=np.float) - kern_ny / 2.0 - zz = np.arange(kern_nz, dtype=np.float) - kern_nz / 2.0 + xx = np.arange(kern_nx, dtype=float) - kern_nx / 2.0 + yy = np.arange(kern_ny, dtype=float) - kern_ny / 2.0 + zz = np.arange(kern_nz, dtype=float) - kern_nz / 2.0 xx, yy, zz = np.meshgrid(xx, yy, zz) # Then evaluate the gaussian function on the kernel grid @@ -314,8 +314,8 @@ def fast_kde_two(x, y, gridsize=(200, 200), extents=None, nocorrelation=False, w inv_cov = np.linalg.inv(cov * scotts_factor**2) # x & y (pixel) coords of the kernel grid, with = <0,0> in center - xx = np.arange(kern_nx, dtype=np.float) - kern_nx / 2.0 - yy = np.arange(kern_ny, dtype=np.float) - kern_ny / 2.0 + xx = np.arange(kern_nx, dtype=float) - kern_nx / 2.0 + yy = np.arange(kern_ny, dtype=float) - kern_ny / 2.0 xx, yy = np.meshgrid(xx, yy) kernel = np.vstack((xx.flatten(), yy.flatten())) From 39bafbced839c8c3fdfa09f71d237c6ed846ce8b Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Wed, 18 Sep 2024 13:52:19 +0100 Subject: [PATCH 3/3] convert to pyproject.toml; add first stab at CI --- .github/workflows/continuousintegration.yaml | 27 +++++++++++ LICENSE | 2 +- MANIFEST.in | 3 -- pyproject.toml | 15 ++++++ setup.py | 49 -------------------- test_requirements.txt | 3 ++ 6 files changed, 46 insertions(+), 53 deletions(-) create mode 100644 .github/workflows/continuousintegration.yaml delete mode 100644 MANIFEST.in create mode 100644 pyproject.toml delete mode 100644 setup.py create mode 100644 test_requirements.txt diff --git a/.github/workflows/continuousintegration.yaml b/.github/workflows/continuousintegration.yaml new file mode 100644 index 0000000..8b183cb --- /dev/null +++ b/.github/workflows/continuousintegration.yaml @@ -0,0 +1,27 @@ +name: Continuous Integration + +on: + - push + - pull_request + +jobs: + test: + runs-on: ${{ matrix.os }} + strategy: + matrix: + os: [ubuntu-latest, windows-latest, macos-latest] + python-version: ['3.10', '3.11', '3.12'] + + steps: + - uses: actions/checkout@v3 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v4 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + pip install -r test_requirements.txt +# - name: Test with pytest +# run: | +# coverage run -m pytest -v -s diff --git a/LICENSE b/LICENSE index bfd531b..818daf0 100644 --- a/LICENSE +++ b/LICENSE @@ -1,4 +1,4 @@ -Copyright (c) 2016, Michael Petersen +Copyright (c) 2016-2024, Michael Petersen All rights reserved. Redistribution and use in source and binary forms, with or without diff --git a/MANIFEST.in b/MANIFEST.in deleted file mode 100644 index 48c9770..0000000 --- a/MANIFEST.in +++ /dev/null @@ -1,3 +0,0 @@ -include README.md LICENSE -include exptool/io/*.dat -include exptool/analysis/*.dat diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..d7fc5fc --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,15 @@ +[build-system] +requires = ["setuptools>=42", "wheel"] +build-backend = "setuptools.build_meta" + +[project] +name = "exptool" +version = "0.4.0" +description = "exp analysis in Python" +authors = [{ name = "Michael Petersen", email = "michael.petersen@roe.ac.uk" }] +license = { text = "New BSD" } +dependencies = [ + "numpy>=1.7", + "scipy", + "matplotlib", +] diff --git a/setup.py b/setup.py deleted file mode 100644 index c1ae5a1..0000000 --- a/setup.py +++ /dev/null @@ -1,49 +0,0 @@ -#https://docs.python.org/2/distutils/setupscript.html - -from setuptools import setup -from distutils.core import Extension -import sys -import sysconfig -import os, os.path -import subprocess -import glob - -# just in case ... -#https://github.com/blog/2104-working-with-submodules - -#accumulate C extension -accumulate_c_src= glob.glob('exptool/basis/accumulate_c/*.c') - -# eventually put numpy header for arrays in here -accumulate_include_dirs= ['exptool/basis/accumulate_c'] - -accumulate_c= Extension('exptool.basis._accumulate_c', - sources=accumulate_c_src, - #libraries=pot_libraries, - include_dirs=accumulate_include_dirs) - #extra_compile_args=extra_compile_args, - #extra_link_args=extra_link_args) - -ext_modules = [] -ext_modules.append(accumulate_c) - - -setup(name='exptool', - version='0.3.2', - description='exp analysis in Python', - author='Michael Petersen', - author_email='michael.petersen@roe.ac.uk', - license='New BSD', - #long_description=long_description, - url='https://github.com/michael-petersen/exptool', - package_dir = {'exptool/': ''}, - packages=['exptool','exptool/analysis','exptool/basis', - 'exptool/io','exptool/orbits','exptool/utils', - 'exptool/observables','exptool/models','exptool/tests', - 'exptool/perturbation'], - package_data={'': ['README.md','LICENSE'],'exptool/tests':['*.dat'],'exptool/tests/data':['*.dat']}, - include_package_data=True, - install_requires=['numpy>=1.7','scipy','matplotlib'], - ext_modules=ext_modules, - classifiers = ["Programming Language :: Python", "Intended Audience :: Science/Research"] - ) diff --git a/test_requirements.txt b/test_requirements.txt new file mode 100644 index 0000000..7b500b3 --- /dev/null +++ b/test_requirements.txt @@ -0,0 +1,3 @@ +numpy==1.26.4 +scipy==1.11.4 +-e . \ No newline at end of file