Improve Performance of Symmetry Constraints via Lookup Table

[AndrewSazonov]

Problem

Applying symmetry constraints is currently slow, primarily due to the use of the Cryspy method:
get_symop_pcentr_multiplicity_letter_site_symmetry_coords_xyz_2(...),
which probably parses the wyckoff.dat file.

Proposed

We should replace it with a prebuilt lookup table stored as a Python file (e.g. space_group_lookup_table.py) containing a dictionary that maps IT_number and IT_coordinate_system_code to preformatted symmetry data.

This dictionary would store values such as:

• crystal_system (e.g. "orthorhombic")
• IT_coordinate_system_code (e.g. "cab")
• symmetry operations (coords_xyz)
• multiplicity, letter, site_symmetry, etc.

Where

• All values except two can be directly extracted from wyckoff.dat or obtained via Cryspy.
• IT_coordinate_system_code (e.g. cab) is not in wyckoff.dat, but can be obtained via Cryspy.
• crystal_system is also not present and needs to be fetched from Cryspy.

Implementation

• Create space_group_lookup_table.py (or similar) under easydiffraction/crystallography/.
• Import the dictionary into crystallography.py.
• Replace the call at line 89:

  result = get_symop_pcentr_multiplicity_letter_site_symmetry_coords_xyz_2(it_number, it_coordinate_system_code)

with a direct lookup from our dictionary. Probably, the following code will need to be adapted.

Example of the Suggested Dictionary:

space_groups = {
    (62, "cab"): {
        "IT_number": 62,
        "setting": 3,
        "IT_coordinate_system_code": "cab",
        "name_H-M_alt": "P b n m",
        "crystal_system": "orthorhombic",
        "Wyckoff_positions": {
            "d": {
                "multiplicity": 8,
                "site_symmetry": "1",
                "coords_xyz": [
                    "(x,y,z)", "(x+1/2,-y+1/2,-z)", "(-x,-y,z+1/2)", "(-x+1/2,y+1/2,-z+1/2)",
                    "(-x,-y,-z)", "(-x+1/2,y+1/2,z)", "(x,y,-z+1/2)", "(x+1/2,-y+1/2,z+1/2)"
                ]
            },
            "c": {
                "multiplicity": 4,
                "site_symmetry": ".m.",
                "coords_xyz": [
                    "(x,y,1/4)", "(x+1/2,-y+1/2,3/4)", "(-x,-y,3/4)", "(-x+1/2,y+1/2,1/4)"
                ]
            },
            "b": {
                "multiplicity": 4,
                "site_symmetry": "-1",
                "coords_xyz": [
                    "(1/2,0,0)", "(0,1/2,0)", "(1/2,0,1/2)", "(0,1/2,1/2)"
                ]
            },
            "a": {
                "multiplicity": 4,
                "site_symmetry": "-1",
                "coords_xyz": [
                    "(0,0,0)", "(1/2,1/2,0)", "(0,0,1/2)", "(1/2,1/2,1/2)"
                ]
            }
        }
    },
    (199, "1"): {
        "IT_number": 199,
        "setting": 1,
        "IT_coordinate_system_code": "1",
        "name_H-M_alt": "I 21 3",
        "crystal_system": "cubic",
        "Wyckoff_positions": {
            "c": {
                "multiplicity": 24,
                "site_symmetry": "1",
                "coords_xyz": [
                    "(x,y,z)", "(-x+1/2,-y,z+1/2)", "(-x,y+1/2,-z+1/2)", "(x+1/2,-y+1/2,-z)",
                    "(z,x,y)", "(z+1/2,-x+1/2,-y)", "(-z+1/2,-x,y+1/2)", "(-z,x+1/2,-y+1/2)",
                    "(y,z,x)", "(-y,z+1/2,-x+1/2)", "(y+1/2,-z+1/2,-x)", "(-y+1/2,-z,x+1/2)"
                ]
            },
            "b": {
                "multiplicity": 12,
                "site_symmetry": "2..",
                "coords_xyz": [
                    "(x,0,1/4)", "(-x+1/2,0,3/4)", "(1/4,x,0)",
                    "(3/4,-x+1/2,0)", "(0,1/4,x)", "(0,3/4,-x+1/2)"
                ]
            },
            "a": {
                "multiplicity": 8,
                "site_symmetry": ".3.",
                "coords_xyz": [
                    "(x,x,x)", "(-x+1/2,-x,x+1/2)", "(-x,x+1/2,-x+1/2)", "(x+1/2,-x+1/2,-x)"
                ]
            }
        }
    }
}

