From 9d9413bd575975c4ceb8b53125dd81edcb5ee956 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Wed, 26 Jul 2023 16:30:41 +0100 Subject: [PATCH 01/36] Draft HDF5 converter --- exptool/io/psp_to_hdf5.py | 124 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 124 insertions(+) create mode 100644 exptool/io/psp_to_hdf5.py diff --git a/exptool/io/psp_to_hdf5.py b/exptool/io/psp_to_hdf5.py new file mode 100644 index 0000000..3d4987e --- /dev/null +++ b/exptool/io/psp_to_hdf5.py @@ -0,0 +1,124 @@ + +import numpy as np +import argparse + +try: + import h5py +except ImportError: + raise ImportError("You will need to 'pip install h5py --user' to use this converter.") + +from . import particle + +class HDFConverter(): + def __init__(self, filename,comp=None,verbose=0): + """the main driver""" + self.filename = filename + + convert_psp_to_hdf5(self.filename) + # if a component is selected, perhaps we should only convert that? + + @staticmethod + def convert_psp_to_hdf5(inputfilename): + """ + """ + # define file + outputfilename = inputfilename + '.h5' + + # steps: open file, get all components + O = particle.Input(inputfilename) + comps = list(O.header.keys()) + + # open the new file for printing + f = h5py.File(outputfilename, 'w') + + # not using a global header, so add some attributes here + f['time'] = O.time + + for comp in comps: + + f.create_group(comp) + + HDFConverter.print_component_header(f,O,comp) + + HDFConverter.make_phasespace(f,inputfilename,comp) + + f.close() + + + @staticmethod + def print_component_header(f,O,comp): + f[comp].create_group('header') + for key in O.header[comp].keys(): + try: # see if there is another dictionary level + for subkey in O.header[comp][key].keys(): + try: # see if there is another dictionary level + for subsubkey in O.header[comp][key][subkey].keys(): # this is the deepest we ever have to go + try: + f['{}/header/{}/{}'.format(comp,key,subkey)].attrs.create(subsubkey,O.header[comp][key][subkey][subsubkey]) + #print('{}/header/{}/{}'.format(comp,key,subkey),subsubkey,O.header[comp][key][subkey][subsubkey]) + except: + f['{}/header/{}'.format(comp,key)].create_group(subkey) + f['{}/header/{}/{}'.format(comp,key,subkey)].attrs.create(subsubkey,O.header[comp][key][subkey][subsubkey]) + #print('CREATE ','{}/header/{}/{}'.format(comp,key,subkey),subsubkey,O.header[comp][key][subkey][subsubkey]) + except: + try: + f['{}/header/{}'.format(comp,key)].attrs.create(subkey,O.header[comp][key][subkey]) + #print('{}/header/{}'.format(comp,key),subkey,O.header[comp][key][subkey]) + except: + f['{}/header'.format(comp)].create_group(key) + f['{}/header/{}'.format(comp,key)].attrs.create(subkey,O.header[comp][key][subkey]) + #print('CREATE ','{}/header/{}'.format(comp,key),subkey,O.header[comp][key][subkey]) + except: + f['{}/header'.format(comp)].attrs.create(key,O.header[comp][key]) + #print('{}/header'.format(comp),key,O.header[comp][key]) + + @staticmethod + def make_phasespace(f,inputfilename,comp): + # get the data + O1 = particle.Input(inputfilename,comp) + # make the phase space + PS = np.array([O1.data['m'],O1.data['x'],O1.data['y'],O1.data['z'],O1.data['vx'],O1.data['vy'],O1.data['vz'],O1.data['potE']]).T + # print the phase space + dset = f[comp].create_dataset('phasespace', data=PS) + + +""" +if __name__ == '__main__': + parser = argparse.ArgumentParser( + prog='PSP2HDF5', + description='Convert OUTPSN files to HDF5 format.')#epilog='Text at the bottom of help') + + parser.add_argument('filename',help='the file to be converted') + + args = parser.parse_args() + + HDFConverter.convert_psp_to_hdf5(args.filename) + + + +# example usage +# to run, do PSP2HDF5 OUT.run0.00000 + +import h5py +outputfilename = 'OUT.run0.00000.h5' +outputfilename = 'OUT.run0.00000.h5' +f = h5py.File(outputfilename, 'r') + +# if the component is called 'halo', then +f['halo/phasespace'] # is the dataset, (Nx8), with (mass,x,y,z,vx,vy,vz,potential) for each particle + +# the header data is saved in a group, +f['halo/header'] +# and then each group may have subgroups or attributes. typically, you will have +f['halo/header/parameters'] +print(f['halo/header/parameters'].attrs.keys()) +print(f['halo/header/parameters'].attrs['nEJwant']) + +f['halo/header/force'] +print(f['halo/header/force'].attrs['id']) +print(f['halo/header/force/parameters'].attrs.keys()) +for key in f['halo/header/force/parameters'].attrs.keys(): + print(key,f['halo/header/force/parameters'].attrs[key]) + + +""" From 544d5fbbad14174014347bbf11bdd7dd5871699f Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Sun, 24 Sep 2023 20:36:54 +0100 Subject: [PATCH 02/36] Improve documentation in the HDF5 format --- exptool/io/psp_to_hdf5.py | 213 ++++++++++++++++++++++++-------------- 1 file changed, 135 insertions(+), 78 deletions(-) diff --git a/exptool/io/psp_to_hdf5.py b/exptool/io/psp_to_hdf5.py index 3d4987e..d05c029 100644 --- a/exptool/io/psp_to_hdf5.py +++ b/exptool/io/psp_to_hdf5.py @@ -1,124 +1,181 @@ +""" +draft conversion from PSP format the HDF5. -import numpy as np + +For each component group, there is a subgroup named 'header', which stores header information related to that component. This information may include various parameters and metadata. The header information is organized into nested groups and attributes within the 'header' subgroup. The structure of the header data may vary depending on the specific PSP file format. + +The main data associated with each component is stored in a dataset named 'phasespace' within the component group. This dataset is an Nx8 array, where N represents the number of particles. Each row of the dataset corresponds to a particle and contains the following information in this order: mass, x-coordinate, y-coordinate, z-coordinate, x-velocity, y-velocity, z-velocity, and potential energy. + + +Example usage: import argparse -try: - import h5py -except ImportError: - raise ImportError("You will need to 'pip install h5py --user' to use this converter.") +if __name__ == '__main__': + parser = argparse.ArgumentParser( + prog='PSP2HDF5', + description='Convert OUTPSN files to HDF5 format.')#epilog='Text at the bottom of help') + + parser.add_argument('filename',help='the file to be converted') + + args = parser.parse_args() + + HDFConverter.convert_psp_to_hdf5(args.filename) + + + +# example usage +# to run, do PSP2HDF5 OUT.run0.00000 + +import h5py +outputfilename = 'OUT.run0.00000.h5' +outputfilename = 'OUT.run0.00000.h5' + +# how is the global header information saved? +# only time is saved +f = h5py.File(outputfilename, 'r') +# if the component is called 'halo', then +f['halo/phasespace'] # is the dataset, (Nx8), with (mass,x,y,z,vx,vy,vz,potential) for each particle + +# the header data is saved in a group, +f['halo/header'] +# and then each group may have subgroups or attributes. typically, you will have +f['halo/header/parameters'] +print(f['halo/header/parameters'].attrs.keys()) +print(f['halo/header/parameters'].attrs['nEJwant']) + +f['halo/header/force'] +print(f['halo/header/force'].attrs['id']) +print(f['halo/header/force/parameters'].attrs.keys()) +for key in f['halo/header/force/parameters'].attrs.keys(): + print(key,f['halo/header/force/parameters'].attrs[key]) + + +""" + + + +import numpy as np +import h5py from . import particle class HDFConverter(): - def __init__(self, filename,comp=None,verbose=0): - """the main driver""" + """ + HDFConverter class for converting custom PSP (Particle Simulation Program) files to HDF5 format. + + This class allows you to convert data from PSP format into HDF5 format, which is a versatile and efficient data storage format. + + Parameters: + filename (str): The name of the PSP input file to be converted. + comp (str): Optional. The specific component to convert. If provided, only the data for the specified component will be converted. + verbose (int): Optional. Verbosity level for printing progress and messages during conversion. + + Attributes: + filename (str): The name of the PSP input file. + """ + + def __init__(self, filename, comp=None, verbose=0): + """ + Initialize the HDFConverter instance. + + Args: + filename (str): The name of the PSP input file to be converted. + comp (str, optional): The specific component to convert. If provided, only the data for the specified component will be converted. + verbose (int, optional): Verbosity level for printing progress and messages during conversion. + """ self.filename = filename + self.comp = comp + self.verbose = verbose - convert_psp_to_hdf5(self.filename) - # if a component is selected, perhaps we should only convert that? + # Start the conversion process + self.convert_psp_to_hdf5() @staticmethod def convert_psp_to_hdf5(inputfilename): """ + Convert a PSP input file to HDF5 format. + + Args: + inputfilename (str): The name of the PSP input file to be converted. """ - # define file + # Define the output file name outputfilename = inputfilename + '.h5' - # steps: open file, get all components + # Open the PSP input file and extract components O = particle.Input(inputfilename) comps = list(O.header.keys()) - # open the new file for printing + # Create a new HDF5 file for storing the converted data f = h5py.File(outputfilename, 'w') - # not using a global header, so add some attributes here + # Store the simulation time as an attribute f['time'] = O.time for comp in comps: - + # Create a group for each component f.create_group(comp) - HDFConverter.print_component_header(f,O,comp) + # Print the header information for the component + HDFConverter.print_component_header(f, O, comp) - HDFConverter.make_phasespace(f,inputfilename,comp) + # Create and store the phase space data for the component + HDFConverter.make_phasespace(f, inputfilename, comp) + # Close the HDF5 file f.close() - @staticmethod - def print_component_header(f,O,comp): + def print_component_header(f, O, comp): + """ + Print header information for a component to an HDF5 file. + + Args: + f (h5py.Group): The HDF5 group to store the header information. + O (particle.Input): The PSP input object. + comp (str): The name of the component. + """ f[comp].create_group('header') for key in O.header[comp].keys(): - try: # see if there is another dictionary level + # Check for nested dictionary levels + try: for subkey in O.header[comp][key].keys(): - try: # see if there is another dictionary level - for subsubkey in O.header[comp][key][subkey].keys(): # this is the deepest we ever have to go + # Check for further nested levels + try: + for subsubkey in O.header[comp][key][subkey].keys(): + # Create attributes for the deepest level try: - f['{}/header/{}/{}'.format(comp,key,subkey)].attrs.create(subsubkey,O.header[comp][key][subkey][subsubkey]) - #print('{}/header/{}/{}'.format(comp,key,subkey),subsubkey,O.header[comp][key][subkey][subsubkey]) + f['{}/header/{}/{}'.format(comp, key, subkey)].attrs.create(subsubkey, O.header[comp][key][subkey][subsubkey]) except: - f['{}/header/{}'.format(comp,key)].create_group(subkey) - f['{}/header/{}/{}'.format(comp,key,subkey)].attrs.create(subsubkey,O.header[comp][key][subkey][subsubkey]) - #print('CREATE ','{}/header/{}/{}'.format(comp,key,subkey),subsubkey,O.header[comp][key][subkey][subsubkey]) + # Create subgroups if necessary + f['{}/header/{}/{}'.format(comp, key, subkey)].create_group(subsubkey) + f['{}/header/{}/{}'.format(comp, key, subkey)].attrs.create(subsubkey, O.header[comp][key][subkey][subsubkey]) except: + # Create attributes for the intermediate level try: - f['{}/header/{}'.format(comp,key)].attrs.create(subkey,O.header[comp][key][subkey]) - #print('{}/header/{}'.format(comp,key),subkey,O.header[comp][key][subkey]) + f['{}/header/{}/{}'.format(comp, key, subkey)].attrs.create(subkey, O.header[comp][key][subkey]) except: + # Create subgroups if necessary f['{}/header'.format(comp)].create_group(key) - f['{}/header/{}'.format(comp,key)].attrs.create(subkey,O.header[comp][key][subkey]) - #print('CREATE ','{}/header/{}'.format(comp,key),subkey,O.header[comp][key][subkey]) + f['{}/header/{}'.format(comp, key)].attrs.create(subkey, O.header[comp][key][subkey]) except: - f['{}/header'.format(comp)].attrs.create(key,O.header[comp][key]) - #print('{}/header'.format(comp),key,O.header[comp][key]) + # Create attributes for the top-level header + f['{}/header'.format(comp)].attrs.create(key, O.header[comp][key]) @staticmethod - def make_phasespace(f,inputfilename,comp): - # get the data - O1 = particle.Input(inputfilename,comp) - # make the phase space - PS = np.array([O1.data['m'],O1.data['x'],O1.data['y'],O1.data['z'],O1.data['vx'],O1.data['vy'],O1.data['vz'],O1.data['potE']]).T - # print the phase space - dset = f[comp].create_dataset('phasespace', data=PS) - - -""" -if __name__ == '__main__': - parser = argparse.ArgumentParser( - prog='PSP2HDF5', - description='Convert OUTPSN files to HDF5 format.')#epilog='Text at the bottom of help') - - parser.add_argument('filename',help='the file to be converted') - - args = parser.parse_args() - - HDFConverter.convert_psp_to_hdf5(args.filename) - - - -# example usage -# to run, do PSP2HDF5 OUT.run0.00000 - -import h5py -outputfilename = 'OUT.run0.00000.h5' -outputfilename = 'OUT.run0.00000.h5' -f = h5py.File(outputfilename, 'r') - -# if the component is called 'halo', then -f['halo/phasespace'] # is the dataset, (Nx8), with (mass,x,y,z,vx,vy,vz,potential) for each particle - -# the header data is saved in a group, -f['halo/header'] -# and then each group may have subgroups or attributes. typically, you will have -f['halo/header/parameters'] -print(f['halo/header/parameters'].attrs.keys()) -print(f['halo/header/parameters'].attrs['nEJwant']) + def make_phasespace(f, inputfilename, comp): + """ + Convert and store phase space data for a component in an HDF5 file. -f['halo/header/force'] -print(f['halo/header/force'].attrs['id']) -print(f['halo/header/force/parameters'].attrs.keys()) -for key in f['halo/header/force/parameters'].attrs.keys(): - print(key,f['halo/header/force/parameters'].attrs[key]) + Args: + f (h5py.Group): The HDF5 group to store the phase space data. + inputfilename (str): The name of the PSP input file. + comp (str): The name of the component. + """ + # Read data from the PSP input file + O1 = particle.Input(inputfilename, comp) + # Create a phase space array + PS = np.array([O1.data['m'], O1.data['x'], O1.data['y'], O1.data['z'], O1.data['vx'], O1.data['vy'], O1.data['vz'], O1.data['potE']]).T -""" + # Store the phase space data as a dataset + dset = f[comp].create_dataset('phasespace', data=PS) From 58a2e0a3b67297c8056057a6853574a64313f5d8 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Sat, 18 May 2024 11:42:11 +0100 Subject: [PATCH 03/36] draft hdf5 implementation for aps arrays --- exptool/analysis/trapping.py | 68 ++++++++++++++++++++++++------------ 1 file changed, 45 insertions(+), 23 deletions(-) diff --git a/exptool/analysis/trapping.py b/exptool/analysis/trapping.py index 158f9dd..8061eb7 100644 --- a/exptool/analysis/trapping.py +++ b/exptool/analysis/trapping.py @@ -50,6 +50,9 @@ import os from scipy import interpolate +# io import +import h5py + # multiprocessing imports import itertools from multiprocessing import Pool, freeze_support @@ -153,10 +156,14 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= else: # limit to the maximum number desired particle_indx = np.arange(0,particle_indx,1) - else: - # assume an array has been passed and accept: could check + + # assume an array has been passed + elif isinstance(particle_indx,np.ndarray): pass + else: + raise ValueError("exptool.ApsFinding.trapping._determin_r_aps: particle_indx must be an integer or an array.") + # sort the particle indices particle_indx = particle_indx[particle_indx.argsort()] @@ -166,17 +173,18 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= # # stamps the output file with the current time. do we like this? # - tstamp = datetime.datetime.strftime(datetime.datetime.now(), '%Y-%m-%d+%H:%M:%S') + #tstamp = datetime.datetime.strftime(datetime.datetime.now(), '%Y-%m-%d+%H:%M:%S') + tstamp = datetime.datetime.strftime(datetime.datetime.now(), '%Y-%m-%d+%H') # create a new file using the particle number, runtag, and time outputfile = out_directory+'RadialAps_N{}_r{}_T{}.dat'.format(total_orbits,runtag,tstamp) - f = open(outputfile,'wb+') + f = h5py.File(outputfile,"w") - # - # print descriptor string - # + # createdescriptor string desc = 'apsfile for '+comp+' in '+out_directory+', norbits='+str(total_orbits)+', threedee='+str(threedee)+', using '+filelist - np.array([desc],dtype='S200').tofile(f) + + # Write the descriptor string as an attribute + f.attrs['description'] = desc aps_dictionary = dict() # make blank dictionary for the aps for i in range(0,total_orbits): aps_dictionary[i] = [] @@ -220,9 +228,9 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= R3 = np.sqrt(X3*X3 + Y3*Y3 + Z3*Z3) else: - R1 = np.sqrt(X1*X1 + Y1*Y1) - R2 = np.sqrt(X2*X2 + Y2*Y2) - R3 = np.sqrt(X3*X3 + Y3*Y3) + R1 = np.linalg.norm([X1,Y1],axis=0) + R2 = np.linalg.norm([X2,Y2],axis=0) + R3 = np.linalg.norm([X3,Y3],axis=0) else: # i!=1 @@ -254,7 +262,7 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= if threedee: R3 = np.sqrt(X3*X3 + Y3*Y3 + Z3*Z3) else: - R3 = np.sqrt(X3*X3 + Y3*Y3) + R3 = np.linalg.norm([X3,Y3],axis=0) # R1 might be the shortest, if particles are added, so only compare up to the length of r1 @@ -276,7 +284,7 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= #orderid = IN[:r1length][aps] if self.verbose > 0: - print('Current time: {4.3f}'.format(tval),end='\r', flush=True) + print('exptool.ApsFinding.trapping._determin_r_aps: Current time: {4.3f}'.format(tval),end='\r', flush=True) # under this convention, the user needs to keep track of the particle index that was input #for j in range(0,len(index_tags)): @@ -286,14 +294,15 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= for j in range(0,len(index_tags)): aps_dictionary[id[j]].append([tval,x[j],y[j],z[j]]) - + # create a tracker for the number of aps per orbit self.napsides = np.zeros([total_orbits,2]) # print a header with the number of orbits - np.array([total_orbits],dtype='i').tofile(f) + f.attrs['total_orbits'] = total_orbits orbits_with_apocentre = 0 + # go back through all the orbits and write to file for j in range(0,total_orbits): orbit_aps_array = np.array(aps_dictionary[particle_indx[j]]) @@ -307,28 +316,41 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= orbits_with_apocentre += 1 naps = len(orbit_aps_array[:,0]) # this might be better as shape - np.array([naps],dtype='i').tofile(f) + #np.array([naps],dtype='i').tofile(f) + + #self.napsides[j,0] = naps + #self.napsides[j,1] = len(orbit_aps_array.reshape(-1,)) + + #np.array( orbit_aps_array.reshape(-1,),dtype='f').tofile(f) - self.napsides[j,0] = naps - self.napsides[j,1] = len(orbit_aps_array.reshape(-1,)) + # create a dataset with the index of the particle as the tag + dataset = f.create_dataset(str(particle_indx[j]), data=orbit_aps_array) + + # create attributes for the dataset + dataset.attrs['naps'] = naps - np.array( orbit_aps_array.reshape(-1,),dtype='f').tofile(f) # no valid turning points: put in a blank else: + # create a dataset with the index of the particle as the tag + dataset = f.create_dataset(str(particle_indx[j]), data=np.array([-1.])) + + # create attributes for the dataset + dataset.attrs['naps'] = 0 + # guard against zero length - np.array([1],dtype='i').tofile(f) + #np.array([1],dtype='i').tofile(f) # indices start at 1 - np.array( np.array(([-1.,-1.,-1.,-1.])).reshape(-1,),dtype='f').tofile(f) + #np.array( np.array(([-1.,-1.,-1.,-1.])).reshape(-1,),dtype='f').tofile(f) f.close() - print('trapping.ApsFinding.determine_r_aps: found {} orbits (out of {}) with valid apocentres.'.format(orbits_with_apocentre,total_orbits)) + print('exptool.trapping.ApsFinding.determine_r_aps: found {} orbits (out of {}) with valid apocentres.'.format(orbits_with_apocentre,total_orbits)) - print('trapping.ApsFinding.determine_r_aps: savefile is {}'.format(outputfile)) + print('exptool.trapping.ApsFinding.determine_r_aps: savefile is {}'.format(outputfile)) if (return_aps): ApsDict = ApsFinding.read_aps_file(self,outputfile) From ce1d4151d06d3ccc56d3de0b4678d90671db1aa3 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Sat, 18 May 2024 11:48:07 +0100 Subject: [PATCH 04/36] bug in changed names --- exptool/analysis/trapping.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exptool/analysis/trapping.py b/exptool/analysis/trapping.py index 8061eb7..b683109 100644 --- a/exptool/analysis/trapping.py +++ b/exptool/analysis/trapping.py @@ -291,7 +291,7 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= # aps_dictionary[orderid[j]].append([tval,x[j],y[j],z[j]]) # under this convention, the id of the orbit is preserved and used as the dictionary key - for j in range(0,len(index_tags)): + for j in range(0,len(id)): aps_dictionary[id[j]].append([tval,x[j],y[j],z[j]]) # create a tracker for the number of aps per orbit From 669b39f5f1afacc4b788933468e7f40fffec46fc Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Sat, 18 May 2024 11:58:54 +0100 Subject: [PATCH 05/36] handle particle_indx properly --- exptool/analysis/trapping.py | 29 ++++++++++++++++------------- 1 file changed, 16 insertions(+), 13 deletions(-) diff --git a/exptool/analysis/trapping.py b/exptool/analysis/trapping.py index b683109..f9499b3 100644 --- a/exptool/analysis/trapping.py +++ b/exptool/analysis/trapping.py @@ -151,8 +151,8 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= if isinstance(particle_indx,int): if particle_indx < 0: # if particle_indx < 0, make the comparison index all particles - Oa = particle.Input(self.SLIST[0],legacy=False,comp=comp,verbose=0) - particle_indx = np.arange(0,Oa.nbodies,1) + Oa = particle.Input(self.SLIST[0],comp=comp,verbose=0) + particle_indx = np.arange(0,Oa.data['id'].size,1) else: # limit to the maximum number desired particle_indx = np.arange(0,particle_indx,1) @@ -186,10 +186,13 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= # Write the descriptor string as an attribute f.attrs['description'] = desc - aps_dictionary = dict() # make blank dictionary for the aps - for i in range(0,total_orbits): aps_dictionary[i] = [] + # make blank dictionary for the aps + aps_dictionary = dict() + # make a blank array for each orbit + for i in particle_indx: aps_dictionary[i] = [] + # loop through files for i in range(1,len(self.SLIST)-1): if i==1: @@ -217,16 +220,16 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= X3 = Oc.data['x'][p3indx];Y3 = Oc.data['y'][p3indx];Z3 = Oc.data['z'][p3indx];I3 = Oc.data['id'][p3indx] else: - X1 = Oa.data['x'];Y1 = Oa.data['y'];Z1 = Oa.data['z'];I1 = Oa.data['id'] - X2 = Ob.data['x'];Y2 = Ob.data['y'];Z2 = Ob.data['z'];I2 = Ob.data['id'] - X3 = Oc.data['x'];Y3 = Oc.data['y'];Z3 = Oc.data['z'];I3 = Oc.data['id'] + X1 = Oa.data['x'][particle_indx];Y1 = Oa.data['y'][particle_indx];Z1 = Oa.data['z'][particle_indx];I1 = Oa.data['id'][particle_indx] + X2 = Ob.data['x'][particle_indx];Y2 = Ob.data['y'][particle_indx];Z2 = Ob.data['z'][particle_indx];I2 = Ob.data['id'][particle_indx] + X3 = Oc.data['x'][particle_indx];Y3 = Oc.data['y'][particle_indx];Z3 = Oc.data['z'][particle_indx];I3 = Oc.data['id'][particle_indx] + # compute radial positions if threedee: - R1 = np.sqrt(X1*X1 + Y1*Y1 + Z1*Z1) - R2 = np.sqrt(X2*X2 + Y2*Y2 + Z2*Z2) - R3 = np.sqrt(X3*X3 + Y3*Y3 + Z3*Z3) - + R1 = np.linalg.norm([X1,Y1,Z1],axis=0) + R2 = np.linalg.norm([X2,Y2,Z2],axis=0) + R3 = np.linalg.norm([X3,Y3,Z3],axis=0) else: R1 = np.linalg.norm([X1,Y1],axis=0) R2 = np.linalg.norm([X2,Y2],axis=0) @@ -257,10 +260,10 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= X3 = Oc.data['x'][p3indx];Y3 = Oc.data['y'][p3indx];Z3 = Oc.data['z'][p3indx];I3 = Oc.data['id'][p3indx] else: - X3 = Oc.data['x'];Y3 = Oc.data['y'];Z3 = Oc.data['z'];I3 = Oc.data['id'] + X3 = Oc.data['x'][particle_indx];Y3 = Oc.data['y'][particle_indx];Z3 = Oc.data['z'][particle_indx];I3 = Oc.data['id'][particle_indx] if threedee: - R3 = np.sqrt(X3*X3 + Y3*Y3 + Z3*Z3) + R3 = np.linalg.norm([X3,Y3,Z3],axis=0) else: R3 = np.linalg.norm([X3,Y3],axis=0) From d73c8e7a7a9d2f160dc8ad4840316fe1d7d93e1a Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Sat, 18 May 2024 12:11:30 +0100 Subject: [PATCH 06/36] revert particle index tracking --- exptool/analysis/trapping.py | 34 +++++++++++----------------------- 1 file changed, 11 insertions(+), 23 deletions(-) diff --git a/exptool/analysis/trapping.py b/exptool/analysis/trapping.py index f9499b3..f8da306 100644 --- a/exptool/analysis/trapping.py +++ b/exptool/analysis/trapping.py @@ -112,7 +112,7 @@ def accept_files(self,filelist,verbose=0): ApsFinding.parse_list(self) - def parse_list(self): + def _parse_list(self): """ parse files from the input list @@ -131,7 +131,7 @@ def parse_list(self): self.SLIST = np.array(s_list) if self.verbose >= 1: - print('ApsFinding.parse_list: Accepted {0:d} files.'.format(len(self.SLIST))) + print('exptool.trapping.ApsFinding.parse_list: Accepted {0:d} files.'.format(len(self.SLIST))) def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory='',threedee=False,return_aps=False,changingindx=False): @@ -144,7 +144,7 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= # take the inputs and identify all files that we will loop through self.slist = filelist - ApsFinding.parse_list(self) + ApsFinding._parse_list(self) # now we have self.SLIST, the parsed list of files we will analyse # first, check type of particle_index @@ -159,7 +159,7 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= # assume an array has been passed elif isinstance(particle_indx,np.ndarray): - pass + changingindx = True else: raise ValueError("exptool.ApsFinding.trapping._determin_r_aps: particle_indx must be an integer or an array.") @@ -220,9 +220,9 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= X3 = Oc.data['x'][p3indx];Y3 = Oc.data['y'][p3indx];Z3 = Oc.data['z'][p3indx];I3 = Oc.data['id'][p3indx] else: - X1 = Oa.data['x'][particle_indx];Y1 = Oa.data['y'][particle_indx];Z1 = Oa.data['z'][particle_indx];I1 = Oa.data['id'][particle_indx] - X2 = Ob.data['x'][particle_indx];Y2 = Ob.data['y'][particle_indx];Z2 = Ob.data['z'][particle_indx];I2 = Ob.data['id'][particle_indx] - X3 = Oc.data['x'][particle_indx];Y3 = Oc.data['y'][particle_indx];Z3 = Oc.data['z'][particle_indx];I3 = Oc.data['id'][particle_indx] + X1 = Oa.data['x'];Y1 = Oa.data['y'];Z1 = Oa.data['z'];I1 = Oa.data['id'] + X2 = Ob.data['x'];Y2 = Ob.data['y'];Z2 = Ob.data['z'];I2 = Ob.data['id'] + X3 = Oc.data['x'];Y3 = Oc.data['y'];Z3 = Oc.data['z'];I3 = Oc.data['id'] # compute radial positions @@ -260,7 +260,7 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= X3 = Oc.data['x'][p3indx];Y3 = Oc.data['y'][p3indx];Z3 = Oc.data['z'][p3indx];I3 = Oc.data['id'][p3indx] else: - X3 = Oc.data['x'][particle_indx];Y3 = Oc.data['y'][particle_indx];Z3 = Oc.data['z'][particle_indx];I3 = Oc.data['id'][particle_indx] + X3 = Oc.data['x'];Y3 = Oc.data['y'];Z3 = Oc.data['z'];I3 = Oc.data['id'] if threedee: R3 = np.linalg.norm([X3,Y3,Z3],axis=0) @@ -289,11 +289,7 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= if self.verbose > 0: print('exptool.ApsFinding.trapping._determin_r_aps: Current time: {4.3f}'.format(tval),end='\r', flush=True) - # under this convention, the user needs to keep track of the particle index that was input - #for j in range(0,len(index_tags)): - # aps_dictionary[orderid[j]].append([tval,x[j],y[j],z[j]]) - - # under this convention, the id of the orbit is preserved and used as the dictionary key + # the id of the orbit is preserved and used as the dictionary key for j in range(0,len(id)): aps_dictionary[id[j]].append([tval,x[j],y[j],z[j]]) @@ -310,21 +306,13 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= orbit_aps_array = np.array(aps_dictionary[particle_indx[j]]) - # print the index to the file - np.array([particle_indx[j]],dtype='i').tofile(f) - # if there are valid turning points: if (len(orbit_aps_array) > 0): orbits_with_apocentre += 1 - naps = len(orbit_aps_array[:,0]) # this might be better as shape - - #np.array([naps],dtype='i').tofile(f) - - #self.napsides[j,0] = naps - #self.napsides[j,1] = len(orbit_aps_array.reshape(-1,)) - #np.array( orbit_aps_array.reshape(-1,),dtype='f').tofile(f) + # count the number of turning points + naps = len(orbit_aps_array[:,0]) # create a dataset with the index of the particle as the tag dataset = f.create_dataset(str(particle_indx[j]), data=orbit_aps_array) From d4c6fbc3b5c4e0e0259ca7bbdac3538572374a4d Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Sat, 18 May 2024 12:19:54 +0100 Subject: [PATCH 07/36] restore finegrained time output --- exptool/analysis/trapping.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/exptool/analysis/trapping.py b/exptool/analysis/trapping.py index f8da306..946b580 100644 --- a/exptool/analysis/trapping.py +++ b/exptool/analysis/trapping.py @@ -173,8 +173,8 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= # # stamps the output file with the current time. do we like this? # - #tstamp = datetime.datetime.strftime(datetime.datetime.now(), '%Y-%m-%d+%H:%M:%S') - tstamp = datetime.datetime.strftime(datetime.datetime.now(), '%Y-%m-%d+%H') + tstamp = datetime.datetime.strftime(datetime.datetime.now(), '%Y-%m-%d+%H:%M:%S') + #tstamp = datetime.datetime.strftime(datetime.datetime.now(), '%Y-%m-%d+%H') # create a new file using the particle number, runtag, and time outputfile = out_directory+'RadialAps_N{}_r{}_T{}.dat'.format(total_orbits,runtag,tstamp) @@ -287,7 +287,7 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= #orderid = IN[:r1length][aps] if self.verbose > 0: - print('exptool.ApsFinding.trapping._determin_r_aps: Current time: {4.3f}'.format(tval),end='\r', flush=True) + print('exptool.ApsFinding.trapping._determin_r_aps: Current time: {0:4.3f}'.format(tval),end='\r', flush=True) # the id of the orbit is preserved and used as the dictionary key for j in range(0,len(id)): From f8a9d9ead5e4466b9cb96f6f9ab0781011d24d7a Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Sat, 18 May 2024 16:59:57 +0100 Subject: [PATCH 08/36] update trapping to allow for coefficient pattern speed calculation --- exptool/analysis/pattern.py | 264 +++++++++++++++++++++++++---------- exptool/analysis/trapping.py | 8 +- 2 files changed, 196 insertions(+), 76 deletions(-) diff --git a/exptool/analysis/pattern.py b/exptool/analysis/pattern.py index ebf0025..83a7a13 100644 --- a/exptool/analysis/pattern.py +++ b/exptool/analysis/pattern.py @@ -1,19 +1,11 @@ - -# 08-29-16: added maximum radius capabilities to bar_fourier_compute - -# 10-25-16: some redundancies noticed (bar_fourier_compute) and should be unified - ''' - pattern.py (part of exptool) tools to find patterns in the global simulation outputs - - - - +MSP 29 Aug 2016 Added maximum radius capabilities to bar_fourier_compute +MSP 25 Oct 2016 Some redundancies noticed (bar_fourier_compute) and should be unified BarTransform @@ -37,10 +29,6 @@ ''' -from __future__ import absolute_import, division, print_function, unicode_literals - - - # general imports import time import numpy as np @@ -51,14 +39,124 @@ # exptool imports -from exptool.io import particle -from exptool.utils import kmeans -from exptool.utils import utils - - - - +from ..io import particle +from ..utils import kmeans +from ..utils import utils + + +class BarFromCoefficients: + """ + Class to calculate and analyze the bar position from Fourier coefficients. + + Parameters + ---------- + times : array-like + Array of time values. + coefs : array-like + Array of complex Fourier coefficients. + unwrap_threshold : float, optional + Threshold for unwrapping the bar position phase (default is -1.). + smooth : bool, optional + Whether to smooth the unwrapped position (default is False). + reverse : bool, optional + Whether to reverse the unwrapping direction (default is False). + adjust : float, optional + Adjustment value used in unwrapping (default is np.pi). + verbose : int, optional + Verbosity level (default is 0). + smth_derivative : int, optional + Smoothing factor for the polynomial derivative calculation (default is 0). + spline_derivative : int, optional + Smoothing factor for the spline derivative calculation (default is 0). + """ + def __init__(self, times, coefs, unwrap_threshold=-1., smooth=False, reverse=False, adjust=np.pi, verbose=0, smth_derivative=0, spline_derivative=0): + self.verbose = verbose + self.time = times + self.cos = np.real(coefs) + self.sin = np.imag(coefs) + self.barposition = np.arctan2(self.sin, self.cos) + # compute the unwrapped position + self._unwrap_position(unwrap_threshold, smooth, reverse, adjust) + # compute the derivative + self._frequency_and_derivative(smth_derivative,spline_derivative) + + def _unwrap_position(self, unwrap_threshold, smooth, reverse, adjust): + """ + Unwrap the bar position to avoid discontinuities. + + Parameters + ---------- + unwrap_threshold : float + Threshold for phase unwrapping. + smooth : bool + Whether to smooth the unwrapped position. + reverse : bool + Whether to reverse the unwrapping direction. + adjust : float + Adjustment value used in unwrapping. + """ + running_number_of_rotations = 0 + number_of_rotations = np.zeros_like(self.barposition) + # start from zero + for i in range(1, len(self.barposition)): + if reverse: + if (self.barposition[i] - self.barposition[i-1]) > -1. * unwrap_threshold: + running_number_of_rotations -= 1 + else: + if (self.barposition[i] - self.barposition[i-1]) < unwrap_threshold: + running_number_of_rotations += 1 + number_of_rotations[i] = running_number_of_rotations + if reverse: + unwrapped_barposition = self.barposition + number_of_rotations * adjust + else: + unwrapped_barposition = self.barposition - number_of_rotations * adjust + self.pos = unwrapped_barposition + + def _frequency_and_derivative(self, smth_derivative,spline_derivative): + """ + Calculate the frequency and derivative of the unwrapped bar position. + + Parameters + ---------- + spline_derivative : int + Smoothing factor for the spline derivative calculation. + """ + # make a numerical derivative estimate + self.deriv = np.zeros_like(self.pos) + for i in range(1, len(self.pos) - 1): + self.deriv[i] = (self.pos[i+1] - self.pos[i-1]) / (2 * (self.time[i] - self.time[i-1])) + + if (smth_derivative): + smth_params = np.polyfit(self.time, self.deriv, smth_derivative) + pos_func = np.poly1d(smth_params) + self.deriv = pos_func(self.time) + # hard set as a cubic spline, + # number is a smoothing factor between knots, see scipy.UnivariateSpline + # + # recommended: 7 for dt=0.002 spacing + if spline_derivative: + spl = UnivariateSpline(self.time, self.pos, k=3, s=spline_derivative) + self.deriv = (spl.derivative())(self.time) + self.dderiv = np.zeros_like(self.deriv) + # + # can also do a second deriv + for indx, timeval in enumerate(self.time): + self.dderiv[indx] = spl.derivatives(timeval)[2] + + def print_bar(self, outfile): + """ + Print the bar position, its derivative, and time to a file. + + Parameters + ---------- + outfile : str + Path to the output file. + """ + with open(outfile, 'w') as f: + for i in range(len(self.time)): + print(self.time[i], self.pos[i], self.deriv[i], file=f) + return None class BarTransform(): @@ -502,77 +600,96 @@ def compute_bar_lag(ParticleInstance,rcut=0.01,verbose=0): +def find_barangle(time, BarInstance, interpolate=True): + """ + Use a bar instance to match the output time to a bar position. + Parameters + ---------- + time : array-like + Array of time values at which to find the bar position. + BarInstance : object + An instance of a class (such as `BarFromCoefficients`) that contains + bar positions and corresponding times. + interpolate : bool, optional + Whether to interpolate the bar position using a spline. If False, + the function finds the closest available bar position (default is True). -def find_barangle(time,BarInstance,interpolate=True): - ''' - # - # use a bar instance to match the output time to a bar position - # - # can take arrays! - # - # but feels like it only goes one direction? - # - ''' - # place in a guard against nan values - BarInstance.pos[BarInstance.pos == np.nan] = 0. + Returns + ------- + indx_barpos : array-like + Array of bar positions corresponding to the input time values. + Notes + ----- + This function can take arrays as input. It currently handles only one + direction of bar position matching and places a guard against NaN values + in the bar positions. + """ + # Place a guard against NaN values + BarInstance.pos[np.isnan(BarInstance.pos)] = 0. + + sord = 0 # Should this be a variable? + + if interpolate: + not_nan = np.where(np.isnan(BarInstance.pos) == False) + bar_func = UnivariateSpline(BarInstance.time[not_nan], -BarInstance.pos[not_nan], s=sord) - sord = 0 # should this be a variable? - # - if (interpolate): - not_nan = np.where(np.isnan(BarInstance.pos)==False) - bar_func = UnivariateSpline(BarInstance.time[not_nan],-BarInstance.pos[not_nan],s=sord) - # try: indx_barpos = np.zeros([len(time)]) - for indx,timeval in enumerate(time): - # - if (interpolate): + for indx, timeval in enumerate(time): + if interpolate: indx_barpos[indx] = bar_func(timeval) - # - # else: - indx_barpos[indx] = -BarInstance.pos[ abs(timeval-BarInstance.time).argmin()] - # - except: - if (interpolate): + indx_barpos[indx] = -BarInstance.pos[np.abs(timeval - BarInstance.time).argmin()] + except TypeError: # Catching a specific exception type is better practice + if interpolate: indx_barpos = bar_func(time) - # else: - indx_barpos = -BarInstance.pos[ abs(time-BarInstance.time).argmin()] - # - return indx_barpos - - + indx_barpos = -BarInstance.pos[np.abs(time - BarInstance.time).argmin()] -def find_barpattern(intime,BarInstance,smth_order=2): - ''' - # - # use a bar instance to match the output time to a bar pattern speed - # - # simple differencing--may want to be careful with this. - # needs a guard for the end points - # - ''' + return indx_barpos - # grab the derivative at whatever smoothing order - BarInstance.frequency_and_derivative(smth_order=smth_order) +def find_barpattern(intime, BarInstance, smth_order=2): + """ + Use a bar instance to match the output time to a bar pattern speed. + + Parameters + ---------- + intime : array-like + Array of time values at which to find the bar pattern speed. + BarInstance : object + An instance of a class (such as `BarFromCoefficients`) that contains + bar positions, times, and their derivatives. + smth_order : int, optional + Smoothing factor for the derivative calculation (default is 2). + + Returns + ------- + barpattern : array-like + Array of bar pattern speeds corresponding to the input time values. + """ + + # Compute the derivative of the bar position at the specified smoothing order + BarInstance._frequency_and_derivative(spline_derivative=smth_order) try: - + # Initialize an array to hold the bar pattern speeds barpattern = np.zeros([len(intime)]) - for indx,timeval in enumerate(intime): - - best_time = abs(timeval-BarInstance.time).argmin() + # Loop over each time value in the input array + for indx, timeval in enumerate(intime): + # Find the index of the closest time in BarInstance.time + best_time = abs(timeval - BarInstance.time).argmin() + # Get the derivative (bar pattern speed) at the closest time barpattern[indx] = BarInstance.deriv[best_time] except: + # Handle the case where intime is a single value + best_time = abs(intime - BarInstance.time).argmin() - best_time = abs(intime-BarInstance.time).argmin() - + # Get the derivative (bar pattern speed) at the closest time barpattern = BarInstance.deriv[best_time] return barpattern @@ -581,7 +698,10 @@ def find_barpattern(intime,BarInstance,smth_order=2): '''Not sure if this is the best place for this - wrote code to make a barfile using fourier analysis to find m=2 phase angle + then pattern speed based on this angle. This is to replace the EOF info if the EOF info is weird. Output file formats should be identical''' -class fourier_barfiles(): +class BarFromFourier(): + + def __init__(self,inputfiles): + def parse_list(self): f = open(self.slist) diff --git a/exptool/analysis/trapping.py b/exptool/analysis/trapping.py index 946b580..6b0d88f 100644 --- a/exptool/analysis/trapping.py +++ b/exptool/analysis/trapping.py @@ -5,6 +5,7 @@ MSP 23 Dec 2017 Break out bar finding algorithms to the more general pattern.py MSP 1 Mar 2019 Work on homogenizing docstrings and general commenting MSP 27 Oct 2021 Enable flexible particle number handling +MSP 18 May 2024 Create HDF5 input/ouput CLASSES: @@ -287,7 +288,7 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= #orderid = IN[:r1length][aps] if self.verbose > 0: - print('exptool.ApsFinding.trapping._determin_r_aps: Current time: {0:4.3f}'.format(tval),end='\r', flush=True) + print('exptool.ApsFinding.trapping._determine_r_aps: Current time: {0:4.3f}'.format(tval),end='\r', flush=True) # the id of the orbit is preserved and used as the dictionary key for j in range(0,len(id)): @@ -1106,9 +1107,8 @@ def do_kmeans_dict(TrappingInstanceDict,BarInstance,\ norb = TrappingInstanceDict['norb'] nfamilies = len(criteria.keys()) if nfamilies == 0: - print('trapping.do_kmeans_dict: no families defined?') - #break - return + return ValueError('exptool.trapping.do_kmeans_dict: no families defined?') + # set up final array trapping_array = np.zeros([nfamilies,norb,len(BarInstance.time)],dtype='i1') From 92e414a4e799acccbac18c44a2b87c912c899aac Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Sat, 18 May 2024 19:43:02 +0100 Subject: [PATCH 09/36] more doc cleanups --- exptool/analysis/pattern.py | 356 +++++++++++++++++++----------------- 1 file changed, 192 insertions(+), 164 deletions(-) diff --git a/exptool/analysis/pattern.py b/exptool/analysis/pattern.py index 83a7a13..2e137ef 100644 --- a/exptool/analysis/pattern.py +++ b/exptool/analysis/pattern.py @@ -159,118 +159,117 @@ def print_bar(self, outfile): return None -class BarTransform(): - ''' - BarTransform : class to do the work to calculate the bar position and transform particles - - on it's own, BarTransform will reset the particles to be in the bar frame (planar transformation) - - inputs - ----------------------- - ParticleInstanceIn : the input PSP instance - bar_angle : (default=None) the known bar angle - rel_bar_angle : (default=0.) the desired rotation angle relative to the bar major axis, counterclockwise (known or computed) - minr : (default=0.) the MINIMUM radius of particles to use to compute the bar angle - maxr : (default=1.) the MAXIMUM radius of particles to use in compute the bar angle - - outputs - ----------------------- - None - (ParticleInstanceIn will be modified to be in the planar bar transformation) - - - helper routines - ----------------------- - calculate_transform_and_return : overwrite the input PSP instance to have the raw positions be transformed - bar_fourier_compute : use m=2 fourier to transform to the bar frame - - ''' - def __init__(self,ParticleInstanceIn,bar_angle=None,rel_bar_angle=0.,minr=0.,maxr=1.): +class BarTransform: + """ + BarTransform: A class to calculate the bar position and transform particles into the bar frame. - ''' - see documentation above + On its own, BarTransform will reset the particles to be in the bar frame (planar transformation). - ''' + Parameters + ---------- + ParticleInstanceIn : object + The input particle instance. + bar_angle : float, optional + The known bar angle. If None, it will be computed (default is None). + rel_bar_angle : float, optional + The desired rotation angle relative to the bar major axis, counterclockwise (default is 0.). + minr : float, optional + The minimum radius of particles to use to compute the bar angle (default is 0.). + maxr : float, optional + The maximum radius of particles to use to compute the bar angle (default is 1.). + + Attributes + ---------- + ParticleInstanceIn : object + The input particle instance. + bar_angle : float + The computed or provided bar angle. + data : dict + Dictionary containing the transformed particle data. + time : float + The time of the particle instance. + filename : str + The filename of the particle instance. + comp : str + The component of the particle instance. + + Methods + ------- + calculate_transform_and_return() + Modify the input particle instance to be in the bar frame. + bar_fourier_compute(posx, posy, minr=0., maxr=1.) + Use m=2 Fourier analysis to compute the bar angle. + """ + def __init__(self, ParticleInstanceIn, bar_angle=None, rel_bar_angle=0., minr=0., maxr=1.): self.ParticleInstanceIn = ParticleInstanceIn - self.bar_angle = bar_angle - self.data = dict() - if self.bar_angle == None: - self.bar_angle = -1.*BarTransform.bar_fourier_compute(self,self.ParticleInstanceIn.data['x'],self.ParticleInstanceIn.data['y'],maxr=maxr) + if self.bar_angle is None: + self.bar_angle = -1. * self.bar_fourier_compute(self.ParticleInstanceIn.data['x'], self.ParticleInstanceIn.data['y'], minr=minr, maxr=maxr) - #-1.*BarTransform.bar_fourier_compute(self,self.ParticleInstanceIn.xpos,self.ParticleInstanceIn.ypos,maxr=maxr) - - # do an arbitary rotation of the particles relative to the bar? + # Apply relative bar angle rotation self.bar_angle += rel_bar_angle + # Perform the transformation self.calculate_transform_and_return() - def calculate_transform_and_return(self): - ''' - calculate_transform_and_return - do the modification of the input PSP instance to be in the bar frame. - - inputs - ---------------------------- - self (BarTransform) - - - ''' - - - transformed_x = self.ParticleInstanceIn.data['x']*np.cos(self.bar_angle) - self.ParticleInstanceIn.data['y']*np.sin(self.bar_angle) - #self.ParticleInstanceIn.xpos*np.cos(self.bar_angle) - self.ParticleInstanceIn.ypos*np.sin(self.bar_angle) - transformed_y = self.ParticleInstanceIn.data['x']*np.sin(self.bar_angle) + self.ParticleInstanceIn.data['y']*np.cos(self.bar_angle) - #self.ParticleInstanceIn.xpos*np.sin(self.bar_angle) + self.ParticleInstanceIn.ypos*np.cos(self.bar_angle) - - transformed_vx = self.ParticleInstanceIn.data['vx']*np.cos(self.bar_angle) - self.ParticleInstanceIn.data['vy']*np.sin(self.bar_angle) - #self.ParticleInstanceIn.xvel*np.cos(self.bar_angle) - self.ParticleInstanceIn.yvel*np.sin(self.bar_angle) - transformed_vy = self.ParticleInstanceIn.data['vx']*np.sin(self.bar_angle) + self.ParticleInstanceIn.data['vy']*np.cos(self.bar_angle) - #self.ParticleInstanceIn.xvel*np.sin(self.bar_angle) + self.ParticleInstanceIn.yvel*np.cos(self.bar_angle) + """ + Modify the input particle instance to be in the bar frame. + """ + # Transform positions + transformed_x = self.ParticleInstanceIn.data['x'] * np.cos(self.bar_angle) - self.ParticleInstanceIn.data['y'] * np.sin(self.bar_angle) + transformed_y = self.ParticleInstanceIn.data['x'] * np.sin(self.bar_angle) + self.ParticleInstanceIn.data['y'] * np.cos(self.bar_angle) + # Transform velocities + transformed_vx = self.ParticleInstanceIn.data['vx'] * np.cos(self.bar_angle) - self.ParticleInstanceIn.data['vy'] * np.sin(self.bar_angle) + transformed_vy = self.ParticleInstanceIn.data['vx'] * np.sin(self.bar_angle) + self.ParticleInstanceIn.data['vy'] * np.cos(self.bar_angle) + # Update the data dictionary self.data['x'] = transformed_x self.data['y'] = transformed_y self.data['z'] = np.copy(self.ParticleInstanceIn.data['z']) - #np.copy(self.ParticleInstanceIn.zpos) # interesting. needs to be a copy for later operations to work! - self.data['vx'] = transformed_vx self.data['vy'] = transformed_vy self.data['vz'] = np.copy(self.ParticleInstanceIn.data['vz']) - #np.copy(self.ParticleInstanceIn.zvel) - self.data['m'] = self.ParticleInstanceIn.data['m'] - #self.ParticleInstanceIn.mass self.data['potE'] = self.ParticleInstanceIn.data['potE'] - #self.ParticleInstanceIn.pote + # Update metadata self.time = self.ParticleInstanceIn.time self.filename = self.ParticleInstanceIn.filename self.comp = self.ParticleInstanceIn.comp + def bar_fourier_compute(self, posx, posy, minr=0., maxr=1.): + """ + Use x and y positions to compute the m=2 Fourier phase angle. - def bar_fourier_compute(self,posx,posy,minr=0.,maxr=1.): - ''' - - use x and y positions to compute the m=2 power, and find phase angle - - TODO: - generalize to transform to any azimuthal order? - - ''' - w = np.where( ( (posx*posx + posy*posy)**0.5 > minr ) & ((posx*posx + posy*posy)**0.5 < maxr ))[0] - - aval = np.sum( np.cos( 2.*np.arctan2(posy[w],posx[w]) ) ) - bval = np.sum( np.sin( 2.*np.arctan2(posy[w],posx[w]) ) ) - - return np.arctan2(bval,aval)/2. + Parameters + ---------- + posx : array-like + x positions of particles. + posy : array-like + y positions of particles. + minr : float, optional + Minimum radius to consider (default is 0.). + maxr : float, optional + Maximum radius to consider (default is 1.). + + Returns + ------- + float + The m=2 phase angle. + """ + radius = np.sqrt(posx**2 + posy**2) + w = np.where((radius > minr) & (radius < maxr))[0] + aval = np.sum(np.cos(2. * np.arctan2(posy[w], posx[w]))) + bval = np.sum(np.sin(2. * np.arctan2(posy[w], posx[w]))) + return np.arctan2(bval, aval) / 2. @@ -566,40 +565,6 @@ def read_bar(self,infile): - - - -def compute_bar_lag(ParticleInstance,rcut=0.01,verbose=0): - ''' - # - # simple fourier method to calculate where the particles are in relation to the bar - # - ''' - R = (ParticleInstance.data['x']*ParticleInstance.data['x'] + ParticleInstance.data['y']*ParticleInstance.data['y'])**0.5 - #(ParticleInstance.xpos*ParticleInstance.xpos + ParticleInstance.ypos*ParticleInstance.ypos)**0.5 - TH = np.arctan2(ParticleInstance.data['y'],ParticleInstance.data['x']) - #np.arctan2(ParticleInstance.ypos,ParticleInstance.xpos) - loR = np.where( R < rcut)[0] - A2 = np.sum(ParticleInstance.mass[loR] * np.cos(2.*TH[loR])) - B2 = np.sum(ParticleInstance.mass[loR] * np.sin(2.*TH[loR])) - bar_angle = 0.5*np.arctan2(B2,A2) - - if (verbose): - print('Position angle is {0:4.3f} . . .'.format(bar_angle)) - - # - # two steps: - # 1. rotate theta so that the bar is aligned at 0,2pi - # 2. fold onto 0,pi to compute the lag - # - tTH = (TH - bar_angle + np.pi/2.) % np.pi # compute lag with bar at pi/2 - # - # verification plot - #plt.scatter( R[0:10000]*np.cos(tTH[0:10000]-np.pi/2.),R[0:10000]*np.sin(tTH[0:10000]-np.pi/2.),color='black',s=0.5) - return tTH - np.pi/2. # retransform to bar at 0 - - - def find_barangle(time, BarInstance, interpolate=True): """ Use a bar instance to match the output time to a bar position. @@ -642,7 +607,7 @@ def find_barangle(time, BarInstance, interpolate=True): indx_barpos[indx] = bar_func(timeval) else: indx_barpos[indx] = -BarInstance.pos[np.abs(timeval - BarInstance.time).argmin()] - except TypeError: # Catching a specific exception type is better practice + except TypeError: if interpolate: indx_barpos = bar_func(time) else: @@ -695,81 +660,144 @@ def find_barpattern(intime, BarInstance, smth_order=2): return barpattern -'''Not sure if this is the best place for this - wrote code to make a barfile using fourier -analysis to find m=2 phase angle + then pattern speed based on this angle. This is to replace -the EOF info if the EOF info is weird. Output file formats should be identical''' -class BarFromFourier(): +class BarFromFourier: + """ + Class to compute the bar pattern speed and phase angle using Fourier analysis. + + This class replaces the EOF information if it appears to be incorrect. The output file + formats are designed to be identical to those generated using EOF information. + + Parameters + ---------- + inputfiles : str + Path to a file containing a list of input files to be processed. + + Attributes + ---------- + slist : str + Path to the input files list. + SLIST : array-like + List of input files parsed from the file. + pos : array-like + Array of bar positions. + deriv : array-like + Array of bar pattern speeds (derivatives). + time : array-like + Array of time steps corresponding to the bar positions and speeds. + """ - def __init__(self,inputfiles): + def __init__(self, inputfiles): + self.slist = inputfiles def parse_list(self): - - f = open(self.slist) - s_list = [] - for line in f: - d = [q for q in line.split()] - s_list.append(d[0]) + """ + Parse the list of input files from the provided file path. + This method reads the file specified in `self.slist` and stores the list + of input files in the attribute `self.SLIST`. + """ + with open(self.slist, 'r') as f: + s_list = [line.split()[0] for line in f] self.SLIST = np.array(s_list) + def bar_fourier_compute(self, posx, posy, maxr=0.5, minr=0.001): + """ + Compute the m=2 Fourier phase angle from particle positions. - def bar_fourier_compute(self,posx,posy,maxr=0.5, minr=.001): - - # - # use x and y positions tom compute the m=2 power, and find phase angle - # - w = np.where( ((posx*posx + posy*posy)**0.5 < maxr) & - ((posx*posx + posy*posy)**0.5 > minr) )[0] + Parameters + ---------- + posx : array-like + x positions of particles. + posy : array-like + y positions of particles. + maxr : float, optional + Maximum radius to consider (default is 0.5). + minr : float, optional + Minimum radius to consider (default is 0.001). + + Returns + ------- + float + The m=2 phase angle. + """ + # Select particles within the specified radius range + radius = np.sqrt(posx**2 + posy**2) + w = np.where((radius < maxr) & (radius > minr))[0] - aval = np.sum( np.cos( 2.*np.arctan2(posy[w],posx[w]) ) ) - bval = np.sum( np.sin( 2.*np.arctan2(posy[w],posx[w]) ) ) + # Compute m=2 Fourier components + aval = np.sum(np.cos(2. * np.arctan2(posy[w], posx[w]))) + bval = np.sum(np.sin(2. * np.arctan2(posy[w], posx[w]))) - return np.arctan2(bval,aval)/2. + return np.arctan2(bval, aval) / 2. def bar_speed(self, filelist, comp='star'): + """ + Compute the bar pattern speed from the list of input files. + + Parameters + ---------- + filelist : array-like + List of input files to process. + comp : str, optional + Component to analyze (default is 'star'). + + Returns + ------- + dict + Dictionary containing arrays of time steps, bar positions, and bar pattern speeds. + """ self.slist = filelist - fourier_barfiles.parse_list(self) - pos = particle.Input(self.SLIST[0],comp=comp,verbose=0) - pos_p1 = particle.Input(self.SLIST[1],comp=comp,verbose=0) + self.parse_list() + + pos = particle.Input(self.SLIST[0], comp=comp, verbose=0) + pos_p1 = particle.Input(self.SLIST[1], comp=comp, verbose=0) first_bar_angle = self.bar_fourier_compute(pos.data['x'], pos.data['y']) - #get time step + + # Calculate the timestep timestep = pos_p1.time - pos.time tt = np.array([]) pp = np.array([]) rot = np.array([]) - for i in range(0,len(self.SLIST)): - #loop through snapshot files in simulation, open file - pos = particle.Input(self.SLIST[i],comp=comp,verbose=0) - #compute bar angle + + for i in range(len(self.SLIST)): + pos = particle.Input(self.SLIST[i], comp=comp, verbose=0) bar_angle = self.bar_fourier_compute(pos.data['x'], pos.data['y']) - #if first time step, old bar angle = current bar angle + if i == 0: old_bar_angle = first_bar_angle - pattern_speed = (old_bar_angle-bar_angle)/(timestep) - #bar_angle > old_bar_angle, if difference is near 180, flip it (took this bit from rachel) - if abs(old_bar_angle-bar_angle)>=(np.pi*3/4): - pattern_speed = (old_bar_angle - (bar_angle + np.pi))/(timestep) + + pattern_speed = (old_bar_angle - bar_angle) / timestep + + if abs(old_bar_angle - bar_angle) >= (np.pi * 3 / 4): + pattern_speed = (old_bar_angle - (bar_angle + np.pi)) / timestep + pp = np.append(pp, bar_angle) rot = np.append(rot, pattern_speed) tt = np.append(tt, pos.time) - #use this bar angle as 'old' angle for next step + old_bar_angle = bar_angle + self.pos = pp self.deriv = rot self.time = tt - return {'time':tt,'pos':pp, 'deriv':rot} - - def print_bar(self,simulation_directory,simulation_name): - - # - # print the barfile to file - # - - - f = open(simulation_directory+simulation_name+'fourier_barpos.dat','w') - for i in range(0,len(self.SLIST)): - print(self.time[i],self.pos[i],self.deriv[i],end="\n",file=f) + + return {'time': tt, 'pos': pp, 'deriv': rot} - f.close() + def print_bar(self, simulation_directory, simulation_name): + """ + Print the bar positions and pattern speeds to a file. - return None \ No newline at end of file + Parameters + ---------- + simulation_directory : str + Directory where the output file will be saved. + simulation_name : str + Base name for the output file. + """ + output_file = simulation_directory + simulation_name + 'fourier_barpos.dat' + + with open(output_file, 'w') as f: + for i in range(len(self.SLIST)): + print(self.time[i], self.pos[i], self.deriv[i], file=f) + + return None From b07d2237f4a18add3b18c088b1ad61070fa0a118 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Sat, 18 May 2024 20:06:43 +0100 Subject: [PATCH 10/36] improve coefficient bar determination --- exptool/analysis/pattern.py | 27 ++++++++++++++++++--------- 1 file changed, 18 insertions(+), 9 deletions(-) diff --git a/exptool/analysis/pattern.py b/exptool/analysis/pattern.py index 2e137ef..c7ca007 100644 --- a/exptool/analysis/pattern.py +++ b/exptool/analysis/pattern.py @@ -16,8 +16,10 @@ -Filtering algorithms for bar determination (e.g. look at better time-series algorithms) -Need a partial pattern calculator for bars that grow and disappear. Perhaps also in eof.py? +-Filter bad values from the pattern speed of the bar somehow +-Combine multiple coefficient series for a better estimate -BASIC USAGE: +BASIC USAGE Examples: # to transform a PSP output to have the bar on the X axis PSPTransform = pattern.BarTransform(PSPInput) @@ -69,7 +71,7 @@ class BarFromCoefficients: spline_derivative : int, optional Smoothing factor for the spline derivative calculation (default is 0). """ - def __init__(self, times, coefs, unwrap_threshold=-1., smooth=False, reverse=False, adjust=np.pi, verbose=0, smth_derivative=0, spline_derivative=0): + def __init__(self, times, coefs, unwrap_threshold=-np.pi/2., smooth=False, reverse=False, adjust=2*np.pi, verbose=0, smth_derivative=0, spline_derivative=0): self.verbose = verbose self.time = times self.cos = np.real(coefs) @@ -97,19 +99,25 @@ def _unwrap_position(self, unwrap_threshold, smooth, reverse, adjust): """ running_number_of_rotations = 0 number_of_rotations = np.zeros_like(self.barposition) - # start from zero + + # which way are we rotating + primarydirection = np.nanmedian(np.ediff1d(B.barposition)) + + # start from the beginning and keep track of number of rotations for i in range(1, len(self.barposition)): - if reverse: - if (self.barposition[i] - self.barposition[i-1]) > -1. * unwrap_threshold: - running_number_of_rotations -= 1 - else: + if (primarydirection > 0): if (self.barposition[i] - self.barposition[i-1]) < unwrap_threshold: running_number_of_rotations += 1 + else: + if (self.barposition[i] - self.barposition[i-1]) > -1. * unwrap_threshold: + running_number_of_rotations += 1 number_of_rotations[i] = running_number_of_rotations - if reverse: + + # now straighten out the zeros + if (primarydirection > 0): unwrapped_barposition = self.barposition + number_of_rotations * adjust else: - unwrapped_barposition = self.barposition - number_of_rotations * adjust + unwrapped_barposition = - self.barposition + number_of_rotations * adjust self.pos = unwrapped_barposition def _frequency_and_derivative(self, smth_derivative,spline_derivative): @@ -160,6 +168,7 @@ def print_bar(self, outfile): + class BarTransform: """ BarTransform: A class to calculate the bar position and transform particles into the bar frame. From 5b8fb60bc8bcc19cbfbbd95c3e4326d69a556d50 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Sat, 18 May 2024 20:07:56 +0100 Subject: [PATCH 11/36] typo in attributes --- exptool/analysis/pattern.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exptool/analysis/pattern.py b/exptool/analysis/pattern.py index c7ca007..a44b936 100644 --- a/exptool/analysis/pattern.py +++ b/exptool/analysis/pattern.py @@ -101,7 +101,7 @@ def _unwrap_position(self, unwrap_threshold, smooth, reverse, adjust): number_of_rotations = np.zeros_like(self.barposition) # which way are we rotating - primarydirection = np.nanmedian(np.ediff1d(B.barposition)) + primarydirection = np.nanmedian(np.ediff1d(self.barposition)) # start from the beginning and keep track of number of rotations for i in range(1, len(self.barposition)): From 8e6a8730d10c12917b07934395db657637e79e65 Mon Sep 17 00:00:00 2001 From: michael-petersen Date: Mon, 20 May 2024 06:46:09 +0100 Subject: [PATCH 12/36] add docs; add Beane 2024 criteria --- exptool/analysis/trapping.py | 491 ++++++++++++++++------------------- exptool/utils/kmeans.py | 303 +++++++++++---------- 2 files changed, 385 insertions(+), 409 deletions(-) diff --git a/exptool/analysis/trapping.py b/exptool/analysis/trapping.py index 6b0d88f..f3f34e0 100644 --- a/exptool/analysis/trapping.py +++ b/exptool/analysis/trapping.py @@ -178,7 +178,7 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= #tstamp = datetime.datetime.strftime(datetime.datetime.now(), '%Y-%m-%d+%H') # create a new file using the particle number, runtag, and time - outputfile = out_directory+'RadialAps_N{}_r{}_T{}.dat'.format(total_orbits,runtag,tstamp) + outputfile = out_directory+'RadialAps_N{}_r{}_T{}.h5'.format(total_orbits,runtag,tstamp) f = h5py.File(outputfile,"w") # createdescriptor string @@ -440,324 +440,272 @@ def read_trapping_file(t_file,tdtype='i1'): -def reduce_aps_dictionary(TrappingInstance,norb): - ''' - sometimes you just don't need all those apsides - ''' +def reduce_aps_dictionary(TrappingInstance, norb): + """ + Reduces the apsides data in a trapping instance to only include a specified number of orbits. + + This function takes a trapping instance dictionary and reduces it to include only the + specified number of orbits (`norb`). It retains the description and the first `norb` + orbits' data. + + Parameters + ---------- + TrappingInstance : dict + A dictionary containing the trapping instance data with multiple orbits. + norb : int + The number of orbits to include in the reduced dictionary. + + Returns + ------- + TrappingInstanceOut : dict + A reduced dictionary containing only the specified number of orbits + from the original trapping instance. + + Examples + -------- + >>> trapping_instance = { + 'desc': 'Sample trapping instance', + 0: {'apside_1': [1, 2], 'apside_2': [3, 4]}, + 1: {'apside_1': [5, 6], 'apside_2': [7, 8]}, + 2: {'apside_1': [9, 10], 'apside_2': [11, 12]} + } + >>> reduced_instance = reduce_aps_dictionary(trapping_instance, 2) + >>> print(reduced_instance) + {'norb': 2, 'desc': 'Sample trapping instance', 0: {'apside_1': [1, 2], 'apside_2': [3, 4]}, 1: {'apside_1': [5, 6], 'apside_2': [7, 8]}} + """ + # Initialize the output dictionary TrappingInstanceOut = {} + + # Add the number of orbits to the output dictionary TrappingInstanceOut['norb'] = norb + + # Add the description to the output dictionary TrappingInstanceOut['desc'] = TrappingInstance['desc'] - for i in range(0,norb): + # Loop through the specified number of orbits and add them to the output dictionary + for i in range(norb): TrappingInstanceOut[i] = TrappingInstance[i] return TrappingInstanceOut - -def evaluate_clusters_polar_legacy(K,maxima=False,rank=False,perc=0.): - ''' - evaluate_clusters_polar - calculate statistics for clusters in polar coordinates - - inputs - ------------- - K : number of clusters - maxima : (boolean, False) if True, use the maximum value from the clusters - rank - perc - - - returns - ------------- - theta_n - clustermean - clusterstd_r - clusterstd_t - - - - ''' +def beane_criteria(K): + """Implement the criteria from Beane et al. (2024) + + works best for polar classifications""" # how many clusters? k = K.K - if (rank) & (perc==0.): - print('evaluate_clusters_polar: Perc must be >0.') - return np.nan,np.nan,np.nan,np.nan + # Compute radii and theta values from clusters + rad_clusters = np.array([np.linalg.norm(K.clusters[i], axis=1) for i in range(k)]) + the_clusters = np.array([np.arctan2(np.abs(K.clusters[i][:, 1]), np.abs(K.clusters[i][:, 0])) for i in range(k)]) + # implement equation A1: the maximum angle from the bar for the clusters + thetadiff = np.max([np.arctan2(np.abs(K.mu[i][1]), np.abs(K.mu[i][0])) for i in range(k)]) - # compute radii and theta values from clusters - rad_clusters = np.array([np.sum(np.array(K.clusters[i])*np.array(K.clusters[i]),axis=1)**0.5 for i in range(0,k)]) + # if thetadiff < pi/8, consider the particle trapped - # for computing the theta values, can decide on a version with - # (legacy) or without (modern) folding + # implement equation A2: + clusterstd = np.sum([np.std(rad_clusters[i]) for i in range(k)]) + clustermean = np.sum([np.mean(rad_clusters[i]) for i in range(k)]) + tightness = clusterstd/clustermean - legacy = True + # if tightness is < 0.22, consider the particle trapped - if legacy: - the_clusters = np.array([np.arctan(np.abs(np.array(K.clusters[i])[:,1])/np.abs(np.array(K.clusters[i])[:,0])) for i in range(0,k)]) - else: - the_clusters = np.array([np.arctan( (np.array(K.clusters[i])[:,1])/np.abs(np.array(K.clusters[i])[:,0])) for i in range(0,k)]) - - if maxima: - # use maxima + return thetadiff,tightness - clustermean = np.max([np.mean(rad_clusters[i]) for i in range(0,k)]) - - if legacy: - theta_n = np.max([np.abs(np.arctan( K.mu[i][1]/K.mu[i][0])) for i in range(0,k)]) - else: - theta_n = np.max([ np.arctan(np.abs(K.mu[i][1]/K.mu[i][0])) for i in range(0,k)]) - - if rank: - # use rank ordered - - organized_rad = np.array([rad_clusters[i][rad_clusters[i].argsort()] for i in range(0,k)]) - organized_the = np.array([the_clusters[i][the_clusters[i].argsort()] for i in range(0,k)]) - - clusterstd_r = np.max([np.percentile(organized_rad[i] - np.mean(rad_clusters[i]),perc) for i in range(0,k)]) - - # I think this needs an absolute value - clusterstd_t = np.max([np.percentile(organized_the[i] - np.mean(the_clusters[i]),perc) for i in range(0,k)]) - - else: - clusterstd_r = np.max([np.std(rad_clusters[i]) for i in range(0,k)]) - - if legacy: - clusterstd_t = np.max([np.std(the_clusters[i]) for i in range(0,k)]) - else: - clusterstd_t = np.max([np.abs(np.max(the_clusters[i])-np.min(the_clusters[i])) for i in range(0,k)]) - - - else: - - # not maxima - - if rank: - # use rank ordered - - organized_rad = np.array([rad_clusters[i][rad_clusters[i].argsort()] for i in range(0,k)]) - organized_the = np.array([the_clusters[i][the_clusters[i].argsort()] for i in range(0,k)]) - - clusterstd_r = np.mean([np.percentile(organized_rad[i] - np.mean(rad_clusters[i]),perc) for i in range(0,k)]) - clusterstd_t = np.mean([np.percentile(organized_the[i] - np.mean(the_clusters[i]),perc) for i in range(0,k)]) - - else: - clusterstd_r = np.mean([np.std(rad_clusters[i]) for i in range(0,k)]) - - if legacy: - clusterstd_t = np.mean([np.std(the_clusters[i]) for i in range(0,k)]) - else: - clusterstd_t = np.mean([np.abs(np.max(the_clusters[i])-np.min(the_clusters[i])) for i in range(0,k)]) - - - - clustermean = np.mean([np.mean(rad_clusters[i]) for i in range(0,k)]) - - # compute the mean of the cluster centers - if legacy: - theta_n = np.mean([np.abs(np.arctan( K.mu[i][1]/K.mu[i][0])) for i in range(0,k)]) - else: - theta_n = np.mean([ np.arctan(np.abs(K.mu[i][1]/K.mu[i][0])) for i in range(0,k)]) - - # return values - return theta_n,clustermean,clusterstd_r,clusterstd_t - - -def evaluate_clusters_polar(K,maxima=False,rank=False,perc=0.): - ''' - evaluate_clusters_polar - calculate statistics for clusters in polar coordinates - - inputs - ------------- - K - maxima - rank - perc - - - returns - ------------- - theta_n - clustermean - clusterstd_r - clusterstd_t +def evaluate_clusters_polar(K, maxima=False, rank=False, perc=0.): + """ + Calculate statistics for clusters in polar coordinates. + This function evaluates the clustering results in polar coordinates (r, theta). + It computes various statistics such as mean, standard deviation, and optionally + ranks and percentiles. - ''' - - # how many clusters? + Parameters + ---------- + K : KMeans + An instance of a K-means clustering result. + maxima : bool, optional + If True, calculate maximum quantities. If False, calculate average quantities. + Default is False. + rank : bool, optional + If True, use rank-ordered statistics. Default is False. + perc : float, optional + Percentage threshold for rank ordering. Default is 0. + + Returns + ------- + theta_n : float + Angle measure in the context of the clusters. + clustermean : float + The mean value of the clusters. + clusterstd_r : float + The standard deviation of the clusters in the radial direction. + clusterstd_t : float + The standard deviation of the clusters in the angular direction. + + Notes + ----- + This function computes the radii and theta values from the clusters, then calculates + either the maximum or average statistics based on the `maxima` parameter. It also + handles rank-ordered statistics if `rank` is True and `perc` is greater than 0. + + Examples + -------- + >>> K = kmeans.KMeans(k=2, X=ApsArray) + >>> K.find_centers() + >>> theta_n, clustermean, clusterstd_r, clusterstd_t = evaluate_clusters_polar(K) + """ + # Number of clusters k = K.K - if (rank) & (perc==0.): - print('evaluate_clusters_polar: Perc must be >0.') - return np.nan,np.nan,np.nan,np.nan + # Check if rank is True but perc is not set + if rank and perc == 0.: + raise SyntaxError('exptool.trapping.evaluate_clusters_polar: Perc must be >0.') - - # compute radii and theta values from clusters - rad_clusters = np.array([np.sum(np.array(K.clusters[i])*np.array(K.clusters[i]),axis=1)**0.5 for i in range(0,k)]) - the_clusters = np.array([np.arctan(np.abs(np.array(K.clusters[i])[:,1])/np.abs(np.array(K.clusters[i])[:,0])) for i in range(0,k)]) + # Compute radii and theta values from clusters + rad_clusters = np.array([np.linalg.norm(K.clusters[i], axis=1) for i in range(k)]) + the_clusters = np.array([np.arctan2(np.abs(K.clusters[i][:, 1]), np.abs(K.clusters[i][:, 0])) for i in range(k)]) if maxima: - # use maxima - - clustermean = np.max([np.mean(rad_clusters[i]) for i in range(0,k)]) - - #theta_n = np.max([abs(np.arctan(K.mu[i][1]/K.mu[i][0])) for i in range(0,k)]) - theta_n = np.max([np.arctan(np.abs(K.mu[i][1]/K.mu[i][0])) for i in range(0,k)]) + # Calculate maximum quantities + clustermean = np.max([np.mean(rad_clusters[i]) for i in range(k)]) + theta_n = np.max([np.arctan2(np.abs(K.mu[i][1]), np.abs(K.mu[i][0])) for i in range(k)]) if rank: - # use rank ordered - - organized_rad = np.array([rad_clusters[i][rad_clusters[i].argsort()] for i in range(0,k)]) - organized_the = np.array([the_clusters[i][the_clusters[i].argsort()] for i in range(0,k)]) - - clusterstd_r = np.max([np.percentile(organized_rad[i] - np.mean(rad_clusters[i]),perc) for i in range(0,k)]) - clusterstd_t = np.max([np.percentile(organized_the[i] - np.mean(the_clusters[i]),perc) for i in range(0,k)]) - + # Use rank-ordered statistics + organized_rad = np.array([np.sort(rad_clusters[i]) for i in range(k)]) + organized_the = np.array([np.sort(the_clusters[i]) for i in range(k)]) + clusterstd_r = np.max([np.percentile(organized_rad[i] - np.mean(rad_clusters[i]), perc) for i in range(k)]) + clusterstd_t = np.max([np.percentile(organized_the[i] - np.mean(the_clusters[i]), perc) for i in range(k)]) else: - clusterstd_r = np.max([np.std(rad_clusters[i]) for i in range(0,k)]) - clusterstd_t = np.max([np.std(the_clusters[i]) for i in range(0,k)]) - + clusterstd_r = np.max([np.std(rad_clusters[i]) for i in range(k)]) + clusterstd_t = np.max([np.std(the_clusters[i]) for i in range(k)]) else: - - # not maxima - + # Calculate average quantities if rank: - # use rank ordered - - organized_rad = np.array([rad_clusters[i][rad_clusters[i].argsort()] for i in range(0,k)]) - organized_the = np.array([the_clusters[i][the_clusters[i].argsort()] for i in range(0,k)]) - - clusterstd_r = np.mean([np.percentile(organized_rad[i] - np.mean(rad_clusters[i]),perc) for i in range(0,k)]) - clusterstd_t = np.mean([np.percentile(organized_the[i] - np.mean(the_clusters[i]),perc) for i in range(0,k)]) - + # Use rank-ordered statistics + organized_rad = np.array([np.sort(rad_clusters[i]) for i in range(k)]) + organized_the = np.array([np.sort(the_clusters[i]) for i in range(k)]) + clusterstd_r = np.mean([np.percentile(organized_rad[i] - np.mean(rad_clusters[i]), perc) for i in range(k)]) + clusterstd_t = np.mean([np.percentile(organized_the[i] - np.mean(the_clusters[i]), perc) for i in range(k)]) else: - clusterstd_r = np.mean([np.std(rad_clusters[i]) for i in range(0,k)]) - clusterstd_t = np.mean([np.std(the_clusters[i]) for i in range(0,k)]) - - clustermean = np.mean([np.mean(rad_clusters[i]) for i in range(0,k)]) - #theta_n = np.mean([abs(np.arctan(K.mu[i][1]/K.mu[i][0])) for i in range(0,k)]) - theta_n = np.mean([np.arctan(np.abs(K.mu[i][1]/K.mu[i][0])) for i in range(0,k)]) - - return theta_n,clustermean,clusterstd_r,clusterstd_t - + clusterstd_r = np.mean([np.std(rad_clusters[i]) for i in range(k)]) + clusterstd_t = np.mean([np.std(the_clusters[i]) for i in range(k)]) + clustermean = np.mean([np.mean(rad_clusters[i]) for i in range(k)]) + theta_n = np.mean([np.arctan2(np.abs(K.mu[i][1]), np.abs(K.mu[i][0])) for i in range(k)]) + return theta_n, clustermean, clusterstd_r, clusterstd_t -def process_kmeans_polar(ApsArray,indx=-1,k=2,maxima=False,rank=False,perc=0.): - ''' - # - # robust kmeans implementation - # - # -can be edited for speed - # -confined to two dimensions - # -computes trapping metrics in polar coordinates - - inputs - ---------- - ApsArray : the array of aps for an individual orbit - indx : a designation of the orbit, for use with multiprocessing - k : the number of clusters - maxima : calculate average (if False) or maximum (if True) quantities - mad : toggle median absolute deviation calculation +def process_kmeans_polar(ApsArray, indx=-1, k=2, maxima=False, rank=False, perc=0.): + """ + Perform robust K-means clustering on apsidal data in polar coordinates. + This function performs K-means clustering on the provided apsidal array, + computes trapping metrics in polar coordinates, and handles potential + edge cases where clusters may have very few points. - returns + Parameters ---------- - theta_n : (see explanation at beginning for definitions) - clustermean : - clusterstd_r : - clusterstd_theta : - kmeans_plus_flag : - - - - - ''' + ApsArray : array-like + The array of apsides for an individual orbit. Each element should contain + the (r, theta) coordinates of an apsis. + indx : int, optional + A designation of the orbit, for use with multiprocessing. Default is -1. + k : int, optional + The number of clusters to form. Default is 2. + maxima : bool, optional + Calculate average (if False) or maximum (if True) quantities. Default is False. + rank : bool, optional + Toggle ranking. Default is False. + perc : float, optional + Percentage threshold for ranking. Default is 0. + + Returns + ------- + theta_n : float + Some angle measure in the context of the clusters. + clustermean : float + The mean value of the clusters. + clusterstd_r : float + The standard deviation of the clusters in the radial direction. + clusterstd_theta : float + The standard deviation of the clusters in the angular direction. + kmeans_plus_flag : int + Indicator flag: 0 for successful basic K-means, 1 for successful K-means++, + 2 for failure in both K-means and K-means++. + + Notes + ----- + This implementation confines the clustering to two dimensions and includes + robustness checks to handle small cluster sizes. In case of failure in + basic K-means, it retries using the K-means++ initialization method. + + Examples + -------- + >>> ApsArray = np.array([[1.0, 0.0], [2.0, 1.0], [1.5, 0.5], [3.0, 1.5]]) + >>> theta_n, clustermean, clusterstd_r, clusterstd_theta, flag = process_kmeans_polar(ApsArray) + """ kmeans_plus_flag = 0 - K = kmeans.KMeans(k,X=ApsArray) + K = kmeans.KMeans(k, X=ApsArray) K.find_centers() - # add an evaluation for if a cluster ends up with only X members, here hard coded to 2 - - - # find the standard deviation of clusters - - # first, check to make sure no single-point clusters were detected - # set rejection threshold + # Minimum cluster size threshold min_cluster_size = 1 - try: + clustersize = np.array([np.array(K.clusters[c]).size / 2. for c in range(k)]) - clustersize = np.array([np.array(K.clusters[c]).size/2. for c in range(0,k)]) - - # eliminate + # Ensure no single-point clusters while np.min(clustersize) <= min_cluster_size: w = np.where(clustersize > min_cluster_size)[0] - new_aps = np.array([np.concatenate([np.array(K.clusters[x])[:,0] for x in w]),\ - np.concatenate([np.array(K.clusters[x])[:,1] for x in w])]).T + new_aps = np.array([np.concatenate([np.array(K.clusters[x])[:, 0] for x in w]), \ + np.concatenate([np.array(K.clusters[x])[:, 1] for x in w])]).T - K = kmeans.KMeans(k,X=new_aps) + K = kmeans.KMeans(k, X=new_aps) K.find_centers() - clustersize = np.array([np.array(K.clusters[c]).size/2. for c in range(0,k)]) + clustersize = np.array([np.array(K.clusters[c]).size / 2. for c in range(k)]) - theta_n,clustermean,clusterstd_r,clusterstd_t = \ - evaluate_clusters_polar(K,maxima=maxima,rank=rank,perc=perc) + theta_n, clustermean, clusterstd_r, clusterstd_theta = \ + evaluate_clusters_polar(K, maxima=maxima, rank=rank, perc=perc) - - - # failure on basic kmeans except: - K = kmeans.KPlusPlus(2,X=ApsArray) + # If basic K-means fails, try K-means++ + K = kmeans.KPlusPlus(k, X=ApsArray) K.init_centers() K.find_centers(method='++') kmeans_plus_flag = 1 try: - - clustersize = np.array([np.array(K.clusters[c]).size/2. for c in range(0,k)]) - + clustersize = np.array([np.array(K.clusters[c]).size / 2. for c in range(k)]) while np.min(clustersize) <= min_cluster_size: w = np.where(clustersize > min_cluster_size)[0] - new_aps = np.array([np.concatenate([np.array(K.clusters[x])[:,0] for x in w]),\ - np.concatenate([np.array(K.clusters[x])[:,1] for x in w])]).T + new_aps = np.array([np.concatenate([np.array(K.clusters[x])[:, 0] for x in w]), \ + np.concatenate([np.array(K.clusters[x])[:, 1] for x in w])]).T - K = kmeans.KPlusPlus(k,X=new_aps) + K = kmeans.KPlusPlus(k, X=new_aps) K.init_centers() K.find_centers(method='++') - clustersize = np.array([np.array(K.clusters[c]).size/2. for c in range(0,k)]) - - - theta_n,clustermean,clusterstd_r,clusterstd_t = \ - evaluate_clusters_polar(K,maxima=maxima,rank=rank,perc=perc) - + clustersize = np.array([np.array(K.clusters[c]).size / 2. for c in range(k)]) + theta_n, clustermean, clusterstd_r, clusterstd_theta = \ + evaluate_clusters_polar(K, maxima=maxima, rank=rank, perc=perc) - # failure mode for advanced kmeans except: - - # - # would like a more intelligent way to diagnose - #if indx >= 0: - # print 'Orbit %i even failed in Kmeans++!!' %indx + # If both methods fail, set all outputs to NaN clusterstd_r = np.nan - clusterstd_t = np.nan + clusterstd_theta = np.nan clustermean = np.nan theta_n = np.nan kmeans_plus_flag = 2 - - - return theta_n,clustermean,clusterstd_r,clusterstd_t,kmeans_plus_flag - - + return theta_n, clustermean, clusterstd_r, clusterstd_theta, kmeans_plus_flag @@ -892,29 +840,43 @@ def process_kmeans(ApsArray,indx=-1,k=2,maxima=False,mad=False): +def transform_aps(ApsArray, BarInstance): + """ + Transform the apsides array into the bar frame of reference. + This function transforms the apsides array, aligning it with the bar frame as determined + by the BarInstance. The transformation is offloaded for clarity. -def transform_aps(ApsArray,BarInstance): - ''' - transform_aps : simple transformation for the aps array, offloaded for clarity. + Parameters + ---------- + ApsArray : np.ndarray + The array of apsides, where each row represents a time step and contains + [time, x_position, y_position]. + BarInstance : object + An instance that contains information about the bar's position and motion. + + Returns + ------- + np.ndarray + The transformed positions in the bar frame. The output array has the same number of rows + as ApsArray and two columns corresponding to the transformed x and y positions. + + Notes + ----- + This transformation assumes that the bar motion is in one direction. + """ - inputs - ------------------ - ApsArray : the array of apsides - BarInstance : + # Find the bar angle positions corresponding to the times in ApsArray + bar_positions = pattern.find_barangle(ApsArray[:, 0], BarInstance) - outputs - ------------------ - X : + # Initialize the output array for transformed positions + X = np.zeros([len(ApsArray[:, 1]), 2]) - stuck in one direction, watch out - ''' - bar_positions = pattern.find_barangle(ApsArray[:,0],BarInstance) - X = np.zeros([len(ApsArray[:,1]),2]) - X[:,0] = ApsArray[:,1]*np.cos(bar_positions) - ApsArray[:,2]*np.sin(bar_positions) - X[:,1] = -ApsArray[:,1]*np.sin(bar_positions) - ApsArray[:,2]*np.cos(bar_positions) - return X + # Apply the transformation to align with the bar frame + X[:, 0] = ApsArray[:, 1] * np.cos(bar_positions) - ApsArray[:, 2] * np.sin(bar_positions) + X[:, 1] = -ApsArray[:, 1] * np.sin(bar_positions) - ApsArray[:, 2] * np.cos(bar_positions) + return X def do_single_kmeans_step(TrappingInstanceDict,BarInstance,desired_time,\ @@ -1476,7 +1438,7 @@ def do_kmeans_multi(TrappingInstanceDict,BarInstance,\ print('Total trapping calculation took {0:3.2f} seconds, or {1:3.2f} milliseconds per orbit.'.format(time.time()-t1, 1.e3*(time.time()-t1)/len(TrappingInstanceDict))) # go through the dictionary of trapping criteria and re-make the arrays - trapped = {} + trapped = dict() for nfam,family in enumerate(np.array(list(criteria.keys()))): @@ -1532,4 +1494,3 @@ def re_form_trapping_arrays(array,array_number): return net_array -#warnings.filterwarnings("ignore",category =RuntimeWarning) diff --git a/exptool/utils/kmeans.py b/exptool/utils/kmeans.py index a4c79b9..d139d7c 100644 --- a/exptool/utils/kmeans.py +++ b/exptool/utils/kmeans.py @@ -1,119 +1,97 @@ -# -# kmeans.py -# -# robust implementation of kmeans based on -# https://datasciencelab.wordpress.com +""" +kmeans.py -import random -import numpy as np +purpose-built, robust implementation of kmeans based on +https://datasciencelab.wordpress.com -import matplotlib.pyplot as plt -import matplotlib +MSP 19 May 2024 Improve documentation; code cleanup +""" -class KMeans(): - ''' - class to implement K-means in Python +import numpy as np +import random +class KMeans: + """ + A class to implement K-means clustering in Python. + """ - ''' - def __init__(self, K, X=None, N=0): - ''' - initialize kmeans - - - inputs - -------------- - self : KMeans class - K : number of clusters - X : array of observations - N : error guard - - - returns - ------------- - self : Kmeans class - - - ''' - + """ + Initialize KMeans. + + Parameters + ---------- + K : int + Number of clusters. + X : array-like, optional + Array of observations. + N : int, optional + Number of points (needed if X is not provided). + + Raises + ------ + Exception + If no data is provided and N is not specified. + """ self.K = K - try: - - tmp = len(X) + if X is not None: self.X = X self.N = len(X) - - except: - + else: if N == 0: - raise Exception("kmeans.KMeans: If no data is provided, \ - a parameter N (number of points) is needed") + raise Exception("kmeans.KMeans: If no data is provided, a parameter N (number of points) is needed") else: self.N = N self.X = self._init_board_gauss(N, K) self.mu = None + self.oldmu = None self.clusters = None self.method = None - def _init_board_gauss(self, N, k): - ''' - _init_board_gauss - initialize the guess points - - - inputs - ------------------- - self - N - k - - - returns - ------------------ - self - X : randomly partitioned clusters - - - ''' - - # number of points to put in each cluster - n = float(N)/k - + def _init_board_gauss(self, N, K): + """ + Initialize the guess points using a Gaussian distribution. + + Parameters + ---------- + N : int + Number of points. + K : int + Number of clusters. + + Returns + ------- + np.array + Randomly partitioned clusters. + """ + n = float(N) / K X = [] - # set up - for i in range(k): - - c = (random.uniform(-1,1), random.uniform(-1,1)) - - s = random.uniform(0.05,0.15) - - # just reflecting--but this means that the clusters are forced to ahve the same number of points - x = [] - while len(x) < n: - - a,b = np.array([np.random.normal(c[0],s),np.random.normal(c[1],s)]) - - # Continue drawing points from the distribution in the range [-1,1] - if abs(a) and abs(b)<1: - x.append([a,b]) - - X.extend(x) + for i in range(K): + c = (random.uniform(-1, 1), random.uniform(-1, 1)) + s = random.uniform(0.05, 0.15) + + cluster_points = [] + while len(cluster_points) < n: + a, b = np.array([np.random.normal(c[0], s), np.random.normal(c[1], s)]) + if abs(a) < 1 and abs(b) < 1: + cluster_points.append([a, b]) + X.extend(cluster_points) X = np.array(X)[:N] return X - def _cluster_points(self): + """ + Assign each point to the nearest cluster center. + """ mu = self.mu - clusters = {} + clusters = {} for x in self.X: - bestmukey = min([(i[0], np.linalg.norm(x-mu[i[0]])) \ - for i in enumerate(mu)], key=lambda t:t[1])[0] + bestmukey = min([(i[0], np.linalg.norm(x - mu[i[0]])) for i in enumerate(mu)], key=lambda t: t[1])[0] try: clusters[bestmukey].append(x) except KeyError: @@ -121,99 +99,136 @@ def _cluster_points(self): self.clusters = clusters def _reevaluate_centers(self): - ''' - _reevaluate_centers - draw new centers based on which center they are closest to - NOTE that this can create asymmetric cluster sizes - - ''' + """ + Compute new cluster centers based on the current cluster assignments. + """ clusters = self.clusters newmu = [] - keys = sorted(self.clusters.keys()) + keys = sorted(clusters.keys()) for k in keys: - newmu.append(np.mean(clusters[k], axis = 0)) + newmu.append(np.mean(clusters[k], axis=0)) self.mu = newmu def _has_converged(self): - ''' - _has_converged - check to see whether clusters change from one step to the next. if not, declare convergence! - - - ''' + """ + Check if the algorithm has converged. + + Returns + ------- + bool + True if the cluster centers do not change, False otherwise. + """ K = len(self.oldmu) - return(set([tuple(a) for a in self.mu]) == \ - set([tuple(a) for a in self.oldmu])\ - and len(set([tuple(a) for a in self.mu])) == K) + return (set([tuple(a) for a in self.mu]) == set([tuple(a) for a in self.oldmu]) and + len(set([tuple(a) for a in self.mu])) == K) - def find_centers(self, method='random',nitermax=1000): - ''' - find_centers - iteratively select new centers - - inputs - ---------------- - - - returns - --------------- - - - ''' - + def find_centers(self, method='random', nitermax=1000): + """ + Find the cluster centers using the specified method. + + Parameters + ---------- + method : str, optional + Method to initialize the cluster centers ('random' or '++', default is 'random'). + nitermax : int, optional + Maximum number of iterations (default is 1000). + + Returns + ------- + None + + Raises + ------ + ValueError + If the specified method is not supported. + """ self.method = method X = self.X K = self.K - - #self.oldmu = random.sample(X, K) - - # draw K samples from the array of clusters - self.oldmu = random.sample(list(X), K) - - if method != '++': + if method == '++': + # Initialize using K-means++ + self.mu = self._init_kmeans_plusplus() + elif method == 'random': # Initialize to K random centers - - # this has a python2/3 compatibility issue - #self.mu = random.sample(X, K) self.mu = random.sample(list(X), K) + else: + raise ValueError("Unsupported method: {}".format(method)) - # put in a guard against iter = 0 - - while (not self._has_converged()) & (iter= r)[0][0] - return(self.X[ind]) + return self.X[ind] def init_centers(self): - #self.mu = random.sample(self.X, 1) + """ + Initialize the cluster centers using K-means++ initialization. + + Returns + ------- + None + """ self.mu = random.sample(list(self.X), 1) while len(self.mu) < self.K: self._dist_from_centers() self.mu.append(self._choose_next_center()) - From b6116073d8281f1b70d85197c8ce1d61d8dff3f5 Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Fri, 12 Dec 2025 10:22:11 +0000 Subject: [PATCH 13/36] Update exptool/io/psp_to_hdf5.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exptool/io/psp_to_hdf5.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exptool/io/psp_to_hdf5.py b/exptool/io/psp_to_hdf5.py index d05c029..d9f5ece 100644 --- a/exptool/io/psp_to_hdf5.py +++ b/exptool/io/psp_to_hdf5.py @@ -1,5 +1,5 @@ """ -draft conversion from PSP format the HDF5. +draft conversion from PSP format to HDF5. For each component group, there is a subgroup named 'header', which stores header information related to that component. This information may include various parameters and metadata. The header information is organized into nested groups and attributes within the 'header' subgroup. The structure of the header data may vary depending on the specific PSP file format. From ca39f4fc32117535465959215d4f361973850bec Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 12 Dec 2025 10:22:33 +0000 Subject: [PATCH 14/36] Initial plan From 7bb0fa8687570736076ec14237fe2e19fad284e2 Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Fri, 12 Dec 2025 10:22:35 +0000 Subject: [PATCH 15/36] Update exptool/io/psp_to_hdf5.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exptool/io/psp_to_hdf5.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/exptool/io/psp_to_hdf5.py b/exptool/io/psp_to_hdf5.py index d9f5ece..50ae028 100644 --- a/exptool/io/psp_to_hdf5.py +++ b/exptool/io/psp_to_hdf5.py @@ -145,19 +145,19 @@ def print_component_header(f, O, comp): # Create attributes for the deepest level try: f['{}/header/{}/{}'.format(comp, key, subkey)].attrs.create(subsubkey, O.header[comp][key][subkey][subsubkey]) - except: + except (KeyError, AttributeError, TypeError): # Create subgroups if necessary f['{}/header/{}/{}'.format(comp, key, subkey)].create_group(subsubkey) f['{}/header/{}/{}'.format(comp, key, subkey)].attrs.create(subsubkey, O.header[comp][key][subkey][subsubkey]) - except: + except (KeyError, AttributeError, TypeError): # Create attributes for the intermediate level try: f['{}/header/{}/{}'.format(comp, key, subkey)].attrs.create(subkey, O.header[comp][key][subkey]) - except: + except (KeyError, AttributeError, TypeError): # Create subgroups if necessary f['{}/header'.format(comp)].create_group(key) f['{}/header/{}'.format(comp, key)].attrs.create(subkey, O.header[comp][key][subkey]) - except: + except (KeyError, AttributeError, TypeError): # Create attributes for the top-level header f['{}/header'.format(comp)].attrs.create(key, O.header[comp][key]) From c9f8c2cb29c349f8db272c792bfa9feea6eb869b Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Fri, 12 Dec 2025 10:22:46 +0000 Subject: [PATCH 16/36] Update exptool/analysis/pattern.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exptool/analysis/pattern.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exptool/analysis/pattern.py b/exptool/analysis/pattern.py index a44b936..78247d3 100644 --- a/exptool/analysis/pattern.py +++ b/exptool/analysis/pattern.py @@ -606,7 +606,7 @@ def find_barangle(time, BarInstance, interpolate=True): sord = 0 # Should this be a variable? if interpolate: - not_nan = np.where(np.isnan(BarInstance.pos) == False) + not_nan = np.where(~np.isnan(BarInstance.pos)) bar_func = UnivariateSpline(BarInstance.time[not_nan], -BarInstance.pos[not_nan], s=sord) try: From ff80e442e3ce0194fdd9056ab532ba5f33ca4012 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 12 Dec 2025 10:23:09 +0000 Subject: [PATCH 17/36] Initial plan From 252cb42b72b6795664d09ece253519c35a1f5113 Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Fri, 12 Dec 2025 10:23:28 +0000 Subject: [PATCH 18/36] Update exptool/analysis/pattern.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exptool/analysis/pattern.py | 1 - 1 file changed, 1 deletion(-) diff --git a/exptool/analysis/pattern.py b/exptool/analysis/pattern.py index 78247d3..ba38843 100644 --- a/exptool/analysis/pattern.py +++ b/exptool/analysis/pattern.py @@ -43,7 +43,6 @@ # exptool imports from ..io import particle from ..utils import kmeans -from ..utils import utils class BarFromCoefficients: From 559f7cce13c57eec9acc61d9291e9386ddb30447 Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Fri, 12 Dec 2025 10:24:06 +0000 Subject: [PATCH 19/36] Update exptool/analysis/pattern.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exptool/analysis/pattern.py | 1 - 1 file changed, 1 deletion(-) diff --git a/exptool/analysis/pattern.py b/exptool/analysis/pattern.py index ba38843..af58432 100644 --- a/exptool/analysis/pattern.py +++ b/exptool/analysis/pattern.py @@ -42,7 +42,6 @@ # exptool imports from ..io import particle -from ..utils import kmeans class BarFromCoefficients: From 09b235ba9877b800e8e360021f34045f92e9b410 Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Fri, 12 Dec 2025 10:24:32 +0000 Subject: [PATCH 20/36] Update exptool/analysis/pattern.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exptool/analysis/pattern.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exptool/analysis/pattern.py b/exptool/analysis/pattern.py index af58432..6d9d57f 100644 --- a/exptool/analysis/pattern.py +++ b/exptool/analysis/pattern.py @@ -657,7 +657,7 @@ def find_barpattern(intime, BarInstance, smth_order=2): # Get the derivative (bar pattern speed) at the closest time barpattern[indx] = BarInstance.deriv[best_time] - except: + except TypeError: # Handle the case where intime is a single value best_time = abs(intime - BarInstance.time).argmin() From b91a7020fe001bf66dd8a51c7caf6fa5d1dc5d1d Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 12 Dec 2025 10:25:54 +0000 Subject: [PATCH 21/36] Pass comp and verbose parameters to convert_psp_to_hdf5 method Co-authored-by: michael-petersen <19195541+michael-petersen@users.noreply.github.com> --- exptool/io/psp_to_hdf5.py | 24 +++++++++++++++++++++--- 1 file changed, 21 insertions(+), 3 deletions(-) diff --git a/exptool/io/psp_to_hdf5.py b/exptool/io/psp_to_hdf5.py index d9f5ece..d53ff78 100644 --- a/exptool/io/psp_to_hdf5.py +++ b/exptool/io/psp_to_hdf5.py @@ -88,22 +88,34 @@ def __init__(self, filename, comp=None, verbose=0): self.verbose = verbose # Start the conversion process - self.convert_psp_to_hdf5() + self.convert_psp_to_hdf5(self.filename, comp=self.comp, verbose=self.verbose) @staticmethod - def convert_psp_to_hdf5(inputfilename): + def convert_psp_to_hdf5(inputfilename, comp=None, verbose=0): """ Convert a PSP input file to HDF5 format. Args: inputfilename (str): The name of the PSP input file to be converted. + comp (str, optional): The specific component to convert. If provided, only the data for the specified component will be converted. + verbose (int, optional): Verbosity level for printing progress and messages during conversion. """ # Define the output file name outputfilename = inputfilename + '.h5' + if verbose > 0: + print(f"Converting {inputfilename} to {outputfilename}") + # Open the PSP input file and extract components O = particle.Input(inputfilename) - comps = list(O.header.keys()) + + # Determine which components to convert + if comp is not None: + comps = [comp] if comp in O.header.keys() else [] + if not comps and verbose > 0: + print(f"Warning: Component '{comp}' not found in file") + else: + comps = list(O.header.keys()) # Create a new HDF5 file for storing the converted data f = h5py.File(outputfilename, 'w') @@ -112,6 +124,9 @@ def convert_psp_to_hdf5(inputfilename): f['time'] = O.time for comp in comps: + if verbose > 0: + print(f"Processing component: {comp}") + # Create a group for each component f.create_group(comp) @@ -123,6 +138,9 @@ def convert_psp_to_hdf5(inputfilename): # Close the HDF5 file f.close() + + if verbose > 0: + print(f"Conversion complete: {outputfilename}") @staticmethod def print_component_header(f, O, comp): From 8fa24f5c8ccbef6cc7d45ed3f72fdf09dc250da6 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 12 Dec 2025 10:27:42 +0000 Subject: [PATCH 22/36] Fix bare except clauses with specific exception types Co-authored-by: michael-petersen <19195541+michael-petersen@users.noreply.github.com> --- exptool/analysis/pattern.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/exptool/analysis/pattern.py b/exptool/analysis/pattern.py index 78247d3..ff855fe 100644 --- a/exptool/analysis/pattern.py +++ b/exptool/analysis/pattern.py @@ -299,7 +299,7 @@ def __init__(self,**kwargs): self.read_bar(kwargs['file']) print('pattern.BarDetermine: BarInstance sucessfully read.') - except: + except (FileNotFoundError, IOError, ValueError): print('pattern.BarDetermine: no compatible bar file found.') @@ -559,7 +559,7 @@ def read_bar(self,infile): pos.append(q[1]) try: deriv.append(q[2]) - except: + except IndexError: pass self.time = np.array(time) @@ -659,7 +659,7 @@ def find_barpattern(intime, BarInstance, smth_order=2): # Get the derivative (bar pattern speed) at the closest time barpattern[indx] = BarInstance.deriv[best_time] - except: + except TypeError: # Handle the case where intime is a single value best_time = abs(intime - BarInstance.time).argmin() From 65efe5a20dba7fb064c47450eb40e143398376d6 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 12 Dec 2025 10:28:52 +0000 Subject: [PATCH 23/36] Fix typo: sucessfully -> successfully Co-authored-by: michael-petersen <19195541+michael-petersen@users.noreply.github.com> --- exptool/analysis/pattern.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exptool/analysis/pattern.py b/exptool/analysis/pattern.py index ff855fe..c479afa 100644 --- a/exptool/analysis/pattern.py +++ b/exptool/analysis/pattern.py @@ -297,7 +297,7 @@ def __init__(self,**kwargs): try: # check to see if bar file has already been created self.read_bar(kwargs['file']) - print('pattern.BarDetermine: BarInstance sucessfully read.') + print('pattern.BarDetermine: BarInstance successfully read.') except (FileNotFoundError, IOError, ValueError): print('pattern.BarDetermine: no compatible bar file found.') From 3989076f0371cfff8d99ef4815c423d613073c1b Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Fri, 12 Dec 2025 11:50:40 +0000 Subject: [PATCH 24/36] Update exptool/analysis/trapping.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exptool/analysis/trapping.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exptool/analysis/trapping.py b/exptool/analysis/trapping.py index f3f34e0..91e8458 100644 --- a/exptool/analysis/trapping.py +++ b/exptool/analysis/trapping.py @@ -163,7 +163,7 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= changingindx = True else: - raise ValueError("exptool.ApsFinding.trapping._determin_r_aps: particle_indx must be an integer or an array.") + raise ValueError("exptool.ApsFinding.trapping._determine_r_aps: particle_indx must be an integer or an array.") # sort the particle indices particle_indx = particle_indx[particle_indx.argsort()] From da25360b0e3d46fbf1b4b29c882a34e239ff68cc Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Fri, 12 Dec 2025 11:50:57 +0000 Subject: [PATCH 25/36] Update exptool/analysis/trapping.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exptool/analysis/trapping.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exptool/analysis/trapping.py b/exptool/analysis/trapping.py index 91e8458..17da5b7 100644 --- a/exptool/analysis/trapping.py +++ b/exptool/analysis/trapping.py @@ -1069,7 +1069,7 @@ def do_kmeans_dict(TrappingInstanceDict,BarInstance,\ norb = TrappingInstanceDict['norb'] nfamilies = len(criteria.keys()) if nfamilies == 0: - return ValueError('exptool.trapping.do_kmeans_dict: no families defined?') + raise ValueError('exptool.trapping.do_kmeans_dict: no families defined?') # set up final array From 00fdae16a4d94d17fb5ecf59ce51cbc9bd86d1c4 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 12 Dec 2025 11:50:59 +0000 Subject: [PATCH 26/36] Initial plan From aefa89c78ee399d096db4c431e8323af3630d152 Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Fri, 12 Dec 2025 11:51:16 +0000 Subject: [PATCH 27/36] Update exptool/analysis/pattern.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exptool/analysis/pattern.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exptool/analysis/pattern.py b/exptool/analysis/pattern.py index 6d9d57f..47d736a 100644 --- a/exptool/analysis/pattern.py +++ b/exptool/analysis/pattern.py @@ -599,7 +599,7 @@ def find_barangle(time, BarInstance, interpolate=True): in the bar positions. """ # Place a guard against NaN values - BarInstance.pos[np.isnan(BarInstance.pos)] = 0. + BarInstance.pos = np.nan_to_num(BarInstance.pos, nan=0.0) sord = 0 # Should this be a variable? From 2c19abddb992b38d2d81b8067db01a291cb432ef Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Fri, 12 Dec 2025 11:51:31 +0000 Subject: [PATCH 28/36] Update exptool/analysis/trapping.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exptool/analysis/trapping.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exptool/analysis/trapping.py b/exptool/analysis/trapping.py index 17da5b7..95591b5 100644 --- a/exptool/analysis/trapping.py +++ b/exptool/analysis/trapping.py @@ -697,7 +697,7 @@ def process_kmeans_polar(ApsArray, indx=-1, k=2, maxima=False, rank=False, perc= theta_n, clustermean, clusterstd_r, clusterstd_theta = \ evaluate_clusters_polar(K, maxima=maxima, rank=rank, perc=perc) - except: + except Exception: # If both methods fail, set all outputs to NaN clusterstd_r = np.nan clusterstd_theta = np.nan From 908c94f8107a291c04b7debc71af2e99418afe9b Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Fri, 12 Dec 2025 11:51:45 +0000 Subject: [PATCH 29/36] Update exptool/io/psp_to_hdf5.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exptool/io/psp_to_hdf5.py | 1 - 1 file changed, 1 deletion(-) diff --git a/exptool/io/psp_to_hdf5.py b/exptool/io/psp_to_hdf5.py index 50ae028..9f49df4 100644 --- a/exptool/io/psp_to_hdf5.py +++ b/exptool/io/psp_to_hdf5.py @@ -28,7 +28,6 @@ import h5py outputfilename = 'OUT.run0.00000.h5' -outputfilename = 'OUT.run0.00000.h5' # how is the global header information saved? # only time is saved From 5a670bdb230240a9df88f7f4869a2c8e112a51a4 Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Fri, 12 Dec 2025 11:51:57 +0000 Subject: [PATCH 30/36] Update exptool/analysis/trapping.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exptool/analysis/trapping.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exptool/analysis/trapping.py b/exptool/analysis/trapping.py index 95591b5..6ec3aad 100644 --- a/exptool/analysis/trapping.py +++ b/exptool/analysis/trapping.py @@ -181,7 +181,7 @@ def determine_r_aps(self,filelist,comp,particle_indx=-1,runtag='',out_directory= outputfile = out_directory+'RadialAps_N{}_r{}_T{}.h5'.format(total_orbits,runtag,tstamp) f = h5py.File(outputfile,"w") - # createdescriptor string + # create descriptor string desc = 'apsfile for '+comp+' in '+out_directory+', norbits='+str(total_orbits)+', threedee='+str(threedee)+', using '+filelist # Write the descriptor string as an attribute From ec9e4683447fa12eda5b4ffd25dfa4fc6f1f576d Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 12 Dec 2025 11:52:40 +0000 Subject: [PATCH 31/36] Initial plan From b968b2f8fbf2a0368785d82288aea5b3b254f26a Mon Sep 17 00:00:00 2001 From: Michael Petersen Date: Fri, 12 Dec 2025 11:53:10 +0000 Subject: [PATCH 32/36] Update exptool/io/psp_to_hdf5.py Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- exptool/io/psp_to_hdf5.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exptool/io/psp_to_hdf5.py b/exptool/io/psp_to_hdf5.py index 9f49df4..e7c879d 100644 --- a/exptool/io/psp_to_hdf5.py +++ b/exptool/io/psp_to_hdf5.py @@ -177,4 +177,4 @@ def make_phasespace(f, inputfilename, comp): PS = np.array([O1.data['m'], O1.data['x'], O1.data['y'], O1.data['z'], O1.data['vx'], O1.data['vy'], O1.data['vz'], O1.data['potE']]).T # Store the phase space data as a dataset - dset = f[comp].create_dataset('phasespace', data=PS) + f[comp].create_dataset('phasespace', data=PS) From 79c874e53827a195c541e3ada0208f960228793a Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 12 Dec 2025 11:56:02 +0000 Subject: [PATCH 33/36] Make convert_psp_to_hdf5 and related methods instance methods Co-authored-by: michael-petersen <19195541+michael-petersen@users.noreply.github.com> --- exptool/io/psp_to_hdf5.py | 27 ++++++++++++--------------- 1 file changed, 12 insertions(+), 15 deletions(-) diff --git a/exptool/io/psp_to_hdf5.py b/exptool/io/psp_to_hdf5.py index 50ae028..7fdcbfb 100644 --- a/exptool/io/psp_to_hdf5.py +++ b/exptool/io/psp_to_hdf5.py @@ -19,7 +19,7 @@ args = parser.parse_args() - HDFConverter.convert_psp_to_hdf5(args.filename) + HDFConverter(args.filename) @@ -90,19 +90,17 @@ def __init__(self, filename, comp=None, verbose=0): # Start the conversion process self.convert_psp_to_hdf5() - @staticmethod - def convert_psp_to_hdf5(inputfilename): + def convert_psp_to_hdf5(self): """ Convert a PSP input file to HDF5 format. - Args: - inputfilename (str): The name of the PSP input file to be converted. + Uses the filename stored in self.filename during initialization. """ # Define the output file name - outputfilename = inputfilename + '.h5' + outputfilename = self.filename + '.h5' # Open the PSP input file and extract components - O = particle.Input(inputfilename) + O = particle.Input(self.filename) comps = list(O.header.keys()) # Create a new HDF5 file for storing the converted data @@ -116,16 +114,15 @@ def convert_psp_to_hdf5(inputfilename): f.create_group(comp) # Print the header information for the component - HDFConverter.print_component_header(f, O, comp) + self.print_component_header(f, O, comp) # Create and store the phase space data for the component - HDFConverter.make_phasespace(f, inputfilename, comp) + self.make_phasespace(f, comp) # Close the HDF5 file f.close() - @staticmethod - def print_component_header(f, O, comp): + def print_component_header(self, f, O, comp): """ Print header information for a component to an HDF5 file. @@ -161,18 +158,18 @@ def print_component_header(f, O, comp): # Create attributes for the top-level header f['{}/header'.format(comp)].attrs.create(key, O.header[comp][key]) - @staticmethod - def make_phasespace(f, inputfilename, comp): + def make_phasespace(self, f, comp): """ Convert and store phase space data for a component in an HDF5 file. Args: f (h5py.Group): The HDF5 group to store the phase space data. - inputfilename (str): The name of the PSP input file. comp (str): The name of the component. + + Uses the filename stored in self.filename during initialization. """ # Read data from the PSP input file - O1 = particle.Input(inputfilename, comp) + O1 = particle.Input(self.filename, comp) # Create a phase space array PS = np.array([O1.data['m'], O1.data['x'], O1.data['y'], O1.data['z'], O1.data['vx'], O1.data['vy'], O1.data['vz'], O1.data['potE']]).T From 06ef64685634b81e664baa9e36b0a65c1f76fd77 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 12 Dec 2025 11:58:07 +0000 Subject: [PATCH 34/36] Fix self.oldmu initialization issue in KMeans.find_centers() Co-authored-by: michael-petersen <19195541+michael-petersen@users.noreply.github.com> --- exptool/utils/kmeans.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/exptool/utils/kmeans.py b/exptool/utils/kmeans.py index d139d7c..36358cb 100644 --- a/exptool/utils/kmeans.py +++ b/exptool/utils/kmeans.py @@ -155,6 +155,9 @@ def find_centers(self, method='random', nitermax=1000): else: raise ValueError("Unsupported method: {}".format(method)) + # Initialize oldmu to mu before the first iteration + self.oldmu = self.mu + iter = 0 while not self._has_converged() and iter < nitermax: self.oldmu = self.mu From f78b0c0d9f81abb865476503f4617d4efac5fb77 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 12 Dec 2025 12:00:39 +0000 Subject: [PATCH 35/36] Improve oldmu initialization to ensure convergence check works properly Co-authored-by: michael-petersen <19195541+michael-petersen@users.noreply.github.com> --- exptool/utils/kmeans.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/exptool/utils/kmeans.py b/exptool/utils/kmeans.py index 36358cb..7f32a71 100644 --- a/exptool/utils/kmeans.py +++ b/exptool/utils/kmeans.py @@ -155,8 +155,9 @@ def find_centers(self, method='random', nitermax=1000): else: raise ValueError("Unsupported method: {}".format(method)) - # Initialize oldmu to mu before the first iteration - self.oldmu = self.mu + # Initialize oldmu to ensure first iteration runs + # Use values guaranteed to be different from initial mu + self.oldmu = [np.array([float('inf')] * len(X[0]))] * K iter = 0 while not self._has_converged() and iter < nitermax: From 94dd250b649ffa4871dea51e6ad8f199f5bd61a7 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Fri, 12 Dec 2025 12:01:56 +0000 Subject: [PATCH 36/36] Use list comprehension to create independent array objects in oldmu Co-authored-by: michael-petersen <19195541+michael-petersen@users.noreply.github.com> --- exptool/utils/kmeans.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/exptool/utils/kmeans.py b/exptool/utils/kmeans.py index 7f32a71..f843bbf 100644 --- a/exptool/utils/kmeans.py +++ b/exptool/utils/kmeans.py @@ -157,7 +157,7 @@ def find_centers(self, method='random', nitermax=1000): # Initialize oldmu to ensure first iteration runs # Use values guaranteed to be different from initial mu - self.oldmu = [np.array([float('inf')] * len(X[0]))] * K + self.oldmu = [np.array([float('inf')] * len(X[0])) for _ in range(K)] iter = 0 while not self._has_converged() and iter < nitermax: