From c66eb5e2446f775ab6dd3501d002fa5459b45ae9 Mon Sep 17 00:00:00 2001 From: jngaravitoc Date: Sat, 11 Nov 2023 13:28:43 -0500 Subject: [PATCH 01/19] basis builder routines --- EXP_tools/basis_builder/__init__.py | 3 + EXP_tools/basis_builder/basis_utils.py | 285 +++++++++++++++++++++++++ EXP_tools/basis_builder/makemodel.py | 99 +++++++++ EXP_tools/basis_builder/profiles.py | 127 +++++++++++ 4 files changed, 514 insertions(+) create mode 100644 EXP_tools/basis_builder/__init__.py create mode 100644 EXP_tools/basis_builder/basis_utils.py create mode 100644 EXP_tools/basis_builder/makemodel.py create mode 100644 EXP_tools/basis_builder/profiles.py diff --git a/EXP_tools/basis_builder/__init__.py b/EXP_tools/basis_builder/__init__.py new file mode 100644 index 0000000..47cd4c4 --- /dev/null +++ b/EXP_tools/basis_builder/__init__.py @@ -0,0 +1,3 @@ +from .makemodel import makemodel +from .profiles import Profiles +from .basis_utils import makebasis diff --git a/EXP_tools/basis_builder/basis_utils.py b/EXP_tools/basis_builder/basis_utils.py new file mode 100644 index 0000000..49f3395 --- /dev/null +++ b/EXP_tools/basis_builder/basis_utils.py @@ -0,0 +1,285 @@ +import os, sys, pickle, pyEXP +import numpy as np +import makemodel + +def make_config(basis_id, numr, rmin, rmax, lmax, nmax, scale, + modelname='', cachename='.slgrid_sph_cache'): + """ + Creates a configuration file required to build a basis model. + + Args: + basis_id (str): The identity of the basis model. + numr (int): The number of radial grid points. + rmin (float): The minimum radius value. + rmax (float): The maximum radius value. + lmax (int): The maximum l value of the basis. + nmax (int): The maximum n value of the basis. + scale (float): Scaling factor for the basis. + modelname (str, optional): Name of the model. Default is an empty string. + cachename (str, optional): Name of the cache file. Default is '.slgrid_sph_cache'. + + Returns: + str: A string representation of the configuration file. + + Raises: + None + """ + + config = '\n---\nid: {:s}\n'.format(basis_id) + config += 'parameters:\n' + config += ' numr: {:d}\n'.format(numr) + config += ' rmin: {:.7f}\n'.format(rmin) + config += ' rmax: {:.3f}\n'.format(rmax) + config += ' Lmax: {:d}\n'.format(lmax) + config += ' nmax: {:d}\n'.format(nmax) + config += ' scale: {:.3f}\n'.format(scale) + config += ' modelname: {}\n'.format(modelname) + config += ' cachename: {}\n'.format(cachename) + config += '...\n' + return config + +def empirical_density_profile(pos, mass, nbins=500, rmin=0, rmax=600, log_space=False): + """ + Computes the number density radial profile assuming all particles have the same mass. + + Args: + pos (ndarray): array of particle positions in cartesian coordinates with shape (n,3). + mass (ndarray): array of particle masses with shape (n,). + nbins (int, optional): number of bins in the radial profile. Default is 500. + rmin (float, optional): minimum radius of the radial profile. Default is 0. + rmax (float, optional): maximum radius of the radial profile. Default is 600. + log_space (bool, optional): whether to use logarithmic binning. Default is False. + + Returns: + tuple: a tuple containing the arrays of radius and density with shapes (nbins,) and (nbins,), respectively. + + Raises: + ValueError: if pos and mass arrays have different lengths or if nbins is not a positive integer. + """ + if len(pos) != len(mass): + raise ValueError("pos and mass arrays must have the same length") + if not isinstance(nbins, int) or nbins <= 0: + raise ValueError("nbins must be a positive integer") + + # Compute radial distances + r_p = np.sqrt(np.sum(pos**2, axis=1)) + + # Compute bin edges and shell volumes + if log_space: + bins = np.logspace(np.log10(rmin), np.log10(rmax), nbins+1) + else: + bins = np.linspace(rmin, rmax, nbins+1) + V_shells = 4/3 * np.pi * (bins[1:]**3 - bins[:-1]**3) + + # Compute density profile + density, _ = np.histogram(r_p, bins=bins, weights=mass) + density /= V_shells + + # Compute bin centers and return profile + radius = 0.5 * (bins[1:] + bins[:-1]) + return radius, density + +def make_exp_basis_table(radius, density, outfile='', plabel='', + verbose=True, physical_units=True, return_values=False): + + """ + Create a table of basis functions for an exponential density profile, compatible with EXP. + + Parameters: + ----------- + radius : array-like + Array of sampled radius points. + density : array-like + Array of density values at radius points. + outfile : str, optional + The name of the output file. If not provided, the file will not be written. + plabel : str, optional + A comment string to add to the output file. + verbose : bool, optional + Whether to print scaling factors and other details during execution. + physical_units : bool, optional + Whether to use physical units (True) or normalized units (False) in the output file. + return_values : bool, optional + Whether to return the radius, density, mass, and potential arrays (True) or not (False). + + Returns: + -------- + If return_values is True: + radius : array-like + The radius values. + density : array-like + The density values. + mass : array-like + The mass enclosed at each radius. + potential : array-like + The potential energy at each radius. + If return_values is False (default): + None + + Notes: + ------ + This function assumes an exponential density profile: + + rho(r) = rho_0 * exp(-r/a) + + where rho_0 and a are constants. + + The table of basis functions is used by EXP to create a density profile. + + Reference: + ---------- + https://gist.github.com/michael-petersen/ec4f20641eedac8f63ec409c9cc65ed7 + """ + + + # make the mass and potential arrays + rvals, dvals = radius, density + mvals = np.zeros(density.size) + pvals = np.zeros(density.size) + pwvals = np.zeros(density.size) + + M = 1. + R = np.nanmax(rvals) + + + # initialise the mass enclosed an potential energy + mvals[0] = 1.e-15 + pwvals[0] = 0. + + # evaluate mass enclosed and potential energy by recursion + for indx in range(1, dvals.size): + mvals[indx] = mvals[indx-1] +\ + 2.0*np.pi*(rvals[indx-1]*rvals[indx-1]*dvals[indx-1] +\ + rvals[indx]*rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); + pwvals[indx] = pwvals[indx-1] + \ + 2.0*np.pi*(rvals[indx-1]*dvals[indx-1] + rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); + + # evaluate potential (see theory document) + pvals = -mvals/(rvals+1.e-10) - (pwvals[dvals.size-1] - pwvals) + + # get the maximum mass and maximum radius + M0 = mvals[dvals.size-1] + R0 = rvals[dvals.size-1] + + # compute scaling factors + Beta = (M/M0) * (R0/R) + Gamma = np.sqrt( (M0*R0)/(M*R) ) * (R0/R) + if verbose: + print("! Scaling: R=",R," M=",M) + + rfac = np.power(Beta,-0.25) * np.power(Gamma,-0.5); + dfac = np.power(Beta,1.5) * Gamma; + mfac = np.power(Beta,0.75) * np.power(Gamma,-0.5); + pfac = Beta; + + if verbose: + print(rfac,dfac,mfac,pfac) + + # save file if desired + if outfile != '': + f = open(outfile,'w') + print('! ', plabel, file=f) + print('! R D M P',file=f) + + print(rvals.size,file=f) + + if physical_units: + for indx in range(0,rvals.size): + print('{0} {1} {2} {3}'.format(rvals[indx],\ + dvals[indx],\ + mvals[indx],\ + pvals[indx]),file=f) + else: + for indx in range(0,rvals.size): + print('{0} {1} {2} {3}'.format(rfac*rvals[indx],\ + dfac*dvals[indx],\ + mfac*mvals[indx],\ + pfac*pvals[indx]),file=f) + + f.close() + + if return_values: + if physical_units: + return rvals, dvals, mvals, pvals + + return rvals*rfac, dfac*dvals, mfac*mvals, pfac*pvals + + +def makebasis(pos, mass, basis_model, config=None, basis_id='sphereSL', time=0, + numr=500, rmin=0.61, rmax=599, lmax=4, nmax=20, scale=22.5, + modelname='dens_table.txt', cachename='.slgrid_sph_cache', add_coef = False, coef_file=''): + """ + Create a BFE expansion for a given set of particle positions and masses. + + Parameters: + pos (numpy.ndarray): The positions of particles. Each row represents one particle, + and each column represents the coordinate of that particle. + mass (numpy.ndarray): The masses of particles. The length of this array should be the same + as the number of particles. + basismodel (): + config (pyEXP.config.Config, optional): A configuration object that specifies the basis set. + If not provided, an empirical density profile will be computed + and a configuration object will be created automatically. + basis_id (str, optional): The type of basis set to be used. Default is 'sphereSL'. + time (float, optional): The time at which the expansion is being computed. Default is 0. + numr (int, optional): The number of radial grid points in the basis set. Default is 200. + rmin (float, optional): The minimum radius of the basis set. Default is 0.61. + rmax (float, optional): The maximum radius of the basis set. Default is 599. + lmax (int, optional): The maximum harmonic order in the basis set. Default is 4. + nmax (int, optional): The maximum number of polynomials in the basis set. Default is 20. + scale (float, optional): The scale of the basis set in physical units. + modelname (str, optional): The name of the file containing the density profile model. + Default is 'dens_table.txt'. + cachename (str, optional): The name of the file that will be used to cache the basis set. + Default is '.slgrid_sph_cache'. + save_dir (str, optional): The name of the file if provided that will be used to save the coef files as .h5. + Default is ''. + Returns: + tuple: A tuple containing the basis and the coefficients of the expansion. + The basis is an instance of pyEXP.basis.Basis, and the coefficients are + an instance of pyEXP.coefs.Coefs. + """ + + if os.path.isfile(modelname) == False: + print("-> File model not found so we are computing one \n") + if empirical == True: + print('-> Computing empirical model') + rad, rho = empirical_density_profile(pos, mass, nbins=numr) + elif empirical == False: + makemodel.hernquist_halo() + + makemodel.Profiles(density_profile) + + R, D, M, P = makemodel.makemodel(rvals=np.logspace(np.log10(rmin), np.log10(rmax), numr), rho, outfile=modelname, return_values=True) + print('-> Computing analytical Hernquist model') + R, D, M, P = makemodel.makemodel(hernquist_halo, 1, [scale], rvals=, pfile=modelname) + print('-> Model computed: rmin={}, rmax={}, numr={}'.format(R[0], R[-1], len(R))) + else: + R, D, M, P = np.loadtxt(modelname, skiprows=3, unpack=True) + # check if config file is passed + if config is None: + print('No config file provided.') + print(f'Computing empirical density') + #rad, rho = empirical_density_profile(pos, mass, nbins=500) + #makemodel_empirical(r_exact, rho, outfile=modelname) + #R = [0.01, 600] + config = make_config(basis_id, numr, R[0], R[-1], lmax, nmax, scale, + modelname, cachename) + + # Construct the basis instances + basis = pyEXP.basis.Basis.factory(config) + + # Prints info from Cache + basis.cacheInfo(cachename) + + #compute coefficients + coef = basis.createFromArray(mass/np.sum(mass), pos, time=time) + coefs = pyEXP.coefs.Coefs.makecoefs(coef, 'dark halo') + coefs.add(coef) + if add_coef == False: + coefs.WriteH5Coefs(coef_file) + elif add_coef == True: + coefs.ExtendH5Coefs(coef_file) + + return basis, coefs + diff --git a/EXP_tools/basis_builder/makemodel.py b/EXP_tools/basis_builder/makemodel.py new file mode 100644 index 0000000..5c1aa31 --- /dev/null +++ b/EXP_tools/basis_builder/makemodel.py @@ -0,0 +1,99 @@ +""" +TODO: +""" + +import numpy as np + + +def makemodel(dvals, rvals, M, unit_physical=False, + pfile='', plabel = '', verbose=True): + """ + Make an EXP-compatible spherical basis function table + + inputs + ------------- + basis_type : (str) + func : (function) the callable functional form of the density + M : (float) the total mass of the model, sets normalisations + funcargs : (list) a list of arguments for the density function. + rvals : (array of floats) radius values to evaluate the density function + = 10.**np.linspace(-2.,4.,2000) + pfile : (string) the name of the output file. If '', will not print file + plabel : (string) comment string + verbose : (boolean) + + outputs + ------------- + R : (array of floats) the radius values + D : (array of floats) the density + M : (array of floats) the mass enclosed + P : (array of floats) the potential + + """ + + R = np.nanmax(rvals) + + # query out the density values + + # make the mass and potential arrays + mvals = np.zeros(dvals.size) + pvals = np.zeros(dvals.size) + pwvals = np.zeros(dvals.size) + + # initialise the mass enclosed an potential energy + mvals[0] = 1.e-15 + pwvals[0] = 0. + + # evaluate mass enclosed and potential energy by recursion + for indx in range(1,dvals.size): + mvals[indx] = mvals[indx-1] +\ + 2.0*np.pi*(rvals[indx-1]*rvals[indx-1]*dvals[indx-1] +\ + rvals[indx]*rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); + pwvals[indx] = pwvals[indx-1] + \ + 2.0*np.pi*(rvals[indx-1]*dvals[indx-1] + rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); + + # evaluate potential (see theory document) + pvals = -mvals/(rvals+1.e-10) - (pwvals[dvals.size-1] - pwvals) + + # get the maximum mass and maximum radius + M0 = mvals[dvals.size-1] + R0 = rvals[dvals.size-1] + + # compute scaling factors + Beta = (M/M0) * (R0/R); + Gamma = np.sqrt((M0*R0)/(M*R)) * (R0/R); + if verbose: + print("! Scaling: R=",R," M=",M) + + rfac = np.power(Beta,-0.25) * np.power(Gamma,-0.5); + dfac = np.power(Beta,1.5) * Gamma; + mfac = np.power(Beta,0.75) * np.power(Gamma,-0.5); + pfac = Beta; + + if verbose: + print(rfac,dfac,mfac,pfac) + + # save file if desired + if pfile != '': + f = open(pfile,'w') + print('! ',plabel,file=f) + print('! R D M P',file=f) + + print(rvals.size,file=f) + + if unit_physical = True: + for indx in range(0,rvals.size): + print('{0} {1} {2} {3}'.format(rvals[indx],\ + dvals[indx],\ + mvals[indx],\ + pvals[indx]),file=f) + else: + + f.close() + + if return_val: + if unit_phy: + return rvals, dvals, mvals, pvals + else: + return rvals*rfac, dfac*dvals, mfac*mvals, pfac*pvals + diff --git a/EXP_tools/basis_builder/profiles.py b/EXP_tools/basis_builder/profiles.py new file mode 100644 index 0000000..fe4f7e3 --- /dev/null +++ b/EXP_tools/basis_builder/profiles.py @@ -0,0 +1,127 @@ +import numpy as np +from scipy import special + +Class Profiles: + def __init__(self, radius, scale_radius, alpha=1., beta=2.0): + """ + Class with density profiles of Dark Matter halos + + inputs + ---------- + r : (float) radius values + rs : (float, default=1.) scale radius + alpha : (float, default=1.) inner halo slope + beta : (float, default=1.e-7) outer halo slope + + notes + ---------- + different combinations are known distributions. + alpha=1,beta=2 is NFW (default) + alpha=1,beta=3 is Hernquist + alpha=2.5, beta=0 is a typical single power law halo + + """ + + self.radius = radius + self.scale_radius = scale_radius + self.alpha = alpha + self.beta = beta + self.ra = self.radius/self.scale_radius + + + def powerhalo(self, rc=0.0): + """ + Generic two power law distribution + + inputs + ---------- + rc : (float, default=0. i.e. no core) core radius + + returns + ---------- + densities evaluated at r + + notes + ---------- + different combinations are known distributions. + alpha=1,beta=2 is NFW + alpha=1,beta=3 is Hernquist + alpha=2.5,beta=0 is a typical single power law halo + + + + """ + return 1./(((self.ra+rc/self.scale_radius)**self.alpha)*((1+self.ra)**self.beta)) + + + def powerhalorolloff(self, rc=0.0): + """ + Generic twopower law distribution with an erf rolloff + + + inputs + ---------- + rc : (float, default=0. i.e. no core) core radius + + returns + ---------- + densities evaluated at r + + notes + ---------- + different combinations are known distributions. + alpha=1,beta=2 is NFW + alpha=1,beta=3 is Hernquist + alpha=2.5,beta=0 is a typical single power law halo + + + """ + dens = 1./(((self.ra+rc/self.scale_radius)**self.alpha)*((1+self.ra)**self.beta)) + rtrunc = 25*self.scale_radius + wtrunc = rtrunc*0.2 + rolloff = 0.5 - 0.5*special.erf((self.radius-rtrunc)/wtrunc) + return dens*rolloff + + + def plummer_density(self): + """ + Plummer density profile + + inputs + --------- + + + returns + ---------- + densities evaluated at r + """ + return ((3.0)/(4*np.pi))*(self.scale_radius**2.)*((self.scale_radius**2 + self.radius**2)**(-2.5)) + + def twopower_density_withrolloff(self, rcen, wcen): + """ + inputs + --------- + + + returns + ---------- + densities evaluated at r + """ + + + prefac = 0.5*(1.-scipy.special.erf((ra-rcen)/wcen)) + return prefac*(ra**-self.alpha)*(1+ra)**(-self.beta+self.alpha) + + def hernquist_halo(self): + """ + TODO: DO we need these profile + Hernquist halo + + inputs + --------- + + returns + ---------- + densities evaluated at r + """ + return 1 / ( 2*np.pi * (self.radius/self.scale_radius) * (1 + self.radius/self.scale_radius)**3) \ No newline at end of file From 2592b876067012cdbda853a06f4d2c40f08a1ded Mon Sep 17 00:00:00 2001 From: jngaravitoc Date: Sat, 11 Nov 2023 13:28:59 -0500 Subject: [PATCH 02/19] visuals routines --- EXP_tools/visuals/visualize.py | 304 +++++++++++++++++++++++++++++++++ 1 file changed, 304 insertions(+) create mode 100644 EXP_tools/visuals/visualize.py diff --git a/EXP_tools/visuals/visualize.py b/EXP_tools/visuals/visualize.py new file mode 100644 index 0000000..8db665d --- /dev/null +++ b/EXP_tools/visuals/visualize.py @@ -0,0 +1,304 @@ +import os, sys, pickle, pyEXP +import numpy as np + +###Field computations for plotting### +def make_basis_plot(basis, savefile=None, nsnap='mean', y=0.92, dpi=200): + """ + Plots the potential of the basis functions for different values of l and n. + + Args: + basis (obj): object containing the basis functions for the simulation + savefile (str, optional): name of the file to save the plot as + nsnap (str, optional): description of the snapshot being plotted + y (float, optional): vertical position of the main title + dpi (int, optional): resolution of the plot in dots per inch + + Returns: + None + + """ + # Set up grid for plotting potential + lrmin, lrmax, rnum = 0.5, 2.7, 100 + halo_grid = basis.getBasis(lrmin, lrmax, rnum) + r = np.linspace(lrmin, lrmax, rnum) + r = np.power(10.0, r) + + # Create subplots and plot potential for each l and n + import matplotlib.pyplot as plt + fig, ax = plt.subplots(4, 5, figsize=(6,6), dpi=dpi, + sharex='col', sharey='row') + plt.subplots_adjust(wspace=0, hspace=0) + ax = ax.flatten() + + for l in range(len(ax)): + ax[l].set_title(f"$\ell = {l}$", y=0.8, fontsize=6) + for n in range(20): + ax[l].semilogx(r, halo_grid[l][n]['potential'], '-', label="n={}".format(n), lw=0.5) + + # Add labels and main title + fig.supylabel('Potential', weight='bold', x=-0.02) + fig.supxlabel('Radius', weight='bold', y=0.02) + fig.suptitle(f'nsnap = {nsnap}', + fontsize=12, + weight='bold', + y=y, + ) + + # Save plot if a filename was provided + if savefile: + plt.savefig(f'{savefile}', bbox_inches='tight') + +def find_field(basis, coefficients, time=0, xyz=(0, 0, 0), property='dens', include_monopole=True): + """ + Finds the value of the specified property of the field at the given position. + + Args: + basis (obj): Object containing the basis functions for the simulation. + coefficients (obj): Object containing the coefficients for the simulation. + time (float, optional): The time at which to evaluate the field. Default is 0. + xyz (tuple or list, optional): The (x, y, z) position at which to evaluate the field. Default is (0, 0, 0). + property (str, optional): The property of the field to evaluate. Can be 'dens', 'pot', or 'force'. Default is 'dens'. + include_monopole (bool, optional): Whether to return the monopole contribution to the property only. Default is True. + + Returns: + float or list: The value of the specified property of the field at the given position. If property is 'force', a list of + three values is returned representing the force vector in (x, y, z) directions. + + Raises: + ValueError: If the property argument is not 'dens', 'pot', or 'force'. + """ + + coefficients.set_coefs(coefficients.getCoefStruct(time)) + dens0, pot0, dens, pot, fx, fy, fz = basis.getFields(xyz[0], xyz[1], xyz[2]) + + if property == 'dens': + if include_monopole: + return dens0 + return dens + dens0 + + elif property == 'pot': + if include_monopole: + return pot0 + return pot + pot0 + + elif property == 'force': + return [fx, fy, fz] + + else: + raise ValueError("Invalid property specified. Possible values are 'dens', 'pot', and 'force'.") + +def spherical_avg_prop(basis, coefficients, time=0, radius=np.linspace(0.1, 600, 100), property='dens'): + """ + Computes the spherically averaged value of the specified property of the field over the given radii. + + Args: + basis (obj): Object containing the basis functions for the simulation. + coefficients (obj): Object containing the coefficients for the simulation. + time (float, optional): The time at which to evaluate the field. Default is 0. + radius (ndarray, optional): An array of radii over which to compute the spherically averaged property. Default is an + array of 100 values logarithmically spaced between 0.1 and 600. + property (str, optional): The property of the field to evaluate. Can be 'dens', 'pot', or 'force'. Default is 'dens'. + + Returns: + ndarray: An array of spherically averaged values of the specified property over the given radii. + + Raises: + ValueError: If the property argument is not 'dens', 'pot', or 'force'. + """ + + coefficients.set_coefs(coefficients.getCoefStruct(time)) + field = [find_field(basis, np.hstack([[rad], [0], [0]]), property=property, include_monopole=True) for rad in radius] + + if property == 'force': + return np.vstack(field), radius + + return np.array(field), radius + +def return_fields_in_grid(basis, coefficients, times=[0], + projection='3D', proj_plane=0, + grid_lim=300, npoints=150,): + + """ + Returns the 2D or 3D field grids as a dict at each time step as a key. + + Args: + basis (obj): object containing the basis functions for the simulation + coefficients (obj): object containing the coefficients for the simulation + times (list): list of the time to compute the field grid. + projection (str): the slice projection to plot. Can be '3D', XY', 'XZ', or 'YZ'. + proj_plane (float, optional): the value of the coordinate that is held constant in the 2D slice projection + npoints (int, optional): the number of grid points in each dimension + grid_limits (float, optional): the limits of the grid in the dimensions. + + Returns: fields along with the specified 2D or 3D grid. The fields are in a + heirarchical dict structure with for each time: + dict_keys(['d', 'd0', 'd1', 'dd', 'fp', 'fr', 'ft', 'p', 'p0', 'p1']) + where + - 'd': Total density at each point. + - 'd0' : Density in the monopole term. + - 'd1' : Density in l>0 terms. + - 'dd': Density difference. + - 'fp': Force in the phi direction. + - 'fr': Force in the radial direction. + - 'ft': Force in the theta direction. + - 'p': Total potential at each point. + - 'p0': Potential in the monopole term. + - 'p1': Potential in l>0 terms. + """ + + assert projection in ['3D', 'XY', 'XZ', 'YZ'], "Invalid value for 'projection'. Must be one of '3D', 'XY', 'XZ', or 'YZ'." + + + + if projection == '3D': + pmin = [-grid_lim, -grid_lim, -grid_lim] + pmax = [grid_lim, grid_lim, grid_lim] + grid = [npoints, npoints, npoints] + + ##create a projection grid along with it + x = np.linspace(-grid_lim, grid_lim, npoints) + xgrid = np.meshgrid(x, x, x) + + else: + x = np.linspace(-grid_lim, grid_lim, npoints) + xgrid = np.meshgrid(x, x) + grid = [npoints, npoints] + + if projection == 'XY': + pmin = [-grid_lim, -grid_lim, proj_plane] + pmax = [grid_lim, grid_lim, proj_plane] + + elif projection == 'XZ': + pmin = [-grid_lim, proj_plane, -grid_lim] + pmax = [grid_lim, proj_plane, grid_lim] + + elif projection == 'YZ': + pmin = [proj_plane, -grid_lim, -grid_lim] + pmax = [proj_plane, grid_lim, grid_lim] + + + field_gen = pyEXP.field.FieldGenerator(times, pmin, pmax, grid) + + return field_gen.volumes(basis, coefs), xgrid + +def slice_fields(basis, coefficients, time=0, + projection='XY', proj_plane=0, npoints=300, + grid_limits=(-300, 300), prop='dens', monopole_only=False): + """ + Plots a slice projection of the fields of a simulation. + + Args: + basis (obj): object containing the basis functions for the simulation + coefficients (obj): object containing the coefficients for the simulation + time (float): the time at which to plot the fields + projection (str): the slice projection to plot. Can be 'XY', 'XZ', or 'YZ'. + proj_plane (float, optional): the value of the coordinate that is held constant in the slice projection + npoints (int, optional): the number of grid points in each dimension + grid_limits (tuple, optional): the limits of the grid in the x and y dimensions, in the form (x_min, x_max) + prop (str, optional): the property to return. Can be 'dens' (density), 'pot' (potential), or 'force' (force). + monopole_only (bool, optional): whether to return the monopole component in the returned property value. + + Returns: + array or list: the property specified by `prop`. If `prop` is 'force', a list of the x, y, and z components of the force is returned. + Also returns the grid used to compute slice fields. + """ + x = np.linspace(grid_limits[0], grid_limits[1], npoints) + xgrid = np.meshgrid(x, x) + xg = xgrid[0].flatten() + yg = xgrid[1].flatten() + + + if projection not in ['XY', 'XZ', 'YZ']: + raise ValueError("Invalid projection specified. Possible values are 'XY', 'XZ', and 'YZ'.") + + N = len(xg) + rho0 = np.zeros_like(xg) + pot0 = np.zeros_like(xg) + rho = np.zeros_like(xg) + pot = np.zeros_like(xg) + basis.set_coefs(coefficients.getCoefStruct(time)) + + for k in range(0, N): + if projection == 'XY': + rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(xg[k], yg[k], proj_plane) + elif projection == 'XZ': + rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(xg[k], proj_plane, yg[k]) + elif projection == 'YZ': + rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(proj_plane, xg[k], yg[k]) + + dens = rho.reshape(npoints, npoints) + pot = pot.reshape(npoints, npoints) + dens0 = rho0.reshape(npoints, npoints) + pot0 = pot0.reshape(npoints, npoints) + + if prop == 'dens': + if monopole_only: + return dens0 + return dens0, dens, xgrid + + if prop == 'pot': + if monopole_only: + return pot0 + return pot0, pot, xgrid + + if prop == 'force': + return [fx.reshape(npoints, npoints), fy.reshape(npoints, npoints), fz.reshape(npoints, npoints)], xgrid + + +def slice_3d_fields(basis, coefficients, time=0, + projection='XY', proj_plane=0, npoints=300, + grid_limits=(-300, 300), prop='dens', monopole_only=False): + """ + Plots a slice projection of the fields of a simulation. + + Args: + basis (obj): object containing the basis functions for the simulation + coefficients (obj): object containing the coefficients for the simulation + time (float): the time at which to plot the fields + projection (str): the slice projection to plot. Can be 'XY', 'XZ', or 'YZ'. + proj_plane (float, optional): the value of the coordinate that is held constant in the slice projection + npoints (int, optional): the number of grid points in each dimension + grid_limits (tuple, optional): the limits of the grid in the x and y dimensions, in the form (x_min, x_max) + prop (str, optional): the property to return. Can be 'dens' (density), 'pot' (potential), or 'force' (force). + monopole_only (bool, optional): whether to return the monopole component in the returned property value. + + Returns: + array or list: the property specified by `prop`. If `prop` is 'force', a list of the x, y, and z components of the force is returned. + Also returns the grid used to compute slice fields. + """ + x = np.linspace(grid_limits[0], grid_limits[1], npoints) + xgrid = np.meshgrid(x, x, x) + xg = xgrid[0].flatten() + yg = xgrid[1].flatten() + zg = xgrid[2].flatten() + + + + N = len(xg) + rho0 = np.zeros_like(xg) + pot0 = np.zeros_like(xg) + rho = np.zeros_like(xg) + pot = np.zeros_like(xg) + basis.set_coefs(coefficients.getCoefStruct(time)) + + for k in range(0, N): + rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(xg[k], yg[k], zg[k]) + + dens = rho.reshape(npoints, npoints, npoints) + pot = pot.reshape(npoints, npoints, npoints) + dens0 = rho0.reshape(npoints, npoints, npoints) + pot0 = pot0.reshape(npoints, npoints, npoints) + + if prop == 'dens': + if monopole_only: + return dens0 + return dens0, dens, xgrid + + if prop == 'pot': + if monopole_only: + return pot0 + return pot0, pot, xgrid + + if prop == 'force': + return [fx.reshape(npoints, npoints, npoints), fy.reshape(npoints, npoints, npoints), fz.reshape(npoints, npoints, npoints)], xgrid + From 420b914274e6386ed2b8c152f29bcefa7bce3df2 Mon Sep 17 00:00:00 2001 From: jngaravitoc Date: Sat, 11 Nov 2023 13:29:42 -0500 Subject: [PATCH 03/19] re-organizing folders --- EXP_tools/makemodel.py | 172 ----------------------- EXP_tools/visualize.py | 304 ----------------------------------------- 2 files changed, 476 deletions(-) delete mode 100644 EXP_tools/makemodel.py delete mode 100644 EXP_tools/visualize.py diff --git a/EXP_tools/makemodel.py b/EXP_tools/makemodel.py deleted file mode 100644 index 31a1ba8..0000000 --- a/EXP_tools/makemodel.py +++ /dev/null @@ -1,172 +0,0 @@ -import numpy as np - -def makemodel(func,M,funcargs,rvals = 10.**np.linspace(-2.,4.,2000),pfile='',plabel = '',verbose=True): - """make an EXP-compatible spherical basis function table - - inputs - ------------- - func : (function) the callable functional form of the density - M : (float) the total mass of the model, sets normalisations - funcargs : (list) a list of arguments for the density function. - rvals : (array of floats) radius values to evaluate the density function - pfile : (string) the name of the output file. If '', will not print file - plabel : (string) comment string - verbose : (boolean) - - outputs - ------------- - R : (array of floats) the radius values - D : (array of floats) the density - M : (array of floats) the mass enclosed - P : (array of floats) the potential - - """ - - R = np.nanmax(rvals) - - # query out the density values - dvals = func(rvals,*funcargs) - - # make the mass and potential arrays - mvals = np.zeros(dvals.size) - pvals = np.zeros(dvals.size) - pwvals = np.zeros(dvals.size) - - # initialise the mass enclosed an potential energy - mvals[0] = 1.e-15 - pwvals[0] = 0. - - # evaluate mass enclosed and potential energy by recursion - for indx in range(1,dvals.size): - mvals[indx] = mvals[indx-1] +\ - 2.0*np.pi*(rvals[indx-1]*rvals[indx-1]*dvals[indx-1] +\ - rvals[indx]*rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); - pwvals[indx] = pwvals[indx-1] + \ - 2.0*np.pi*(rvals[indx-1]*dvals[indx-1] + rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); - - # evaluate potential (see theory document) - pvals = -mvals/(rvals+1.e-10) - (pwvals[dvals.size-1] - pwvals) - - # get the maximum mass and maximum radius - M0 = mvals[dvals.size-1] - R0 = rvals[dvals.size-1] - - # compute scaling factors - Beta = (M/M0) * (R0/R); - Gamma = np.sqrt((M0*R0)/(M*R)) * (R0/R); - if verbose: - print("! Scaling: R=",R," M=",M) - - rfac = np.power(Beta,-0.25) * np.power(Gamma,-0.5); - dfac = np.power(Beta,1.5) * Gamma; - mfac = np.power(Beta,0.75) * np.power(Gamma,-0.5); - pfac = Beta; - - if verbose: - print(rfac,dfac,mfac,pfac) - - # save file if desired - if pfile != '': - f = open(pfile,'w') - print('! ',plabel,file=f) - print('! R D M P',file=f) - - print(rvals.size,file=f) - - for indx in range(0,rvals.size): - print('{0} {1} {2} {3}'.format( rfac*rvals[indx],\ - dfac*dvals[indx],\ - mfac*mvals[indx],\ - pfac*pvals[indx]),file=f) - - f.close() - - return rvals*rfac,dfac*dvals,mfac*mvals,pfac*pvals - - - -def powerhalo(r,rs=1.,rc=0.,alpha=1.,beta=1.e-7): - """return generic twopower law distribution - - inputs - ---------- - r : (float) radius values - rs : (float, default=1.) scale radius - rc : (float, default=0. i.e. no core) core radius - alpha : (float, default=1.) inner halo slope - beta : (float, default=1.e-7) outer halo slope - - returns - ---------- - densities evaluated at r - - notes - ---------- - different combinations are known distributions. - alpha=1,beta=2 is NFW - alpha=1,beta=3 is Hernquist - alpha=2.5,beta=0 is a typical single power law halo - - - """ - ra = r/rs - return 1./(((ra+rc/rs)**alpha)*((1+ra)**beta)) - - -def powerhalorolloff(r,rs=1.,rc=0.,alpha=1.,beta=1.e-7): - """return generic twopower law distribution with an erf rolloff - - inputs - ---------- - r : (float) radius values - rs : (float, default=1.) scale radius - rc : (float, default=0. i.e. no core) core radius - alpha : (float, default=1.) inner halo slope - beta : (float, default=1.e-7) outer halo slope - - returns - ---------- - densities evaluated at r - - notes - ---------- - different combinations are known distributions. - alpha=1,beta=2 is NFW - alpha=1,beta=3 is Hernquist - alpha=2.5,beta=0 is a typical single power law halo - - - """ - ra = r/rs - dens = 1./(((ra+rc/rs)**alpha)*((1+ra)**beta)) - rtrunc = 25*rs - wtrunc = rtrunc*0.2 - rolloff = 0.5 - 0.5*special.erf((r-rtrunc)/wtrunc) - return dens*rolloff - - -def plummer_density(radius,scale_radius=1.0,mass=1.0,astronomicalG=1.0): - """basic plummer density profile""" - return ((3.0*mass)/(4*np.pi))*(scale_radius**2.)*((scale_radius**2 + radius**2)**(-2.5)) - -def twopower_density_withrolloff(r,a,alpha,beta,rcen,wcen): - """a twopower density profile""" - ra = r/a - prefac = 0.5*(1.-scipy.special.erf((ra-rcen)/wcen)) - return prefac*(ra**-alpha)*(1+ra)**(-beta+alpha) - -def hernquist_halo(r, a): - return 1 / ( 2*np.pi * (r/a) * (1 + r/a)**3) - - -""" -plummer_b = 1.0 -R,D,M,P = makemodel(plummer_density,1.,[plummer_b],rvals = 10.**np.linspace(-3.,1.,2000)) - -alpha=1 -beta=2 -concentration = 15. -rs = 1/concentration -rc = 0.0001 -R,D,M,P = makemodel(powerhalo,1.,[rs,rc,alpha,beta],rvals = 10.**np.linspace(-5.,0.3,2000)pfile='NFWModelc{}.txt'.format(concentration),plabel = 'NFWModelc{}.txt'.format(concentration)) -""" \ No newline at end of file diff --git a/EXP_tools/visualize.py b/EXP_tools/visualize.py deleted file mode 100644 index 8db665d..0000000 --- a/EXP_tools/visualize.py +++ /dev/null @@ -1,304 +0,0 @@ -import os, sys, pickle, pyEXP -import numpy as np - -###Field computations for plotting### -def make_basis_plot(basis, savefile=None, nsnap='mean', y=0.92, dpi=200): - """ - Plots the potential of the basis functions for different values of l and n. - - Args: - basis (obj): object containing the basis functions for the simulation - savefile (str, optional): name of the file to save the plot as - nsnap (str, optional): description of the snapshot being plotted - y (float, optional): vertical position of the main title - dpi (int, optional): resolution of the plot in dots per inch - - Returns: - None - - """ - # Set up grid for plotting potential - lrmin, lrmax, rnum = 0.5, 2.7, 100 - halo_grid = basis.getBasis(lrmin, lrmax, rnum) - r = np.linspace(lrmin, lrmax, rnum) - r = np.power(10.0, r) - - # Create subplots and plot potential for each l and n - import matplotlib.pyplot as plt - fig, ax = plt.subplots(4, 5, figsize=(6,6), dpi=dpi, - sharex='col', sharey='row') - plt.subplots_adjust(wspace=0, hspace=0) - ax = ax.flatten() - - for l in range(len(ax)): - ax[l].set_title(f"$\ell = {l}$", y=0.8, fontsize=6) - for n in range(20): - ax[l].semilogx(r, halo_grid[l][n]['potential'], '-', label="n={}".format(n), lw=0.5) - - # Add labels and main title - fig.supylabel('Potential', weight='bold', x=-0.02) - fig.supxlabel('Radius', weight='bold', y=0.02) - fig.suptitle(f'nsnap = {nsnap}', - fontsize=12, - weight='bold', - y=y, - ) - - # Save plot if a filename was provided - if savefile: - plt.savefig(f'{savefile}', bbox_inches='tight') - -def find_field(basis, coefficients, time=0, xyz=(0, 0, 0), property='dens', include_monopole=True): - """ - Finds the value of the specified property of the field at the given position. - - Args: - basis (obj): Object containing the basis functions for the simulation. - coefficients (obj): Object containing the coefficients for the simulation. - time (float, optional): The time at which to evaluate the field. Default is 0. - xyz (tuple or list, optional): The (x, y, z) position at which to evaluate the field. Default is (0, 0, 0). - property (str, optional): The property of the field to evaluate. Can be 'dens', 'pot', or 'force'. Default is 'dens'. - include_monopole (bool, optional): Whether to return the monopole contribution to the property only. Default is True. - - Returns: - float or list: The value of the specified property of the field at the given position. If property is 'force', a list of - three values is returned representing the force vector in (x, y, z) directions. - - Raises: - ValueError: If the property argument is not 'dens', 'pot', or 'force'. - """ - - coefficients.set_coefs(coefficients.getCoefStruct(time)) - dens0, pot0, dens, pot, fx, fy, fz = basis.getFields(xyz[0], xyz[1], xyz[2]) - - if property == 'dens': - if include_monopole: - return dens0 - return dens + dens0 - - elif property == 'pot': - if include_monopole: - return pot0 - return pot + pot0 - - elif property == 'force': - return [fx, fy, fz] - - else: - raise ValueError("Invalid property specified. Possible values are 'dens', 'pot', and 'force'.") - -def spherical_avg_prop(basis, coefficients, time=0, radius=np.linspace(0.1, 600, 100), property='dens'): - """ - Computes the spherically averaged value of the specified property of the field over the given radii. - - Args: - basis (obj): Object containing the basis functions for the simulation. - coefficients (obj): Object containing the coefficients for the simulation. - time (float, optional): The time at which to evaluate the field. Default is 0. - radius (ndarray, optional): An array of radii over which to compute the spherically averaged property. Default is an - array of 100 values logarithmically spaced between 0.1 and 600. - property (str, optional): The property of the field to evaluate. Can be 'dens', 'pot', or 'force'. Default is 'dens'. - - Returns: - ndarray: An array of spherically averaged values of the specified property over the given radii. - - Raises: - ValueError: If the property argument is not 'dens', 'pot', or 'force'. - """ - - coefficients.set_coefs(coefficients.getCoefStruct(time)) - field = [find_field(basis, np.hstack([[rad], [0], [0]]), property=property, include_monopole=True) for rad in radius] - - if property == 'force': - return np.vstack(field), radius - - return np.array(field), radius - -def return_fields_in_grid(basis, coefficients, times=[0], - projection='3D', proj_plane=0, - grid_lim=300, npoints=150,): - - """ - Returns the 2D or 3D field grids as a dict at each time step as a key. - - Args: - basis (obj): object containing the basis functions for the simulation - coefficients (obj): object containing the coefficients for the simulation - times (list): list of the time to compute the field grid. - projection (str): the slice projection to plot. Can be '3D', XY', 'XZ', or 'YZ'. - proj_plane (float, optional): the value of the coordinate that is held constant in the 2D slice projection - npoints (int, optional): the number of grid points in each dimension - grid_limits (float, optional): the limits of the grid in the dimensions. - - Returns: fields along with the specified 2D or 3D grid. The fields are in a - heirarchical dict structure with for each time: - dict_keys(['d', 'd0', 'd1', 'dd', 'fp', 'fr', 'ft', 'p', 'p0', 'p1']) - where - - 'd': Total density at each point. - - 'd0' : Density in the monopole term. - - 'd1' : Density in l>0 terms. - - 'dd': Density difference. - - 'fp': Force in the phi direction. - - 'fr': Force in the radial direction. - - 'ft': Force in the theta direction. - - 'p': Total potential at each point. - - 'p0': Potential in the monopole term. - - 'p1': Potential in l>0 terms. - """ - - assert projection in ['3D', 'XY', 'XZ', 'YZ'], "Invalid value for 'projection'. Must be one of '3D', 'XY', 'XZ', or 'YZ'." - - - - if projection == '3D': - pmin = [-grid_lim, -grid_lim, -grid_lim] - pmax = [grid_lim, grid_lim, grid_lim] - grid = [npoints, npoints, npoints] - - ##create a projection grid along with it - x = np.linspace(-grid_lim, grid_lim, npoints) - xgrid = np.meshgrid(x, x, x) - - else: - x = np.linspace(-grid_lim, grid_lim, npoints) - xgrid = np.meshgrid(x, x) - grid = [npoints, npoints] - - if projection == 'XY': - pmin = [-grid_lim, -grid_lim, proj_plane] - pmax = [grid_lim, grid_lim, proj_plane] - - elif projection == 'XZ': - pmin = [-grid_lim, proj_plane, -grid_lim] - pmax = [grid_lim, proj_plane, grid_lim] - - elif projection == 'YZ': - pmin = [proj_plane, -grid_lim, -grid_lim] - pmax = [proj_plane, grid_lim, grid_lim] - - - field_gen = pyEXP.field.FieldGenerator(times, pmin, pmax, grid) - - return field_gen.volumes(basis, coefs), xgrid - -def slice_fields(basis, coefficients, time=0, - projection='XY', proj_plane=0, npoints=300, - grid_limits=(-300, 300), prop='dens', monopole_only=False): - """ - Plots a slice projection of the fields of a simulation. - - Args: - basis (obj): object containing the basis functions for the simulation - coefficients (obj): object containing the coefficients for the simulation - time (float): the time at which to plot the fields - projection (str): the slice projection to plot. Can be 'XY', 'XZ', or 'YZ'. - proj_plane (float, optional): the value of the coordinate that is held constant in the slice projection - npoints (int, optional): the number of grid points in each dimension - grid_limits (tuple, optional): the limits of the grid in the x and y dimensions, in the form (x_min, x_max) - prop (str, optional): the property to return. Can be 'dens' (density), 'pot' (potential), or 'force' (force). - monopole_only (bool, optional): whether to return the monopole component in the returned property value. - - Returns: - array or list: the property specified by `prop`. If `prop` is 'force', a list of the x, y, and z components of the force is returned. - Also returns the grid used to compute slice fields. - """ - x = np.linspace(grid_limits[0], grid_limits[1], npoints) - xgrid = np.meshgrid(x, x) - xg = xgrid[0].flatten() - yg = xgrid[1].flatten() - - - if projection not in ['XY', 'XZ', 'YZ']: - raise ValueError("Invalid projection specified. Possible values are 'XY', 'XZ', and 'YZ'.") - - N = len(xg) - rho0 = np.zeros_like(xg) - pot0 = np.zeros_like(xg) - rho = np.zeros_like(xg) - pot = np.zeros_like(xg) - basis.set_coefs(coefficients.getCoefStruct(time)) - - for k in range(0, N): - if projection == 'XY': - rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(xg[k], yg[k], proj_plane) - elif projection == 'XZ': - rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(xg[k], proj_plane, yg[k]) - elif projection == 'YZ': - rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(proj_plane, xg[k], yg[k]) - - dens = rho.reshape(npoints, npoints) - pot = pot.reshape(npoints, npoints) - dens0 = rho0.reshape(npoints, npoints) - pot0 = pot0.reshape(npoints, npoints) - - if prop == 'dens': - if monopole_only: - return dens0 - return dens0, dens, xgrid - - if prop == 'pot': - if monopole_only: - return pot0 - return pot0, pot, xgrid - - if prop == 'force': - return [fx.reshape(npoints, npoints), fy.reshape(npoints, npoints), fz.reshape(npoints, npoints)], xgrid - - -def slice_3d_fields(basis, coefficients, time=0, - projection='XY', proj_plane=0, npoints=300, - grid_limits=(-300, 300), prop='dens', monopole_only=False): - """ - Plots a slice projection of the fields of a simulation. - - Args: - basis (obj): object containing the basis functions for the simulation - coefficients (obj): object containing the coefficients for the simulation - time (float): the time at which to plot the fields - projection (str): the slice projection to plot. Can be 'XY', 'XZ', or 'YZ'. - proj_plane (float, optional): the value of the coordinate that is held constant in the slice projection - npoints (int, optional): the number of grid points in each dimension - grid_limits (tuple, optional): the limits of the grid in the x and y dimensions, in the form (x_min, x_max) - prop (str, optional): the property to return. Can be 'dens' (density), 'pot' (potential), or 'force' (force). - monopole_only (bool, optional): whether to return the monopole component in the returned property value. - - Returns: - array or list: the property specified by `prop`. If `prop` is 'force', a list of the x, y, and z components of the force is returned. - Also returns the grid used to compute slice fields. - """ - x = np.linspace(grid_limits[0], grid_limits[1], npoints) - xgrid = np.meshgrid(x, x, x) - xg = xgrid[0].flatten() - yg = xgrid[1].flatten() - zg = xgrid[2].flatten() - - - - N = len(xg) - rho0 = np.zeros_like(xg) - pot0 = np.zeros_like(xg) - rho = np.zeros_like(xg) - pot = np.zeros_like(xg) - basis.set_coefs(coefficients.getCoefStruct(time)) - - for k in range(0, N): - rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(xg[k], yg[k], zg[k]) - - dens = rho.reshape(npoints, npoints, npoints) - pot = pot.reshape(npoints, npoints, npoints) - dens0 = rho0.reshape(npoints, npoints, npoints) - pot0 = pot0.reshape(npoints, npoints, npoints) - - if prop == 'dens': - if monopole_only: - return dens0 - return dens0, dens, xgrid - - if prop == 'pot': - if monopole_only: - return pot0 - return pot0, pot, xgrid - - if prop == 'force': - return [fx.reshape(npoints, npoints, npoints), fy.reshape(npoints, npoints, npoints), fz.reshape(npoints, npoints, npoints)], xgrid - From 2ca031bd6c2856dd5a77430bfbf4de91ca8bb641 Mon Sep 17 00:00:00 2001 From: jngaravitoc Date: Sat, 11 Nov 2023 13:31:06 -0500 Subject: [PATCH 04/19] rm old basis utils --- EXP_tools/basis_utils.py | 281 --------------------------------------- 1 file changed, 281 deletions(-) delete mode 100644 EXP_tools/basis_utils.py diff --git a/EXP_tools/basis_utils.py b/EXP_tools/basis_utils.py deleted file mode 100644 index f1d70d5..0000000 --- a/EXP_tools/basis_utils.py +++ /dev/null @@ -1,281 +0,0 @@ -import os, sys, pickle, pyEXP -import numpy as np -import makemodel -from makemodel import hernquist_halo - -def make_config(basis_id, numr, rmin, rmax, lmax, nmax, scale, - modelname='', cachename='.slgrid_sph_cache'): - """ - Creates a configuration file required to build a basis model. - - Args: - basis_id (str): The identity of the basis model. - numr (int): The number of radial grid points. - rmin (float): The minimum radius value. - rmax (float): The maximum radius value. - lmax (int): The maximum l value of the basis. - nmax (int): The maximum n value of the basis. - scale (float): Scaling factor for the basis. - modelname (str, optional): Name of the model. Default is an empty string. - cachename (str, optional): Name of the cache file. Default is '.slgrid_sph_cache'. - - Returns: - str: A string representation of the configuration file. - - Raises: - None - """ - - config = '\n---\nid: {:s}\n'.format(basis_id) - config += 'parameters:\n' - config += ' numr: {:d}\n'.format(numr) - config += ' rmin: {:.7f}\n'.format(rmin) - config += ' rmax: {:.3f}\n'.format(rmax) - config += ' Lmax: {:d}\n'.format(lmax) - config += ' nmax: {:d}\n'.format(nmax) - config += ' scale: {:.3f}\n'.format(scale) - config += ' modelname: {}\n'.format(modelname) - config += ' cachename: {}\n'.format(cachename) - config += '...\n' - return config - -def empirical_density_profile(pos, mass, nbins=500, rmin=0, rmax=600, log_space=False): - """ - Computes the number density radial profile assuming all particles have the same mass. - - Args: - pos (ndarray): array of particle positions in cartesian coordinates with shape (n,3). - mass (ndarray): array of particle masses with shape (n,). - nbins (int, optional): number of bins in the radial profile. Default is 500. - rmin (float, optional): minimum radius of the radial profile. Default is 0. - rmax (float, optional): maximum radius of the radial profile. Default is 600. - log_space (bool, optional): whether to use logarithmic binning. Default is False. - - Returns: - tuple: a tuple containing the arrays of radius and density with shapes (nbins,) and (nbins,), respectively. - - Raises: - ValueError: if pos and mass arrays have different lengths or if nbins is not a positive integer. - """ - if len(pos) != len(mass): - raise ValueError("pos and mass arrays must have the same length") - if not isinstance(nbins, int) or nbins <= 0: - raise ValueError("nbins must be a positive integer") - - # Compute radial distances - r_p = np.sqrt(np.sum(pos**2, axis=1)) - - # Compute bin edges and shell volumes - if log_space: - bins = np.logspace(np.log10(rmin), np.log10(rmax), nbins+1) - else: - bins = np.linspace(rmin, rmax, nbins+1) - V_shells = 4/3 * np.pi * (bins[1:]**3 - bins[:-1]**3) - - # Compute density profile - density, _ = np.histogram(r_p, bins=bins, weights=mass) - density /= V_shells - - # Compute bin centers and return profile - radius = 0.5 * (bins[1:] + bins[:-1]) - return radius, density - -def make_exp_basis_table(radius, density, outfile='', plabel='', - verbose=True, physical_units=True, return_values=False): - - """ - Create a table of basis functions for an exponential density profile, compatible with EXP. - - Parameters: - ----------- - radius : array-like - Array of sampled radius points. - density : array-like - Array of density values at radius points. - outfile : str, optional - The name of the output file. If not provided, the file will not be written. - plabel : str, optional - A comment string to add to the output file. - verbose : bool, optional - Whether to print scaling factors and other details during execution. - physical_units : bool, optional - Whether to use physical units (True) or normalized units (False) in the output file. - return_values : bool, optional - Whether to return the radius, density, mass, and potential arrays (True) or not (False). - - Returns: - -------- - If return_values is True: - radius : array-like - The radius values. - density : array-like - The density values. - mass : array-like - The mass enclosed at each radius. - potential : array-like - The potential energy at each radius. - If return_values is False (default): - None - - Notes: - ------ - This function assumes an exponential density profile: - - rho(r) = rho_0 * exp(-r/a) - - where rho_0 and a are constants. - - The table of basis functions is used by EXP to create a density profile. - - Reference: - ---------- - https://gist.github.com/michael-petersen/ec4f20641eedac8f63ec409c9cc65ed7 - """ - - - # make the mass and potential arrays - rvals, dvals = radius, density - mvals = np.zeros(density.size) - pvals = np.zeros(density.size) - pwvals = np.zeros(density.size) - - M = 1. - R = np.nanmax(rvals) - - - # initialise the mass enclosed an potential energy - mvals[0] = 1.e-15 - pwvals[0] = 0. - - # evaluate mass enclosed and potential energy by recursion - for indx in range(1, dvals.size): - mvals[indx] = mvals[indx-1] +\ - 2.0*np.pi*(rvals[indx-1]*rvals[indx-1]*dvals[indx-1] +\ - rvals[indx]*rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); - pwvals[indx] = pwvals[indx-1] + \ - 2.0*np.pi*(rvals[indx-1]*dvals[indx-1] + rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); - - # evaluate potential (see theory document) - pvals = -mvals/(rvals+1.e-10) - (pwvals[dvals.size-1] - pwvals) - - # get the maximum mass and maximum radius - M0 = mvals[dvals.size-1] - R0 = rvals[dvals.size-1] - - # compute scaling factors - Beta = (M/M0) * (R0/R) - Gamma = np.sqrt( (M0*R0)/(M*R) ) * (R0/R) - if verbose: - print("! Scaling: R=",R," M=",M) - - rfac = np.power(Beta,-0.25) * np.power(Gamma,-0.5); - dfac = np.power(Beta,1.5) * Gamma; - mfac = np.power(Beta,0.75) * np.power(Gamma,-0.5); - pfac = Beta; - - if verbose: - print(rfac,dfac,mfac,pfac) - - # save file if desired - if outfile != '': - f = open(outfile,'w') - print('! ', plabel, file=f) - print('! R D M P',file=f) - - print(rvals.size,file=f) - - if physical_units: - for indx in range(0,rvals.size): - print('{0} {1} {2} {3}'.format(rvals[indx],\ - dvals[indx],\ - mvals[indx],\ - pvals[indx]),file=f) - else: - for indx in range(0,rvals.size): - print('{0} {1} {2} {3}'.format(rfac*rvals[indx],\ - dfac*dvals[indx],\ - mfac*mvals[indx],\ - pfac*pvals[indx]),file=f) - - f.close() - - if return_values: - if physical_units: - return rvals, dvals, mvals, pvals - - return rvals*rfac, dfac*dvals, mfac*mvals, pfac*pvals - - -def make_a_BFE(pos, mass, config=None, basis_id='sphereSL', time=0, - numr=500, rmin=0.61, rmax=599, lmax=4, nmax=20, scale=22.5, - modelname='dens_table.txt', cachename='.slgrid_sph_cache', add_coef = False, coef_file='', empirical=False): - """ - Create a BFE expansion for a given set of particle positions and masses. - - Parameters: - pos (numpy.ndarray): The positions of particles. Each row represents one particle, - and each column represents the coordinate of that particle. - mass (numpy.ndarray): The masses of particles. The length of this array should be the same - as the number of particles. - config (pyEXP.config.Config, optional): A configuration object that specifies the basis set. - If not provided, an empirical density profile will be computed - and a configuration object will be created automatically. - basis_id (str, optional): The type of basis set to be used. Default is 'sphereSL'. - time (float, optional): The time at which the expansion is being computed. Default is 0. - numr (int, optional): The number of radial grid points in the basis set. Default is 200. - rmin (float, optional): The minimum radius of the basis set. Default is 0.61. - rmax (float, optional): The maximum radius of the basis set. Default is 599. - lmax (int, optional): The maximum harmonic order in the basis set. Default is 4. - nmax (int, optional): The maximum number of polynomials in the basis set. Default is 20. - scale (float, optional): The scale of the basis set in physical units. - modelname (str, optional): The name of the file containing the density profile model. - Default is 'dens_table.txt'. - cachename (str, optional): The name of the file that will be used to cache the basis set. - Default is '.slgrid_sph_cache'. - save_dir (str, optional): The name of the file if provided that will be used to save the coef files as .h5. - Default is ''. - Returns: - tuple: A tuple containing the basis and the coefficients of the expansion. - The basis is an instance of pyEXP.basis.Basis, and the coefficients are - an instance of pyEXP.coefs.Coefs. - """ - - if os.path.isfile(modelname) == False: - print("-> File model not found so we are computing one \n") - if empirical == True: - print('-> Computing empirical model') - rad, rho = empirical_density_profile(pos, mass, nbins=numr) - R, D, M, P = makemodel_empirical(r_exact, rho, outfile=modelname, return_values=True) - elif empirical == False: - print('-> Computing analytical Hernquist model') - R, D, M, P = makemodel.makemodel(hernquist_halo, 1, [scale], rvals=np.logspace(np.log10(rmin), np.log10(rmax), numr), pfile=modelname) - print('-> Model computed: rmin={}, rmax={}, numr={}'.format(R[0], R[-1], len(R))) - else: - R, D, M, P = np.loadtxt(modelname, skiprows=3, unpack=True) - # check if config file is passed - if config is None: - print('No config file provided.') - print(f'Computing empirical density') - #rad, rho = empirical_density_profile(pos, mass, nbins=500) - #makemodel_empirical(r_exact, rho, outfile=modelname) - #R = [0.01, 600] - config = make_config(basis_id, numr, R[0], R[-1], lmax, nmax, scale, - modelname, cachename) - - # Construct the basis instances - basis = pyEXP.basis.Basis.factory(config) - - # Prints info from Cache - basis.cacheInfo(cachename) - - #compute coefficients - coef = basis.createFromArray(mass/np.sum(mass), pos, time=time) - coefs = pyEXP.coefs.Coefs.makecoefs(coef, 'dark halo') - coefs.add(coef) - if add_coef == False: - coefs.WriteH5Coefs(coef_file) - elif add_coef == True: - coefs.ExtendH5Coefs(coef_file) - - return basis, coefs - From 3a7658733a060460fc94be358b2d2bf4518a1a43 Mon Sep 17 00:00:00 2001 From: jngaravitoc Date: Sat, 11 Nov 2023 15:47:58 -0500 Subject: [PATCH 05/19] rename --- EXP_tools/__init__.py | 3 - EXP_tools/basis_builder/__init__.py | 3 - EXP_tools/basis_builder/basis_utils.py | 285 ----------------------- EXP_tools/basis_builder/makemodel.py | 99 -------- EXP_tools/basis_builder/profiles.py | 127 ----------- EXP_tools/visuals/visualize.py | 304 ------------------------- 6 files changed, 821 deletions(-) delete mode 100755 EXP_tools/__init__.py delete mode 100644 EXP_tools/basis_builder/__init__.py delete mode 100644 EXP_tools/basis_builder/basis_utils.py delete mode 100644 EXP_tools/basis_builder/makemodel.py delete mode 100644 EXP_tools/basis_builder/profiles.py delete mode 100644 EXP_tools/visuals/visualize.py diff --git a/EXP_tools/__init__.py b/EXP_tools/__init__.py deleted file mode 100755 index 11c9ca6..0000000 --- a/EXP_tools/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -from . import visualize -from . import makemodel -from . import basis_utils diff --git a/EXP_tools/basis_builder/__init__.py b/EXP_tools/basis_builder/__init__.py deleted file mode 100644 index 47cd4c4..0000000 --- a/EXP_tools/basis_builder/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -from .makemodel import makemodel -from .profiles import Profiles -from .basis_utils import makebasis diff --git a/EXP_tools/basis_builder/basis_utils.py b/EXP_tools/basis_builder/basis_utils.py deleted file mode 100644 index 49f3395..0000000 --- a/EXP_tools/basis_builder/basis_utils.py +++ /dev/null @@ -1,285 +0,0 @@ -import os, sys, pickle, pyEXP -import numpy as np -import makemodel - -def make_config(basis_id, numr, rmin, rmax, lmax, nmax, scale, - modelname='', cachename='.slgrid_sph_cache'): - """ - Creates a configuration file required to build a basis model. - - Args: - basis_id (str): The identity of the basis model. - numr (int): The number of radial grid points. - rmin (float): The minimum radius value. - rmax (float): The maximum radius value. - lmax (int): The maximum l value of the basis. - nmax (int): The maximum n value of the basis. - scale (float): Scaling factor for the basis. - modelname (str, optional): Name of the model. Default is an empty string. - cachename (str, optional): Name of the cache file. Default is '.slgrid_sph_cache'. - - Returns: - str: A string representation of the configuration file. - - Raises: - None - """ - - config = '\n---\nid: {:s}\n'.format(basis_id) - config += 'parameters:\n' - config += ' numr: {:d}\n'.format(numr) - config += ' rmin: {:.7f}\n'.format(rmin) - config += ' rmax: {:.3f}\n'.format(rmax) - config += ' Lmax: {:d}\n'.format(lmax) - config += ' nmax: {:d}\n'.format(nmax) - config += ' scale: {:.3f}\n'.format(scale) - config += ' modelname: {}\n'.format(modelname) - config += ' cachename: {}\n'.format(cachename) - config += '...\n' - return config - -def empirical_density_profile(pos, mass, nbins=500, rmin=0, rmax=600, log_space=False): - """ - Computes the number density radial profile assuming all particles have the same mass. - - Args: - pos (ndarray): array of particle positions in cartesian coordinates with shape (n,3). - mass (ndarray): array of particle masses with shape (n,). - nbins (int, optional): number of bins in the radial profile. Default is 500. - rmin (float, optional): minimum radius of the radial profile. Default is 0. - rmax (float, optional): maximum radius of the radial profile. Default is 600. - log_space (bool, optional): whether to use logarithmic binning. Default is False. - - Returns: - tuple: a tuple containing the arrays of radius and density with shapes (nbins,) and (nbins,), respectively. - - Raises: - ValueError: if pos and mass arrays have different lengths or if nbins is not a positive integer. - """ - if len(pos) != len(mass): - raise ValueError("pos and mass arrays must have the same length") - if not isinstance(nbins, int) or nbins <= 0: - raise ValueError("nbins must be a positive integer") - - # Compute radial distances - r_p = np.sqrt(np.sum(pos**2, axis=1)) - - # Compute bin edges and shell volumes - if log_space: - bins = np.logspace(np.log10(rmin), np.log10(rmax), nbins+1) - else: - bins = np.linspace(rmin, rmax, nbins+1) - V_shells = 4/3 * np.pi * (bins[1:]**3 - bins[:-1]**3) - - # Compute density profile - density, _ = np.histogram(r_p, bins=bins, weights=mass) - density /= V_shells - - # Compute bin centers and return profile - radius = 0.5 * (bins[1:] + bins[:-1]) - return radius, density - -def make_exp_basis_table(radius, density, outfile='', plabel='', - verbose=True, physical_units=True, return_values=False): - - """ - Create a table of basis functions for an exponential density profile, compatible with EXP. - - Parameters: - ----------- - radius : array-like - Array of sampled radius points. - density : array-like - Array of density values at radius points. - outfile : str, optional - The name of the output file. If not provided, the file will not be written. - plabel : str, optional - A comment string to add to the output file. - verbose : bool, optional - Whether to print scaling factors and other details during execution. - physical_units : bool, optional - Whether to use physical units (True) or normalized units (False) in the output file. - return_values : bool, optional - Whether to return the radius, density, mass, and potential arrays (True) or not (False). - - Returns: - -------- - If return_values is True: - radius : array-like - The radius values. - density : array-like - The density values. - mass : array-like - The mass enclosed at each radius. - potential : array-like - The potential energy at each radius. - If return_values is False (default): - None - - Notes: - ------ - This function assumes an exponential density profile: - - rho(r) = rho_0 * exp(-r/a) - - where rho_0 and a are constants. - - The table of basis functions is used by EXP to create a density profile. - - Reference: - ---------- - https://gist.github.com/michael-petersen/ec4f20641eedac8f63ec409c9cc65ed7 - """ - - - # make the mass and potential arrays - rvals, dvals = radius, density - mvals = np.zeros(density.size) - pvals = np.zeros(density.size) - pwvals = np.zeros(density.size) - - M = 1. - R = np.nanmax(rvals) - - - # initialise the mass enclosed an potential energy - mvals[0] = 1.e-15 - pwvals[0] = 0. - - # evaluate mass enclosed and potential energy by recursion - for indx in range(1, dvals.size): - mvals[indx] = mvals[indx-1] +\ - 2.0*np.pi*(rvals[indx-1]*rvals[indx-1]*dvals[indx-1] +\ - rvals[indx]*rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); - pwvals[indx] = pwvals[indx-1] + \ - 2.0*np.pi*(rvals[indx-1]*dvals[indx-1] + rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); - - # evaluate potential (see theory document) - pvals = -mvals/(rvals+1.e-10) - (pwvals[dvals.size-1] - pwvals) - - # get the maximum mass and maximum radius - M0 = mvals[dvals.size-1] - R0 = rvals[dvals.size-1] - - # compute scaling factors - Beta = (M/M0) * (R0/R) - Gamma = np.sqrt( (M0*R0)/(M*R) ) * (R0/R) - if verbose: - print("! Scaling: R=",R," M=",M) - - rfac = np.power(Beta,-0.25) * np.power(Gamma,-0.5); - dfac = np.power(Beta,1.5) * Gamma; - mfac = np.power(Beta,0.75) * np.power(Gamma,-0.5); - pfac = Beta; - - if verbose: - print(rfac,dfac,mfac,pfac) - - # save file if desired - if outfile != '': - f = open(outfile,'w') - print('! ', plabel, file=f) - print('! R D M P',file=f) - - print(rvals.size,file=f) - - if physical_units: - for indx in range(0,rvals.size): - print('{0} {1} {2} {3}'.format(rvals[indx],\ - dvals[indx],\ - mvals[indx],\ - pvals[indx]),file=f) - else: - for indx in range(0,rvals.size): - print('{0} {1} {2} {3}'.format(rfac*rvals[indx],\ - dfac*dvals[indx],\ - mfac*mvals[indx],\ - pfac*pvals[indx]),file=f) - - f.close() - - if return_values: - if physical_units: - return rvals, dvals, mvals, pvals - - return rvals*rfac, dfac*dvals, mfac*mvals, pfac*pvals - - -def makebasis(pos, mass, basis_model, config=None, basis_id='sphereSL', time=0, - numr=500, rmin=0.61, rmax=599, lmax=4, nmax=20, scale=22.5, - modelname='dens_table.txt', cachename='.slgrid_sph_cache', add_coef = False, coef_file=''): - """ - Create a BFE expansion for a given set of particle positions and masses. - - Parameters: - pos (numpy.ndarray): The positions of particles. Each row represents one particle, - and each column represents the coordinate of that particle. - mass (numpy.ndarray): The masses of particles. The length of this array should be the same - as the number of particles. - basismodel (): - config (pyEXP.config.Config, optional): A configuration object that specifies the basis set. - If not provided, an empirical density profile will be computed - and a configuration object will be created automatically. - basis_id (str, optional): The type of basis set to be used. Default is 'sphereSL'. - time (float, optional): The time at which the expansion is being computed. Default is 0. - numr (int, optional): The number of radial grid points in the basis set. Default is 200. - rmin (float, optional): The minimum radius of the basis set. Default is 0.61. - rmax (float, optional): The maximum radius of the basis set. Default is 599. - lmax (int, optional): The maximum harmonic order in the basis set. Default is 4. - nmax (int, optional): The maximum number of polynomials in the basis set. Default is 20. - scale (float, optional): The scale of the basis set in physical units. - modelname (str, optional): The name of the file containing the density profile model. - Default is 'dens_table.txt'. - cachename (str, optional): The name of the file that will be used to cache the basis set. - Default is '.slgrid_sph_cache'. - save_dir (str, optional): The name of the file if provided that will be used to save the coef files as .h5. - Default is ''. - Returns: - tuple: A tuple containing the basis and the coefficients of the expansion. - The basis is an instance of pyEXP.basis.Basis, and the coefficients are - an instance of pyEXP.coefs.Coefs. - """ - - if os.path.isfile(modelname) == False: - print("-> File model not found so we are computing one \n") - if empirical == True: - print('-> Computing empirical model') - rad, rho = empirical_density_profile(pos, mass, nbins=numr) - elif empirical == False: - makemodel.hernquist_halo() - - makemodel.Profiles(density_profile) - - R, D, M, P = makemodel.makemodel(rvals=np.logspace(np.log10(rmin), np.log10(rmax), numr), rho, outfile=modelname, return_values=True) - print('-> Computing analytical Hernquist model') - R, D, M, P = makemodel.makemodel(hernquist_halo, 1, [scale], rvals=, pfile=modelname) - print('-> Model computed: rmin={}, rmax={}, numr={}'.format(R[0], R[-1], len(R))) - else: - R, D, M, P = np.loadtxt(modelname, skiprows=3, unpack=True) - # check if config file is passed - if config is None: - print('No config file provided.') - print(f'Computing empirical density') - #rad, rho = empirical_density_profile(pos, mass, nbins=500) - #makemodel_empirical(r_exact, rho, outfile=modelname) - #R = [0.01, 600] - config = make_config(basis_id, numr, R[0], R[-1], lmax, nmax, scale, - modelname, cachename) - - # Construct the basis instances - basis = pyEXP.basis.Basis.factory(config) - - # Prints info from Cache - basis.cacheInfo(cachename) - - #compute coefficients - coef = basis.createFromArray(mass/np.sum(mass), pos, time=time) - coefs = pyEXP.coefs.Coefs.makecoefs(coef, 'dark halo') - coefs.add(coef) - if add_coef == False: - coefs.WriteH5Coefs(coef_file) - elif add_coef == True: - coefs.ExtendH5Coefs(coef_file) - - return basis, coefs - diff --git a/EXP_tools/basis_builder/makemodel.py b/EXP_tools/basis_builder/makemodel.py deleted file mode 100644 index 5c1aa31..0000000 --- a/EXP_tools/basis_builder/makemodel.py +++ /dev/null @@ -1,99 +0,0 @@ -""" -TODO: -""" - -import numpy as np - - -def makemodel(dvals, rvals, M, unit_physical=False, - pfile='', plabel = '', verbose=True): - """ - Make an EXP-compatible spherical basis function table - - inputs - ------------- - basis_type : (str) - func : (function) the callable functional form of the density - M : (float) the total mass of the model, sets normalisations - funcargs : (list) a list of arguments for the density function. - rvals : (array of floats) radius values to evaluate the density function - = 10.**np.linspace(-2.,4.,2000) - pfile : (string) the name of the output file. If '', will not print file - plabel : (string) comment string - verbose : (boolean) - - outputs - ------------- - R : (array of floats) the radius values - D : (array of floats) the density - M : (array of floats) the mass enclosed - P : (array of floats) the potential - - """ - - R = np.nanmax(rvals) - - # query out the density values - - # make the mass and potential arrays - mvals = np.zeros(dvals.size) - pvals = np.zeros(dvals.size) - pwvals = np.zeros(dvals.size) - - # initialise the mass enclosed an potential energy - mvals[0] = 1.e-15 - pwvals[0] = 0. - - # evaluate mass enclosed and potential energy by recursion - for indx in range(1,dvals.size): - mvals[indx] = mvals[indx-1] +\ - 2.0*np.pi*(rvals[indx-1]*rvals[indx-1]*dvals[indx-1] +\ - rvals[indx]*rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); - pwvals[indx] = pwvals[indx-1] + \ - 2.0*np.pi*(rvals[indx-1]*dvals[indx-1] + rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); - - # evaluate potential (see theory document) - pvals = -mvals/(rvals+1.e-10) - (pwvals[dvals.size-1] - pwvals) - - # get the maximum mass and maximum radius - M0 = mvals[dvals.size-1] - R0 = rvals[dvals.size-1] - - # compute scaling factors - Beta = (M/M0) * (R0/R); - Gamma = np.sqrt((M0*R0)/(M*R)) * (R0/R); - if verbose: - print("! Scaling: R=",R," M=",M) - - rfac = np.power(Beta,-0.25) * np.power(Gamma,-0.5); - dfac = np.power(Beta,1.5) * Gamma; - mfac = np.power(Beta,0.75) * np.power(Gamma,-0.5); - pfac = Beta; - - if verbose: - print(rfac,dfac,mfac,pfac) - - # save file if desired - if pfile != '': - f = open(pfile,'w') - print('! ',plabel,file=f) - print('! R D M P',file=f) - - print(rvals.size,file=f) - - if unit_physical = True: - for indx in range(0,rvals.size): - print('{0} {1} {2} {3}'.format(rvals[indx],\ - dvals[indx],\ - mvals[indx],\ - pvals[indx]),file=f) - else: - - f.close() - - if return_val: - if unit_phy: - return rvals, dvals, mvals, pvals - else: - return rvals*rfac, dfac*dvals, mfac*mvals, pfac*pvals - diff --git a/EXP_tools/basis_builder/profiles.py b/EXP_tools/basis_builder/profiles.py deleted file mode 100644 index fe4f7e3..0000000 --- a/EXP_tools/basis_builder/profiles.py +++ /dev/null @@ -1,127 +0,0 @@ -import numpy as np -from scipy import special - -Class Profiles: - def __init__(self, radius, scale_radius, alpha=1., beta=2.0): - """ - Class with density profiles of Dark Matter halos - - inputs - ---------- - r : (float) radius values - rs : (float, default=1.) scale radius - alpha : (float, default=1.) inner halo slope - beta : (float, default=1.e-7) outer halo slope - - notes - ---------- - different combinations are known distributions. - alpha=1,beta=2 is NFW (default) - alpha=1,beta=3 is Hernquist - alpha=2.5, beta=0 is a typical single power law halo - - """ - - self.radius = radius - self.scale_radius = scale_radius - self.alpha = alpha - self.beta = beta - self.ra = self.radius/self.scale_radius - - - def powerhalo(self, rc=0.0): - """ - Generic two power law distribution - - inputs - ---------- - rc : (float, default=0. i.e. no core) core radius - - returns - ---------- - densities evaluated at r - - notes - ---------- - different combinations are known distributions. - alpha=1,beta=2 is NFW - alpha=1,beta=3 is Hernquist - alpha=2.5,beta=0 is a typical single power law halo - - - - """ - return 1./(((self.ra+rc/self.scale_radius)**self.alpha)*((1+self.ra)**self.beta)) - - - def powerhalorolloff(self, rc=0.0): - """ - Generic twopower law distribution with an erf rolloff - - - inputs - ---------- - rc : (float, default=0. i.e. no core) core radius - - returns - ---------- - densities evaluated at r - - notes - ---------- - different combinations are known distributions. - alpha=1,beta=2 is NFW - alpha=1,beta=3 is Hernquist - alpha=2.5,beta=0 is a typical single power law halo - - - """ - dens = 1./(((self.ra+rc/self.scale_radius)**self.alpha)*((1+self.ra)**self.beta)) - rtrunc = 25*self.scale_radius - wtrunc = rtrunc*0.2 - rolloff = 0.5 - 0.5*special.erf((self.radius-rtrunc)/wtrunc) - return dens*rolloff - - - def plummer_density(self): - """ - Plummer density profile - - inputs - --------- - - - returns - ---------- - densities evaluated at r - """ - return ((3.0)/(4*np.pi))*(self.scale_radius**2.)*((self.scale_radius**2 + self.radius**2)**(-2.5)) - - def twopower_density_withrolloff(self, rcen, wcen): - """ - inputs - --------- - - - returns - ---------- - densities evaluated at r - """ - - - prefac = 0.5*(1.-scipy.special.erf((ra-rcen)/wcen)) - return prefac*(ra**-self.alpha)*(1+ra)**(-self.beta+self.alpha) - - def hernquist_halo(self): - """ - TODO: DO we need these profile - Hernquist halo - - inputs - --------- - - returns - ---------- - densities evaluated at r - """ - return 1 / ( 2*np.pi * (self.radius/self.scale_radius) * (1 + self.radius/self.scale_radius)**3) \ No newline at end of file diff --git a/EXP_tools/visuals/visualize.py b/EXP_tools/visuals/visualize.py deleted file mode 100644 index 8db665d..0000000 --- a/EXP_tools/visuals/visualize.py +++ /dev/null @@ -1,304 +0,0 @@ -import os, sys, pickle, pyEXP -import numpy as np - -###Field computations for plotting### -def make_basis_plot(basis, savefile=None, nsnap='mean', y=0.92, dpi=200): - """ - Plots the potential of the basis functions for different values of l and n. - - Args: - basis (obj): object containing the basis functions for the simulation - savefile (str, optional): name of the file to save the plot as - nsnap (str, optional): description of the snapshot being plotted - y (float, optional): vertical position of the main title - dpi (int, optional): resolution of the plot in dots per inch - - Returns: - None - - """ - # Set up grid for plotting potential - lrmin, lrmax, rnum = 0.5, 2.7, 100 - halo_grid = basis.getBasis(lrmin, lrmax, rnum) - r = np.linspace(lrmin, lrmax, rnum) - r = np.power(10.0, r) - - # Create subplots and plot potential for each l and n - import matplotlib.pyplot as plt - fig, ax = plt.subplots(4, 5, figsize=(6,6), dpi=dpi, - sharex='col', sharey='row') - plt.subplots_adjust(wspace=0, hspace=0) - ax = ax.flatten() - - for l in range(len(ax)): - ax[l].set_title(f"$\ell = {l}$", y=0.8, fontsize=6) - for n in range(20): - ax[l].semilogx(r, halo_grid[l][n]['potential'], '-', label="n={}".format(n), lw=0.5) - - # Add labels and main title - fig.supylabel('Potential', weight='bold', x=-0.02) - fig.supxlabel('Radius', weight='bold', y=0.02) - fig.suptitle(f'nsnap = {nsnap}', - fontsize=12, - weight='bold', - y=y, - ) - - # Save plot if a filename was provided - if savefile: - plt.savefig(f'{savefile}', bbox_inches='tight') - -def find_field(basis, coefficients, time=0, xyz=(0, 0, 0), property='dens', include_monopole=True): - """ - Finds the value of the specified property of the field at the given position. - - Args: - basis (obj): Object containing the basis functions for the simulation. - coefficients (obj): Object containing the coefficients for the simulation. - time (float, optional): The time at which to evaluate the field. Default is 0. - xyz (tuple or list, optional): The (x, y, z) position at which to evaluate the field. Default is (0, 0, 0). - property (str, optional): The property of the field to evaluate. Can be 'dens', 'pot', or 'force'. Default is 'dens'. - include_monopole (bool, optional): Whether to return the monopole contribution to the property only. Default is True. - - Returns: - float or list: The value of the specified property of the field at the given position. If property is 'force', a list of - three values is returned representing the force vector in (x, y, z) directions. - - Raises: - ValueError: If the property argument is not 'dens', 'pot', or 'force'. - """ - - coefficients.set_coefs(coefficients.getCoefStruct(time)) - dens0, pot0, dens, pot, fx, fy, fz = basis.getFields(xyz[0], xyz[1], xyz[2]) - - if property == 'dens': - if include_monopole: - return dens0 - return dens + dens0 - - elif property == 'pot': - if include_monopole: - return pot0 - return pot + pot0 - - elif property == 'force': - return [fx, fy, fz] - - else: - raise ValueError("Invalid property specified. Possible values are 'dens', 'pot', and 'force'.") - -def spherical_avg_prop(basis, coefficients, time=0, radius=np.linspace(0.1, 600, 100), property='dens'): - """ - Computes the spherically averaged value of the specified property of the field over the given radii. - - Args: - basis (obj): Object containing the basis functions for the simulation. - coefficients (obj): Object containing the coefficients for the simulation. - time (float, optional): The time at which to evaluate the field. Default is 0. - radius (ndarray, optional): An array of radii over which to compute the spherically averaged property. Default is an - array of 100 values logarithmically spaced between 0.1 and 600. - property (str, optional): The property of the field to evaluate. Can be 'dens', 'pot', or 'force'. Default is 'dens'. - - Returns: - ndarray: An array of spherically averaged values of the specified property over the given radii. - - Raises: - ValueError: If the property argument is not 'dens', 'pot', or 'force'. - """ - - coefficients.set_coefs(coefficients.getCoefStruct(time)) - field = [find_field(basis, np.hstack([[rad], [0], [0]]), property=property, include_monopole=True) for rad in radius] - - if property == 'force': - return np.vstack(field), radius - - return np.array(field), radius - -def return_fields_in_grid(basis, coefficients, times=[0], - projection='3D', proj_plane=0, - grid_lim=300, npoints=150,): - - """ - Returns the 2D or 3D field grids as a dict at each time step as a key. - - Args: - basis (obj): object containing the basis functions for the simulation - coefficients (obj): object containing the coefficients for the simulation - times (list): list of the time to compute the field grid. - projection (str): the slice projection to plot. Can be '3D', XY', 'XZ', or 'YZ'. - proj_plane (float, optional): the value of the coordinate that is held constant in the 2D slice projection - npoints (int, optional): the number of grid points in each dimension - grid_limits (float, optional): the limits of the grid in the dimensions. - - Returns: fields along with the specified 2D or 3D grid. The fields are in a - heirarchical dict structure with for each time: - dict_keys(['d', 'd0', 'd1', 'dd', 'fp', 'fr', 'ft', 'p', 'p0', 'p1']) - where - - 'd': Total density at each point. - - 'd0' : Density in the monopole term. - - 'd1' : Density in l>0 terms. - - 'dd': Density difference. - - 'fp': Force in the phi direction. - - 'fr': Force in the radial direction. - - 'ft': Force in the theta direction. - - 'p': Total potential at each point. - - 'p0': Potential in the monopole term. - - 'p1': Potential in l>0 terms. - """ - - assert projection in ['3D', 'XY', 'XZ', 'YZ'], "Invalid value for 'projection'. Must be one of '3D', 'XY', 'XZ', or 'YZ'." - - - - if projection == '3D': - pmin = [-grid_lim, -grid_lim, -grid_lim] - pmax = [grid_lim, grid_lim, grid_lim] - grid = [npoints, npoints, npoints] - - ##create a projection grid along with it - x = np.linspace(-grid_lim, grid_lim, npoints) - xgrid = np.meshgrid(x, x, x) - - else: - x = np.linspace(-grid_lim, grid_lim, npoints) - xgrid = np.meshgrid(x, x) - grid = [npoints, npoints] - - if projection == 'XY': - pmin = [-grid_lim, -grid_lim, proj_plane] - pmax = [grid_lim, grid_lim, proj_plane] - - elif projection == 'XZ': - pmin = [-grid_lim, proj_plane, -grid_lim] - pmax = [grid_lim, proj_plane, grid_lim] - - elif projection == 'YZ': - pmin = [proj_plane, -grid_lim, -grid_lim] - pmax = [proj_plane, grid_lim, grid_lim] - - - field_gen = pyEXP.field.FieldGenerator(times, pmin, pmax, grid) - - return field_gen.volumes(basis, coefs), xgrid - -def slice_fields(basis, coefficients, time=0, - projection='XY', proj_plane=0, npoints=300, - grid_limits=(-300, 300), prop='dens', monopole_only=False): - """ - Plots a slice projection of the fields of a simulation. - - Args: - basis (obj): object containing the basis functions for the simulation - coefficients (obj): object containing the coefficients for the simulation - time (float): the time at which to plot the fields - projection (str): the slice projection to plot. Can be 'XY', 'XZ', or 'YZ'. - proj_plane (float, optional): the value of the coordinate that is held constant in the slice projection - npoints (int, optional): the number of grid points in each dimension - grid_limits (tuple, optional): the limits of the grid in the x and y dimensions, in the form (x_min, x_max) - prop (str, optional): the property to return. Can be 'dens' (density), 'pot' (potential), or 'force' (force). - monopole_only (bool, optional): whether to return the monopole component in the returned property value. - - Returns: - array or list: the property specified by `prop`. If `prop` is 'force', a list of the x, y, and z components of the force is returned. - Also returns the grid used to compute slice fields. - """ - x = np.linspace(grid_limits[0], grid_limits[1], npoints) - xgrid = np.meshgrid(x, x) - xg = xgrid[0].flatten() - yg = xgrid[1].flatten() - - - if projection not in ['XY', 'XZ', 'YZ']: - raise ValueError("Invalid projection specified. Possible values are 'XY', 'XZ', and 'YZ'.") - - N = len(xg) - rho0 = np.zeros_like(xg) - pot0 = np.zeros_like(xg) - rho = np.zeros_like(xg) - pot = np.zeros_like(xg) - basis.set_coefs(coefficients.getCoefStruct(time)) - - for k in range(0, N): - if projection == 'XY': - rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(xg[k], yg[k], proj_plane) - elif projection == 'XZ': - rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(xg[k], proj_plane, yg[k]) - elif projection == 'YZ': - rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(proj_plane, xg[k], yg[k]) - - dens = rho.reshape(npoints, npoints) - pot = pot.reshape(npoints, npoints) - dens0 = rho0.reshape(npoints, npoints) - pot0 = pot0.reshape(npoints, npoints) - - if prop == 'dens': - if monopole_only: - return dens0 - return dens0, dens, xgrid - - if prop == 'pot': - if monopole_only: - return pot0 - return pot0, pot, xgrid - - if prop == 'force': - return [fx.reshape(npoints, npoints), fy.reshape(npoints, npoints), fz.reshape(npoints, npoints)], xgrid - - -def slice_3d_fields(basis, coefficients, time=0, - projection='XY', proj_plane=0, npoints=300, - grid_limits=(-300, 300), prop='dens', monopole_only=False): - """ - Plots a slice projection of the fields of a simulation. - - Args: - basis (obj): object containing the basis functions for the simulation - coefficients (obj): object containing the coefficients for the simulation - time (float): the time at which to plot the fields - projection (str): the slice projection to plot. Can be 'XY', 'XZ', or 'YZ'. - proj_plane (float, optional): the value of the coordinate that is held constant in the slice projection - npoints (int, optional): the number of grid points in each dimension - grid_limits (tuple, optional): the limits of the grid in the x and y dimensions, in the form (x_min, x_max) - prop (str, optional): the property to return. Can be 'dens' (density), 'pot' (potential), or 'force' (force). - monopole_only (bool, optional): whether to return the monopole component in the returned property value. - - Returns: - array or list: the property specified by `prop`. If `prop` is 'force', a list of the x, y, and z components of the force is returned. - Also returns the grid used to compute slice fields. - """ - x = np.linspace(grid_limits[0], grid_limits[1], npoints) - xgrid = np.meshgrid(x, x, x) - xg = xgrid[0].flatten() - yg = xgrid[1].flatten() - zg = xgrid[2].flatten() - - - - N = len(xg) - rho0 = np.zeros_like(xg) - pot0 = np.zeros_like(xg) - rho = np.zeros_like(xg) - pot = np.zeros_like(xg) - basis.set_coefs(coefficients.getCoefStruct(time)) - - for k in range(0, N): - rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(xg[k], yg[k], zg[k]) - - dens = rho.reshape(npoints, npoints, npoints) - pot = pot.reshape(npoints, npoints, npoints) - dens0 = rho0.reshape(npoints, npoints, npoints) - pot0 = pot0.reshape(npoints, npoints, npoints) - - if prop == 'dens': - if monopole_only: - return dens0 - return dens0, dens, xgrid - - if prop == 'pot': - if monopole_only: - return pot0 - return pot0, pot, xgrid - - if prop == 'force': - return [fx.reshape(npoints, npoints, npoints), fy.reshape(npoints, npoints, npoints), fz.reshape(npoints, npoints, npoints)], xgrid - From 5a0525f8d1ab298cae45ce192a6bd0a2e559a1d2 Mon Sep 17 00:00:00 2001 From: jngaravitoc Date: Sat, 11 Nov 2023 15:48:09 -0500 Subject: [PATCH 06/19] EXPtools --- EXPtools/__init__.py | 2 + EXPtools/basis_builder/__init__.py | 3 + EXPtools/basis_builder/basis_utils.py | 285 ++++++++++++++++++++++++ EXPtools/basis_builder/makemodel.py | 99 +++++++++ EXPtools/basis_builder/profiles.py | 127 +++++++++++ EXPtools/visuals/visualize.py | 304 ++++++++++++++++++++++++++ 6 files changed, 820 insertions(+) create mode 100755 EXPtools/__init__.py create mode 100644 EXPtools/basis_builder/__init__.py create mode 100644 EXPtools/basis_builder/basis_utils.py create mode 100644 EXPtools/basis_builder/makemodel.py create mode 100644 EXPtools/basis_builder/profiles.py create mode 100644 EXPtools/visuals/visualize.py diff --git a/EXPtools/__init__.py b/EXPtools/__init__.py new file mode 100755 index 0000000..372e92d --- /dev/null +++ b/EXPtools/__init__.py @@ -0,0 +1,2 @@ +from . import visuals +from . import basis_builder diff --git a/EXPtools/basis_builder/__init__.py b/EXPtools/basis_builder/__init__.py new file mode 100644 index 0000000..47cd4c4 --- /dev/null +++ b/EXPtools/basis_builder/__init__.py @@ -0,0 +1,3 @@ +from .makemodel import makemodel +from .profiles import Profiles +from .basis_utils import makebasis diff --git a/EXPtools/basis_builder/basis_utils.py b/EXPtools/basis_builder/basis_utils.py new file mode 100644 index 0000000..49f3395 --- /dev/null +++ b/EXPtools/basis_builder/basis_utils.py @@ -0,0 +1,285 @@ +import os, sys, pickle, pyEXP +import numpy as np +import makemodel + +def make_config(basis_id, numr, rmin, rmax, lmax, nmax, scale, + modelname='', cachename='.slgrid_sph_cache'): + """ + Creates a configuration file required to build a basis model. + + Args: + basis_id (str): The identity of the basis model. + numr (int): The number of radial grid points. + rmin (float): The minimum radius value. + rmax (float): The maximum radius value. + lmax (int): The maximum l value of the basis. + nmax (int): The maximum n value of the basis. + scale (float): Scaling factor for the basis. + modelname (str, optional): Name of the model. Default is an empty string. + cachename (str, optional): Name of the cache file. Default is '.slgrid_sph_cache'. + + Returns: + str: A string representation of the configuration file. + + Raises: + None + """ + + config = '\n---\nid: {:s}\n'.format(basis_id) + config += 'parameters:\n' + config += ' numr: {:d}\n'.format(numr) + config += ' rmin: {:.7f}\n'.format(rmin) + config += ' rmax: {:.3f}\n'.format(rmax) + config += ' Lmax: {:d}\n'.format(lmax) + config += ' nmax: {:d}\n'.format(nmax) + config += ' scale: {:.3f}\n'.format(scale) + config += ' modelname: {}\n'.format(modelname) + config += ' cachename: {}\n'.format(cachename) + config += '...\n' + return config + +def empirical_density_profile(pos, mass, nbins=500, rmin=0, rmax=600, log_space=False): + """ + Computes the number density radial profile assuming all particles have the same mass. + + Args: + pos (ndarray): array of particle positions in cartesian coordinates with shape (n,3). + mass (ndarray): array of particle masses with shape (n,). + nbins (int, optional): number of bins in the radial profile. Default is 500. + rmin (float, optional): minimum radius of the radial profile. Default is 0. + rmax (float, optional): maximum radius of the radial profile. Default is 600. + log_space (bool, optional): whether to use logarithmic binning. Default is False. + + Returns: + tuple: a tuple containing the arrays of radius and density with shapes (nbins,) and (nbins,), respectively. + + Raises: + ValueError: if pos and mass arrays have different lengths or if nbins is not a positive integer. + """ + if len(pos) != len(mass): + raise ValueError("pos and mass arrays must have the same length") + if not isinstance(nbins, int) or nbins <= 0: + raise ValueError("nbins must be a positive integer") + + # Compute radial distances + r_p = np.sqrt(np.sum(pos**2, axis=1)) + + # Compute bin edges and shell volumes + if log_space: + bins = np.logspace(np.log10(rmin), np.log10(rmax), nbins+1) + else: + bins = np.linspace(rmin, rmax, nbins+1) + V_shells = 4/3 * np.pi * (bins[1:]**3 - bins[:-1]**3) + + # Compute density profile + density, _ = np.histogram(r_p, bins=bins, weights=mass) + density /= V_shells + + # Compute bin centers and return profile + radius = 0.5 * (bins[1:] + bins[:-1]) + return radius, density + +def make_exp_basis_table(radius, density, outfile='', plabel='', + verbose=True, physical_units=True, return_values=False): + + """ + Create a table of basis functions for an exponential density profile, compatible with EXP. + + Parameters: + ----------- + radius : array-like + Array of sampled radius points. + density : array-like + Array of density values at radius points. + outfile : str, optional + The name of the output file. If not provided, the file will not be written. + plabel : str, optional + A comment string to add to the output file. + verbose : bool, optional + Whether to print scaling factors and other details during execution. + physical_units : bool, optional + Whether to use physical units (True) or normalized units (False) in the output file. + return_values : bool, optional + Whether to return the radius, density, mass, and potential arrays (True) or not (False). + + Returns: + -------- + If return_values is True: + radius : array-like + The radius values. + density : array-like + The density values. + mass : array-like + The mass enclosed at each radius. + potential : array-like + The potential energy at each radius. + If return_values is False (default): + None + + Notes: + ------ + This function assumes an exponential density profile: + + rho(r) = rho_0 * exp(-r/a) + + where rho_0 and a are constants. + + The table of basis functions is used by EXP to create a density profile. + + Reference: + ---------- + https://gist.github.com/michael-petersen/ec4f20641eedac8f63ec409c9cc65ed7 + """ + + + # make the mass and potential arrays + rvals, dvals = radius, density + mvals = np.zeros(density.size) + pvals = np.zeros(density.size) + pwvals = np.zeros(density.size) + + M = 1. + R = np.nanmax(rvals) + + + # initialise the mass enclosed an potential energy + mvals[0] = 1.e-15 + pwvals[0] = 0. + + # evaluate mass enclosed and potential energy by recursion + for indx in range(1, dvals.size): + mvals[indx] = mvals[indx-1] +\ + 2.0*np.pi*(rvals[indx-1]*rvals[indx-1]*dvals[indx-1] +\ + rvals[indx]*rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); + pwvals[indx] = pwvals[indx-1] + \ + 2.0*np.pi*(rvals[indx-1]*dvals[indx-1] + rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); + + # evaluate potential (see theory document) + pvals = -mvals/(rvals+1.e-10) - (pwvals[dvals.size-1] - pwvals) + + # get the maximum mass and maximum radius + M0 = mvals[dvals.size-1] + R0 = rvals[dvals.size-1] + + # compute scaling factors + Beta = (M/M0) * (R0/R) + Gamma = np.sqrt( (M0*R0)/(M*R) ) * (R0/R) + if verbose: + print("! Scaling: R=",R," M=",M) + + rfac = np.power(Beta,-0.25) * np.power(Gamma,-0.5); + dfac = np.power(Beta,1.5) * Gamma; + mfac = np.power(Beta,0.75) * np.power(Gamma,-0.5); + pfac = Beta; + + if verbose: + print(rfac,dfac,mfac,pfac) + + # save file if desired + if outfile != '': + f = open(outfile,'w') + print('! ', plabel, file=f) + print('! R D M P',file=f) + + print(rvals.size,file=f) + + if physical_units: + for indx in range(0,rvals.size): + print('{0} {1} {2} {3}'.format(rvals[indx],\ + dvals[indx],\ + mvals[indx],\ + pvals[indx]),file=f) + else: + for indx in range(0,rvals.size): + print('{0} {1} {2} {3}'.format(rfac*rvals[indx],\ + dfac*dvals[indx],\ + mfac*mvals[indx],\ + pfac*pvals[indx]),file=f) + + f.close() + + if return_values: + if physical_units: + return rvals, dvals, mvals, pvals + + return rvals*rfac, dfac*dvals, mfac*mvals, pfac*pvals + + +def makebasis(pos, mass, basis_model, config=None, basis_id='sphereSL', time=0, + numr=500, rmin=0.61, rmax=599, lmax=4, nmax=20, scale=22.5, + modelname='dens_table.txt', cachename='.slgrid_sph_cache', add_coef = False, coef_file=''): + """ + Create a BFE expansion for a given set of particle positions and masses. + + Parameters: + pos (numpy.ndarray): The positions of particles. Each row represents one particle, + and each column represents the coordinate of that particle. + mass (numpy.ndarray): The masses of particles. The length of this array should be the same + as the number of particles. + basismodel (): + config (pyEXP.config.Config, optional): A configuration object that specifies the basis set. + If not provided, an empirical density profile will be computed + and a configuration object will be created automatically. + basis_id (str, optional): The type of basis set to be used. Default is 'sphereSL'. + time (float, optional): The time at which the expansion is being computed. Default is 0. + numr (int, optional): The number of radial grid points in the basis set. Default is 200. + rmin (float, optional): The minimum radius of the basis set. Default is 0.61. + rmax (float, optional): The maximum radius of the basis set. Default is 599. + lmax (int, optional): The maximum harmonic order in the basis set. Default is 4. + nmax (int, optional): The maximum number of polynomials in the basis set. Default is 20. + scale (float, optional): The scale of the basis set in physical units. + modelname (str, optional): The name of the file containing the density profile model. + Default is 'dens_table.txt'. + cachename (str, optional): The name of the file that will be used to cache the basis set. + Default is '.slgrid_sph_cache'. + save_dir (str, optional): The name of the file if provided that will be used to save the coef files as .h5. + Default is ''. + Returns: + tuple: A tuple containing the basis and the coefficients of the expansion. + The basis is an instance of pyEXP.basis.Basis, and the coefficients are + an instance of pyEXP.coefs.Coefs. + """ + + if os.path.isfile(modelname) == False: + print("-> File model not found so we are computing one \n") + if empirical == True: + print('-> Computing empirical model') + rad, rho = empirical_density_profile(pos, mass, nbins=numr) + elif empirical == False: + makemodel.hernquist_halo() + + makemodel.Profiles(density_profile) + + R, D, M, P = makemodel.makemodel(rvals=np.logspace(np.log10(rmin), np.log10(rmax), numr), rho, outfile=modelname, return_values=True) + print('-> Computing analytical Hernquist model') + R, D, M, P = makemodel.makemodel(hernquist_halo, 1, [scale], rvals=, pfile=modelname) + print('-> Model computed: rmin={}, rmax={}, numr={}'.format(R[0], R[-1], len(R))) + else: + R, D, M, P = np.loadtxt(modelname, skiprows=3, unpack=True) + # check if config file is passed + if config is None: + print('No config file provided.') + print(f'Computing empirical density') + #rad, rho = empirical_density_profile(pos, mass, nbins=500) + #makemodel_empirical(r_exact, rho, outfile=modelname) + #R = [0.01, 600] + config = make_config(basis_id, numr, R[0], R[-1], lmax, nmax, scale, + modelname, cachename) + + # Construct the basis instances + basis = pyEXP.basis.Basis.factory(config) + + # Prints info from Cache + basis.cacheInfo(cachename) + + #compute coefficients + coef = basis.createFromArray(mass/np.sum(mass), pos, time=time) + coefs = pyEXP.coefs.Coefs.makecoefs(coef, 'dark halo') + coefs.add(coef) + if add_coef == False: + coefs.WriteH5Coefs(coef_file) + elif add_coef == True: + coefs.ExtendH5Coefs(coef_file) + + return basis, coefs + diff --git a/EXPtools/basis_builder/makemodel.py b/EXPtools/basis_builder/makemodel.py new file mode 100644 index 0000000..5c1aa31 --- /dev/null +++ b/EXPtools/basis_builder/makemodel.py @@ -0,0 +1,99 @@ +""" +TODO: +""" + +import numpy as np + + +def makemodel(dvals, rvals, M, unit_physical=False, + pfile='', plabel = '', verbose=True): + """ + Make an EXP-compatible spherical basis function table + + inputs + ------------- + basis_type : (str) + func : (function) the callable functional form of the density + M : (float) the total mass of the model, sets normalisations + funcargs : (list) a list of arguments for the density function. + rvals : (array of floats) radius values to evaluate the density function + = 10.**np.linspace(-2.,4.,2000) + pfile : (string) the name of the output file. If '', will not print file + plabel : (string) comment string + verbose : (boolean) + + outputs + ------------- + R : (array of floats) the radius values + D : (array of floats) the density + M : (array of floats) the mass enclosed + P : (array of floats) the potential + + """ + + R = np.nanmax(rvals) + + # query out the density values + + # make the mass and potential arrays + mvals = np.zeros(dvals.size) + pvals = np.zeros(dvals.size) + pwvals = np.zeros(dvals.size) + + # initialise the mass enclosed an potential energy + mvals[0] = 1.e-15 + pwvals[0] = 0. + + # evaluate mass enclosed and potential energy by recursion + for indx in range(1,dvals.size): + mvals[indx] = mvals[indx-1] +\ + 2.0*np.pi*(rvals[indx-1]*rvals[indx-1]*dvals[indx-1] +\ + rvals[indx]*rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); + pwvals[indx] = pwvals[indx-1] + \ + 2.0*np.pi*(rvals[indx-1]*dvals[indx-1] + rvals[indx]*dvals[indx])*(rvals[indx] - rvals[indx-1]); + + # evaluate potential (see theory document) + pvals = -mvals/(rvals+1.e-10) - (pwvals[dvals.size-1] - pwvals) + + # get the maximum mass and maximum radius + M0 = mvals[dvals.size-1] + R0 = rvals[dvals.size-1] + + # compute scaling factors + Beta = (M/M0) * (R0/R); + Gamma = np.sqrt((M0*R0)/(M*R)) * (R0/R); + if verbose: + print("! Scaling: R=",R," M=",M) + + rfac = np.power(Beta,-0.25) * np.power(Gamma,-0.5); + dfac = np.power(Beta,1.5) * Gamma; + mfac = np.power(Beta,0.75) * np.power(Gamma,-0.5); + pfac = Beta; + + if verbose: + print(rfac,dfac,mfac,pfac) + + # save file if desired + if pfile != '': + f = open(pfile,'w') + print('! ',plabel,file=f) + print('! R D M P',file=f) + + print(rvals.size,file=f) + + if unit_physical = True: + for indx in range(0,rvals.size): + print('{0} {1} {2} {3}'.format(rvals[indx],\ + dvals[indx],\ + mvals[indx],\ + pvals[indx]),file=f) + else: + + f.close() + + if return_val: + if unit_phy: + return rvals, dvals, mvals, pvals + else: + return rvals*rfac, dfac*dvals, mfac*mvals, pfac*pvals + diff --git a/EXPtools/basis_builder/profiles.py b/EXPtools/basis_builder/profiles.py new file mode 100644 index 0000000..fe4f7e3 --- /dev/null +++ b/EXPtools/basis_builder/profiles.py @@ -0,0 +1,127 @@ +import numpy as np +from scipy import special + +Class Profiles: + def __init__(self, radius, scale_radius, alpha=1., beta=2.0): + """ + Class with density profiles of Dark Matter halos + + inputs + ---------- + r : (float) radius values + rs : (float, default=1.) scale radius + alpha : (float, default=1.) inner halo slope + beta : (float, default=1.e-7) outer halo slope + + notes + ---------- + different combinations are known distributions. + alpha=1,beta=2 is NFW (default) + alpha=1,beta=3 is Hernquist + alpha=2.5, beta=0 is a typical single power law halo + + """ + + self.radius = radius + self.scale_radius = scale_radius + self.alpha = alpha + self.beta = beta + self.ra = self.radius/self.scale_radius + + + def powerhalo(self, rc=0.0): + """ + Generic two power law distribution + + inputs + ---------- + rc : (float, default=0. i.e. no core) core radius + + returns + ---------- + densities evaluated at r + + notes + ---------- + different combinations are known distributions. + alpha=1,beta=2 is NFW + alpha=1,beta=3 is Hernquist + alpha=2.5,beta=0 is a typical single power law halo + + + + """ + return 1./(((self.ra+rc/self.scale_radius)**self.alpha)*((1+self.ra)**self.beta)) + + + def powerhalorolloff(self, rc=0.0): + """ + Generic twopower law distribution with an erf rolloff + + + inputs + ---------- + rc : (float, default=0. i.e. no core) core radius + + returns + ---------- + densities evaluated at r + + notes + ---------- + different combinations are known distributions. + alpha=1,beta=2 is NFW + alpha=1,beta=3 is Hernquist + alpha=2.5,beta=0 is a typical single power law halo + + + """ + dens = 1./(((self.ra+rc/self.scale_radius)**self.alpha)*((1+self.ra)**self.beta)) + rtrunc = 25*self.scale_radius + wtrunc = rtrunc*0.2 + rolloff = 0.5 - 0.5*special.erf((self.radius-rtrunc)/wtrunc) + return dens*rolloff + + + def plummer_density(self): + """ + Plummer density profile + + inputs + --------- + + + returns + ---------- + densities evaluated at r + """ + return ((3.0)/(4*np.pi))*(self.scale_radius**2.)*((self.scale_radius**2 + self.radius**2)**(-2.5)) + + def twopower_density_withrolloff(self, rcen, wcen): + """ + inputs + --------- + + + returns + ---------- + densities evaluated at r + """ + + + prefac = 0.5*(1.-scipy.special.erf((ra-rcen)/wcen)) + return prefac*(ra**-self.alpha)*(1+ra)**(-self.beta+self.alpha) + + def hernquist_halo(self): + """ + TODO: DO we need these profile + Hernquist halo + + inputs + --------- + + returns + ---------- + densities evaluated at r + """ + return 1 / ( 2*np.pi * (self.radius/self.scale_radius) * (1 + self.radius/self.scale_radius)**3) \ No newline at end of file diff --git a/EXPtools/visuals/visualize.py b/EXPtools/visuals/visualize.py new file mode 100644 index 0000000..8db665d --- /dev/null +++ b/EXPtools/visuals/visualize.py @@ -0,0 +1,304 @@ +import os, sys, pickle, pyEXP +import numpy as np + +###Field computations for plotting### +def make_basis_plot(basis, savefile=None, nsnap='mean', y=0.92, dpi=200): + """ + Plots the potential of the basis functions for different values of l and n. + + Args: + basis (obj): object containing the basis functions for the simulation + savefile (str, optional): name of the file to save the plot as + nsnap (str, optional): description of the snapshot being plotted + y (float, optional): vertical position of the main title + dpi (int, optional): resolution of the plot in dots per inch + + Returns: + None + + """ + # Set up grid for plotting potential + lrmin, lrmax, rnum = 0.5, 2.7, 100 + halo_grid = basis.getBasis(lrmin, lrmax, rnum) + r = np.linspace(lrmin, lrmax, rnum) + r = np.power(10.0, r) + + # Create subplots and plot potential for each l and n + import matplotlib.pyplot as plt + fig, ax = plt.subplots(4, 5, figsize=(6,6), dpi=dpi, + sharex='col', sharey='row') + plt.subplots_adjust(wspace=0, hspace=0) + ax = ax.flatten() + + for l in range(len(ax)): + ax[l].set_title(f"$\ell = {l}$", y=0.8, fontsize=6) + for n in range(20): + ax[l].semilogx(r, halo_grid[l][n]['potential'], '-', label="n={}".format(n), lw=0.5) + + # Add labels and main title + fig.supylabel('Potential', weight='bold', x=-0.02) + fig.supxlabel('Radius', weight='bold', y=0.02) + fig.suptitle(f'nsnap = {nsnap}', + fontsize=12, + weight='bold', + y=y, + ) + + # Save plot if a filename was provided + if savefile: + plt.savefig(f'{savefile}', bbox_inches='tight') + +def find_field(basis, coefficients, time=0, xyz=(0, 0, 0), property='dens', include_monopole=True): + """ + Finds the value of the specified property of the field at the given position. + + Args: + basis (obj): Object containing the basis functions for the simulation. + coefficients (obj): Object containing the coefficients for the simulation. + time (float, optional): The time at which to evaluate the field. Default is 0. + xyz (tuple or list, optional): The (x, y, z) position at which to evaluate the field. Default is (0, 0, 0). + property (str, optional): The property of the field to evaluate. Can be 'dens', 'pot', or 'force'. Default is 'dens'. + include_monopole (bool, optional): Whether to return the monopole contribution to the property only. Default is True. + + Returns: + float or list: The value of the specified property of the field at the given position. If property is 'force', a list of + three values is returned representing the force vector in (x, y, z) directions. + + Raises: + ValueError: If the property argument is not 'dens', 'pot', or 'force'. + """ + + coefficients.set_coefs(coefficients.getCoefStruct(time)) + dens0, pot0, dens, pot, fx, fy, fz = basis.getFields(xyz[0], xyz[1], xyz[2]) + + if property == 'dens': + if include_monopole: + return dens0 + return dens + dens0 + + elif property == 'pot': + if include_monopole: + return pot0 + return pot + pot0 + + elif property == 'force': + return [fx, fy, fz] + + else: + raise ValueError("Invalid property specified. Possible values are 'dens', 'pot', and 'force'.") + +def spherical_avg_prop(basis, coefficients, time=0, radius=np.linspace(0.1, 600, 100), property='dens'): + """ + Computes the spherically averaged value of the specified property of the field over the given radii. + + Args: + basis (obj): Object containing the basis functions for the simulation. + coefficients (obj): Object containing the coefficients for the simulation. + time (float, optional): The time at which to evaluate the field. Default is 0. + radius (ndarray, optional): An array of radii over which to compute the spherically averaged property. Default is an + array of 100 values logarithmically spaced between 0.1 and 600. + property (str, optional): The property of the field to evaluate. Can be 'dens', 'pot', or 'force'. Default is 'dens'. + + Returns: + ndarray: An array of spherically averaged values of the specified property over the given radii. + + Raises: + ValueError: If the property argument is not 'dens', 'pot', or 'force'. + """ + + coefficients.set_coefs(coefficients.getCoefStruct(time)) + field = [find_field(basis, np.hstack([[rad], [0], [0]]), property=property, include_monopole=True) for rad in radius] + + if property == 'force': + return np.vstack(field), radius + + return np.array(field), radius + +def return_fields_in_grid(basis, coefficients, times=[0], + projection='3D', proj_plane=0, + grid_lim=300, npoints=150,): + + """ + Returns the 2D or 3D field grids as a dict at each time step as a key. + + Args: + basis (obj): object containing the basis functions for the simulation + coefficients (obj): object containing the coefficients for the simulation + times (list): list of the time to compute the field grid. + projection (str): the slice projection to plot. Can be '3D', XY', 'XZ', or 'YZ'. + proj_plane (float, optional): the value of the coordinate that is held constant in the 2D slice projection + npoints (int, optional): the number of grid points in each dimension + grid_limits (float, optional): the limits of the grid in the dimensions. + + Returns: fields along with the specified 2D or 3D grid. The fields are in a + heirarchical dict structure with for each time: + dict_keys(['d', 'd0', 'd1', 'dd', 'fp', 'fr', 'ft', 'p', 'p0', 'p1']) + where + - 'd': Total density at each point. + - 'd0' : Density in the monopole term. + - 'd1' : Density in l>0 terms. + - 'dd': Density difference. + - 'fp': Force in the phi direction. + - 'fr': Force in the radial direction. + - 'ft': Force in the theta direction. + - 'p': Total potential at each point. + - 'p0': Potential in the monopole term. + - 'p1': Potential in l>0 terms. + """ + + assert projection in ['3D', 'XY', 'XZ', 'YZ'], "Invalid value for 'projection'. Must be one of '3D', 'XY', 'XZ', or 'YZ'." + + + + if projection == '3D': + pmin = [-grid_lim, -grid_lim, -grid_lim] + pmax = [grid_lim, grid_lim, grid_lim] + grid = [npoints, npoints, npoints] + + ##create a projection grid along with it + x = np.linspace(-grid_lim, grid_lim, npoints) + xgrid = np.meshgrid(x, x, x) + + else: + x = np.linspace(-grid_lim, grid_lim, npoints) + xgrid = np.meshgrid(x, x) + grid = [npoints, npoints] + + if projection == 'XY': + pmin = [-grid_lim, -grid_lim, proj_plane] + pmax = [grid_lim, grid_lim, proj_plane] + + elif projection == 'XZ': + pmin = [-grid_lim, proj_plane, -grid_lim] + pmax = [grid_lim, proj_plane, grid_lim] + + elif projection == 'YZ': + pmin = [proj_plane, -grid_lim, -grid_lim] + pmax = [proj_plane, grid_lim, grid_lim] + + + field_gen = pyEXP.field.FieldGenerator(times, pmin, pmax, grid) + + return field_gen.volumes(basis, coefs), xgrid + +def slice_fields(basis, coefficients, time=0, + projection='XY', proj_plane=0, npoints=300, + grid_limits=(-300, 300), prop='dens', monopole_only=False): + """ + Plots a slice projection of the fields of a simulation. + + Args: + basis (obj): object containing the basis functions for the simulation + coefficients (obj): object containing the coefficients for the simulation + time (float): the time at which to plot the fields + projection (str): the slice projection to plot. Can be 'XY', 'XZ', or 'YZ'. + proj_plane (float, optional): the value of the coordinate that is held constant in the slice projection + npoints (int, optional): the number of grid points in each dimension + grid_limits (tuple, optional): the limits of the grid in the x and y dimensions, in the form (x_min, x_max) + prop (str, optional): the property to return. Can be 'dens' (density), 'pot' (potential), or 'force' (force). + monopole_only (bool, optional): whether to return the monopole component in the returned property value. + + Returns: + array or list: the property specified by `prop`. If `prop` is 'force', a list of the x, y, and z components of the force is returned. + Also returns the grid used to compute slice fields. + """ + x = np.linspace(grid_limits[0], grid_limits[1], npoints) + xgrid = np.meshgrid(x, x) + xg = xgrid[0].flatten() + yg = xgrid[1].flatten() + + + if projection not in ['XY', 'XZ', 'YZ']: + raise ValueError("Invalid projection specified. Possible values are 'XY', 'XZ', and 'YZ'.") + + N = len(xg) + rho0 = np.zeros_like(xg) + pot0 = np.zeros_like(xg) + rho = np.zeros_like(xg) + pot = np.zeros_like(xg) + basis.set_coefs(coefficients.getCoefStruct(time)) + + for k in range(0, N): + if projection == 'XY': + rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(xg[k], yg[k], proj_plane) + elif projection == 'XZ': + rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(xg[k], proj_plane, yg[k]) + elif projection == 'YZ': + rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(proj_plane, xg[k], yg[k]) + + dens = rho.reshape(npoints, npoints) + pot = pot.reshape(npoints, npoints) + dens0 = rho0.reshape(npoints, npoints) + pot0 = pot0.reshape(npoints, npoints) + + if prop == 'dens': + if monopole_only: + return dens0 + return dens0, dens, xgrid + + if prop == 'pot': + if monopole_only: + return pot0 + return pot0, pot, xgrid + + if prop == 'force': + return [fx.reshape(npoints, npoints), fy.reshape(npoints, npoints), fz.reshape(npoints, npoints)], xgrid + + +def slice_3d_fields(basis, coefficients, time=0, + projection='XY', proj_plane=0, npoints=300, + grid_limits=(-300, 300), prop='dens', monopole_only=False): + """ + Plots a slice projection of the fields of a simulation. + + Args: + basis (obj): object containing the basis functions for the simulation + coefficients (obj): object containing the coefficients for the simulation + time (float): the time at which to plot the fields + projection (str): the slice projection to plot. Can be 'XY', 'XZ', or 'YZ'. + proj_plane (float, optional): the value of the coordinate that is held constant in the slice projection + npoints (int, optional): the number of grid points in each dimension + grid_limits (tuple, optional): the limits of the grid in the x and y dimensions, in the form (x_min, x_max) + prop (str, optional): the property to return. Can be 'dens' (density), 'pot' (potential), or 'force' (force). + monopole_only (bool, optional): whether to return the monopole component in the returned property value. + + Returns: + array or list: the property specified by `prop`. If `prop` is 'force', a list of the x, y, and z components of the force is returned. + Also returns the grid used to compute slice fields. + """ + x = np.linspace(grid_limits[0], grid_limits[1], npoints) + xgrid = np.meshgrid(x, x, x) + xg = xgrid[0].flatten() + yg = xgrid[1].flatten() + zg = xgrid[2].flatten() + + + + N = len(xg) + rho0 = np.zeros_like(xg) + pot0 = np.zeros_like(xg) + rho = np.zeros_like(xg) + pot = np.zeros_like(xg) + basis.set_coefs(coefficients.getCoefStruct(time)) + + for k in range(0, N): + rho0[k], pot0[k], rho[k], pot[k], fx, fy, fz = basis.getFields(xg[k], yg[k], zg[k]) + + dens = rho.reshape(npoints, npoints, npoints) + pot = pot.reshape(npoints, npoints, npoints) + dens0 = rho0.reshape(npoints, npoints, npoints) + pot0 = pot0.reshape(npoints, npoints, npoints) + + if prop == 'dens': + if monopole_only: + return dens0 + return dens0, dens, xgrid + + if prop == 'pot': + if monopole_only: + return pot0 + return pot0, pot, xgrid + + if prop == 'force': + return [fx.reshape(npoints, npoints, npoints), fy.reshape(npoints, npoints, npoints), fz.reshape(npoints, npoints, npoints)], xgrid + From 9abc3528259c3fa1ed5947cea551e24fb94bcd18 Mon Sep 17 00:00:00 2001 From: jngaravitoc Date: Sat, 11 Nov 2023 15:56:54 -0500 Subject: [PATCH 07/19] tests --- tests/test_basis_utils.py | 4 ++++ tests/test_profiles.py | 6 ++++++ 2 files changed, 10 insertions(+) create mode 100644 tests/test_basis_utils.py create mode 100644 tests/test_profiles.py diff --git a/tests/test_basis_utils.py b/tests/test_basis_utils.py new file mode 100644 index 0000000..401831b --- /dev/null +++ b/tests/test_basis_utils.py @@ -0,0 +1,4 @@ +import numpy as np +from EXPtools import basis_utils + +# test integration between functions! diff --git a/tests/test_profiles.py b/tests/test_profiles.py new file mode 100644 index 0000000..6128a23 --- /dev/null +++ b/tests/test_profiles.py @@ -0,0 +1,6 @@ +import numpy as np +from EXPtools import profiles + +# Test output of functions +# Test document test functions + From f34d452878c7a763eeca1135b6d8dd3117c9d315 Mon Sep 17 00:00:00 2001 From: jngaravitoc Date: Sat, 11 Nov 2023 16:58:04 -0500 Subject: [PATCH 08/19] adding the first test --- EXPtools/basis_builder/__init__.py | 2 ++ EXPtools/basis_builder/basis_utils.py | 7 ++++--- EXPtools/basis_builder/makemodel.py | 17 ++++++++++------ EXPtools/basis_builder/profiles.py | 4 ++-- setup.py | 2 +- tests/test_basis_utils.py | 2 +- tests/test_profiles.py | 28 ++++++++++++++++++++++++++- 7 files changed, 48 insertions(+), 14 deletions(-) diff --git a/EXPtools/basis_builder/__init__.py b/EXPtools/basis_builder/__init__.py index 47cd4c4..d07fca5 100644 --- a/EXPtools/basis_builder/__init__.py +++ b/EXPtools/basis_builder/__init__.py @@ -1,3 +1,5 @@ +from . import makemodel +from . import profiles from .makemodel import makemodel from .profiles import Profiles from .basis_utils import makebasis diff --git a/EXPtools/basis_builder/basis_utils.py b/EXPtools/basis_builder/basis_utils.py index 49f3395..0afd0d0 100644 --- a/EXPtools/basis_builder/basis_utils.py +++ b/EXPtools/basis_builder/basis_utils.py @@ -1,6 +1,6 @@ import os, sys, pickle, pyEXP import numpy as np -import makemodel +from . import makemodel def make_config(basis_id, numr, rmin, rmax, lmax, nmax, scale, modelname='', cachename='.slgrid_sph_cache'): @@ -239,7 +239,7 @@ def makebasis(pos, mass, basis_model, config=None, basis_id='sphereSL', time=0, The basis is an instance of pyEXP.basis.Basis, and the coefficients are an instance of pyEXP.coefs.Coefs. """ - + """ if os.path.isfile(modelname) == False: print("-> File model not found so we are computing one \n") if empirical == True: @@ -282,4 +282,5 @@ def makebasis(pos, mass, basis_model, config=None, basis_id='sphereSL', time=0, coefs.ExtendH5Coefs(coef_file) return basis, coefs - + """ + return None diff --git a/EXPtools/basis_builder/makemodel.py b/EXPtools/basis_builder/makemodel.py index 5c1aa31..ddee600 100644 --- a/EXPtools/basis_builder/makemodel.py +++ b/EXPtools/basis_builder/makemodel.py @@ -81,14 +81,19 @@ def makemodel(dvals, rvals, M, unit_physical=False, print(rvals.size,file=f) - if unit_physical = True: + if unit_physical == True: for indx in range(0,rvals.size): - print('{0} {1} {2} {3}'.format(rvals[indx],\ - dvals[indx],\ - mvals[indx],\ - pvals[indx]),file=f) + print('{0} {1} {2} {3}'.format(rvals[indx], + dvals[indx], + mvals[indx], + pvals[indx]),file=f) else: - + for indx in range(0,rvals.size): + print('{0} {1} {2} {3}'.format(rvals[indx], + dvals[indx]*dfac, + mvals[indx]*mfac, + pvals[indx]*pfac),file=f) + f.close() if return_val: diff --git a/EXPtools/basis_builder/profiles.py b/EXPtools/basis_builder/profiles.py index fe4f7e3..efb18f8 100644 --- a/EXPtools/basis_builder/profiles.py +++ b/EXPtools/basis_builder/profiles.py @@ -1,7 +1,7 @@ import numpy as np from scipy import special -Class Profiles: +class Profiles: def __init__(self, radius, scale_radius, alpha=1., beta=2.0): """ Class with density profiles of Dark Matter halos @@ -124,4 +124,4 @@ def hernquist_halo(self): ---------- densities evaluated at r """ - return 1 / ( 2*np.pi * (self.radius/self.scale_radius) * (1 + self.radius/self.scale_radius)**3) \ No newline at end of file + return 1 / ( 2*np.pi * (self.radius/self.scale_radius) * (1 + self.radius/self.scale_radius)**3) diff --git a/setup.py b/setup.py index c81664a..64d6e2b 100755 --- a/setup.py +++ b/setup.py @@ -6,5 +6,5 @@ author="B-BFE collaboration", author_email="ngaravito@flatironinstitute.org", description="Analysis tools for EXP and other BFE packages", - packages=["EXP_tools"], + packages=["EXPtools", "EXPtools/basis_builder/", "EXPtools/visuals"], ) diff --git a/tests/test_basis_utils.py b/tests/test_basis_utils.py index 401831b..0c1c514 100644 --- a/tests/test_basis_utils.py +++ b/tests/test_basis_utils.py @@ -1,4 +1,4 @@ import numpy as np -from EXPtools import basis_utils +from EXPtools import basis_builder # test integration between functions! diff --git a/tests/test_profiles.py b/tests/test_profiles.py index 6128a23..2fc671f 100644 --- a/tests/test_profiles.py +++ b/tests/test_profiles.py @@ -1,6 +1,32 @@ import numpy as np -from EXPtools import profiles +from EXPtools.basis_builder import profiles # Test output of functions # Test document test functions +def test_profile_class(): + """ + Density profiles tests + """ + r = np.logspace(0, 3, 100) + scale_length = 1 + + alpha1 = 1 + beta2 = 2 + alpha2 = 2.5 + beta3 = 3 + beta0 = 0 + + # NFW + NFW = profiles.Profiles(r, scale_length, alpha1, beta2) + rho_NFW = NFW.powerhalo() + + # Hernquist + Hernquist = profiles.Profiles(r, scale_length, alpha1, beta3) + rho_Hern = Hernquist.powerhalo() + + # Single power law + SPL = profiles.Profiles(r, scale_length, alpha2, beta0) + rho_SPL = SPL.powerhalo() + + return None From 1d2e303b8829c7cf75c58af891175b0770c92ee4 Mon Sep 17 00:00:00 2001 From: NGC Date: Sat, 11 Nov 2023 17:06:11 -0500 Subject: [PATCH 09/19] Create python-package.yml --- .github/workflows/python-package.yml | 40 ++++++++++++++++++++++++++++ 1 file changed, 40 insertions(+) create mode 100644 .github/workflows/python-package.yml diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml new file mode 100644 index 0000000..fd04aae --- /dev/null +++ b/.github/workflows/python-package.yml @@ -0,0 +1,40 @@ +# This workflow will install Python dependencies, run tests and lint with a variety of Python versions +# For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-python + +name: Python package + +on: + push: + branches: [ "main" ] + pull_request: + branches: [ "main" ] + +jobs: + build: + + runs-on: ubuntu-latest + strategy: + fail-fast: false + matrix: + python-version: ["3.9", "3.10", "3.11"] + + steps: + - uses: actions/checkout@v3 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v3 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + python -m pip install flake8 pytest + if [ -f requirements.txt ]; then pip install -r requirements.txt; fi + - name: Lint with flake8 + run: | + # stop the build if there are Python syntax errors or undefined names + flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics + # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide + flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics + - name: Test with pytest + run: | + pytest From 8126c9c877be90a1138733d744d19274a87cba36 Mon Sep 17 00:00:00 2001 From: jngaravitoc Date: Sat, 11 Nov 2023 17:09:44 -0500 Subject: [PATCH 10/19] scipy special --- EXPtools/basis_builder/profiles.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/EXPtools/basis_builder/profiles.py b/EXPtools/basis_builder/profiles.py index efb18f8..dc858be 100644 --- a/EXPtools/basis_builder/profiles.py +++ b/EXPtools/basis_builder/profiles.py @@ -109,12 +109,12 @@ def twopower_density_withrolloff(self, rcen, wcen): """ - prefac = 0.5*(1.-scipy.special.erf((ra-rcen)/wcen)) + prefac = 0.5*(1.-special.erf((ra-rcen)/wcen)) return prefac*(ra**-self.alpha)*(1+ra)**(-self.beta+self.alpha) def hernquist_halo(self): """ - TODO: DO we need these profile + TODO: Do we need these profile Hernquist halo inputs From 67222b9f992faebd886115173b0f0a1d2a8e83d5 Mon Sep 17 00:00:00 2001 From: jngaravitoc Date: Sat, 11 Nov 2023 17:12:50 -0500 Subject: [PATCH 11/19] physical units variable name --- EXPtools/basis_builder/makemodel.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/EXPtools/basis_builder/makemodel.py b/EXPtools/basis_builder/makemodel.py index ddee600..6489ce4 100644 --- a/EXPtools/basis_builder/makemodel.py +++ b/EXPtools/basis_builder/makemodel.py @@ -5,7 +5,7 @@ import numpy as np -def makemodel(dvals, rvals, M, unit_physical=False, +def makemodel(dvals, rvals, M, unit_physical=False, return_val=False, pfile='', plabel = '', verbose=True): """ Make an EXP-compatible spherical basis function table @@ -96,8 +96,8 @@ def makemodel(dvals, rvals, M, unit_physical=False, f.close() - if return_val: - if unit_phy: + if return_val==True: + if unit_physical==True: return rvals, dvals, mvals, pvals else: return rvals*rfac, dfac*dvals, mfac*mvals, pfac*pvals From 84d9b61f3585a1c79f97e67883b68f9fb9ab7ba5 Mon Sep 17 00:00:00 2001 From: jngaravitoc Date: Sat, 11 Nov 2023 17:14:17 -0500 Subject: [PATCH 12/19] ra -> self.ra --- EXPtools/basis_builder/profiles.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/EXPtools/basis_builder/profiles.py b/EXPtools/basis_builder/profiles.py index dc858be..cc7a6d2 100644 --- a/EXPtools/basis_builder/profiles.py +++ b/EXPtools/basis_builder/profiles.py @@ -109,8 +109,8 @@ def twopower_density_withrolloff(self, rcen, wcen): """ - prefac = 0.5*(1.-special.erf((ra-rcen)/wcen)) - return prefac*(ra**-self.alpha)*(1+ra)**(-self.beta+self.alpha) + prefac = 0.5*(1.-special.erf((self.ra-rcen)/wcen)) + return prefac*(self.ra**-self.alpha)*(1+self.ra)**(-self.beta+self.alpha) def hernquist_halo(self): """ From 4d86db188994917880eebfed8066140ca6d6b4da Mon Sep 17 00:00:00 2001 From: jngaravitoc Date: Sat, 11 Nov 2023 17:16:27 -0500 Subject: [PATCH 13/19] coefs -> coefficients --- EXPtools/visuals/visualize.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/EXPtools/visuals/visualize.py b/EXPtools/visuals/visualize.py index 8db665d..f0aae3f 100644 --- a/EXPtools/visuals/visualize.py +++ b/EXPtools/visuals/visualize.py @@ -179,7 +179,7 @@ def return_fields_in_grid(basis, coefficients, times=[0], field_gen = pyEXP.field.FieldGenerator(times, pmin, pmax, grid) - return field_gen.volumes(basis, coefs), xgrid + return field_gen.volumes(basis, coefficients), xgrid def slice_fields(basis, coefficients, time=0, projection='XY', proj_plane=0, npoints=300, From f305ca57af10c84afdef1f643cb8f62e6f8d08c3 Mon Sep 17 00:00:00 2001 From: NGC Date: Sat, 11 Nov 2023 17:19:00 -0500 Subject: [PATCH 14/19] removeing pytest in CI --- .github/workflows/python-package.yml | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/.github/workflows/python-package.yml b/.github/workflows/python-package.yml index fd04aae..841659e 100644 --- a/.github/workflows/python-package.yml +++ b/.github/workflows/python-package.yml @@ -35,6 +35,7 @@ jobs: flake8 . --count --select=E9,F63,F7,F82 --show-source --statistics # exit-zero treats all errors as warnings. The GitHub editor is 127 chars wide flake8 . --count --exit-zero --max-complexity=10 --max-line-length=127 --statistics - - name: Test with pytest - run: | - pytest + # Tests not working at the moment -- pyEXP requirement! + #- name: Test with pytest + # run: | + # pytest From c8eaee5b3769b6123d80f3060ecf34f4b9c1e958 Mon Sep 17 00:00:00 2001 From: NGC Date: Sat, 11 Nov 2023 17:21:33 -0500 Subject: [PATCH 15/19] Update README.md --- README.md | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/README.md b/README.md index 7930a4f..a41c507 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,9 @@ # EXP_tools This is a shared repo with all the tools developed for cosmological+external simulations + +To install EXPtools run: + +```bash +python -m pip install . +``` +[![Python package](https://github.com/jngaravitoc/EXP_tools/actions/workflows/python-package.yml/badge.svg)](https://github.com/jngaravitoc/EXP_tools/actions/workflows/python-package.yml) From fb5861ea00eba0208174808e41f3839b59b8c956 Mon Sep 17 00:00:00 2001 From: NGC Date: Sat, 11 Nov 2023 17:25:22 -0500 Subject: [PATCH 16/19] Update README.md --- README.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/README.md b/README.md index a41c507..a1e84df 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ -# EXP_tools -This is a shared repo with all the tools developed for cosmological+external simulations +# EXPtools +EXPtools is a Python library that contains analysis routines for [pyEXP](https://github.com/EXP-code/EXP). EXPtools was designed with the aim of analyze cosmological and idealized N-body simulations -To install EXPtools run: +### Installation: ```bash python -m pip install . From da31a3005de809919f07f00fdf31d546916998f9 Mon Sep 17 00:00:00 2001 From: NGC Date: Sat, 11 Nov 2023 21:53:55 -0500 Subject: [PATCH 17/19] Create LICENSE --- LICENSE | 339 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 339 insertions(+) create mode 100644 LICENSE diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..d159169 --- /dev/null +++ b/LICENSE @@ -0,0 +1,339 @@ + GNU GENERAL PUBLIC LICENSE + Version 2, June 1991 + + Copyright (C) 1989, 1991 Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +License is intended to guarantee your freedom to share and change free +software--to make sure the software is free for all its users. This +General Public License applies to most of the Free Software +Foundation's software and to any other program whose authors commit to +using it. (Some other Free Software Foundation software is covered by +the GNU Lesser General Public License instead.) You can apply it to +your programs, too. + + When we speak of free software, we are referring to freedom, not +price. Our General Public Licenses are designed to make sure that you +have the freedom to distribute copies of free software (and charge for +this service if you wish), that you receive source code or can get it +if you want it, that you can change the software or use pieces of it +in new free programs; and that you know you can do these things. + + To protect your rights, we need to make restrictions that forbid +anyone to deny you these rights or to ask you to surrender the rights. +These restrictions translate to certain responsibilities for you if you +distribute copies of the software, or if you modify it. + + For example, if you distribute copies of such a program, whether +gratis or for a fee, you must give the recipients all the rights that +you have. You must make sure that they, too, receive or can get the +source code. And you must show them these terms so they know their +rights. + + We protect your rights with two steps: (1) copyright the software, and +(2) offer you this license which gives you legal permission to copy, +distribute and/or modify the software. + + Also, for each author's protection and ours, we want to make certain +that everyone understands that there is no warranty for this free +software. If the software is modified by someone else and passed on, we +want its recipients to know that what they have is not the original, so +that any problems introduced by others will not reflect on the original +authors' reputations. + + Finally, any free program is threatened constantly by software +patents. We wish to avoid the danger that redistributors of a free +program will individually obtain patent licenses, in effect making the +program proprietary. To prevent this, we have made it clear that any +patent must be licensed for everyone's free use or not licensed at all. + + The precise terms and conditions for copying, distribution and +modification follow. + + GNU GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License applies to any program or other work which contains +a notice placed by the copyright holder saying it may be distributed +under the terms of this General Public License. The "Program", below, +refers to any such program or work, and a "work based on the Program" +means either the Program or any derivative work under copyright law: +that is to say, a work containing the Program or a portion of it, +either verbatim or with modifications and/or translated into another +language. (Hereinafter, translation is included without limitation in +the term "modification".) Each licensee is addressed as "you". + +Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running the Program is not restricted, and the output from the Program +is covered only if its contents constitute a work based on the +Program (independent of having been made by running the Program). +Whether that is true depends on what the Program does. + + 1. You may copy and distribute verbatim copies of the Program's +source code as you receive it, in any medium, provided that you +conspicuously and appropriately publish on each copy an appropriate +copyright notice and disclaimer of warranty; keep intact all the +notices that refer to this License and to the absence of any warranty; +and give any other recipients of the Program a copy of this License +along with the Program. + +You may charge a fee for the physical act of transferring a copy, and +you may at your option offer warranty protection in exchange for a fee. + + 2. You may modify your copy or copies of the Program or any portion +of it, thus forming a work based on the Program, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) You must cause the modified files to carry prominent notices + stating that you changed the files and the date of any change. + + b) You must cause any work that you distribute or publish, that in + whole or in part contains or is derived from the Program or any + part thereof, to be licensed as a whole at no charge to all third + parties under the terms of this License. + + c) If the modified program normally reads commands interactively + when run, you must cause it, when started running for such + interactive use in the most ordinary way, to print or display an + announcement including an appropriate copyright notice and a + notice that there is no warranty (or else, saying that you provide + a warranty) and that users may redistribute the program under + these conditions, and telling the user how to view a copy of this + License. (Exception: if the Program itself is interactive but + does not normally print such an announcement, your work based on + the Program is not required to print an announcement.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Program, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Program, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Program. + +In addition, mere aggregation of another work not based on the Program +with the Program (or with a work based on the Program) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may copy and distribute the Program (or a work based on it, +under Section 2) in object code or executable form under the terms of +Sections 1 and 2 above provided that you also do one of the following: + + a) Accompany it with the complete corresponding machine-readable + source code, which must be distributed under the terms of Sections + 1 and 2 above on a medium customarily used for software interchange; or, + + b) Accompany it with a written offer, valid for at least three + years, to give any third party, for a charge no more than your + cost of physically performing source distribution, a complete + machine-readable copy of the corresponding source code, to be + distributed under the terms of Sections 1 and 2 above on a medium + customarily used for software interchange; or, + + c) Accompany it with the information you received as to the offer + to distribute corresponding source code. (This alternative is + allowed only for noncommercial distribution and only if you + received the program in object code or executable form with such + an offer, in accord with Subsection b above.) + +The source code for a work means the preferred form of the work for +making modifications to it. For an executable work, complete source +code means all the source code for all modules it contains, plus any +associated interface definition files, plus the scripts used to +control compilation and installation of the executable. However, as a +special exception, the source code distributed need not include +anything that is normally distributed (in either source or binary +form) with the major components (compiler, kernel, and so on) of the +operating system on which the executable runs, unless that component +itself accompanies the executable. + +If distribution of executable or object code is made by offering +access to copy from a designated place, then offering equivalent +access to copy the source code from the same place counts as +distribution of the source code, even though third parties are not +compelled to copy the source along with the object code. + + 4. You may not copy, modify, sublicense, or distribute the Program +except as expressly provided under this License. Any attempt +otherwise to copy, modify, sublicense or distribute the Program is +void, and will automatically terminate your rights under this License. +However, parties who have received copies, or rights, from you under +this License will not have their licenses terminated so long as such +parties remain in full compliance. + + 5. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Program or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Program (or any work based on the +Program), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Program or works based on it. + + 6. Each time you redistribute the Program (or any work based on the +Program), the recipient automatically receives a license from the +original licensor to copy, distribute or modify the Program subject to +these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties to +this License. + + 7. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Program at all. For example, if a patent +license would not permit royalty-free redistribution of the Program by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Program. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system, which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 8. If the distribution and/or use of the Program is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Program under this License +may add an explicit geographical distribution limitation excluding +those countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 9. The Free Software Foundation may publish revised and/or new versions +of the General Public License from time to time. Such new versions will +be similar in spirit to the present version, but may differ in detail to +address new problems or concerns. + +Each version is given a distinguishing version number. If the Program +specifies a version number of this License which applies to it and "any +later version", you have the option of following the terms and conditions +either of that version or of any later version published by the Free +Software Foundation. If the Program does not specify a version number of +this License, you may choose any version ever published by the Free Software +Foundation. + + 10. If you wish to incorporate parts of the Program into other free +programs whose distribution conditions are different, write to the author +to ask for permission. For software which is copyrighted by the Free +Software Foundation, write to the Free Software Foundation; we sometimes +make exceptions for this. Our decision will be guided by the two goals +of preserving the free status of all derivatives of our free software and +of promoting the sharing and reuse of software generally. + + NO WARRANTY + + 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY +FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN +OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES +PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY OF ANY KIND, EITHER EXPRESSED +OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF +MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. THE ENTIRE RISK AS +TO THE QUALITY AND PERFORMANCE OF THE PROGRAM IS WITH YOU. SHOULD THE +PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING, +REPAIR OR CORRECTION. + + 12. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING +WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY AND/OR +REDISTRIBUTE THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, +INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING +OUT OF THE USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED +TO LOSS OF DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY +YOU OR THIRD PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER +PROGRAMS), EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE +POSSIBILITY OF SUCH DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Programs + + If you develop a new program, and you want it to be of the greatest +possible use to the public, the best way to achieve this is to make it +free software which everyone can redistribute and change under these terms. + + To do so, attach the following notices to the program. It is safest +to attach them to the start of each source file to most effectively +convey the exclusion of warranty; and each file should have at least +the "copyright" line and a pointer to where the full notice is found. + + + Copyright (C) + + This program is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +Also add information on how to contact you by electronic and paper mail. + +If the program is interactive, make it output a short notice like this +when it starts in an interactive mode: + + Gnomovision version 69, Copyright (C) year name of author + Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'. + This is free software, and you are welcome to redistribute it + under certain conditions; type `show c' for details. + +The hypothetical commands `show w' and `show c' should show the appropriate +parts of the General Public License. Of course, the commands you use may +be called something other than `show w' and `show c'; they could even be +mouse-clicks or menu items--whatever suits your program. + +You should also get your employer (if you work as a programmer) or your +school, if any, to sign a "copyright disclaimer" for the program, if +necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the program + `Gnomovision' (which makes passes at compilers) written by James Hacker. + + , 1 April 1989 + Ty Coon, President of Vice + +This General Public License does not permit incorporating your program into +proprietary programs. If your program is a subroutine library, you may +consider it more useful to permit linking proprietary applications with the +library. If this is what you want to do, use the GNU Lesser General +Public License instead of this License. From c92245535886a9f7fcb81f6bf955d44affe36164 Mon Sep 17 00:00:00 2001 From: jngaravitoc Date: Mon, 13 Nov 2023 12:04:10 -0500 Subject: [PATCH 18/19] docs --- docs/Makefile | 20 ++++++++++++++++++++ docs/make.bat | 35 +++++++++++++++++++++++++++++++++++ docs/source/conf.py | 28 ++++++++++++++++++++++++++++ docs/source/index.rst | 20 ++++++++++++++++++++ 4 files changed, 103 insertions(+) create mode 100644 docs/Makefile create mode 100644 docs/make.bat create mode 100644 docs/source/conf.py create mode 100644 docs/source/index.rst diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..d0c3cbf --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = source +BUILDDIR = build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 index 0000000..747ffb7 --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=source +set BUILDDIR=build + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.https://www.sphinx-doc.org/ + exit /b 1 +) + +if "%1" == "" goto help + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/docs/source/conf.py b/docs/source/conf.py new file mode 100644 index 0000000..a9c81a6 --- /dev/null +++ b/docs/source/conf.py @@ -0,0 +1,28 @@ +# Configuration file for the Sphinx documentation builder. +# +# For the full list of built-in configuration values, see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Project information ----------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information + +project = 'EXPtools' +copyright = '2023, Nicolas Garavito-Camargo' +author = 'Nicolas Garavito-Camargo' +release = '0.1' + +# -- General configuration --------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration + +extensions = [] + +templates_path = ['_templates'] +exclude_patterns = [] + + + +# -- Options for HTML output ------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output + +html_theme = 'alabaster' +html_static_path = ['_static'] diff --git a/docs/source/index.rst b/docs/source/index.rst new file mode 100644 index 0000000..e3e81bf --- /dev/null +++ b/docs/source/index.rst @@ -0,0 +1,20 @@ +.. EXPtools documentation master file, created by + sphinx-quickstart on Sun Nov 12 12:37:16 2023. + You can adapt this file completely to your liking, but it should at least + contain the root `toctree` directive. + +Welcome to EXPtools's documentation! +==================================== + +.. toctree:: + :maxdepth: 2 + :caption: Contents: + + + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` From 5f363d8a96027fba353dff42cddb90b76ce20104 Mon Sep 17 00:00:00 2001 From: NGC Date: Wed, 15 Nov 2023 22:58:11 -0500 Subject: [PATCH 19/19] Create EXP build.yml --- .github/workflows/build.yml | 94 +++++++++++++++++++++++++++++++++++++ 1 file changed, 94 insertions(+) create mode 100644 .github/workflows/build.yml diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml new file mode 100644 index 0000000..92dd368 --- /dev/null +++ b/.github/workflows/build.yml @@ -0,0 +1,94 @@ +name: "Test Ubuntu Build" + +on: + push: + branches: + - main + pull_request: + +jobs: + pyexp: + strategy: + matrix: + os: [ubuntu-latest] + cc: [gcc, mpicc] + + name: "Test pyEXP Build" + runs-on: ${{ matrix.os }} + steps: + - name: Checkout + uses: actions/checkout@v3 + + - name: Install core dependencies - ubuntu + if: runner.os == 'Linux' + run: | + sudo apt-get update + sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev + + - name: Setup submodule and build + run: | + git submodule update --init --recursive + mkdir -p build/install + + - name: Compile pyEXP - Linux + if: runner.os == 'Linux' + env: + CC: ${{ matrix.cc }} + working-directory: ./build + run: >- + cmake + -DENABLE_NBODY=NO + -DENABLE_PYEXP=YES + -DCMAKE_BUILD_TYPE=Release + -DEigen3_DIR=/usr/include/eigen3/share/eigen3/cmake + -DCMAKE_INSTALL_PREFIX=./install + -Wno-dev + .. + + - name: Make + working-directory: ./build + run: make -j 2 + + # ----------------------------------------------------------------------------------- + + exp: + strategy: + matrix: + os: [ubuntu-latest] + cc: [gcc, mpicc] + + name: "Test Full EXP Build" + runs-on: ${{ matrix.os }} + steps: + - name: Checkout + uses: actions/checkout@v3 + + - name: Install core dependencies - ubuntu + if: runner.os == 'Linux' + run: | + sudo apt-get update + sudo apt-get install -y build-essential libeigen3-dev libfftw3-dev libhdf5-dev libopenmpi-dev + + - name: Setup submodule and build + run: | + git submodule update --init --recursive + mkdir -p build/install + + - name: Compile Full EXP - Linux + if: runner.os == 'Linux' + env: + CC: ${{ matrix.cc }} + working-directory: ./build + run: >- + cmake + -DENABLE_NBODY=YES + -DENABLE_PYEXP=NO + -DCMAKE_BUILD_TYPE=Release + -DEigen3_DIR=/usr/include/eigen3/share/eigen3/cmake + -DCMAKE_INSTALL_PREFIX=./install + -Wno-dev + .. + + - name: Make + working-directory: ./build + run: make -j 2