[elindgren]

I've started working on this issue. Currently, I'm using the following snippet to read the table REFERENCE_TABLE_IT_COORDINATE_SYSTEM_CODE_NAME_HM_EXTENDED

# Parse table from CrysPy with all coordinate names and system codes
# Also includes the shortened HM name

data = {}

for entry in REFERENCE_TABLE_IT_COORDINATE_SYSTEM_CODE_NAME_HM_EXTENDED:
    IT_number, IT_coordinate_system_code, hm_name = entry
    crystal_system = get_crystal_system_by_it_number(IT_number)

    # This logic is from the cryspy internals, and is the internal conversion
    # between system_code and choice (equivalent to our `setting`?)
    if IT_coordinate_system_code is None:
        choice = 1
    elif "2" in IT_coordinate_system_code:
        choice = 2
    elif "h" in IT_coordinate_system_code:
        if crystal_system.startswith("trigonal"):
            choice = 2
        else:  # hexagonal
            choice = 1
    elif "3" in IT_coordinate_system_code:
        choice = 3
    else:
        choice = 1
    
    try:
        symop_2, p_centr_new, l_multiplicity, l_letter, l_site_symmetry, l_coord_xyz_2 = get_symop_pcentr_multiplicity_letter_site_symmetry_coords_xyz_2(IT_number, IT_coordinate_system_code)
    except (TypeError, UnboundLocalError, ValueError) as e:
        print(f"Space group {IT_number}-{IT_coordinate_system_code} failed with", type(e))
        continue

    key = (IT_number, IT_coordinate_system_code)
    data[key] = {
            "IT_number": IT_number,
            "setting": choice,
            "IT_coordinate_system_code": IT_coordinate_system_code,
            "name_H-M_alt": hm_name,
            "crystal_system": crystal_system,
            "Wyckoff_positions": {}
    }

    for mult, letter, site_symmetry, xyz in zip(l_multiplicity, l_letter, l_site_symmetry, l_coord_xyz_2):
        data[key]['Wyckoff_positions'][letter] = {
                "multiplicity": mult,
                "site_symmetry": site_symmetry,
                "coords_xyz": xyz
    }

print(len(data.keys()))  # 613 of 711 entries in the REFERENCE_TABLE are captured.

output_file = 'space_group_lookup_table.py'

# Write dictionary to file with clean formatting
with open(output_file, 'w') as f:
    f.write("space_group_lookup_table = {\n")
    for k, v in data.items():
        f.write(f"  {repr(k)}: {repr(v)},\n")
    f.write("}\n")


from space_group_lookup_table import space_group_lookup_table

print(space_group_lookup_table[(62, "cab")])

@AndrewSazonov I have two questions/issues here:

1. What is the "setting" parameter in your example above? Is that equivalent to choice in my code/in cryspy? I don't think so, since in your example you have setting=3 for 62-'cab', but according to cryspy the choice variable is 1 in that case.
2. I can only successfully query get_symop_pcentr_multiplicity_letter_site_symmetry_coords_xyz_2 with 613 of the 711 rows in the REFERENCE_TABLE_IT_COORDINATE_SYSTEM_CODE_NAME_HM_EXTENDED. The other 98 rows fail with either a TypeError, UnboundLocalError or ValueError in subroutines of the function. I assume this is an upstream issue in cryspy. Do you have any insight into this?

[elindgren]

A very quick and rough benchmark indicates that the memo:ed table is about 4 orders of magnitude faster. I haven't been able to figure out why cryspy is so slow, but it does perform quite a lot of calculations via subroutines for each call to get_symop_pcentr_multiplicity_letter_site_symmetry_coords_xyz_2.

Cryspy call time (s):          0.004170659780502319
Memo:ed table call time (s):   6.532669067382813e-08

Benchmark script:

from cryspy.A_functions_base.function_2_space_group import (
    get_symop_pcentr_multiplicity_letter_site_symmetry_coords_xyz_2,
)

from space_group_lookup_table import space_group_lookup_table

from time import time

N = 1000

start = time()

for i in range(N):
    symop_2, p_centr_new, l_multiplicity, l_letter, l_site_symmetry, l_coord_xyz_2 = get_symop_pcentr_multiplicity_letter_site_symmetry_coords_xyz_2(62, "1")
time_cryspy = time() - start
print('Cryspy call time (s): '.ljust(30), time_cryspy / N)

key = (62, "cab")
start = time()
for i in range(N):
    res = space_group_lookup_table[key]
time_tabulated = time() - start
print('Memo:ed table call time (s): '.ljust(30), time_tabulated / N)

[AndrewSazonov]

What is the "setting" parameter in your example above? Is that equivalent to choice in my code/in cryspy? I don't think so, since in your example you have setting=3 for 62-'cab', but according to cryspy the choice variable is 1 in that case.

The setting simply represents the index or sequence number of a given IT_coordinate_system_code for that space group.

In that example, space group number 62 has six different coordinate system settings (from 1 to 6), each corresponding to a unique IT_coordinate_system_code. So when we say setting=3 for space group 62, we mean it corresponds to the third listed coordinate system (which in this case is "cab").

Here is a part of wyckoff.dat, where the second number on the first line indicates the setting:

62  3  5  1  4   P b n m
8	d	1
(x,y,z)	(x+1/2,-y+1/2,-z)	(-x,-y,z+1/2)	(-x+1/2,y+1/2,-z+1/2)
(-x,-y,-z)	(-x+1/2,y+1/2,z)	(x,y,-z+1/2)	(x+1/2,-y+1/2,z+1/2)
4	c	.m.
(x,y,1/4)	(x+1/2,-y+1/2,3/4)	(-x,-y,3/4)	(-x+1/2,y+1/2,1/4)
4	b	-1
(1/2,0,0)	(0,1/2,0)	(1/2,0,1/2)	(0,1/2,1/2)
4	a	-1
(0,0,0)	(1/2,1/2,0)	(0,0,1/2)	(1/2,1/2,1/2)

[elindgren]

I wrote a quick profiling script that just creates a SampleModel and applies the symmetry constraints. This operation takes from 0.1s -> 0.007s. It seems like the next biggest bottleneck is the call to sympify.

Profiling results on develop

Profiling results on the feature branch

Profiling script

from easydiffraction import SampleModel
import cProfile
import pstats

def main():
    model = SampleModel('hs')
    model.space_group.name_h_m = 'R -3 m'
    model.space_group.it_coordinate_system_code = 'h'
    model.cell.length_a = 6.9
    model.cell.length_c = 14.1
    model.atom_sites.add('Zn', 'Zn', 0, 0, 0.5, wyckoff_letter='b', b_iso=0.5)
    model.atom_sites.add('Cu', 'Cu', 0.5, 0, 0, wyckoff_letter='e', b_iso=0.5)
    model.atom_sites.add('O', 'O', 0.21, -0.21, 0.06, wyckoff_letter='h', b_iso=0.5)
    model.atom_sites.add('Cl', 'Cl', 0, 0, 0.197, wyckoff_letter='c', b_iso=0.5)
    model.atom_sites.add('H', '2H', 0.13, -0.13, 0.08, wyckoff_letter='h', b_iso=0.5)
    model.apply_symmetry_constraints()

if __name__ == "__main__":
    with cProfile.Profile() as pr:
        main()
    stats = pstats.Stats(pr)
    stats.dump_stats("profile_output.prof")

[elindgren]

I also benchmarked a more thorough example, which actually performs a fit. The tabulated values helps a lot, the fit takes 42s on develop but only ~8s on the feature branch. On develop, a whopping 35s is taken up by cryspy computing the symmetry operations. On the feature branch, symmetry operations still take up 2s of the runtime, because of sympy.

Profiling results on develop

Profiling results on the feature branch

Benchmarking script

import cProfile
import pstats
from easydiffraction import (
    Project,
    SampleModel,
    Experiment
)

def prebench():
    # Create and configure sample model
    model = SampleModel('hs')
    model.space_group.name_h_m = 'R -3 m'
    model.space_group.it_coordinate_system_code = 'h'
    model.cell.length_a = 6.9
    model.cell.length_c = 14.1
    model.atom_sites.add('Zn', 'Zn', 0, 0, 0.5, wyckoff_letter='b', b_iso=0.5)
    model.atom_sites.add('Cu', 'Cu', 0.5, 0, 0, wyckoff_letter='e', b_iso=0.5)
    model.atom_sites.add('O', 'O', 0.21, -0.21, 0.06, wyckoff_letter='h', b_iso=0.5)
    model.atom_sites.add('Cl', 'Cl', 0, 0, 0.197, wyckoff_letter='c', b_iso=0.5)
    model.atom_sites.add('H', '2H', 0.13, -0.13, 0.08, wyckoff_letter='h', b_iso=0.5)
    model.show_as_cif()
    model.apply_symmetry_constraints()
    model.show_as_cif()

    # Create and configure experiment
    expt = Experiment('hrpt', data_path='examples/data/hrpt_hs.xye')
    expt.instrument.setup_wavelength = 1.89
    expt.instrument.calib_twotheta_offset = 0.0
    expt.peak.broad_gauss_u = 0.1
    expt.peak.broad_gauss_v = -0.2
    expt.peak.broad_gauss_w = 0.2
    expt.peak.broad_lorentz_x = 0.0
    expt.peak.broad_lorentz_y = 0
    expt.linked_phases.add('hs', scale=0.5)
    expt.background.add(x=4.4196, y=500)
    expt.background.add(x=6.6207, y=500)
    expt.background.add(x=10.4918, y=500)
    expt.background.add(x=15.4634, y=500)
    expt.background.add(x=45.6041, y=500)
    expt.background.add(x=74.6844, y=500)
    expt.background.add(x=103.4187, y=500)
    expt.background.add(x=121.6311, y=500)
    expt.background.add(x=159.4116, y=500)
    expt.show_as_cif()

    # Create project and add sample model and experiments
    project = Project()
    project.sample_models.add(model)
    project.experiments.add(expt)

    # Set calculator, minimizer and fit mode
    project.analysis.current_calculator = 'cryspy'
    project.analysis.current_minimizer = 'lmfit (leastsq)'

    # Compare measured and calculated patterns
    project.analysis.show_meas_vs_calc_chart('hrpt', 48, 51)

    # ------------ 1st fitting ------------

    # Define free parameters
    model.cell.length_a.free = True
    model.cell.length_c.free = True
    expt.linked_phases['hs'].scale.free = True
    expt.instrument.calib_twotheta_offset.free = True
    project.analysis.show_free_params()
    return project

def main(project):
    # Start fitting
    project.analysis.fit()

if __name__ == "__main__":
    project = prebench()
    with cProfile.Profile() as pr:
        main(project)
    stats = pstats.Stats(pr)
    stats.dump_stats("profile_output.prof")