diff --git a/.github/workflows/python-package-mac.yml b/.github/workflows/python-package-mac.yml index a951e1b..ea390a6 100644 --- a/.github/workflows/python-package-mac.yml +++ b/.github/workflows/python-package-mac.yml @@ -12,8 +12,8 @@ jobs: runs-on: macos-latest strategy: matrix: - python-version: ["3.9", "3.10", "3.11", "3.12"] - xcode-version: [latest-stable, 15] + python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"] + xcode-version: [latest-stable] steps: - uses: actions/checkout@v3 diff --git a/.github/workflows/python-package-ubuntu.yml b/.github/workflows/python-package-ubuntu.yml index 52ebe41..5bad931 100644 --- a/.github/workflows/python-package-ubuntu.yml +++ b/.github/workflows/python-package-ubuntu.yml @@ -16,7 +16,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - python-version: ["3.9", "3.10", "3.11", "3.12"] + python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"] steps: - uses: actions/checkout@v3 diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml index e6de737..35b281a 100644 --- a/.github/workflows/release.yml +++ b/.github/workflows/release.yml @@ -122,7 +122,7 @@ jobs: CIBW_BUILD_FRONTEND_WINDOWS: "build; args: -C setup-args=-Duse_openmp=disabled" #CIBW_BEFORE_BUILD_WINDOWS: '"C:\Program Files\Microsoft Visual Studio\2022\Enterprise\Common7\Tools\VsDevCmd.bat"' CIBW_BUILD_VERBOSITY: 1 - CIBW_BUILD: cp39-* cp310-* cp311-* cp312-* + CIBW_BUILD: cp310-* cp311-* cp312-* cp313-* cp314-* # Do not build for pypy and muslinux CIBW_SKIP: pp* *-musllinux_* CIBW_ARCHS: ${{ matrix.cibw_archs }} diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index a375ce4..86305b8 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -1,6 +1,11 @@ Change-log ########## +FreeSAS v2026.2.0 06/02/2026 +============================ +- Fix some tests +- Supports Python 3.10-3.14 + FreeSAS v2024.9.0 11/09/2021 ============================ - Rg+Dmax inferance from a Dense Neural Network (free_dnn) @@ -34,4 +39,4 @@ FreeSAS v0.6.0: 17/10/2017 FreeSAS v0.4.0: 17/12/2015 ========================== - First public version - - built around bead-model alignment + - built around bead-model alignment diff --git a/doc/source/installation.rst b/doc/source/installation.rst index 8361847..b46f2f8 100644 --- a/doc/source/installation.rst +++ b/doc/source/installation.rst @@ -3,6 +3,6 @@ Installation As most Python packages:: - pip install https://github.com/kif/freesas/archive/master.zip + pip install https://github.com/silx-kit/freesas/archive/main.zip -This requires Python3.6+ and all packages described in `requirements.txt`. \ No newline at end of file +This requires Python3.10+ and all packages described in `requirements.txt`. \ No newline at end of file diff --git a/doc/source/test.rst b/doc/source/test.rst index 40a805d..2f61d8b 100644 --- a/doc/source/test.rst +++ b/doc/source/test.rst @@ -3,7 +3,6 @@ Tests FreeSAS contains a set of non-regression tests:: - python3 setup.py build python3 run-test.py A connection to internet is needed to download the test material. diff --git a/src/freesas/__init__.py b/src/freesas/__init__.py index 489625b..80b7c1f 100644 --- a/src/freesas/__init__.py +++ b/src/freesas/__init__.py @@ -22,13 +22,11 @@ # THE SOFTWARE. # # ###########################################################################*/ -""" - -""" +""" """ __authors__ = ["Jérôme Kieffer"] __license__ = "MIT" -__date__ = "04/12/2023" +__date__ = "06/02/2026" import os as _os import logging as _logging @@ -40,16 +38,15 @@ try: from .version import __date__ as date # noqa - from .version import ( + from .version import ( # noqa version, version_info, hexversion, strictversion, dated_version, - citation - ) # noqa + citation, + ) except ImportError: raise RuntimeError( - "Do NOT use %s from its sources: build it and use the built version" - % project + "Do NOT use %s from its sources: build it and use the built version" % project ) diff --git a/src/freesas/align.py b/src/freesas/align.py index c4fec7d..84aa075 100644 --- a/src/freesas/align.py +++ b/src/freesas/align.py @@ -6,12 +6,14 @@ import sys import numpy import matplotlib + # matplotlib.use('Agg') import matplotlib.pyplot as plt from freesas.model import SASModel import itertools from scipy.optimize import fmin import logging + logging.basicConfig(level=logging.INFO) logger = logging.getLogger("log_freesas") @@ -47,7 +49,10 @@ def assign_models(self, molecule=None): model.canonical_parameters() self.sasmodels.append(model) if len(self.inputfiles) != len(self.sasmodels): - logger.error("Problem of assignment\n%s models for %s files" % (len(self.sasmodels), len(self.inputfiles))) + logger.error( + "Problem of assignment\n%s models for %s files" + % (len(self.sasmodels), len(self.inputfiles)) + ) elif len(molecule) != 0: model = SASModel() @@ -61,11 +66,11 @@ def assign_models(self, molecule=None): def rcalculation(self): """ - Calculation the maximal value for the R-factors, which is the mean of all the R-factors of + Calculation the maximal value for the R-factors, which is the mean of all the R-factors of inputs plus 2 times the standard deviation. R-factors are saved in the attribute self.rfactors, 1d array, and in percentage. - :return rmax: maximal value for the R-factor + :return rmax: maximal value for the R-factor """ if len(self.sasmodels) == 0: self.assign_models() @@ -123,12 +128,20 @@ def rfactorplot(self, filename=None, save=False): xticks = 1 + numpy.arange(dammif_files) fig = plt.figure(figsize=(7.5, 10)) - labels = [os.path.splitext(os.path.basename(self.inputfiles[i]))[0] for i in range(dammif_files)] + labels = [ + os.path.splitext(os.path.basename(self.inputfiles[i]))[0] + for i in range(dammif_files) + ] ax2 = fig.add_subplot(1, 1, 1) ax2.set_title("Selection of dammif models based on R factor") ax2.bar(xticks - 0.5, R) - ax2.plot([0.5, dammif_files + 0.5], [Rmax, Rmax], "-r", label="R$_{max}$ = %.3f" % Rmax) + ax2.plot( + [0.5, dammif_files + 0.5], + [Rmax, Rmax], + "-r", + label="R$_{max}$ = %.3f" % Rmax, + ) ax2.set_ylabel("R factor in percent") ax2.set_xticks(xticks) ax2.set_xticklabels(labels, rotation=90) @@ -137,7 +150,16 @@ def rfactorplot(self, filename=None, save=False): bbox_props = dict(fc="pink", ec="r", lw=1) for i in range(dammif_files): if not self.validmodels[i]: - ax2.text(i + 0.95, Rmax / 2, "Discarded", ha="center", va="center", rotation=90, size=10, bbox=bbox_props) + ax2.text( + i + 0.95, + Rmax / 2, + "Discarded", + ha="center", + va="center", + rotation=90, + size=10, + bbox=bbox_props, + ) logger.info("model %s discarded, Rfactor > Rmax" % self.inputfiles[i]) if save: @@ -186,7 +208,10 @@ def assign_models(self): model.canonical_parameters() self.models.append(model) if len(self.inputfiles) != len(self.models): - logger.error("Problem of assignment\n%s models for %s files" % (len(self.models), len(self.inputfiles))) + logger.error( + "Problem of assignment\n%s models for %s files" + % (len(self.models), len(self.inputfiles)) + ) return self.models @@ -195,12 +220,20 @@ def optimize(self, reference, molecule, symmetry): Use scipy.optimize to optimize transformation parameters to minimize NSD :param reference: SASmodel - :param molecule: SASmodel + :param molecule: SASmodel :param symmetry: 3-list of +/-1 :return p: transformation parameters optimized :return dist: NSD after optimization """ - p, dist, niter, nfuncalls, warmflag = fmin(reference.dist_after_movement, molecule.can_param, args=(molecule, symmetry), ftol=1e-4, maxiter=200, full_output=True, disp=False) + p, dist, niter, nfuncalls, warmflag = fmin( + reference.dist_after_movement, + molecule.can_param, + args=(molecule, symmetry), + ftol=1e-4, + maxiter=200, + full_output=True, + disp=False, + ) if niter == 200: logger.debug("convergence not reached") else: @@ -307,12 +340,17 @@ def plotNSDarray(self, rmax=None, filename=None, save=False): dammif_files = len(self.inputfiles) valid_models = self.validmodels - labels = [os.path.splitext(os.path.basename(self.outputfiles[i]))[0] for i in range(dammif_files)] - mask2d = (numpy.outer(valid_models, valid_models)) + labels = [ + os.path.splitext(os.path.basename(self.outputfiles[i]))[0] + for i in range(dammif_files) + ] + mask2d = numpy.outer(valid_models, valid_models) tableNSD = self.arrayNSD * mask2d maskedNSD = numpy.ma.masked_array(tableNSD, mask=numpy.logical_not(mask2d)) - data = valid_models * (tableNSD.sum(axis=-1) / (valid_models.sum() - 1)) # mean for the valid models, excluding itself - + data = valid_models * ( + tableNSD.sum(axis=-1) / (valid_models.sum() - 1) + ) # mean for the valid models, excluding itself + fig = plt.figure(figsize=(15, 10)) xticks = 1 + numpy.arange(dammif_files) ax1 = fig.add_subplot(1, 2, 1) @@ -324,16 +362,36 @@ def plotNSDarray(self, rmax=None, filename=None, save=False): for j in range(dammif_files): nsd = maskedNSD[i, j] if not maskedNSD.mask[i, j]: - ax1.text(i, j, "%.2f" % nsd, ha="center", va="center", size=12 * 8 // dammif_files) - ax1.text(j, i, "%.2f" % nsd, ha="center", va="center", size=12 * 8 // dammif_files) + ax1.text( + i, + j, + "%.2f" % nsd, + ha="center", + va="center", + size=12 * 8 // dammif_files, + ) + ax1.text( + j, + i, + "%.2f" % nsd, + ha="center", + va="center", + size=12 * 8 // dammif_files, + ) if i != j: lnsd.append(nsd) lnsd = numpy.array(lnsd) nsd_max = lnsd.mean() + lnsd.std() # threshold for nsd mean - ax1.imshow(maskedNSD, interpolation="nearest", origin="upper", cmap="YlOrRd", norm=matplotlib.colors.Normalize(vmin=min(lnsd))) - ax1.set_title(u"NSD correlation table") + ax1.imshow( + maskedNSD, + interpolation="nearest", + origin="upper", + cmap="YlOrRd", + norm=matplotlib.colors.Normalize(vmin=min(lnsd)), + ) + ax1.set_title("NSD correlation table") ax1.set_xticks(range(dammif_files)) ax1.set_xticklabels(labels, rotation=90) ax1.set_xlim(-0.5, dammif_files - 0.5) @@ -343,26 +401,68 @@ def plotNSDarray(self, rmax=None, filename=None, save=False): # second subplot : the NSD mean for each model ax2.bar(xticks - 0.5, data) - ax2.plot([0.5, dammif_files + 0.5], [nsd_max, nsd_max], "-r", label=u"NSD$_{max}$ = %.2f" % nsd_max) - ax2.set_title(u"NSD between any model and all others") + ax2.plot( + [0.5, dammif_files + 0.5], + [nsd_max, nsd_max], + "-r", + label="NSD$_{max}$ = %.2f" % nsd_max, + ) + ax2.set_title("NSD between any model and all others") ax2.set_ylabel("Normalized Spatial Discrepancy") ax2.set_xticks(xticks) ax2.set_xticklabels(labels, rotation=90) bbox_props = dict(fc="cyan", ec="b", lw=1) - ax2.text(self.reference + 0.95, data[self.reference] / 2, "Reference", ha="center", va="center", rotation=90, size=10, bbox=bbox_props) + ax2.text( + self.reference + 0.95, + data[self.reference] / 2, + "Reference", + ha="center", + va="center", + rotation=90, + size=10, + bbox=bbox_props, + ) ax2.legend(loc=8) bbox_props = dict(fc="pink", ec="r", lw=1) valid_number = 0 for i in range(dammif_files): if data[i] > nsd_max: - ax2.text(i + 0.95, data[self.reference] / 2, "Discarded", ha="center", va="center", rotation=90, size=10, bbox=bbox_props) + ax2.text( + i + 0.95, + data[self.reference] / 2, + "Discarded", + ha="center", + va="center", + rotation=90, + size=10, + bbox=bbox_props, + ) logger.debug("model %s discarded, nsd > nsd_max" % self.inputfiles[i]) elif not valid_models[i]: if rmax: - ax2.text(i + 0.95, data[self.reference] / 2, "Discarded, Rfactor = %s > Rmax = %s" % (100.0 * self.models[i].rfactor, rmax), ha="center", va="center", rotation=90, size=10, bbox=bbox_props) + ax2.text( + i + 0.95, + data[self.reference] / 2, + "Discarded, Rfactor = %s > Rmax = %s" + % (100.0 * self.models[i].rfactor, rmax), + ha="center", + va="center", + rotation=90, + size=10, + bbox=bbox_props, + ) else: - ax2.text(i + 0.95, data[self.reference] / 2, "Discarded", ha="center", va="center", rotation=90, size=10, bbox=bbox_props) + ax2.text( + i + 0.95, + data[self.reference] / 2, + "Discarded", + ha="center", + va="center", + rotation=90, + size=10, + bbox=bbox_props, + ) else: if valid_models[i] == 1.0: valid_number += 1 @@ -379,7 +479,7 @@ def find_reference(self): """ Find the reference model among the models aligned. The reference model is the one with lower average NSD with other models. - + :return ref_number: position of the reference model in the list self.models """ if self.arrayNSD is None: @@ -404,7 +504,7 @@ def alignment_reference(self, ref_number=None): """ if self.reference is None and ref_number is None: self.find_reference() - + ref_number = self.reference models = self.models reference = models[ref_number] @@ -416,8 +516,12 @@ def alignment_reference(self, ref_number=None): symmetry, p = self.alignment_sym(reference, molecule) if not self.slow: p, dist = self.optimize(reference, molecule, symmetry) - molecule.atoms = molecule.transform(p, symmetry) # molecule sent on its canonical position - molecule.atoms = molecule.transform(reference.can_param, [1, 1, 1], reverse=True) # molecule sent on reference position + molecule.atoms = molecule.transform( + p, symmetry + ) # molecule sent on its canonical position + molecule.atoms = molecule.transform( + reference.can_param, [1, 1, 1], reverse=True + ) # molecule sent on reference position molecule.save(self.outputfiles[i]) reference.save(self.outputfiles[ref_number]) return 0 @@ -439,7 +543,9 @@ def alignment_2models(self, save=True): p, dist = self.optimize(reference, molecule, symmetry) molecule.atoms = molecule.transform(p, symmetry) - molecule.atoms = molecule.transform(reference.can_param, [1, 1, 1], reverse=True) + molecule.atoms = molecule.transform( + reference.can_param, [1, 1, 1], reverse=True + ) if self.slow: dist = reference.dist(molecule, reference.atoms, molecule.atoms) if save: diff --git a/src/freesas/app/auto_gpa.py b/src/freesas/app/auto_gpa.py index fe02a90..666afc9 100644 --- a/src/freesas/app/auto_gpa.py +++ b/src/freesas/app/auto_gpa.py @@ -54,9 +54,7 @@ def build_parser() -> GuinierParser: the autorg algorithm originately part of the ATSAS suite. As this tool used a different theory, some results may differ """ - return GuinierParser( - prog="free_gpa", description=description, epilog=epilog - ) + return GuinierParser(prog="free_gpa", description=description, epilog=epilog) def main() -> None: diff --git a/src/freesas/app/auto_guinier.py b/src/freesas/app/auto_guinier.py index 8f918cc..b435c97 100644 --- a/src/freesas/app/auto_guinier.py +++ b/src/freesas/app/auto_guinier.py @@ -54,9 +54,7 @@ def build_parser() -> GuinierParser: the autorg algorithm originately part of the ATSAS suite. As this tool used a different theory, some results may differ """ - return GuinierParser( - prog="free_guinier", description=description, epilog=epilog - ) + return GuinierParser(prog="free_guinier", description=description, epilog=epilog) def main() -> None: diff --git a/src/freesas/app/autorg.py b/src/freesas/app/autorg.py index 2e9ef0c..6190b74 100644 --- a/src/freesas/app/autorg.py +++ b/src/freesas/app/autorg.py @@ -54,9 +54,7 @@ def build_parser() -> GuinierParser: the autorg algorithm originately part of the ATSAS suite. As this is reverse engineered, some constants and results may differ """ - return GuinierParser( - prog="free_rg", description=description, epilog=epilog - ) + return GuinierParser(prog="free_rg", description=description, epilog=epilog) def main() -> None: diff --git a/src/freesas/app/bift.py b/src/freesas/app/bift.py index 2a540b0..a77e00f 100644 --- a/src/freesas/app/bift.py +++ b/src/freesas/app/bift.py @@ -26,8 +26,8 @@ __author__ = "Jérôme Kieffer" __license__ = "MIT" -__copyright__ = "2017, ESRF" -__date__ = "13/10/2020" +__copyright__ = "2017-2026, ESRF" +__date__ = "06/02/2026" import sys import logging @@ -68,9 +68,7 @@ def build_parser() -> SASParser: It aims at being a drop in replacement for datgnom of the ATSAS suite. """ - parser = SASParser( - prog="free_bift", description=description, epilog=epilog - ) + parser = SASParser(prog="free_bift", description=description, epilog=epilog) parser.add_file_argument(help_text="I(q) files to convert into p(r)") parser.add_output_filename_argument() parser.add_q_unit_argument() @@ -118,7 +116,7 @@ def main(): for afile in files: try: data = load_scattering_data(afile) - except: + except Exception: logger.error("Unable to parse file %s", afile) else: if args.unit == "Å": diff --git a/src/freesas/app/cormap.py b/src/freesas/app/cormap.py index abbcbd7..82e5a7f 100644 --- a/src/freesas/app/cormap.py +++ b/src/freesas/app/cormap.py @@ -4,29 +4,27 @@ __author__ = "Jérôme Kieffer" __license__ = "MIT" __copyright__ = "2015, ESRF" -__date__ = "12/07/2024" +__date__ = "06/02/2026" import os import logging import glob +import platform from itertools import combinations from collections import namedtuple -import numpy -logging.basicConfig(level=logging.INFO) -logger = logging.getLogger("cormap") -import freesas from freesas.cormap import gof from freesas.sasio import load_scattering_data from freesas.sas_argparser import SASParser +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger("cormap") + datum = namedtuple("datum", ["index", "filename", "data"]) -import platform operatingSystem = platform.system() - def parse(): """Parse input and return list of files. :return: list of input files @@ -55,7 +53,7 @@ def parse(): def compare(lstfiles): res = [ "Pair-wise Correlation Map", - "" " C Pr(>C)", + " C Pr(>C)", ] data = [] for i, f in enumerate(lstfiles): @@ -69,14 +67,10 @@ def compare(lstfiles): data.append(d) for a, b in combinations(data, 2): r = gof(a.data, b.data) - res.append( - "%6i vs. %6i %6i %8.6f" % (a.index, b.index, r.c, r.P) - ) + res.append("%6i vs. %6i %6i %8.6f" % (a.index, b.index, r.c, r.P)) res.append("") for a in data: - res.append( - "%6i %8f + %8f * %s" % (a.index, 0.0, 1.0, a.filename) - ) + res.append("%6i %8f + %8f * %s" % (a.index, 0.0, 1.0, a.filename)) res.append("") print(os.linesep.join(res)) return res diff --git a/src/freesas/app/dnn.py b/src/freesas/app/dnn.py index c99d942..fb58cf7 100644 --- a/src/freesas/app/dnn.py +++ b/src/freesas/app/dnn.py @@ -63,7 +63,6 @@ def build_parser() -> SASParser: return parser - def main() -> None: """Entry point for free_gpa app""" parser = build_parser() diff --git a/src/freesas/app/extract_ascii.py b/src/freesas/app/extract_ascii.py index 2ffba3b..d859e76 100644 --- a/src/freesas/app/extract_ascii.py +++ b/src/freesas/app/extract_ascii.py @@ -26,8 +26,8 @@ __author__ = "Jérôme Kieffer" __license__ = "MIT" -__copyright__ = "2020-2024, ESRF" -__date__ = "22/04/2025" +__copyright__ = "2020-2026, ESRF" +__date__ = "06/02/2026" import io import os @@ -57,16 +57,13 @@ def parse(): - """Parse input and return list of files. :return: list of input files """ description = "Extract the SAXS data from a Nexus files as a 3 column ascii (q, I, err). Metadata are exported in the headers as needed." epilog = """extract_ascii.py allows you to export the data in inverse nm or inverse A with possible intensity scaling. """ - parser = SASParser( - prog="extract-ascii.py", description=description, epilog=epilog - ) + parser = SASParser(prog="extract-ascii.py", description=description, epilog=epilog) # Commented option need to be implemented # parser.add_argument("-o", "--output", action='store', help="Output filename, by default the same with .dat extension", default=None, type=str) # parser.add_argument("-u", "--unit", action='store', help="Unit for q: inverse nm or Angstrom?", default="nm", type=str) @@ -90,7 +87,7 @@ def parse(): def extract_averaged(filename): - "return some infomations extracted from a HDF5 file " + "return some infomations extracted from a HDF5 file" results = OrderedDict() results["filename"] = filename # Missing: comment normalization @@ -105,13 +102,15 @@ def extract_averaged(filename): # if program == "" default = entry_grp.attrs["default"] if posixpath.split(default)[-1] == "hplc": - default = posixpath.join(posixpath.split(default)[0],"results") + default = posixpath.join(posixpath.split(default)[0], "results") # print(default) nxdata_grp = nxsr.h5[default] signal = nxdata_grp.attrs["signal"] axis = nxdata_grp.attrs["axes"] - if not isinstance(axis, (str,bytes)): - logger.error(f"There are several curves if the dataset {default} from file {filename}, please use the option --all to extract them all") + if not isinstance(axis, (str, bytes)): + logger.error( + f"There are several curves if the dataset {default} from file {filename}, please use the option --all to extract them all" + ) sys.exit(1) results["I"] = nxdata_grp[signal][()] results["q"] = nxdata_grp[axis][()] @@ -120,15 +119,16 @@ def extract_averaged(filename): axis + "_" + nxdata_grp[axis].attrs["units"] ) integration_grp = nxdata_grp.parent - results["geometry"] = json.loads( - integration_grp["configuration/data"][()] - ) - results["polarization"] = integration_grp["configuration/polarization_factor"][()] + results["geometry"] = json.loads(integration_grp["configuration/data"][()]) + results["polarization"] = integration_grp["configuration/polarization_factor"][ + () + ] instrument_grps = nxsr.get_class(entry_grp, class_type="NXinstrument") if instrument_grps: - detector_grp = nxsr.get_class(instrument_grps[0], - class_type="NXdetector")[0] + detector_grp = nxsr.get_class(instrument_grps[0], class_type="NXdetector")[ + 0 + ] results["mask"] = detector_grp["pixel_mask"].attrs["filename"] sample_grp = nxsr.get_class(entry_grp, class_type="NXsample")[0] results["sample"] = posixpath.split(sample_grp.name)[-1] @@ -137,12 +137,14 @@ def extract_averaged(filename): results["exposure temperature"] = sample_grp["temperature"][()] results["concentration"] = sample_grp["concentration"][()] if "2_correlation_mapping" in entry_grp: - results["to_merge"] = entry_grp["2_correlation_mapping/results/to_merge"][()] + results["to_merge"] = entry_grp["2_correlation_mapping/results/to_merge"][ + () + ] return results def extract_all(filename): - """return some infomations extracted from a HDF5 file for all individual frames. + """return some infomations extracted from a HDF5 file for all individual frames. Supports HPLC and SC and freshly integrated blocks of frames""" res = [] results = OrderedDict() @@ -153,7 +155,11 @@ def extract_all(filename): entry_grp = nxsr.h5[default] else: entry_grp = nxsr.get_entries()[0] - program = entry_grp.get("program_name")[()].decode() if "program_name" in entry_grp else None + program = ( + entry_grp.get("program_name")[()].decode() + if "program_name" in entry_grp + else None + ) if program == "bm29.hplc": target = "1_chromatogram" @@ -168,11 +174,13 @@ def extract_all(filename): nxdata_grp = nxsr.h5[target] signal = nxdata_grp.attrs["signal"] axis = nxdata_grp.attrs["axes"][1] - I = nxdata_grp[signal][()] + intensity = nxdata_grp[signal][()] results["q"] = nxdata_grp[axis][()] std = nxdata_grp["errors"][()] try: - results["unit"] = pyFAI.units.to_unit(axis + "_" + nxdata_grp[axis].attrs["units"]) + results["unit"] = pyFAI.units.to_unit( + axis + "_" + nxdata_grp[axis].attrs["units"] + ) except KeyError: logger.warning("Unable to parse radial units") integration_grp = nxdata_grp.parent @@ -181,7 +189,9 @@ def extract_all(filename): else: logger.warning("Unable to parse AzimuthalIntegrator configuration") if "configuration/polarization_factor" in integration_grp: - results["polarization"] = integration_grp["configuration/polarization_factor"][()] + results["polarization"] = integration_grp[ + "configuration/polarization_factor" + ][()] instrument_grp = nxsr.get_class(entry_grp, class_type="NXinstrument") if instrument_grp: detector_grp = nxsr.get_class(instrument_grp[0], class_type="NXdetector") @@ -201,7 +211,7 @@ def extract_all(filename): results["concentration"] = sample_grp["concentration"][()] # if "2_correlation_mapping" in entry_grp: # results["to_merge"] = entry_grp["2_correlation_mapping/results/to_merge"][()] - for i, s in zip(I, std): + for i, s in zip(intensity, std): r = copy.copy(results) r["I"] = i r["std"] = s @@ -257,28 +267,18 @@ def write_ascii(results, output=None, hdr="#", linesep=os.linesep): headers.append(hdr + " " + results["comments"]) else: headers.append(hdr) - headers.append( - hdr + " Sample c= %s mg/ml" % results.get("concentration", -1) - ) + headers.append(hdr + " Sample c= %s mg/ml" % results.get("concentration", -1)) headers += [hdr, hdr + " Sample environment:"] if "geometry" in results: - headers.append( - hdr + " Detector = %s" % results["geometry"]["detector"] - ) - headers.append( - hdr + " SampleDistance = %s" % results["geometry"]["dist"] - ) - headers.append( - hdr + " WaveLength = %s" % results["geometry"]["wavelength"] - ) + headers.append(hdr + " Detector = %s" % results["geometry"]["detector"]) + headers.append(hdr + " SampleDistance = %s" % results["geometry"]["dist"]) + headers.append(hdr + " WaveLength = %s" % results["geometry"]["wavelength"]) headers.append(hdr) if "comments" in results: headers.append(hdr + " title = %s" % results["comment"]) if "to_merge" in results: headers.append( - hdr - + " Frames merged: " - + " ".join([str(i) for i in results["to_merge"]]) + hdr + " Frames merged: " + " ".join([str(i) for i in results["to_merge"]]) ) if "normalization" in results: headers.append(hdr + " Normalization = %s" % results["normalization"]) @@ -302,8 +302,7 @@ def write_ascii(results, output=None, hdr="#", linesep=os.linesep): if "storage temperature" in results: headers.append( hdr - + " Storage Temperature (degrees C): %s" - % results["storage temperature"] + + " Storage Temperature (degrees C): %s" % results["storage temperature"] ) if "exposure temperature" in results: headers.append( @@ -312,9 +311,7 @@ def write_ascii(results, output=None, hdr="#", linesep=os.linesep): % results["exposure temperature"] ) - headers.append( - hdr + " Concentration: %s" % results.get("concentration", -1) - ) + headers.append(hdr + " Concentration: %s" % results.get("concentration", -1)) if "buffer" in results: headers.append(hdr + " Buffer: %s" % results["buffer"]) headers.append(hdr + " Code: %s" % results.get("sample", "")) @@ -326,15 +323,13 @@ def write(headers, file_): if "std" in results: data = [ - "%14.6e\t%14.6e\t%14.6e" % (q, I, std) - for q, I, std in zip( - results["q"], results["I"], results["std"] - ) + "%14.6e\t%14.6e\t%14.6e" % (q, intensity, std) + for q, intensity, std in zip(results["q"], results["I"], results["std"]) ] else: data = [ - "%14.6e\t%14.6e\t" % (q, I) - for q, I in zip(results["q"], results["I"]) + "%14.6e\t%14.6e\t" % (q, intensity) + for q, intensity in zip(results["q"], results["I"]) ] data.append("") file_.writelines(linesep.join(data)) diff --git a/src/freesas/app/plot_sas.py b/src/freesas/app/plot_sas.py index 6378b30..9b5bb62 100644 --- a/src/freesas/app/plot_sas.py +++ b/src/freesas/app/plot_sas.py @@ -76,9 +76,7 @@ def parse(): description = "Generate typical sas plots with matplotlib" epilog = """freesas is an open-source implementation of a bunch of small angle scattering algorithms. """ - parser = SASParser( - prog="freesas.py", description=description, epilog=epilog - ) + parser = SASParser(prog="freesas.py", description=description, epilog=epilog) parser.add_file_argument(help_text="dat files to plot") parser.add_output_filename_argument() parser.add_output_data_format("jpeg", "svg", "png", "pdf") diff --git a/src/freesas/app/supycomb.py b/src/freesas/app/supycomb.py index 2c4f688..53e8343 100644 --- a/src/freesas/app/supycomb.py +++ b/src/freesas/app/supycomb.py @@ -16,7 +16,6 @@ def parse(): - """Parse input and return list of files. :return: list of args """ @@ -113,9 +112,7 @@ def main(): logger.info("%s and %s aligned" % (args.file[0], args.file[1])) logger.info("NSD after optimized alignment = %.2f" % dist) else: - align.outputfiles = [ - "model-%02i.pdb" % (i + 1) for i in range(input_len) - ] + align.outputfiles = ["model-%02i.pdb" % (i + 1) for i in range(input_len)] selection.inputfiles = args.file selection.models_selection() selection.rfactorplot(save=save) @@ -124,9 +121,7 @@ def main(): align.makeNSDarray() align.alignment_reference() - logger.info( - "valid models aligned on the model %s" % (align.reference + 1) - ) + logger.info("valid models aligned on the model %s" % (align.reference + 1)) align.plotNSDarray(rmax=round(selection.rmax, 4), save=save) if not save and input_len > 2: diff --git a/src/freesas/autorg.py b/src/freesas/autorg.py index 4284f94..3255462 100644 --- a/src/freesas/autorg.py +++ b/src/freesas/autorg.py @@ -1,30 +1,36 @@ # -*- coding: utf-8 -*- """Functions for calculating the radius of gyration and forward scattering intensity.""" -__authors__ = ["Jerome Kieffer"] +__authors__ = ["Jérôme Kieffer"] __license__ = "MIT" -__copyright__ = "2020, ESRF" -__date__ = "05/06/2020" +__copyright__ = "2020-2026, ESRF" +__date__ = "06/02/2026" import logging import numpy from scipy.optimize import curve_fit -from ._autorg import ( # pylint: disable=E0401 +from ._autorg import ( # noqa RG_RESULT, - autoRg, - AutoGuinier, - linear_fit, - FIT_RESULT, guinier, NoGuinierRegionError, DTYPE, InsufficientDataError, + autoRg, + AutoGuinier, + linear_fit, + FIT_RESULT, ) logger = logging.getLogger(__name__) +def _gpa(q2, Rg, I0): + """Function to be fitted""" + x_prime = q2 * Rg * Rg + return I0 / Rg * numpy.sqrt(x_prime) * numpy.exp(-x_prime / 3.0) + + def auto_gpa(data, Rg_min=1.0, qRg_max=1.3, qRg_min=0.5): """ Uses the GPA theory to guess quickly Rg, the @@ -48,28 +54,26 @@ def auto_gpa(data, Rg_min=1.0, qRg_max=1.3, qRg_min=0.5): """ def curate_data(data): - q = data.T[0] - I = data.T[1] - err = data.T[2] + q, intensity, err = data.T[:3] - start0 = numpy.argmax(I) + start0 = numpy.argmax(intensity) stop0 = numpy.where(q > qRg_max / Rg_min)[0][0] range0 = slice(start0, stop0) q = q[range0] - I = I[range0] + intensity = intensity[range0] err = err[range0] - q2 = q ** 2 - lnI = numpy.log(I) - I2_over_sigma2 = err ** 2 / I ** 2 + q2 = q**2 + lnI = numpy.log(intensity) + I2_over_sigma2 = err**2 / intensity**2 - y = I * q + y = intensity * q p1 = numpy.argmax(y) # Those are guess from the max position: Rg = (1.5 / q2[p1]) ** 0.5 - I0 = I[p1] * numpy.exp(q2[p1] * Rg ** 2 / 3.0) + I0 = intensity[p1] * numpy.exp(q2[p1] * Rg**2 / 3.0) # Let's cut-down the guinier region from 0.5-1.3 in qRg try: @@ -83,30 +87,20 @@ def curate_data(data): range1 = slice(start1, stop1) q1 = q[range1] - I1 = I[range1] + I1 = intensity[range1] return q1, I1, Rg, I0, q2, lnI, I2_over_sigma2, start0 q1, I1, Rg, I0, q2, lnI, I2_over_sigma2, start0 = curate_data(data) if len(q1) < 3: reduced_data = numpy.delete(data, start0, axis=0) - q1, I1, Rg, I0, q2, lnI, I2_over_sigma2, start0 = curate_data( - reduced_data - ) + q1, I1, Rg, I0, q2, lnI, I2_over_sigma2, start0 = curate_data(reduced_data) x = q1 * q1 y = I1 * q1 - f = ( - lambda x, Rg, I0: I0 - / Rg - * numpy.sqrt(x * Rg * Rg) - * numpy.exp(-x * Rg * Rg / 3.0) - ) - res = curve_fit(f, x, y, [Rg, I0]) - logger.debug( - "GPA upgrade Rg %s-> %s and I0 %s -> %s", Rg, res[0][0], I0, res[0][1] - ) + res = curve_fit(_gpa, x, y, [Rg, I0]) + logger.debug("GPA upgrade Rg %s-> %s and I0 %s -> %s", Rg, res[0][0], I0, res[0][1]) Rg, I0 = res[0] sigma_Rg, sigma_I0 = numpy.sqrt(numpy.diag(res[1])) end = numpy.where(data.T[0] > qRg_max / Rg)[0][0] @@ -117,9 +111,7 @@ def curate_data(data): quality = guinier.calc_quality( Rg, sigma_Rg, data.T[0, start], data.T[0, end], aggregation, qRg_max ) - return RG_RESULT( - Rg, sigma_Rg, I0, sigma_I0, start, end, quality, aggregation - ) + return RG_RESULT(Rg, sigma_Rg, I0, sigma_I0, start, end, quality, aggregation) def auto_guinier(data, Rg_min=1.0, qRg_max=1.3, relax=1.2): @@ -160,9 +152,7 @@ def auto_guinier(data, Rg_min=1.0, qRg_max=1.3, relax=1.2): data, q_ary, i_ary, sigma_ary, Rg_min, qRg_max, relax ) if start0 < 0: - raise InsufficientDataError( - "Minimum region size is %s" % guinier.min_size - ) + raise InsufficientDataError("Minimum region size is %s" % guinier.min_size) guinier.guinier_space( start0, stop0, q_ary, i_ary, sigma_ary, q2_ary, lnI_ary, wg_ary ) @@ -171,9 +161,7 @@ def auto_guinier(data, Rg_min=1.0, qRg_max=1.3, relax=1.2): q2_ary, lnI_ary, wg_ary, start0, stop0, Rg_min, qRg_max, relax ) - cnt, relaxed, qRg_max, aslope_max = guinier.count_valid( - fits, qRg_max, relax - ) + cnt, relaxed, qRg_max, aslope_max = guinier.count_valid(fits, qRg_max, relax) # valid_fits = fits[fits[:, 9] < qRg_max] if cnt == 0: raise NoGuinierRegionError(qRg_max) @@ -182,9 +170,7 @@ def auto_guinier(data, Rg_min=1.0, qRg_max=1.3, relax=1.2): start, stop = guinier.find_region(fits, qRg_max) # Now average out the - Rg_avg, Rg_std, I0_avg, I0_std, good = guinier.average_values( - fits, start, stop - ) + Rg_avg, Rg_std, I0_avg, I0_std, good = guinier.average_values(fits, start, stop) aggregated = guinier.check_aggregation( q2_ary, lnI_ary, wg_ary, start0, stop, Rg=Rg_avg, threshold=False @@ -192,7 +178,5 @@ def auto_guinier(data, Rg_min=1.0, qRg_max=1.3, relax=1.2): quality = guinier.calc_quality( Rg_avg, Rg_std, q_ary[start], q_ary[stop], aggregated, qRg_max ) - result = RG_RESULT( - Rg_avg, Rg_std, I0_avg, I0_std, start, stop, quality, aggregated - ) + result = RG_RESULT(Rg_avg, Rg_std, I0_avg, I0_std, start, stop, quality, aggregated) return result diff --git a/src/freesas/average.py b/src/freesas/average.py index 3c428c4..32c0e2a 100644 --- a/src/freesas/average.py +++ b/src/freesas/average.py @@ -10,6 +10,7 @@ class Grid: """ This class is used to create a grid which include all the input models """ + def __init__(self, inputfiles): """ :param inputfiles: list of pdb files needed for averaging @@ -21,19 +22,19 @@ def __init__(self, inputfiles): self.coordknots = [] def __repr__(self): - return "Grid with %i knots"%self.nbknots + return "Grid with %i knots" % self.nbknots def spatial_extent(self): """ Calculate the maximal extent of input models - + :return self.size: 6-list with x,y,z max and then x,y,z min """ atoms = [] models_fineness = [] for files in self.inputs: m = SASModel(files) - if len(atoms)==0: + if len(atoms) == 0: atoms = m.atoms else: atoms = numpy.append(atoms, m.atoms, axis=0) @@ -42,19 +43,26 @@ def spatial_extent(self): coordmin = atoms.min(axis=0) - mean_fineness coordmax = atoms.max(axis=0) + mean_fineness - self.size = [coordmax[0],coordmax[1],coordmax[2],coordmin[0],coordmin[1],coordmin[2]] + self.size = [ + coordmax[0], + coordmax[1], + coordmax[2], + coordmin[0], + coordmin[1], + coordmin[2], + ] return self.size def calc_radius(self, nbknots=None): """ - Calculate the radius of each point of a hexagonal close-packed grid, + Calculate the radius of each point of a hexagonal close-packed grid, knowing the total volume and the number of knots in this grid. :param nbknots: number of knots wanted for the grid :return radius: the radius of each knot of the grid """ - if len(self.size)==0: + if len(self.size) == 0: self.spatial_extent() nbknots = nbknots if nbknots is not None else 5000 size = self.size @@ -63,8 +71,8 @@ def calc_radius(self, nbknots=None): dz = size[2] - size[5] volume = dx * dy * dz - density = numpy.pi / (3*2**0.5) - radius = ((3 /( 4 * numpy.pi)) * density * volume / nbknots)**(1.0/3) + density = numpy.pi / (3 * 2**0.5) + radius = ((3 / (4 * numpy.pi)) * density * volume / nbknots) ** (1.0 / 3) self.radius = radius return radius @@ -76,13 +84,13 @@ def make_grid(self): :return knots: 2d-array, coordinates of each dot of the grid. Saved as self.coordknots. """ - if len(self.size)==0: + if len(self.size) == 0: self.spatial_extent() if self.radius is None: self.calc_radius() radius = self.radius - a = numpy.sqrt(2.0)*radius + a = numpy.sqrt(2.0) * radius xmax = self.size[0] xmin = self.size[3] @@ -98,7 +106,7 @@ def make_grid(self): xlist = [] ylist = [] zlist = [] - knots = numpy.empty((1,4), dtype="float") + knots = numpy.empty((1, 4), dtype="float") while (zmin + z) <= zmax: zlist.append(z) z += a @@ -111,24 +119,32 @@ def make_grid(self): for i in range(len(zlist)): z = zlist[i] - if i % 2 ==0: + if i % 2 == 0: for j in range(len(xlist)): x = xlist[j] if j % 2 == 0: for y in ylist[0:-1:2]: - knots = numpy.append(knots, [[xmin+x, ymin+y, zmin+z, 0.0]], axis=0) + knots = numpy.append( + knots, [[xmin + x, ymin + y, zmin + z, 0.0]], axis=0 + ) else: for y in ylist[1:-1:2]: - knots = numpy.append(knots, [[xmin+x, ymin+y, zmin+z, 0.0]], axis=0) + knots = numpy.append( + knots, [[xmin + x, ymin + y, zmin + z, 0.0]], axis=0 + ) else: for j in range(len(xlist)): x = xlist[j] if j % 2 == 0: for y in ylist[1:-1:2]: - knots = numpy.append(knots, [[xmin+x, ymin+y, zmin+z, 0.0]], axis=0) + knots = numpy.append( + knots, [[xmin + x, ymin + y, zmin + z, 0.0]], axis=0 + ) else: for y in ylist[0:-1:2]: - knots = numpy.append(knots, [[xmin+x, ymin+y, zmin+z, 0.0]], axis=0) + knots = numpy.append( + knots, [[xmin + x, ymin + y, zmin + z, 0.0]], axis=0 + ) knots = numpy.delete(knots, 0, axis=0) self.nbknots = knots.shape[0] @@ -137,10 +153,11 @@ def make_grid(self): return knots -class AverModels(): +class AverModels: """ Provides tools to create an averaged models using several aligned dummy atom models """ + def __init__(self, inputfiles, grid): """ :param inputfiles: list of pdb files of aligned models @@ -154,7 +171,7 @@ def __init__(self, inputfiles, grid): self.grid = grid def __repr__(self): - return "Average SAS model with %i atoms"%len(self.atoms) + return "Average SAS model with %i atoms" % len(self.atoms) def read_files(self, reference=None): """ @@ -169,7 +186,7 @@ def read_files(self, reference=None): models = [] models.append(SASModel(inputfiles[ref])) for i in range(len(inputfiles)): - if i==ref: + if i == ref: continue else: models.append(SASModel(inputfiles[i])) @@ -203,7 +220,7 @@ def assign_occupancy(self): """ For each point of the grid, total occupancy and contribution factor are computed and saved. The grid is then ordered with decreasing value of occupancy. - The fourth column of the array correspond to the occupancy of the point and the fifth to + The fourth column of the array correspond to the occupancy of the point and the fifth to the contribution for this point. :return sortedgrid: 2d-array, coordinates of each point of the grid @@ -229,25 +246,25 @@ def make_header(self): Create the layout of the pdb file for the averaged model. """ header = [] - header.append("Number of files averaged : %s\n"%len(self.inputfiles)) + header.append("Number of files averaged : %s\n" % len(self.inputfiles)) for i in self.inputfiles: header.append(i + "\n") - header.append("Total number of dots in the grid : %s\n"%self.grid.shape[0]) + header.append("Total number of dots in the grid : %s\n" % self.grid.shape[0]) decade = 1 for i in range(self.grid.shape[0]): line = "ATOM CA ASP 1 20.00 2 201\n" - line = line[:7] + "%4.i"%(i + 1) + line[11:] + line = line[:7] + "%4.i" % (i + 1) + line[11:] if not (i + 1) % 10: decade += 1 - line = line[:21] + "%4.i"%decade + line[25:] + line = line[:21] + "%4.i" % decade + line[25:] header.append(line) self.header = header return header def save_aver(self, filename): """ - Save the position of each occupied dot of the grid, its occupancy and its contribution + Save the position of each occupied dot of the grid, its occupancy and its contribution in a pdb file. :param filename: name of the pdb file to write @@ -264,7 +281,9 @@ def save_aver(self, filename): coord = "%8.3f%8.3f%8.3f" % tuple(self.grid[nr, 0:3]) occ = "%6.2f" % self.grid[nr, 3] contrib = "%2.f" % self.grid[nr, 4] - line = line[:30] + coord + occ + line[60:66] + contrib + line[68:] + line = ( + line[:30] + coord + occ + line[60:66] + contrib + line[68:] + ) else: line = "" nr += 1 diff --git a/src/freesas/bift.py b/src/freesas/bift.py index 058914c..4e932f1 100644 --- a/src/freesas/bift.py +++ b/src/freesas/bift.py @@ -11,23 +11,31 @@ Many thanks to Pierre Paleo for the auto-alpha guess """ -__authors__ = ["Jerome Kieffer", "Jesse Hopkins"] +__authors__ = ["Jérôme Kieffer", "Jesse Hopkins"] __license__ = "MIT" -__copyright__ = "2020, ESRF" -__date__ = "04/12/2023" +__copyright__ = "2020-2026, ESRF" +__date__ = "06/02/2026" import logging -logger = logging.getLogger(__name__) -from math import log, ceil -import numpy +from math import log from scipy.optimize import minimize from ._bift import BIFT -from .autorg import auto_gpa, autoRg, auto_guinier, NoGuinierRegionError +from .autorg import auto_guinier, NoGuinierRegionError + + +logger = logging.getLogger(__name__) -def auto_bift(data, Dmax=None, alpha=None, npt=100, - start_point=None, end_point=None, - scan_size=11, Dmax_over_Rg=3): +def auto_bift( + data, + Dmax=None, + alpha=None, + npt=100, + start_point=None, + end_point=None, + scan_size=11, + Dmax_over_Rg=3, +): """Calculates the inverse Fourier tranform of the data using an optimisation of the evidence :param data: 2D array with q, I(q), δI(q). q can be in 1/nm or 1/A, it imposes the unit for r & Dmax @@ -44,36 +52,46 @@ def auto_bift(data, Dmax=None, alpha=None, npt=100, assert data.shape[1] == 3 # enforce q, I, err use_wisdom = False data = data[slice(start_point, end_point)] - q, I, err = data.T + q, intenisty, err = data.T npt = min(npt, q.size) # no chance for oversampling ! - bo = BIFT(q, I, err) # this is the bift object + bo = BIFT(q, intenisty, err) # this is the bift object if Dmax is None: # Try to get a reasonable guess from Rg try: Guinier = auto_guinier(data) - except: + except Exception: logger.error("Guinier analysis failed !") raise -# print(Guinier) + # print(Guinier) if Guinier.Rg <= 0: raise NoGuinierRegionError Dmax = bo.set_Guinier(Guinier, Dmax_over_Rg) if alpha is None: alpha_max = bo.guess_alpha_max(npt) # First scan on alpha: - key = bo.grid_scan(Dmax, Dmax, 1, - 1.0 / alpha_max, alpha_max, scan_size, npt) + key = bo.grid_scan(Dmax, Dmax, 1, 1.0 / alpha_max, alpha_max, scan_size, npt) Dmax, alpha = key[:2] # Then scan on Dmax: - key = bo.grid_scan(max(Dmax / 2, Dmax * (Dmax_over_Rg - 1) / Dmax_over_Rg), Dmax * (Dmax_over_Rg + 1) / Dmax_over_Rg, scan_size, - alpha, alpha, 1, npt) + key = bo.grid_scan( + max(Dmax / 2, Dmax * (Dmax_over_Rg - 1) / Dmax_over_Rg), + Dmax * (Dmax_over_Rg + 1) / Dmax_over_Rg, + scan_size, + alpha, + alpha, + 1, + npt, + ) Dmax, alpha = key[:2] if bo.evidence_cache[key].converged: bo.update_wisdom() use_wisdom = True # Optimization using Bayesian operator: - logger.info("Start search at Dmax=%.2f alpha=%.2f use wisdom=%s", Dmax, alpha, use_wisdom) - res = minimize(bo.opti_evidence, (Dmax, log(alpha)), args=(npt, use_wisdom), method="powell") + logger.info( + "Start search at Dmax=%.2f alpha=%.2f use wisdom=%s", Dmax, alpha, use_wisdom + ) + res = minimize( + bo.opti_evidence, (Dmax, log(alpha)), args=(npt, use_wisdom), method="powell" + ) logger.info("Result of optimisation:\n %s", res) return bo diff --git a/src/freesas/containers.py b/src/freesas/containers.py index 91fdea8..8355ac5 100644 --- a/src/freesas/containers.py +++ b/src/freesas/containers.py @@ -26,13 +26,13 @@ """ Set of namedtuples defined a bit everywhere """ + __authors__ = ["Jérôme Kieffer"] __license__ = "MIT" __copyright__ = "2020 ESRF" __date__ = "13/10/2020" from collections import namedtuple -from os import linesep import numpy # Used in AutoRg @@ -43,7 +43,7 @@ def _RG_RESULT_repr(self): - return f"Rg={self.Rg:6.4f}(±{self.sigma_Rg:6.4f}) I0={self.I0:6.4f}(±{self.sigma_I0:6.4f}) [{self.start_point}-{self.end_point}] {100.0*self.quality:5.2f}% {'aggregated' if self.aggregated>0.1 else ''}" + return f"Rg={self.Rg:6.4f}(±{self.sigma_Rg:6.4f}) I0={self.I0:6.4f}(±{self.sigma_I0:6.4f}) [{self.start_point}-{self.end_point}] {100.0 * self.quality:5.2f}% {'aggregated' if self.aggregated > 0.1 else ''}" RG_RESULT.__repr__ = _RG_RESULT_repr diff --git a/src/freesas/cormap.py b/src/freesas/cormap.py index d2dbca1..517fd53 100644 --- a/src/freesas/cormap.py +++ b/src/freesas/cormap.py @@ -5,7 +5,6 @@ import numpy from math import log from .containers import GOF - from ._cormap import measure_longest @@ -21,14 +20,14 @@ def __init__(self): def A(self, n, c): """Calculate A(number_of_toss, length_of_longest_run) - + :param n: number of coin toss in the experiment, an integer - :param c: length of the longest run of + :param c: length of the longest run of :return: The A parameter used in the formula - + """ if n <= c: - return 2 ** n + return 2**n elif (n, c) in self.knowledge: return self.knowledge[(n, c)] else: @@ -41,32 +40,32 @@ def A(self, n, c): def B(self, n, c): """Calculate B(number_of_toss, length_of_longest_run) to have either a run of Heads either a run of Tails - + :param n: number of coin toss in the experiment, an integer - :param c: length of the longest run of + :param c: length of the longest run of :return: The B parameter used in the formula """ return 2 * self.A(n - 1, c - 1) def __call__(self, n, c): - """Calculate the probability for the longest run of heads to exceed the observed length - + """Calculate the probability for the longest run of heads to exceed the observed length + :param n: number of coin toss in the experiment, an integer - :param c: length of the longest run of heads, an integer + :param c: length of the longest run of heads, an integer :return: The probablility of having c subsequent heads in a n toss of fair coin """ if c >= n: return 0 - delta = 2 ** n - self.A(n, c) + delta = 2**n - self.A(n, c) if delta <= 0: return 0 return 2.0 ** (log(delta, 2) - n) def probaHeadOrTail(self, n, c): - """Calculate the probability of a longest run of head or tails to occur - + """Calculate the probability of a longest run of head or tails to occur + :param n: number of coin toss in the experiment, an integer - :param c: length of the longest run of heads or tails, an integer + :param c: length of the longest run of heads or tails, an integer :return: The probablility of having c subsequent heads or tails in a n toss of fair coin """ if c > n: @@ -79,17 +78,17 @@ def probaHeadOrTail(self, n, c): return min(2.0 ** (log(delta, 2.0) - n), 1.0) def probaLongerRun(self, n, c): - """Calculate the probability for the longest run of heads or tails to exceed the observed length - + """Calculate the probability for the longest run of heads or tails to exceed the observed length + :param n: number of coin toss in the experiment, an integer - :param c: length of thee observed run of heads or tails, an integer + :param c: length of thee observed run of heads or tails, an integer :return: The probablility of having more than c subsequent heads or tails in a n toss of fair coin """ if c > n: return 0 if c == 0: return 0 - delta = (2 ** n) - self.B(n, c) + delta = (2**n) - self.B(n, c) if delta <= 0: return 0 return min(2.0 ** (log(delta, 2.0) - n), 1.0) @@ -99,11 +98,11 @@ def probaLongerRun(self, n, c): def gof(data1, data2): - """Calculate the probability for a couple of dataset to be equivalent - + """Calculate the probability for a couple of dataset to be equivalent + Implementation according to: http://www.nature.com/nmeth/journal/v12/n5/full/nmeth.3358.html - + :param data1: numpy array :param data2: numpy array :return: probablility for the 2 data to be equivalent diff --git a/src/freesas/decorators.py b/src/freesas/decorators.py index 6bdb517..c135392 100644 --- a/src/freesas/decorators.py +++ b/src/freesas/decorators.py @@ -32,7 +32,7 @@ __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" __date__ = "27/04/2020" __status__ = "development" -__docformat__ = 'restructuredtext' +__docformat__ = "restructuredtext" import sys import time @@ -44,8 +44,8 @@ def timeit(func): def wrapper(*arg, **kw): - '''This is the docstring of timeit: - a decorator that logs the execution time''' + """This is the docstring of timeit: + a decorator that logs the execution time""" t1 = time.perf_counter() res = func(*arg, **kw) t2 = time.perf_counter() diff --git a/src/freesas/dnn.py b/src/freesas/dnn.py index 48a34c0..a3d61d4 100644 --- a/src/freesas/dnn.py +++ b/src/freesas/dnn.py @@ -8,79 +8,88 @@ # Activation functions -def tanh(x): +def tanh(x): return np.tanh(x) -def relu(x): + +def relu(x): return np.maximum(0, x) -def sigmoid(x): + +def sigmoid(x): return 1 / (1 + np.exp(-x)) -def linear(x): + +def linear(x): return x + # Mapping activation names to functions activation_functions = { - 'tanh': tanh, - 'relu': relu, - 'sigmoid': sigmoid, - 'linear': linear + "tanh": tanh, + "relu": relu, + "sigmoid": sigmoid, + "linear": linear, } + # Parse config.json def parse_config(config_path): config = {} if isinstance(config_path, io.IOBase): config = json.load(config_path) elif os.path.exists(config_path): - with open(config_path, 'r') as f: + with open(config_path, "r") as f: config = json.load(f) else: - raise RuntimeError(f"config_path type {type(config_path)} not handled, got {config_path}") + raise RuntimeError( + f"config_path type {type(config_path)} not handled, got {config_path}" + ) layer_dims = [] activations = [] - - for layer in config.get('config',{}).get('layers', {}): - if layer['class_name'] == 'InputLayer': - layer_dims.append(layer['config']['batch_shape'][1]) - elif layer['class_name'] == 'Dense': - layer_dims.append(layer['config']['units']) - activation = layer['config']['activation'] + + for layer in config.get("config", {}).get("layers", {}): + if layer["class_name"] == "InputLayer": + layer_dims.append(layer["config"]["batch_shape"][1]) + elif layer["class_name"] == "Dense": + layer_dims.append(layer["config"]["units"]) + activation = layer["config"]["activation"] activations.append(activation_functions[activation]) - + return layer_dims, activations + # Load weights from model.weights.h5 def load_weights(weights_path, layer_dims): - with h5py.File(weights_path, 'r') as f: + with h5py.File(weights_path, "r") as f: params = [] for i in range(len(layer_dims) - 1): - layer_name = f'dense_{i}' if i > 0 else 'dense' - W = f[f'layers/{layer_name}/vars/0'][:] - b = f[f'layers/{layer_name}/vars/1'][:] + layer_name = f"dense_{i}" if i > 0 else "dense" + W = f[f"layers/{layer_name}/vars/0"][:] + b = f[f"layers/{layer_name}/vars/1"][:] params.extend([W, b]) return params + # Forward propagation def forward_propagation(X, params, activations): A = X L = len(params) // 2 - for l in range(L): - W, b = params[2*l], params[2*l+1] + for i in range(L): + W, b = params[2 * i], params[2 * i + 1] Z = np.dot(A, W) + b - A = activations[l](Z) + A = activations[i](Z) return A class DenseLayer: def __init__(self, weights, bias, activation): - """ Constructor - + """Constructor + :param weight: 2d array of size input_size x output_size :param bias: 1d array of size output_size :param activation: name of the activation function""" - assert weights.ndim == 2 + assert weights.ndim == 2 self.weights = weights assert bias.ndim == 1 assert weights.shape[1] == bias.size @@ -106,76 +115,88 @@ def forward(self, ary): """Forward Propagate the array and return the result of activation(Ax+b)""" Z = np.dot(ary, self.weights) + self.bias return self.activation(Z) + __call__ = forward + class DNN: def __init__(self, *args): """Constructor - + :param args: list of dense layers """ self.list = args - def infer(self,ary): + + def infer(self, ary): """Infer the neural network""" tmp = ary for layer in self.list: tmp = layer(tmp) return tmp + __call__ = infer -def preprocess(q, I): + +def preprocess(q, intensity): """ Preprocess the input data and infer Rg and Dmax using the DNN. - + :param dnn: An instance of the DNN class :param q: 1D array of q values - :param I: 1D array of intensity values + :param intensity: 1D array of intensity values :param dI: 1D array of the intensity error values :param q_interp: 1D array of the interpolated q values :return: Rg and Dmax """ - + # Define q_interp as a regularly spaced array in the range 0-4 nm^-1 with 1024 points q_interp = np.linspace(0, 4, 1024) - # Normalize I by Imax - Imax = I.max() - I_normalized = I / Imax + # Normalize intensity by Imax + Imax = intensity.max() + I_normalized = intensity / Imax - # Interpolate I over q_interp + # Interpolate intensity over q_interp I_interp = np.interp(q_interp, q, I_normalized, left=1, right=0) return I_interp def parse_keras_file(keras_file): - with zipfile.ZipFile(keras_file, 'r') as z: - with z.open('config.json') as config_file: - #config = parse_config(io.TextIOWrapper(config_file, 'tf-8')) + with zipfile.ZipFile(keras_file, "r") as z: + with z.open("config.json") as config_file: + # config = parse_config(io.TextIOWrapper(config_file, 'tf-8')) config = parse_config(config_file) - with z.open('model.weights.h5') as weights_file: + with z.open("model.weights.h5") as weights_file: weights = load_weights(io.BytesIO(weights_file.read()), config[0]) - #weights = load_weights(weights_file, config[0]) + # weights = load_weights(weights_file, config[0]) return config, weights class KerasDNN: def __init__(self, keras_file): config, weights = parse_keras_file(keras_file) - self.dnn = DNN(*[DenseLayer(weights[2*i], weights[2*i+1], a) for i,a in enumerate(config[1])]) - - def infer(self, q, I): - """Infer the neural network with q/I + self.dnn = DNN( + *[ + DenseLayer(weights[2 * i], weights[2 * i + 1], a) + for i, a in enumerate(config[1]) + ] + ) + + def infer(self, q, intensity): + """Infer the neural network with q/intensity :param q: 1D array, in inverse nm - :param I: 1D array, same size as q + :param intensity: 1D array, same size as q :return: result of the neural network. """ - Iprep = preprocess(q,I) + Iprep = preprocess(q, intensity) output = self.dnn.infer(Iprep) # Extract Rg and Dmax from the output Rg = output[0] # Assuming Rg is the first element Dmax = output[1] # Assuming Dmax is the second element - + return Rg, Dmax + __call__ = infer -Rg_Dmax = KerasDNN(resource_filename("keras_models/Rg+Dmax.keras")) \ No newline at end of file + +Rg_Dmax = KerasDNN(resource_filename("keras_models/Rg+Dmax.keras")) diff --git a/src/freesas/fitting.py b/src/freesas/fitting.py index ebebc2f..37ab0a1 100644 --- a/src/freesas/fitting.py +++ b/src/freesas/fitting.py @@ -5,7 +5,7 @@ __contact__ = "martha.brennich@googlemail.com" __license__ = "MIT" __copyright__ = "European Synchrotron Radiation Facility, Grenoble, France" -__date__ = "11/09/2024" +__date__ = "06/02/2026" __status__ = "development" __docformat__ = "restructuredtext" @@ -83,9 +83,7 @@ def get_linesep(output_destination: IO[str]) -> str: return "\n" -def get_guinier_header( - linesep: str, output_format: Optional[str] = None -) -> str: +def get_guinier_header(linesep: str, output_format: Optional[str] = None) -> str: """Return appropriate header line for selected output format :param output_format: output format from string parser :param linesep: correct linesep for chosen destination @@ -110,9 +108,8 @@ def get_guinier_header( else: return "" -def get_dnn_header( - linesep: str, output_format: Optional[str] = None -) -> str: + +def get_dnn_header(linesep: str, output_format: Optional[str] = None) -> str: """Return appropriate header line for selected output format :param output_format: output format from string parser :param linesep: correct linesep for chosen destination @@ -183,6 +180,7 @@ def rg_result_to_output_line( else: return f"{afile} {rg_result}{linesep}" + def dnn_result_to_output_line( dnn_result: tuple, afile: Path, @@ -280,6 +278,8 @@ def run_guinier_fit( ) output_destination.write(res) output_destination.flush() + + def run_dnn( parser: SASParser, logger: logging.Logger, @@ -289,7 +289,7 @@ def run_dnn( :param parser: a function that returns the output of argparse.parse() :param logger: a Logger """ - from .dnn import Rg_Dmax # heavy import + from .dnn import Rg_Dmax # heavy import args = parser.parse_args() set_logging_level(args.verbose) @@ -317,9 +317,9 @@ def run_dnn( else: if args.unit == "Å": data = convert_inverse_angstrom_to_nanometer(data) - q, I = data.T[:2] + q, intensity = data.T[:2] try: - dnn_result = Rg_Dmax(q, I) + dnn_result = Rg_Dmax(q, intensity) except ( InsufficientDataError, NoGuinierRegionError, diff --git a/src/freesas/invariants.py b/src/freesas/invariants.py index 423c208..5cabc89 100644 --- a/src/freesas/invariants.py +++ b/src/freesas/invariants.py @@ -3,7 +3,7 @@ # Project: freesas # https://github.com/kif/freesas # -# Copyright (C) 2020-2024 European Synchrotron Radiation Facility, Grenoble, France +# Copyright (C) 2020-2026 European Synchrotron Radiation Facility, Grenoble, France # # Permission is hereby granted, free of charge, to any person obtaining a copy # of this software and associated documentation files (the "Software"), to deal @@ -26,16 +26,16 @@ This module is mainly about the calculation of the Rambo-Tainer invariant described in: -https://dx.doi.org/10.1038%2Fnature12070 +https://dx.doi.org/10.1038%2Fnature12070 Some formula taken from Putnam et al, 2007, Table 1 in the review """ -__authors__ = ["Martha E. Brennich", "J. Kieffer"] + +__authors__ = ["Martha E. Brennich", "Jérôme Kieffer"] __license__ = "MIT" -__date__ = "31/05/2024" +__date__ = "06/02/2026" import logging -logger = logging.getLogger(__name__) import numpy from .containers import RT_RESULT @@ -45,53 +45,56 @@ from numpy import trapz as trapezoid # numpy1 +logger = logging.getLogger(__name__) + + def extrapolate(data, guinier): """Extrapolate SAS data according to the Guinier fit until q=0 - Uncertainties are extrapolated (linearly) from the Guinier region - - :param data: SAS data in q,I,dI format + Uncertainties are extrapolated (linearly) from the Guinier region + + :param data: SAS data in q, I, dI format :param guinier: result of a Guinier fit - :return: extrapolated SAS data + :return: extrapolated SAS data """ - + dq = data[1, 0] - data[0, 0] qmin = data[guinier.start_point, 0] q_low = numpy.arange(0, qmin, dq) - # Extrapolate I from Guinier approximation: + # Extrapolate I from Guinier approximation: I_low = guinier.I0 * numpy.exp(-(q_low**2 * guinier.Rg**2) / 3.0) # Extrapolate dI from Guinier region: - range_ = slice(guinier.start_point, guinier.end_point+1) + range_ = slice(guinier.start_point, guinier.end_point + 1) slope, intercept = numpy.polyfit(data[range_, 0], data[range_, 2], deg=1) - dI_low = abs(q_low*slope + intercept) + dI_low = abs(q_low * slope + intercept) # Now wrap-up data_low = numpy.vstack((q_low, I_low, dI_low)).T - return numpy.concatenate((data_low, data[guinier.start_point:])) + return numpy.concatenate((data_low, data[guinier.start_point :])) def calc_Porod(data, guinier): """Calculate the particle volume according to Porod's formula: - + V = 2*π²I₀²/(sum_q I(q)q² dq) - + Formula from Putnam's review, 2007, table 1 Intensities are extrapolated to q=0 using Guinier fit. - - :param data: SAS data in q,I,dI format + + :param data: SAS data in q, I, dI format :param Guinier: result of a Guinier fit (instance of RT_RESULT) :return: Volume calculated according to Porrod's formula - """ - q, I, dI = extrapolate(data, guinier).T - - denom = trapezoid(I*q**2, q) - volume = 2*numpy.pi**2*guinier.I0 / denom + """ + q, intensity, dI = extrapolate(data, guinier).T + + denom = trapezoid(intensity * q**2, q) + volume = 2 * numpy.pi**2 * guinier.I0 / denom return volume def calc_Vc(data, Rg, dRg, I0, dI0, imin): """Calculates the Rambo-Tainer invariant Vc, including extrapolation to q=0 - - :param data: SAS data in q,I,dI format, cropped to maximal q that should be used for calculation (normally 2 nm-1) + + :param data: SAS data in q, I, dI format, cropped to maximal q that should be used for calculation (normally 2 nm-1) :param Rg,dRg,I0,dI0: results from Guinier approximation/autorg :param imin: minimal index of the Guinier range, below that index data will be extrapolated by the Guinier approximation :returns: Vc and an error estimate based on non-correlated error propagation @@ -101,7 +104,19 @@ def calc_Vc(data, Rg, dRg, I0, dI0, imin): qlow = numpy.arange(0, qmin, dq) lowqint = trapezoid((qlow * I0 * numpy.exp(-(qlow * qlow * Rg * Rg) / 3.0)), qlow) - dlowqint = trapezoid(qlow * numpy.sqrt((numpy.exp(-(qlow * qlow * Rg * Rg) / 3.0) * dI0) ** 2 + ((I0 * 2.0 * (qlow * qlow) * Rg / 3.0) * numpy.exp(-(qlow * qlow * Rg * Rg) / 3.0) * dRg) ** 2), qlow) + dlowqint = trapezoid( + qlow + * numpy.sqrt( + (numpy.exp(-(qlow * qlow * Rg * Rg) / 3.0) * dI0) ** 2 + + ( + (I0 * 2.0 * (qlow * qlow) * Rg / 3.0) + * numpy.exp(-(qlow * qlow * Rg * Rg) / 3.0) + * dRg + ) + ** 2 + ), + qlow, + ) vabs = trapezoid(data[imin:, 0] * data[imin:, 1], data[imin:, 0]) dvabs = trapezoid(data[imin:, 0] * data[imin:, 2], data[imin:, 0]) vc = I0 / (lowqint + vabs) @@ -109,12 +124,11 @@ def calc_Vc(data, Rg, dRg, I0, dI0, imin): return (vc, dvc) -def calc_Rambo_Tainer(data, - guinier, qmax=2.0): +def calc_Rambo_Tainer(data, guinier, qmax=2.0): """calculates the invariants Vc and Qr from the Rambo & Tainer 2013 Paper, also the the mass estimate based on Qr for proteins - - :param data: data in q,I,dI format, q in nm^-1 + + :param data: data in q, I, dI format, q in nm^-1 :param guinier: RG_RESULT instance with result from the Guinier fit :param qmax: maximum q-value for the calculation in nm^-1 @return: dict with Vc, Qr and mass plus errors @@ -123,13 +137,24 @@ def calc_Rambo_Tainer(data, power_prot = 1.0 imax = abs(data[:, 0] - qmax).argmin() - if (imax <= guinier.start_point) or (guinier.start_point < 0): # unlikely but can happened - logger.error("Guinier region start too late for Rambo_Tainer invariants calculation") + if (imax <= guinier.start_point) or ( + guinier.start_point < 0 + ): # unlikely but can happened + logger.error( + "Guinier region start too late for Rambo_Tainer invariants calculation" + ) return None - vc = calc_Vc(data[:imax, :], guinier.Rg, guinier.sigma_Rg, guinier.I0, guinier.sigma_I0, guinier.start_point) + vc = calc_Vc( + data[:imax, :], + guinier.Rg, + guinier.sigma_Rg, + guinier.I0, + guinier.sigma_I0, + guinier.start_point, + ) qr = vc[0] ** 2 / (guinier.Rg) - mass = scale_prot * qr ** power_prot + mass = scale_prot * qr**power_prot dqr = qr * (guinier.sigma_Rg / guinier.Rg + 2 * ((vc[1]) / (vc[0]))) dmass = mass * dqr / qr diff --git a/src/freesas/model.py b/src/freesas/model.py index 89dff9d..584eca1 100644 --- a/src/freesas/model.py +++ b/src/freesas/model.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # coding: utf-8 -__author__ = "Guillaume" +__author__ = "Guillaume Bonamis" __license__ = "MIT" __copyright__ = "2015, ESRF" @@ -9,6 +9,7 @@ from math import sqrt import threading import numpy + try: from . import _distance except ImportError: @@ -16,6 +17,11 @@ from . import transformations +def delta_kron(i, j): + "Delta Kronecker function" + return 1 if i == j else 0 + + def delta_expand(vec1, vec2): """Create a 2d array with the difference vec1[i]-vec2[j] @@ -43,7 +49,9 @@ def __init__(self, molecule=None): if isinstance(molecule, (str, bytes)) and os.path.exists(molecule): self.read(molecule) else: - self.atoms = molecule if molecule is not None else [] # initial coordinates of each dummy atoms of the molecule, fourth column full of one for the transformation matrix + self.atoms = ( + molecule if molecule is not None else [] + ) # initial coordinates of each dummy atoms of the molecule, fourth column full of one for the transformation matrix self.header = "" # header of the PDB file self.rfactor = None self.radius = 1.0 # unused at the moment @@ -76,12 +84,16 @@ def read(self, filename): y = float(line[38:46]) z = float(line[46:54]) atoms.append([x, y, z]) - if line.startswith("REMARK 265 Final R-factor"): # very dependent of the pdb file format ! + if line.startswith( + "REMARK 265 Final R-factor" + ): # very dependent of the pdb file format ! self.rfactor = float(line[43:56]) header.append(line) self.header = header atom3 = numpy.array(atoms) - self.atoms = numpy.append(atom3, numpy.ones((atom3.shape[0], 1), dtype="float"), axis=1) + self.atoms = numpy.append( + atom3, numpy.ones((atom3.shape[0], 1), dtype="float"), axis=1 + ) def save(self, filename): """ @@ -95,7 +107,11 @@ def save(self, filename): for line in self.header: if line.startswith("ATOM"): if nr < self.atoms.shape[0]: - line = line[:30] + "%8.3f%8.3f%8.3f" % tuple(self.atoms[nr]) + line[54:] + line = ( + line[:30] + + "%8.3f%8.3f%8.3f" % tuple(self.atoms[nr]) + + line[54:] + ) else: line = "" nr += 1 @@ -122,10 +138,12 @@ def inertiatensor(self): mol = self.atoms[:, 0:3] - self.com self.inertensor = numpy.empty((3, 3), dtype="float") - delta_kron = lambda i, j: 1 if i == j else 0 + for i in range(3): for j in range(i, 3): - self.inertensor[i, j] = self.inertensor[j, i] = (delta_kron(i, j) * (mol ** 2).sum(axis=1) - (mol[:, i] * mol[:, j])).sum() / mol.shape[0] + self.inertensor[i, j] = self.inertensor[j, i] = ( + delta_kron(i, j) * (mol**2).sum(axis=1) - (mol[:, i] * mol[:, j]) + ).sum() / mol.shape[0] return self.inertensor def canonical_translate(self): @@ -162,7 +180,10 @@ def canonical_rotate(self): self.enantiomer = [1, 1, 1] else: self.enantiomer = [-1, -1, -1] - mirror = numpy.array([[-1, 0, 0, 0], [0, -1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]], dtype="float") + mirror = numpy.array( + [[-1, 0, 0, 0], [0, -1, 0, 0], [0, 0, -1, 0], [0, 0, 0, 1]], + dtype="float", + ) rot = numpy.dot(mirror, rot) return rot @@ -194,7 +215,11 @@ def calc_invariants(self, use_cython=True): else: size = self.atoms.shape[0] - D = delta_expand(self.atoms[:, 0], self.atoms[:, 0]) ** 2 + delta_expand(self.atoms[:, 1], self.atoms[:, 1]) ** 2 + delta_expand(self.atoms[:, 2], self.atoms[:, 2]) ** 2 + D = ( + delta_expand(self.atoms[:, 0], self.atoms[:, 0]) ** 2 + + delta_expand(self.atoms[:, 1], self.atoms[:, 1]) ** 2 + + delta_expand(self.atoms[:, 2], self.atoms[:, 2]) ** 2 + ) Rg = sqrt(D.sum() / 2.0) / size Dmax = sqrt(D.max()) d12 = (D.max() * numpy.eye(size) + D).min(axis=0).mean() @@ -235,7 +260,9 @@ def dist(self, other, molecule1, molecule2, use_cython=True): :return D: NSD between the 2 molecules, in their position molecule1 and molecule2 """ if _distance and use_cython: - return _distance.calc_distance(molecule1, molecule2, self.fineness, other.fineness) + return _distance.calc_distance( + molecule1, molecule2, self.fineness, other.fineness + ) else: mol1 = molecule1[:, 0:3] @@ -255,9 +282,21 @@ def dist(self, other, molecule1, molecule2, use_cython=True): mol2y.shape = mol2.shape[0], 1 mol2z.shape = mol2.shape[0], 1 - d2 = delta_expand(mol1x, mol2x) ** 2 + delta_expand(mol1y, mol2y) ** 2 + delta_expand(mol1z, mol2z) ** 2 - - D = (0.5 * ((1. / ((mol1.shape[0]) * other.fineness * other.fineness)) * (d2.min(axis=1).sum()) + (1. / ((mol2.shape[0]) * self.fineness * self.fineness)) * (d2.min(axis=0)).sum())) ** 0.5 + d2 = ( + delta_expand(mol1x, mol2x) ** 2 + + delta_expand(mol1y, mol2y) ** 2 + + delta_expand(mol1z, mol2z) ** 2 + ) + + D = ( + 0.5 + * ( + (1.0 / ((mol1.shape[0]) * other.fineness * other.fineness)) + * (d2.min(axis=1).sum()) + + (1.0 / ((mol2.shape[0]) * self.fineness * self.fineness)) + * (d2.min(axis=0)).sum() + ) + ) ** 0.5 return D def transform(self, param, symmetry, reverse=None): @@ -270,10 +309,18 @@ def transform(self, param, symmetry, reverse=None): """ mol = self.atoms - sym = numpy.array([[symmetry[0], 0, 0, 0], [0, symmetry[1], 0, 0], [0, 0, symmetry[2], 0], [0, 0, 0, 1]], dtype="float") + sym = numpy.array( + [ + [symmetry[0], 0, 0, 0], + [0, symmetry[1], 0, 0], + [0, 0, symmetry[2], 0], + [0, 0, 0, 1], + ], + dtype="float", + ) if not reverse: vect = numpy.array([param[0:3]]) - angles = (param[3:6]) + angles = param[3:6] translat1 = transformations.translation_matrix(vect) rotation = transformations.euler_matrix(*angles) @@ -305,9 +352,13 @@ def dist_after_movement(self, param, other, symmetry): self.canonical_parameters() can_param1 = self.can_param - molref_can = self.transform(can_param1, [1, 1, 1]) # molecule reference put on its canonical position + molref_can = self.transform( + can_param1, [1, 1, 1] + ) # molecule reference put on its canonical position - mol2_moved = other.transform(param, symmetry) # movement selected applied to mol2 + mol2_moved = other.transform( + param, symmetry + ) # movement selected applied to mol2 distance = self.dist(other, molref_can, mol2_moved) return distance diff --git a/src/freesas/nexus_parser.py b/src/freesas/nexus_parser.py index f9f55cf..7fb2565 100644 --- a/src/freesas/nexus_parser.py +++ b/src/freesas/nexus_parser.py @@ -1,9 +1,10 @@ __author__ = "Jérôme Kieffer" __license__ = "MIT" __copyright__ = "2017, ESRF" -__date__ = "20/01/2025" +__date__ = "06/02/2026" -import sys, os +import sys +import os import zipfile import posixpath import logging @@ -11,6 +12,7 @@ from silx.io.nxdata import NXdata from dataclasses import dataclass import numpy + logger = logging.getLogger(__name__) try: @@ -18,7 +20,7 @@ except ImportError: logger.error("H5py is mandatory to parse HDF5 files") h5py = None - + @dataclass class IntegratedPattern: @@ -27,7 +29,7 @@ class IntegratedPattern: point: Union[float, int, None] radial: numpy.ndarray intensity: numpy.ndarray - intensity_errors: Union[numpy.ndarray, None]=None + intensity_errors: Union[numpy.ndarray, None] = None radial_name: str = "" radial_units: str = "" intensity_name: str = "" @@ -39,10 +41,10 @@ def __repr__(self): line += " \t uncertainties" res = [line] if self.intensity_errors is None: - for q,i,s in zip(self.radial, self.intensity): + for q, i, s in zip(self.radial, self.intensity): res.append(f"{q} \t {i}") else: - for q,i,s in zip(self.radial, self.intensity, self.intensity_errors): + for q, i, s in zip(self.radial, self.intensity, self.intensity_errors): res.append(f"{q} \t {i} \t {s}") return os.linesep.join(res) @@ -75,17 +77,16 @@ def read_nexus_integrated_patterns(group): if nxdata.axes[-1] is None: radial = numpy.arange(nxdata.signal.shape[1]) radial_units = "" - radial_name = "" + radial_name = "" else: axis_dataset = nxdata.axes[-1] radial = axis_dataset[()] - radial_name = axis_dataset.name.split("/")[-1] + radial_name = axis_dataset.name.split("/")[-1] radial_units = axis_dataset.attrs.get("units", "") intensities = numpy.atleast_2d(nxdata.signal) intensity_name = nxdata.signal.name.split("/")[-1] intensity_units = nxdata.signal.attrs.get("units", "") - if nxdata.errors is None: errors = [None] * intensities.shape[0] @@ -95,19 +96,31 @@ def read_nexus_integrated_patterns(group): if (len(points), len(radial)) != intensities.shape: raise RuntimeError("Shape mismatch between axes and signal") - return [IntegratedPattern( - point, radial, intensity, intensity_errors, radial_name, radial_units, intensity_name, intensity_units) for point, intensity, intensity_errors in zip(points, intensities, errors)] + return [ + IntegratedPattern( + point, + radial, + intensity, + intensity_errors, + radial_name, + radial_units, + intensity_name, + intensity_units, + ) + for point, intensity, intensity_errors in zip(points, intensities, errors) + ] class Tree: def __init__(self, root=None): self.root = root or {} self.skip = set() + def visit_item(self, name, obj): if name in self.skip: - return + return node = self.root - path = [i.replace(" ","_") for i in name.split("/")] + path = [i.replace(" ", "_") for i in name.split("/")] for i in path[:-1]: if i not in node: node[i] = {} @@ -119,57 +132,62 @@ def visit_item(self, name, obj): except (KeyError, OSError) as err: print(f"{type(err).__name__}: {err} while readding {path}") for key in obj: - self.skip.add(posixpath.join(name,key)) + self.skip.add(posixpath.join(name, key)) if isinstance(obj[key], h5py.Group): for sub in obj[key]: - self.skip.add(posixpath.join(name,key, sub)) + self.skip.add(posixpath.join(name, key, sub)) else: node[path[-1]] = {} if isinstance(obj, h5py.Dataset): if len(obj.shape) <= 1: node[path[-1]] = obj[()] - + def save(self, filename): with zipfile.ZipFile(filename, "w") as z: + def write(path, name, obj): new_path = posixpath.join(path, name) if isinstance(obj, dict): - if sys.version_info>=(3,11): z.mkdir(new_path) + if sys.version_info >= (3, 11): + z.mkdir(new_path) for key, value in obj.items(): write(new_path, key, value) elif isinstance(obj, numpy.ndarray): if obj.ndim == 1: z.writestr(new_path, os.linesep.join(str(i) for i in obj)) else: - z.writestr(new_path, str(obj)) + z.writestr(new_path, str(obj)) elif isinstance(obj, list): - if sys.version_info>=(3,11): z.mkdir(new_path) - if len(obj)==1: - fname = new_path+"/biosaxs.dat" - z.writestr(fname, str(obj[0])) + if sys.version_info >= (3, 11): + z.mkdir(new_path) + if len(obj) == 1: + fname = new_path + "/biosaxs.dat" + z.writestr(fname, str(obj[0])) else: - for i,j in enumerate(obj): - fname = new_path+f"/bioxaxs_{i:04d}.dat" + for i, j in enumerate(obj): + fname = new_path + f"/bioxaxs_{i:04d}.dat" z.writestr(fname, str(j)) elif isinstance(obj, (int, float, numpy.number, bool, numpy.bool)): - z.writestr(new_path, str(obj)) + z.writestr(new_path, str(obj)) elif isinstance(obj, (str, bytes)): - z.writestr(new_path, obj) + z.writestr(new_path, obj) else: print(f"skip {new_path} for {obj} of type {obj.__class__.__mro__}") - + root = self.root for key, value in root.items(): write("", key, value) + def get(self, path): node = self.root for i in path.split("/"): node = node[i] return node - + + def convert_nexus2zip(nexusfile, outfile=None): - """ Convert a nexus-file, as produced by BM29 beamline into a zip file - + """Convert a nexus-file, as produced by BM29 beamline into a zip file + :param nexusfile: string with the path of the input file :param outfile: name of the output file, unless, just replace the extension with ".zip" :return: nothing, maybe an error code ? @@ -177,6 +195,5 @@ def convert_nexus2zip(nexusfile, outfile=None): tree = Tree() with h5py.File(nexusfile, "r") as h: h.visititems(tree.visit_item) - outfile = outfile or (os.path.splitext(nexusfile)[0]+".h5") + outfile = outfile or (os.path.splitext(nexusfile)[0] + ".h5") tree.save(outfile) - diff --git a/src/freesas/plot.py b/src/freesas/plot.py index a166584..32290f6 100644 --- a/src/freesas/plot.py +++ b/src/freesas/plot.py @@ -3,17 +3,17 @@ Functions to generating graphs related to SAS. """ -__authors__ = ["Jerome Kieffer"] +__authors__ = ["Jérôme Kieffer"] __license__ = "MIT" -__copyright__ = "2020, ESRF" -__date__ = "15/09/2022" +__copyright__ = "2020-2026, ESRF" +__date__ = "06/02/2026" import logging - -logger = logging.getLogger(__name__) import numpy from matplotlib.pyplot import subplots +logger = logging.getLogger(__name__) + def scatter_plot( data, @@ -30,7 +30,7 @@ def scatter_plot( """ Generate a scattering plot I = f(q) in semi_log_y. - :param data: data read from an ASCII file, 3 column (q,I,err) + :param data: data read from an ASCII file, 3 column (q, I, err) :param filename: name of the file where the cuve should be saved :param img_format: image format :param unit: Unit name for Rg and 1/q @@ -49,10 +49,10 @@ def scatter_plot( assert data.ndim == 2 assert data.shape[1] >= 2 q = data.T[0] - I = data.T[1] + intensity = data.T[1] try: err = data.T[2] - except: + except Exception: err = None if ax: fig = ax.figure @@ -84,23 +84,23 @@ def scatter_plot( if err is not None: ax.errorbar( q, - I, - numpy.maximum(0,err), + intensity, + numpy.maximum(0, err), label=label_exp, capsize=0, color=exp_color, ecolor=err_color, ) else: - ax.plot(q, I, label=label_exp, color="blue") + ax.plot(q, intensity, label=label_exp, color="blue") else: q_guinier = q[first_point:last_point] I_guinier = I0 * numpy.exp(-((q_guinier * rg) ** 2) / 3) if err is not None: ax.errorbar( q, - I, - numpy.maximum(0,err), + intensity, + numpy.maximum(0, err), label=label_exp, capsize=0, color=exp_color, @@ -108,7 +108,7 @@ def scatter_plot( alpha=0.5, ) else: - ax.plot(q, I, label=label_exp, color=exp_color, alpha=0.5) + ax.plot(q, intensity, label=label_exp, color=exp_color, alpha=0.5) label_guinier += ": $R_g=$%.2f %s, $I_0=$%.2f" % (rg, unit, I0) ax.plot( q_guinier, @@ -151,10 +151,10 @@ def scatter_plot( crv, lab = ax.get_legend_handles_labels() ordered_lab = [] ordered_crv = [] - for l in [label_exp, label_guinier, label_ift]: + for lbl in [label_exp, label_guinier, label_ift]: try: - idx = lab.index(l) - except: + idx = lab.index(lbl) + except Exception: continue ordered_lab.append(lab[idx]) ordered_crv.append(crv[idx]) @@ -183,7 +183,7 @@ def kratky_plot( """ Generate a Kratky plot q²Rg²I/I₀ = f(q·Rg) - :param data: data read from an ASCII file, 3 column (q,I,err) + :param data: data read from an ASCII file, 3 column (q, I, err) :param guinier: output of autoRg :param filename: name of the file where the cuve should be saved :param img_format: image format @@ -195,10 +195,10 @@ def kratky_plot( assert data.ndim == 2 assert data.shape[1] >= 2 q = data.T[0] - I = data.T[1] + intensity = data.T[1] try: err = data.T[2] - except: + except Exception: err = None if ax: fig = ax.figure @@ -208,10 +208,10 @@ def kratky_plot( I0 = guinier.I0 xdata = q * Rg - ydata = xdata * xdata * I / I0 + ydata = xdata * xdata * intensity / I0 if err is not None: - dy = xdata * xdata * numpy.maximum(0,err) / abs(I0) - dplot = ax.errorbar( + dy = xdata * xdata * numpy.maximum(0, err) / abs(I0) + ax.errorbar( xdata, ydata, dy, @@ -221,7 +221,7 @@ def kratky_plot( ecolor="lightblue", ) else: - dplot = ax.plot(xdata, ydata, label=label, color="blue") + ax.plot(xdata, ydata, label=label, color="blue") ax.set_ylabel("$(qR_{g})^2 I/I_{0}$", fontsize=fontsize) ax.set_xlabel("$qR_{g}$", fontsize=fontsize) ax.legend(loc=1) @@ -269,7 +269,7 @@ def guinier_plot( """ Generate a guinier plot: ln(I) = f(q²) - :param data: data read from an ASCII file, 3 column (q,I,err) + :param data: data read from an ASCII file, 3 column (q, I, err) :param guinier: A RG_RESULT object from AutoRg :param filename: name of the file where the cuve should be saved :param img_format: image format @@ -279,9 +279,9 @@ def guinier_plot( assert data.ndim == 2 assert data.shape[1] >= 2 - q, I, err = data.T[:3] + q, intensity, err = data.T[:3] - mask = (I > 0) & numpy.isfinite(I) & (q > 0) & numpy.isfinite(q) + mask = (intensity > 0) & numpy.isfinite(intensity) & (q > 0) & numpy.isfinite(q) if err is not None: mask &= (err > 0.0) & numpy.isfinite(err) mask = mask.astype(bool) @@ -295,14 +295,14 @@ def guinier_plot( mask[end:] = False q2 = q[mask] ** 2 - logI = numpy.log(I[mask]) + logI = numpy.log(intensity[mask]) if ax: fig = ax.figure else: fig, ax = subplots(figsize=(12, 10)) if err is not None: - dlogI = err[mask] / I[mask] + dlogI = err[mask] / intensity[mask] ax.errorbar( q2, logI, @@ -324,13 +324,11 @@ def guinier_plot( # ax.plot(q2[first_point:last_point], logI[first_point:last_point], marker='D', markersize=5, label="guinier region") xmin = q[first_point] ** 2 xmax = q[last_point] ** 2 - ymax = numpy.log(I[first_point]) - ymin = numpy.log(I[last_point]) + ymax = numpy.log(intensity[first_point]) + ymin = numpy.log(intensity[last_point]) dy = (ymax - ymin) / 2.0 ax.vlines(xmin, ymin=ymin, ymax=ymax + dy, color="0.75", linewidth=1.0) - ax.vlines( - xmax, ymin=ymin - dy, ymax=ymin + dy, color="0.75", linewidth=1.0 - ) + ax.vlines(xmax, ymin=ymin - dy, ymax=ymin + dy, color="0.75", linewidth=1.0) ax.annotate( "$(qR_{g})_{min}$=%.1f" % (Rg * q[first_point]), (xmin, ymax + dy), @@ -509,14 +507,17 @@ def plot_all( fig.savefig(filename) return fig -def hplc_plot(hplc, - fractions = None, - title="Chromatogram", - filename=None, - img_format="png", - ax=None, - labelsize=None, - fontsize=None,): + +def hplc_plot( + hplc, + fractions=None, + title="Chromatogram", + filename=None, + img_format="png", + ax=None, + labelsize=None, + fontsize=None, +): """ Generate an HPLC plot I=f(t) @@ -531,32 +532,33 @@ def hplc_plot(hplc, fig = ax.figure else: fig, ax = subplots(figsize=(12, 10)) - data = [sum(i) if hasattr(i, '__iter__') else i for i in hplc] - ax.plot(data, label = "Chromatogram") + data = [sum(i) if hasattr(i, "__iter__") else i for i in hplc] + ax.plot(data, label="Chromatogram") ax.set_xlabel("Elution (frame index)", fontsize=fontsize) ax.set_ylabel("Summed intensities", fontsize=fontsize) ax.set_title(title) - + ax.tick_params(axis="x", labelsize=labelsize) ax.tick_params(axis="y", labelsize=labelsize) if fractions is not None and len(fractions): fractions.sort() - l = len(data) - idx = list(range(l)) - for start,stop in fractions: - start = int(min(l-1, max(0, start))) - stop = int(min(l-1, max(0, stop))) - ax.plot(idx[start:stop+1], - data[start:stop+1], - label=f"Fraction {start}-{stop}", - linewidth=10, - alpha=0.5) + nbdata = len(data) + idx = list(range(nbdata)) + for start, stop in fractions: + start = int(min(nbdata - 1, max(0, start))) + stop = int(min(nbdata - 1, max(0, stop))) + ax.plot( + idx[start : stop + 1], + data[start : stop + 1], + label=f"Fraction {start}-{stop}", + linewidth=10, + alpha=0.5, + ) ax.legend() - + if filename: if img_format: fig.savefig(filename, format=img_format) else: fig.savefig(filename) return fig - diff --git a/src/freesas/resources/__init__.py b/src/freesas/resources/__init__.py index 5ddcc71..d9a4c26 100644 --- a/src/freesas/resources/__init__.py +++ b/src/freesas/resources/__init__.py @@ -68,21 +68,22 @@ # importlib_resources is useful when this package is stored in a zip # When importlib.resources is not available, the resources dir defaults to the # directory containing this module. -if sys.version_info >= (3,9): +if sys.version_info >= (3, 9): import importlib.resources as importlib_resources else: try: - import importlib_resources + import importlib_resources except ImportError: logger.info("Unable to import importlib_resources") logger.debug("Backtrace", exc_info=True) importlib_resources = None if importlib_resources is not None: - import atexit - from contextlib import ExitStack - file_manager = ExitStack() - atexit.register(file_manager.close) + import atexit + from contextlib import ExitStack + + file_manager = ExitStack() + atexit.register(file_manager.close) # For packaging purpose, patch this variable to use an alternative directory @@ -96,11 +97,11 @@ # cx_Freeze forzen support # See http://cx-freeze.readthedocs.io/en/latest/faq.html#using-data-files -if getattr(sys, 'frozen', False): +if getattr(sys, "frozen", False): # Running in a frozen application: # We expect resources to be located either in a pyFAI/resources/ dir # relative to the executable or within this package. - _dir = os.path.join(os.path.dirname(sys.executable), 'freesas', 'resources') + _dir = os.path.join(os.path.dirname(sys.executable), "freesas", "resources") if os.path.isdir(_dir): _RESOURCES_DIR = _dir @@ -124,10 +125,11 @@ def resource_filename(resource): # return os.path.join(_RESOURCES_DOC_DIR, *resource.split('/')[1:]) if _RESOURCES_DIR is not None: # if set, use this directory - return os.path.join(_RESOURCES_DIR, *resource.split('/')) + return os.path.join(_RESOURCES_DIR, *resource.split("/")) elif importlib_resources is None: # Fallback if pkg_resources is not available - return os.path.join(os.path.abspath(os.path.dirname(__file__)), - *resource.split('/')) + return os.path.join( + os.path.abspath(os.path.dirname(__file__)), *resource.split("/") + ) else: # Preferred way to get resources as it supports zipfile package ref = importlib_resources.files(__name__) / resource path = file_manager.enter_context(importlib_resources.as_file(ref)) @@ -143,7 +145,6 @@ def silx_integration(): if _integrated: return import silx.resources - silx.resources.register_resource_directory("freesas", - __name__, - _RESOURCES_DIR) + + silx.resources.register_resource_directory("freesas", __name__, _RESOURCES_DIR) _integrated = True diff --git a/src/freesas/sas_argparser.py b/src/freesas/sas_argparser.py index 50b2ef8..7a3fd1c 100644 --- a/src/freesas/sas_argparser.py +++ b/src/freesas/sas_argparser.py @@ -29,7 +29,6 @@ def parse_unit(unit_input: str) -> str: class SASParser: - """ Wrapper class for argparse ArgumentParser that provides predefined argument. """ @@ -69,11 +68,11 @@ def __init__(self, prog: str, description: str, epilog: str, **kwargs): self.add_argument("-V", "--version", action="version", version=version) def parse_args(self, *args, **kwargs): - """ Wrapper for argparse parse_args() """ + """Wrapper for argparse parse_args()""" return self.parser.parse_args(*args, **kwargs) def add_argument(self, *args, **kwargs): - """ Wrapper for argparse add_argument() """ + """Wrapper for argparse add_argument()""" self.parser.add_argument(*args, **kwargs) def add_file_argument(self, help_text: str): @@ -100,7 +99,7 @@ def add_q_unit_argument(self): ) def add_output_filename_argument(self): - """ Add default argument for specifying output format. """ + """Add default argument for specifying output format.""" self.add_argument( "-o", "--output", @@ -111,7 +110,7 @@ def add_output_filename_argument(self): ) def add_output_data_format(self, *formats: str, default: str = None): - """ Add default argument for specifying output format. """ + """Add default argument for specifying output format.""" help_string = "Output format: " + ", ".join(formats) self.add_argument( "-f", @@ -150,16 +149,14 @@ def __init__(self, prog: str, description: str, epilog: str, **kwargs): ) self.parser.add_file_argument(help_text=file_help_text) self.parser.add_output_filename_argument() - self.parser.add_output_data_format( - "native", "csv", "ssf", default="native" - ) + self.parser.add_output_data_format("native", "csv", "ssf", default="native") self.parser.add_q_unit_argument() self.usage = self.parser.usage def parse_args(self, *args, **kwargs): - """ Wrapper for SASParser parse_args() """ + """Wrapper for SASParser parse_args()""" return self.parser.parse_args(*args, **kwargs) def add_argument(self, *args, **kwargs): - """ Wrapper for SASParser add_argument() """ + """Wrapper for SASParser add_argument()""" self.parser.add_argument(*args, **kwargs) diff --git a/src/freesas/sasio.py b/src/freesas/sasio.py index cb7084a..8f39e7d 100644 --- a/src/freesas/sasio.py +++ b/src/freesas/sasio.py @@ -12,6 +12,7 @@ """ Contains helper functions for loading SAS data from differents sources. """ + __authors__ = ["Martha Brennich"] __contact__ = "martha.brennich@googlemail.com" __license__ = "MIT" @@ -37,9 +38,9 @@ def load_scattering_data(filename: PathType) -> ndarray: """ try: data = loadtxt(filename) - except OSError as err: + except OSError: raise OSError("File could not be read.") - except ValueError as err: + except ValueError: text = None if isinstance(filename, (io.StringIO, io.BytesIO)): filename.seek(0) @@ -54,15 +55,11 @@ def load_scattering_data(filename: PathType) -> ndarray: try: data = parse_ascii_data(text, number_of_columns=3) except ValueError: - raise ValueError( - "File does not seem to be " "in the format q, I, err. " - ) + raise ValueError("File does not seem to be in the format q, I, err. ") return data -def parse_ascii_data( - input_file_text: List[str], number_of_columns: int -) -> ndarray: +def parse_ascii_data(input_file_text: List[str], number_of_columns: int) -> ndarray: """ Parse data from an ascii file into an N column numpy array diff --git a/src/freesas/test/__init__.py b/src/freesas/test/__init__.py index dbcadf1..ff7e07c 100644 --- a/src/freesas/test/__init__.py +++ b/src/freesas/test/__init__.py @@ -25,5 +25,5 @@ def run_tests(): run = run_tests -if __name__ == '__main__': +if __name__ == "__main__": sys.exit(run_tests()) diff --git a/src/freesas/test/mock_open_38.py b/src/freesas/test/mock_open_38.py index cf85a66..a61b188 100644 --- a/src/freesas/test/mock_open_38.py +++ b/src/freesas/test/mock_open_38.py @@ -1,7 +1,7 @@ """ This is the Python 3.8 implementation of mock_open taken from https://github.com/python/cpython/blob/3.8/Lib/unittest/mock.py -Hence: +Hence: "Copyright (c) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 Python Software Foundation; All Rights Reserved" @@ -12,8 +12,9 @@ file_spec = None -#sentinel = _Sentinel() -#DEFAULT = sentinel.DEFAULT +# sentinel = _Sentinel() +# DEFAULT = sentinel.DEFAULT + def _to_stream(read_data): if isinstance(read_data, bytes): @@ -22,7 +23,7 @@ def _to_stream(read_data): return io.StringIO(read_data) -def mock_open(mock=None, read_data=''): +def mock_open(mock=None, read_data=""): """ A helper function to create a mock to replace the use of `open`. It works for `open` called directly or used as a context manager. @@ -65,10 +66,11 @@ def _next_side_effect(): global file_spec if file_spec is None: import _io + file_spec = list(set(dir(_io.TextIOWrapper)).union(set(dir(_io.BytesIO)))) if mock is None: - mock = MagicMock(name='open', spec=open) + mock = MagicMock(name="open", spec=open) handle = MagicMock(spec=file_spec) handle.__enter__.return_value = handle diff --git a/src/freesas/test/test_align.py b/src/freesas/test/test_align.py index 3e40eb5..9b6ebd5 100644 --- a/src/freesas/test/test_align.py +++ b/src/freesas/test/test_align.py @@ -7,10 +7,10 @@ import numpy import unittest from .utilstest import get_datafile -from ..model import SASModel from ..align import AlignModels from ..transformations import translation_matrix, euler_matrix import logging + logging.basicConfig(level=logging.INFO) logger = logging.getLogger("AlignModels_test") @@ -80,7 +80,11 @@ def test_usefull_alignment(self): dist_before = mol1.dist(mol2, mol1.atoms, mol2.atoms) symmetry, par = align.alignment_sym(mol1, mol2) dist_after = mol1.dist_after_movement(par, mol2, symmetry) - self.assertGreaterEqual(dist_before, dist_after, "increase of distance after alignment %s<%s" % (dist_before, dist_after)) + self.assertGreaterEqual( + dist_before, + dist_after, + "increase of distance after alignment %s<%s" % (dist_before, dist_after), + ) def test_optimisation_align(self): inputfiles = [self.testfile1, self.testfile2] @@ -94,7 +98,12 @@ def test_optimisation_align(self): align.slow = True sym, p = align.alignment_sym(mol1, mol2) dist_after = mol1.dist_after_movement(p, mol2, sym) - self.assertGreaterEqual(dist_before, dist_after, "increase of distance after optimized alignment %s<%s" % (dist_before, dist_after)) + self.assertGreaterEqual( + dist_before, + dist_after, + "increase of distance after optimized alignment %s<%s" + % (dist_before, dist_after), + ) def test_alignment_intruder(self): intruder = numpy.random.randint(0, 8) @@ -112,12 +121,16 @@ def test_alignment_intruder(self): if table.sum() == 0: logger.error("there is no intruders") - averNSD = ((table.sum(axis=-1)) / (align.validmodels.sum() - 1)) + averNSD = (table.sum(axis=-1)) / (align.validmodels.sum() - 1) num_intr = averNSD.argmax() if not num_intr and num_intr != 0: logger.error("cannot find the intruder") - self.assertEqual(num_intr, intruder, msg="not find the good intruder, %s!=%s" % (num_intr, intruder)) + self.assertEqual( + num_intr, + intruder, + msg="not find the good intruder, %s!=%s" % (num_intr, intruder), + ) def test_reference(self): inputfiles = [self.testfile1] * 8 @@ -146,6 +159,7 @@ def suite(): testSuite.addTest(TestAlign("test_reference")) return testSuite -if __name__ == '__main__': + +if __name__ == "__main__": runner = unittest.TextTestRunner() runner.run(suite()) diff --git a/src/freesas/test/test_autorg.py b/src/freesas/test/test_autorg.py index 0de7f4c..d26c08f 100644 --- a/src/freesas/test/test_autorg.py +++ b/src/freesas/test/test_autorg.py @@ -23,9 +23,9 @@ # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. -__authors__ = ["J. Kieffer"] +__authors__ = ["Jérôme Kieffer"] __license__ = "MIT" -__date__ = "29/11/2023" +__date__ = "06/02/2026" import logging import unittest @@ -59,9 +59,9 @@ def create_synthetic_data(R0=4, I0=100): q = numpy.linspace(0, 10, size) qr = numpy.outer(q, r / pi) T = (4 * pi * (r[-1] - r[0]) / npt) * numpy.sinc(qr) - I = T.dot(p) - err = numpy.sqrt(I) - return numpy.vstack((q, I, err)).T[1:] + intensity = T.dot(p) + err = numpy.sqrt(intensity) + return numpy.vstack((q, intensity, err)).T[1:] class TestAutoRg(unittest.TestCase): @@ -116,9 +116,7 @@ def test_synthetic(self): data = create_synthetic_data(R0=R0, I0=I0) Rg = autoRg(data) logger.info("auto_rg %s", Rg) - self.assertAlmostEqual( - R0 * sqrt(3 / 5), Rg.Rg, 0, "Rg matches for a sphere" - ) + self.assertAlmostEqual(R0 * sqrt(3 / 5), Rg.Rg, 0, "Rg matches for a sphere") self.assertGreater( R0 * sqrt(3 / 5), Rg.Rg - Rg.sigma_Rg, @@ -158,16 +156,12 @@ def test_synthetic(self): "Rg in range matches for a sphere", ) self.assertAlmostEqual(I0, guinier.I0, 0, "I0 matches for a sphere") - self.assertGreater( - I0, guinier.I0 - sigma_I0, "I0 matches for a sphere" - ) + self.assertGreater(I0, guinier.I0 - sigma_I0, "I0 matches for a sphere") self.assertLess(I0, guinier.I0 + sigma_I0, "I0 matches for a sphere") # Check RT invarients... rt = calc_Rambo_Tainer(data, guinier) - self.assertIsNotNone( - rt, "Rambo-Tainer invariants are actually calculated" - ) + self.assertIsNotNone(rt, "Rambo-Tainer invariants are actually calculated") def test_auto_gpa_with_outlier(self): if "outlier_position" not in self.extra_arg: @@ -266,7 +260,6 @@ def test_linspace(self): ) def test_random(self): - """ Tests that our linear regression implementation gives the same results as scipy.stats for random data @@ -288,7 +281,7 @@ def test_random(self): ) self.assertAlmostEqual( fit_result.R2, - ref.rvalue ** 2, + ref.rvalue**2, 5, "R² value matcheswihtin 4(?) digits", ) @@ -378,7 +371,7 @@ def test_curate_synthetic_data(self): def test_curate_synthetic_data_with_negative_points(self): if "negative_point_index" not in self.extra_arg: raise unittest.SkipTest("No negative point index to test") - + """Test that if one of the first three points is negative, all date before it gets ignored.""" negative_point_index = self.extra_arg["negative_point_index"] @@ -415,11 +408,10 @@ def test_curate_synthetic_data_with_negative_points(self): ) self.assertTrue( - data[offsets[data_range[1]] - 1, 1] - > data[negative_point_index + 1, 1] / 10 + data[offsets[data_range[1]] - 1, 1] > data[negative_point_index + 1, 1] / 10 and data[offsets[data_range[1]] + 1, 1] < data[negative_point_index + 1, 1] / 10, - msg=f"curated data for artificial data ends at approx. I(point after negaitve point)/10 if negative point at {negative_point_index + 1}", + msg=f"curated data for artificial data ends at approx. I(point after negative point)/10 if negative point at {negative_point_index + 1}", ) @@ -430,9 +422,7 @@ def suite(): testSuite.addTest(TestAutoRg("test_synthetic")) for outlier_position in range(3): testSuite.addTest( - TestAutoRg( - "test_auto_gpa_with_outlier", outlier_position=outlier_position - ) + TestAutoRg("test_auto_gpa_with_outlier", outlier_position=outlier_position) ) testSuite.addTest(TestFit("test_linear_fit_static")) testSuite.addTest(TestFit("test_linspace")) diff --git a/src/freesas/test/test_average.py b/src/freesas/test/test_average.py index 6498637..33873c0 100644 --- a/src/freesas/test/test_average.py +++ b/src/freesas/test/test_average.py @@ -12,6 +12,7 @@ from ..average import Grid, AverModels import logging + logging.basicConfig(level=logging.INFO) logger = logging.getLogger("AlignModels_test") @@ -45,13 +46,20 @@ def test_knots(self): grid.calc_radius(nbknots) grid.make_grid() gap = (1.0 * (grid.nbknots - nbknots) / nbknots) * 100 - self.assertGreater(threshold, gap, msg="final number of knots too different of wanted number: %s != %s" % (nbknots, grid.nbknots)) + self.assertGreater( + threshold, + gap, + msg="final number of knots too different of wanted number: %s != %s" + % (nbknots, grid.nbknots), + ) def test_makegrid(self): grid = self.grid lattice = grid.make_grid() m = SASModel(lattice) - self.assertAlmostEqual(m.fineness, 2 * grid.radius, 10, msg="grid do not have the computed radius") + self.assertAlmostEqual( + m.fineness, 2 * grid.radius, 10, msg="grid do not have the computed radius" + ) def test_read(self): inputfiles = self.inputfiles @@ -70,7 +78,9 @@ def test_occupancy(self): average.grid = occ_grid assert occ_grid.shape[-1] == 5, "problem in grid shape" diff = occ_grid[:-1, 3] - occ_grid[1:, 3] - self.assertTrue(diff.max() >= 0.0, msg="grid is not properly sorted with occupancy") + self.assertTrue( + diff.max() >= 0.0, msg="grid is not properly sorted with occupancy" + ) def suite(): @@ -83,6 +93,6 @@ def suite(): return testSuite -if __name__ == '__main__': +if __name__ == "__main__": runner = unittest.TextTestRunner() runner.run(suite()) diff --git a/src/freesas/test/test_bift.py b/src/freesas/test/test_bift.py index 36c3043..e6918fd 100644 --- a/src/freesas/test/test_bift.py +++ b/src/freesas/test/test_bift.py @@ -23,28 +23,32 @@ # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN # THE SOFTWARE. -__authors__ = ["J. Kieffer"] +__authors__ = ["Jérôme Kieffer"] __license__ = "MIT" -__date__ = "31/05/2024" +__date__ = "06/02/2026" +import logging +import time import numpy import unittest -from .utilstest import get_datafile from ..bift import auto_bift -from .._bift import BIFT, distribution_parabola, distribution_sphere, \ - ensure_edges_zero, smooth_density -import logging -logger = logging.getLogger(__name__) -import time +from .._bift import ( + BIFT, + distribution_parabola, + distribution_sphere, + ensure_edges_zero, + smooth_density, +) try: from numpy import trapezoid # numpy 2 except ImportError: from numpy import trapz as trapezoid # numpy1 +logger = logging.getLogger(__name__) -class TestBIFT(unittest.TestCase): +class TestBIFT(unittest.TestCase): DMAX = 10 NPT = 100 SIZE = 1000 @@ -57,13 +61,15 @@ def setUpClass(cls): cls.p = -cls.r * (cls.r - cls.DMAX) # Nice parabola q = numpy.linspace(0, 8 * cls.DMAX / 3, cls.SIZE + 1) sincqr = numpy.sinc(numpy.outer(q, cls.r / numpy.pi)) - I = 4 * numpy.pi * (cls.p * sincqr).sum(axis=-1) * dr - err = numpy.sqrt(I) - cls.I0 = I[0] + intensity = 4 * numpy.pi * (cls.p * sincqr).sum(axis=-1) * dr + err = numpy.sqrt(intensity) + cls.I0 = intensity[0] cls.q = q[1:] - cls.I = I[1:] + cls.I = intensity[1:] cls.err = err[1:] - cls.Rg = numpy.sqrt(0.5 * trapezoid(cls.p * cls.r ** 2, cls.r) / trapezoid(cls.p, cls.r)) + cls.Rg = numpy.sqrt( + 0.5 * trapezoid(cls.p * cls.r**2, cls.r) / trapezoid(cls.p, cls.r) + ) print(cls.Rg) @classmethod @@ -76,9 +82,9 @@ def test_autobift(self): t0 = time.perf_counter() bo = auto_bift(data) key, value, valid = bo.get_best() -# print("key is ", key) + # print("key is ", key) stats = bo.calc_stats() -# print("stat is ", stats) + # print("stat is ", stats) logger.info("Auto_bift time: %s", time.perf_counter() - t0) self.assertAlmostEqual(self.DMAX / key.Dmax, 1, 1, "DMax is correct") self.assertAlmostEqual(self.I0 / stats.I0_avg, 1, 1, "I0 is correct") @@ -97,31 +103,68 @@ def test_BIFT(self): def test_disributions(self): pp = numpy.asarray(distribution_parabola(self.I0, self.DMAX, self.NPT)) ps = numpy.asarray(distribution_sphere(self.I0, self.DMAX, self.NPT)) - self.assertAlmostEqual(trapezoid(ps, self.r) * 4 * numpy.pi / self.I0, 1, 3, "Distribution for a sphere looks OK") - self.assertAlmostEqual(trapezoid(pp, self.r) * 4 * numpy.pi / self.I0, 1, 3, "Distribution for a parabola looks OK") + self.assertAlmostEqual( + trapezoid(ps, self.r) * 4 * numpy.pi / self.I0, + 1, + 3, + "Distribution for a sphere looks OK", + ) + self.assertAlmostEqual( + trapezoid(pp, self.r) * 4 * numpy.pi / self.I0, + 1, + 3, + "Distribution for a parabola looks OK", + ) self.assertTrue(numpy.allclose(pp, self.p, 1e-4), "distribution matches") def test_fixEdges(self): ones = numpy.ones(self.NPT) ensure_edges_zero(ones) - self.assertAlmostEqual(ones[0], 0, msg="1st point set to 0") - self.assertAlmostEqual(ones[-1], 0, msg="last point set to 0") - self.assertTrue(numpy.allclose(ones[1:-1], numpy.ones(self.NPT-2), 1e-7), msg="non-edge points unchanged") + self.assertAlmostEqual(ones[0], 0, msg="1st point set to 0") + self.assertAlmostEqual(ones[-1], 0, msg="last point set to 0") + self.assertTrue( + numpy.allclose(ones[1:-1], numpy.ones(self.NPT - 2), 1e-7), + msg="non-edge points unchanged", + ) def test_smoothing(self): ones = numpy.ones(self.NPT) empty = numpy.empty(self.NPT) - smooth_density(ones,empty) - self.assertTrue(numpy.allclose(ones, empty, 1e-7), msg="flat array smoothed into flat array") + smooth_density(ones, empty) + self.assertTrue( + numpy.allclose(ones, empty, 1e-7), msg="flat array smoothed into flat array" + ) random = numpy.random.rand(self.NPT) - smooth = numpy.empty(self.NPT) - smooth_density(random,smooth) - self.assertAlmostEqual(random[0], smooth[0], msg="first points of random array and smoothed random array match") - self.assertAlmostEqual(random[-1], smooth[-1], msg="last points of random array and smoothed random array match") - self.assertTrue(smooth[1]>=min(smooth[0], smooth[2]) and smooth[1]<=max(smooth[0], smooth[2]), msg="second point of random smoothed array between 1st and 3rd") - self.assertTrue(smooth[-2]>=min(smooth[-1], smooth[-3]) and smooth[-2]<= max(smooth[-1], smooth[-3]), msg="second to last point of random smoothed array between 3rd to last and last") - sign = numpy.sign(random[1:-3] - smooth[2:-2]) * numpy.sign(smooth[2:-2] - random[3:-1]) - self.assertTrue(numpy.allclose(sign, numpy.ones(self.NPT-4), 1e-7), msg="central points of random array and smoothed random array alternate") + smooth = numpy.empty(self.NPT) + smooth_density(random, smooth) + self.assertAlmostEqual( + random[0], + smooth[0], + msg="first points of random array and smoothed random array match", + ) + self.assertAlmostEqual( + random[-1], + smooth[-1], + msg="last points of random array and smoothed random array match", + ) + self.assertTrue( + smooth[1] >= min(smooth[0], smooth[2]) + and smooth[1] <= max(smooth[0], smooth[2]), + msg="second point of random smoothed array between 1st and 3rd", + ) + self.assertTrue( + smooth[-2] >= min(smooth[-1], smooth[-3]) + and smooth[-2] <= max(smooth[-1], smooth[-3]), + msg="second to last point of random smoothed array between 3rd to last and last", + ) + sign = numpy.sign(random[1:-3] - smooth[2:-2]) * numpy.sign( + smooth[2:-2] - random[3:-1] + ) + self.assertTrue( + numpy.allclose(sign, numpy.ones(self.NPT - 4), 1e-7), + msg="central points of random array and smoothed random array alternate", + ) + def suite(): testSuite = unittest.TestSuite() @@ -132,6 +175,6 @@ def suite(): return testSuite -if __name__ == '__main__': +if __name__ == "__main__": runner = unittest.TextTestRunner() runner.run(suite()) diff --git a/src/freesas/test/test_cormap.py b/src/freesas/test/test_cormap.py index e4f56ae..c026631 100644 --- a/src/freesas/test/test_cormap.py +++ b/src/freesas/test/test_cormap.py @@ -1,21 +1,21 @@ #!/usr/bin/python # coding: utf-8 -__author__ = "Jerome" +__author__ = "Jérôme Kieffer" __license__ = "MIT" -__copyright__ = "2017, ESRF" +__copyright__ = "2017-2026, ESRF" import numpy import unittest - import logging +from .. import cormap + + logging.basicConfig(level=logging.INFO) logger = logging.getLogger("test_cormap") -from .. import cormap class TestCormap(unittest.TestCase): - def test_longest(self): size = 1000 target = 50 @@ -30,21 +30,25 @@ def test_longest(self): self.assertEqual(res, size, msg="computed size is correct: negative") data[:] = 0 - data[start: start + target] = 1.0 + data[start : start + target] = 1.0 res = cormap.measure_longest(data) self.assertEqual(res, target, msg="computed size is correct: positive/zero") data = numpy.zeros(size, dtype="float32") - data[start: start + target] = -1.0 + data[start : start + target] = -1.0 res = cormap.measure_longest(data) self.assertEqual(res, target, msg="computed size is correct: negative/zero") - data = numpy.fromfunction(lambda n:(-1) ** n, (size,)) - data[start: start + target] = 1.0 + data = numpy.fromfunction(lambda n: (-1) ** n, (size,)) + data[start : start + target] = 1.0 res = cormap.measure_longest(data) - self.assertEqual(res, target + 1, msg="computed size is correct: positive/alternating") - data = numpy.fromfunction(lambda n:(-1) ** n, (size,)) - data[start: start + target] = -1.0 + self.assertEqual( + res, target + 1, msg="computed size is correct: positive/alternating" + ) + data = numpy.fromfunction(lambda n: (-1) ** n, (size,)) + data[start : start + target] = -1.0 res = cormap.measure_longest(data) - self.assertEqual(res, target + 1, msg="computed size is correct: negative/alternating") + self.assertEqual( + res, target + 1, msg="computed size is correct: negative/alternating" + ) def test_stats(self): self.assertEqual(cormap.LROH.A(10, 0), 1) @@ -82,6 +86,6 @@ def suite(): return testSuite -if __name__ == '__main__': +if __name__ == "__main__": runner = unittest.TextTestRunner() runner.run(suite()) diff --git a/src/freesas/test/test_distance.py b/src/freesas/test/test_distance.py index 9287fa6..dea2aa5 100644 --- a/src/freesas/test/test_distance.py +++ b/src/freesas/test/test_distance.py @@ -11,6 +11,7 @@ from .utilstest import get_datafile from ..model import SASModel import logging + logging.basicConfig(level=logging.INFO) logger = logging.getLogger("cdistance_test") @@ -24,7 +25,9 @@ def test_invariants(self): m.read(self.testfile1) f_np, r_np, d_np = m.calc_invariants(False) f_cy, r_cy, d_cy = m.calc_invariants(True) - self.assertAlmostEqual(f_np, f_cy, 10, "fineness is the same %s!=%s" % (f_np, f_cy)) + self.assertAlmostEqual( + f_np, f_cy, 10, "fineness is the same %s!=%s" % (f_np, f_cy) + ) self.assertAlmostEqual(r_np, r_cy, 10, "Rg is the same %s!=%s" % (r_np, r_cy)) self.assertAlmostEqual(d_np, d_cy, 10, "Dmax is the same %s!=%s" % (d_np, d_cy)) @@ -35,7 +38,9 @@ def test_distance(self): n.read(self.testfile2) f_np = m.dist(n, m.atoms, n.atoms, False) f_cy = m.dist(n, m.atoms, n.atoms, True) - self.assertAlmostEqual(f_np, f_cy, 10, "distance is the same %s!=%s" % (f_np, f_cy)) + self.assertAlmostEqual( + f_np, f_cy, 10, "distance is the same %s!=%s" % (f_np, f_cy) + ) def test_same(self): m = SASModel() @@ -56,6 +61,7 @@ def suite(): testSuite.addTest(TestDistance("test_same")) return testSuite -if __name__ == '__main__': + +if __name__ == "__main__": runner = unittest.TextTestRunner() runner.run(suite()) diff --git a/src/freesas/test/test_dnn.py b/src/freesas/test/test_dnn.py index 046c936..f404619 100644 --- a/src/freesas/test/test_dnn.py +++ b/src/freesas/test/test_dnn.py @@ -25,67 +25,80 @@ __authors__ = ["Jérôme Kieffer"] __license__ = "MIT" -__date__ = "03/07/2024" +__date__ = "06/02/2026" import unittest import logging -import os -import io import numpy as np from .utilstest import get_datafile -from ..resources import resource_filename from ..sasio import load_scattering_data -from ..dnn import * +from ..dnn import ( + DNN, + DenseLayer, + forward_propagation, + preprocess, + tanh, + relu, + sigmoid, + linear, +) + logger = logging.getLogger(__name__) -class TestDNN(unittest.TestCase): +class TestDNN(unittest.TestCase): def test_activation_functions(self): """ Test for the activation functions """ x = np.array([-1, 0, 1]) - + # Test tanh expected_tanh = np.tanh(x) self.assertTrue(np.allclose(tanh(x), expected_tanh), msg="tanh function failed") - + # Test relu expected_relu = np.maximum(0, x) self.assertTrue(np.allclose(relu(x), expected_relu), msg="relu function failed") - + # Test sigmoid expected_sigmoid = 1 / (1 + np.exp(-x)) - self.assertTrue(np.allclose(sigmoid(x), expected_sigmoid), msg="sigmoid function failed") - + self.assertTrue( + np.allclose(sigmoid(x), expected_sigmoid), msg="sigmoid function failed" + ) + # Test linear expected_linear = x - self.assertTrue(np.allclose(linear(x), expected_linear), msg="linear function failed") - + self.assertTrue( + np.allclose(linear(x), expected_linear), msg="linear function failed" + ) + logger.info("test_activation_functions ran successfully") - def test_preprocess(self): """ - Test for the preprocessing function + Test for the preprocessing function """ datfile = get_datafile("bsa_005_sub.dat") data = load_scattering_data(datfile) - q, I, sigma = data.T - Iprep = preprocess(q, I) + q, intensity, sigma = data.T + Iprep = preprocess(q, intensity) self.assertEqual(Iprep.max(), 1, msg="range 0-1") self.assertEqual(Iprep.shape, (1024,), msg="size 1024") - - def test_forward_propagation(self): """ Test for the forward_propagation function """ try: X = np.random.rand(1, 10) - params = [np.random.rand(10, 20), np.random.rand(20), np.random.rand(20, 10), np.random.rand(10)] + params = [ + np.random.rand(10, 20), + np.random.rand(20), + np.random.rand(20, 10), + np.random.rand(10), + ] activations = [np.tanh, np.tanh] output = forward_propagation(X, params, activations) self.assertEqual(output.shape, (1, 10)) @@ -94,15 +107,14 @@ def test_forward_propagation(self): logger.error(f"test_forward_propagation failed: {e}") raise - def test_DenseLayer(self): """ Test for the DenseLayer class """ - try : + try: weights = np.random.rand(10, 20) bias = np.random.rand(20) - layer = DenseLayer(weights, bias, 'tanh') + layer = DenseLayer(weights, bias, "tanh") self.assertEqual(layer.input_size, 10) self.assertEqual(layer.output_size, 20) output = layer.forward(np.random.rand(1, 10)) @@ -112,15 +124,15 @@ def test_DenseLayer(self): logger.error(f"test_DenseLayer failed: {e}") raise - - def test_DNN(self): """ Test for the DNN class """ try: - layers = [DenseLayer(np.random.rand(10, 20), np.random.rand(20), 'tanh'), - DenseLayer(np.random.rand(20, 10), np.random.rand(10), 'tanh')] + layers = [ + DenseLayer(np.random.rand(10, 20), np.random.rand(20), "tanh"), + DenseLayer(np.random.rand(20, 10), np.random.rand(10), "tanh"), + ] dnn = DNN(*layers) output = dnn.infer(np.random.rand(1, 10)) self.assertEqual(output.shape, (1, 10)) @@ -138,6 +150,6 @@ def suite(): return test_suite -if __name__ == '__main__': +if __name__ == "__main__": runner = unittest.TextTestRunner() runner.run(suite()) diff --git a/src/freesas/test/test_fitting.py b/src/freesas/test/test_fitting.py index cb9d651..b9a24d8 100644 --- a/src/freesas/test/test_fitting.py +++ b/src/freesas/test/test_fitting.py @@ -96,7 +96,6 @@ def wrapped(*args, **kwargs): def build_mock_for_load_scattering_with_Errors(erronous_file: dict): - """Create mock for loading of data from a file. The resulting function will raise an error, for files for which an error is provided in errenous_file, @@ -106,9 +105,7 @@ def mock_for_load_scattering(file: pathlib.Path): if file.name in erronous_file: raise erronous_file[file.name] else: - return numpy.array( - [[1.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]] - ) + return numpy.array([[1.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]]) return mock_for_load_scattering @@ -117,7 +114,6 @@ class TestFitting(unittest.TestCase): def test_set_logging_level_does_not_change_logging_level_if_input_lower_1( self, ): - """ Test that the logging level only gets changed if the requested level is > 0. """ @@ -141,7 +137,6 @@ def test_set_logging_level_does_not_change_logging_level_if_input_lower_1( def test_set_logging_level_sets_logging_to_INFO_if_input_is_1( self, ): - """ Test that the logging level gets changed to INFO if verbosity is 1. """ @@ -162,7 +157,6 @@ def test_set_logging_level_sets_logging_to_INFO_if_input_is_1( def test_set_logging_level_sets_logging_to_DEBUG_if_input_is_2_or_more( self, ): - """ Test that the logging level gets changed to DEBUG if verbosity is 2 or larger. """ @@ -189,7 +183,6 @@ def test_set_logging_level_sets_logging_to_DEBUG_if_input_is_2_or_more( @patch.dict("sys.modules", {"nt": MagicMock()}) def test_get_linesep_returns_rn_if_output_is_stdout_on_windows(self): - """ Test that get_linesep() returns \r\n if output destination is sys.stdout on Windows. """ @@ -206,7 +199,6 @@ def test_get_linesep_returns_rn_if_output_is_stdout_on_windows(self): def test_get_linesep_returns_n_if_output_is_stdout_on_posix( self, ): - """ Test that get_linesep() returns \n if output destination is sys.stdout on Posix. Only should run on posix. @@ -215,7 +207,6 @@ def test_get_linesep_returns_n_if_output_is_stdout_on_posix( @patch.dict("sys.modules", {"nt": MagicMock()}) def test_get_linesep_returns_n_if_output_is_filestream_on_windows(self): - """ Test that get_linesep() returns \n if output destination is a filestream on Windows. """ @@ -232,7 +223,6 @@ def test_get_linesep_returns_n_if_output_is_filestream_on_windows(self): def test_get_linesep_returns_n_if_output_is_filestream_on_posix( self, ): - """ Test that get_linesep() returns \n if output destination is filestream on Posix. Only should run on posix. @@ -243,7 +233,6 @@ def test_get_linesep_returns_n_if_output_is_filestream_on_posix( def test_get_output_destination_with_path_input_returns_writable_io( self, ): - """Test that by calling get_output_destination with a Path as input we obtain write access to the file of Path.""" mocked_open = mock_open() @@ -258,7 +247,6 @@ def test_get_output_destination_with_path_input_returns_writable_io( def test_get_output_destination_without_input_returns_stdout( self, ): - """Test that by calling get_output_destination without input we obtain sys.stdout.""" with get_output_destination() as destination: @@ -271,7 +259,6 @@ def test_get_output_destination_without_input_returns_stdout( def test_closing_get_output_destination_does_not_close_stdout( self, ): - """Test that get_output_destination() can be safely used without closing sys.stdout.""" with get_output_destination() as _: @@ -288,7 +275,6 @@ def test_closing_get_output_destination_does_not_close_stdout( def test_get_guinier_header_for_csv( self, ): - """Test that by calling get_guinier_header with input csv we get the correct line.""" header = get_guinier_header("linesep", "csv") @@ -301,7 +287,6 @@ def test_get_guinier_header_for_csv( def test_get_guinier_header_for_ssv( self, ): - """Test that by calling get_guinier_header with input ssv we get an empty string.""" header = get_guinier_header("linesep", "ssv") @@ -314,7 +299,6 @@ def test_get_guinier_header_for_ssv( def test_get_guinier_header_for_native( self, ): - """Test that by calling get_guinier_header with input native we get an empty string.""" header = get_guinier_header("linesep", "native") @@ -327,7 +311,6 @@ def test_get_guinier_header_for_native( def test_get_guinier_header_without_input_format( self, ): - """Test that by calling get_guinier_header without input format we get an empty string.""" header = get_guinier_header("linesep", None) @@ -339,7 +322,6 @@ def test_get_guinier_header_without_input_format( @unittest.skipIf(platform.system() == "Windows", "Only POSIX") def test_collect_files_only_returns_existing_files(self): - """Test that collect_files discards strings that do not match an existing file.""" def os_stat_mock(path, **_): @@ -367,7 +349,6 @@ def os_stat_mock(path, **_): @patch("platform.system", MagicMock(return_value="Windows")) def test_collect_files_globs_on_windows(self): - """Test that collect_files globs on Windows if no existent files provided.""" def os_stat_mock(path): @@ -378,9 +359,7 @@ def os_stat_mock(path): mocked_stat = MagicMock(side_effect=os_stat_mock) mocked_glob = MagicMock( - side_effect=[ - (p for p in [pathlib.Path("pathA"), pathlib.Path("pathB")]) - ] + side_effect=[(p for p in [pathlib.Path("pathA"), pathlib.Path("pathB")])] ) with patch("os.stat", mocked_stat): with patch.object(pathlib.Path, "glob", mocked_glob): @@ -397,11 +376,12 @@ def os_stat_mock(path): reload_os_and_fitting() def test_rg_result_line_csv(self): - """Test the formatting of a csv result line for a Guinier fit.""" test_result = RG_RESULT(3.1, 0.1, 103, 2.5, 13, 207, 50.1, 0.05) - expected_line = "test.file,3.1000,0.1000,103.0000,2.5000, 13,207,50.1000,0.0500lineend" + expected_line = ( + "test.file,3.1000,0.1000,103.0000,2.5000, 13,207,50.1000,0.0500lineend" + ) obtained_line = rg_result_to_output_line( rg_result=test_result, afile=pathlib.Path("test.file"), @@ -413,11 +393,12 @@ def test_rg_result_line_csv(self): ) def test_rg_result_line_ssv(self): - """Test the formatting of a ssv result line for a Guinier fit.""" test_result = RG_RESULT(3.1, 0.1, 103, 2.5, 13, 207, 50.1, 0.05) - expected_line = "3.1000 0.1000 103.0000 2.5000 13 207 50.1000 0.0500 test.filelineend" + expected_line = ( + "3.1000 0.1000 103.0000 2.5000 13 207 50.1000 0.0500 test.filelineend" + ) obtained_line = rg_result_to_output_line( rg_result=test_result, afile=pathlib.Path("test.file"), @@ -445,7 +426,6 @@ def test_rg_result_line_native(self): ) def test_rg_result_line_no_format(self): - """Test the formatting of a native result line for a Guinier fit.""" test_result = RG_RESULT(3.1, 0.1, 103, 2.5, 13, 207, 50.1, 0.05) @@ -499,7 +479,9 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: parser=dummy_parser, logger=logger, ) - expected_output = "test Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" + expected_output = ( + "test Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" + ) self.assertEqual( output_catcher.getvalue(), expected_output, @@ -528,7 +510,6 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: MagicMock(return_value="linesep"), ) def test_run_guinier_fit_iterates_over_files(self): - """Test that run_guinier_fit calls the provided fit function for each provided file.""" @counted @@ -550,8 +531,7 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: logger=logger, ) expected_output = ( - "test Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" - * 2 + "test Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" * 2 ) self.assertEqual( output_catcher.getvalue(), @@ -605,7 +585,9 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: parser=dummy_parser, logger=logger, ) - expected_stdout_output = "test2 Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" + expected_stdout_output = ( + "test2 Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" + ) self.assertEqual( output_catcher_stdout.getvalue(), expected_stdout_output, @@ -636,7 +618,6 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: @patch_collect_files @patch_linesep def test_run_guinier_outputs_error_if_file_not_parsable(self): - """Test that run_guinier_fit outputs an error if data loading raises ValueError and continues to the next file.""" @@ -660,7 +641,9 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: parser=dummy_parser, logger=logger, ) - expected_stdout_output = "test2 Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" + expected_stdout_output = ( + "test2 Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" + ) self.assertEqual( output_catcher_stdout.getvalue(), expected_stdout_output, @@ -683,12 +666,8 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: "freesas.fitting.load_scattering_data", MagicMock( side_effect=[ - numpy.array( - [[0.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]] - ), - numpy.array( - [[2.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]] - ), + numpy.array([[0.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]]), + numpy.array([[2.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]]), ] ), ) @@ -697,7 +676,6 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: def test_run_guinier_outputs_error_if_fitting_raises_insufficient_data_error( self, ): - """Test that run_guinier_fit outputs an error if fitting raises InsufficientDataError and continues to the next file.""" @@ -723,7 +701,9 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: parser=dummy_parser, logger=logger, ) - expected_stdout_output = "test2 Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" + expected_stdout_output = ( + "test2 Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" + ) self.assertEqual( output_catcher_stdout.getvalue(), expected_stdout_output, @@ -744,12 +724,8 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: "freesas.fitting.load_scattering_data", MagicMock( side_effect=[ - numpy.array( - [[0.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]] - ), - numpy.array( - [[2.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]] - ), + numpy.array([[0.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]]), + numpy.array([[2.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]]), ] ), ) @@ -758,7 +734,6 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: def test_run_guinier_outputs_error_if_fitting_raises_no_guinier_region_error( self, ): - """Test that run_guinier_fit outputs an error if fitting raises NoGuinierRegionError and continues to the next file.""" @@ -784,7 +759,9 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: parser=dummy_parser, logger=logger, ) - expected_stdout_output = "test2 Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" + expected_stdout_output = ( + "test2 Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" + ) self.assertEqual( output_catcher_stdout.getvalue(), expected_stdout_output, @@ -805,12 +782,8 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: "freesas.fitting.load_scattering_data", MagicMock( side_effect=[ - numpy.array( - [[0.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]] - ), - numpy.array( - [[2.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]] - ), + numpy.array([[0.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]]), + numpy.array([[2.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]]), ] ), ) @@ -819,7 +792,6 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: def test_run_guinier_outputs_error_if_fitting_raises_value_error( self, ): - """Test that run_guinier_fit outputs an error if fitting raises ValueError and continues to the next file.""" @@ -845,7 +817,9 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: parser=dummy_parser, logger=logger, ) - expected_stdout_output = "test2 Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" + expected_stdout_output = ( + "test2 Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" + ) self.assertEqual( output_catcher_stdout.getvalue(), expected_stdout_output, @@ -866,12 +840,8 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: "freesas.fitting.load_scattering_data", MagicMock( side_effect=[ - numpy.array( - [[0.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]] - ), - numpy.array( - [[2.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]] - ), + numpy.array([[0.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]]), + numpy.array([[2.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]]), ] ), ) @@ -880,7 +850,6 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: def test_run_guinier_outputs_error_if_fitting_raises_index_error( self, ): - """Test that run_guinier_fit outputs an error if fitting raises IndexError and continues to the next file.""" @@ -906,7 +875,9 @@ def dummy_fit_function(input_data: numpy.ndarray) -> RG_RESULT: parser=dummy_parser, logger=logger, ) - expected_stdout_output = "test2 Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" + expected_stdout_output = ( + "test2 Rg=3.1000(±0.1000) I0=103.0000(±2.5000) [13-207] 5010.00% linesep" + ) self.assertEqual( output_catcher_stdout.getvalue(), expected_stdout_output, diff --git a/src/freesas/test/test_model.py b/src/freesas/test/test_model.py index 2bd924c..c95b2b0 100644 --- a/src/freesas/test/test_model.py +++ b/src/freesas/test/test_model.py @@ -13,6 +13,7 @@ from ..model import SASModel from ..transformations import translation_from_matrix, euler_from_matrix import logging + logging.basicConfig(level=logging.INFO) logger = logging.getLogger("SASModel_test") @@ -49,8 +50,10 @@ def test_same(self): m = SASModel() m.read(self.testfile) m.save(self.outfile) - with open(self.testfile) as f: infile=f.read() - with open(self.outfile) as f: outfile=f.read() + with open(self.testfile) as f: + infile = f.read() + with open(self.outfile) as f: + outfile = f.read() self.assertEqual(infile, outfile, msg="file content is the same") def test_rfactor(self): @@ -58,7 +61,11 @@ def test_rfactor(self): m.read(self.testfile) n = SASModel() n.read(self.testfile) - self.assertEqual(m.rfactor, n.rfactor, msg="R-factor is not the same %s != %s" % (m.rfactor, n.rfactor)) + self.assertEqual( + m.rfactor, + n.rfactor, + msg="R-factor is not the same %s != %s" % (m.rfactor, n.rfactor), + ) def test_init(self): m = SASModel() @@ -72,17 +79,29 @@ def test_centroid(self): m = assign_random_mol() m.centroid() if len(m.com) != 3: - logger.error("center of mass has not been saved correctly : length of COM position vector = %s!=3" % (len(m.com))) + logger.error( + "center of mass has not been saved correctly : length of COM position vector = %s!=3" + % (len(m.com)) + ) mol_centered = m.atoms[:, 0:3] - m.com center = mol_centered.mean(axis=0) norm = (center * center).sum() - self.assertAlmostEqual(norm, 0, 12, msg="molecule is not centered : norm of the COM position vector %s!=0" % (norm)) + self.assertAlmostEqual( + norm, + 0, + 12, + msg="molecule is not centered : norm of the COM position vector %s!=0" + % (norm), + ) def test_inertia_tensor(self): m = assign_random_mol() m.inertiatensor() tensor = m.inertensor - assert tensor.shape == (3, 3), "inertia tensor has not been saved correctly : shape of inertia matrix = %s" % (tensor.shape) + assert tensor.shape == (3, 3), ( + "inertia tensor has not been saved correctly : shape of inertia matrix = %s" + % (tensor.shape) + ) def test_canonical_translate(self): m = assign_random_mol() @@ -92,7 +111,9 @@ def test_canonical_translate(self): com = m.com com_componants = [com[0], com[1], com[2]] trans_vect = [-trans[0, -1], -trans[1, -1], -trans[2, -1]] - self.assertEqual(com_componants, trans_vect, msg="do not translate on canonical position") + self.assertEqual( + com_componants, trans_vect, msg="do not translate on canonical position" + ) def test_canonical_rotate(self): m = assign_random_mol() @@ -102,7 +123,9 @@ def test_canonical_rotate(self): if not m.enantiomer: logger.error("enantiomer has not been selected") det = numpy.linalg.det(rot) - self.assertAlmostEqual(det, 1, 10, msg="rotation matrix determinant is not 1: %s" % (det)) + self.assertAlmostEqual( + det, 1, 10, msg="rotation matrix determinant is not 1: %s" % (det) + ) def test_canonical_parameters(self): m = assign_random_mol() @@ -112,8 +135,17 @@ def test_canonical_parameters(self): logger.error("canonical parameters has not been saved properly") com_trans = translation_from_matrix(m.canonical_translate()) euler_rot = euler_from_matrix(m.canonical_rotate()) - out_param = [com_trans[0], com_trans[1], com_trans[2], euler_rot[0], euler_rot[1], euler_rot[2]] - self.assertEqual(can_param, out_param, msg="canonical parameters are not the good ones") + out_param = [ + com_trans[0], + com_trans[1], + com_trans[2], + euler_rot[0], + euler_rot[1], + euler_rot[2], + ] + self.assertEqual( + can_param, out_param, msg="canonical parameters are not the good ones" + ) def test_dist(self): m = assign_random_mol() @@ -135,7 +167,9 @@ def test_can_transform(self): tensor = m.inertensor diag = numpy.eye(3) matrix = tensor - tensor * diag - self.assertAlmostEqual(abs(com).sum(), 0, 10, msg="molecule not on its center of mass") + self.assertAlmostEqual( + abs(com).sum(), 0, 10, msg="molecule not on its center of mass" + ) self.assertAlmostEqual(abs(matrix).sum(), 0, 10, "inertia moments unaligned ") def test_dist_move(self): @@ -147,7 +181,9 @@ def test_dist_move(self): logger.error("molecules are different") p0 = m.can_param dist_after_mvt = m.dist_after_movement(p0, n, [1, 1, 1]) - self.assertEqual(dist_after_mvt, 0, msg="NSD different of 0: %s!=0" % (dist_after_mvt)) + self.assertEqual( + dist_after_mvt, 0, msg="NSD different of 0: %s!=0" % (dist_after_mvt) + ) def test_reverse_transform(self): m = assign_random_mol() @@ -156,7 +192,9 @@ def test_reverse_transform(self): m.atoms = m.transform(m.can_param, [1, 1, 1], reverse=None) m.atoms = m.transform(m.can_param, [1, 1, 1], reverse=True) dist = m.dist(n, m.atoms, n.atoms) - self.assertAlmostEqual(dist, 0.0, 10, msg="pb with reverse transformation : %s != 0.0" % dist) + self.assertAlmostEqual( + dist, 0.0, 10, msg="pb with reverse transformation : %s != 0.0" % dist + ) def suite(): @@ -176,6 +214,6 @@ def suite(): return testSuite -if __name__ == '__main__': +if __name__ == "__main__": runner = unittest.TextTestRunner() runner.run(suite()) diff --git a/src/freesas/test/test_resources.py b/src/freesas/test/test_resources.py index c660a8f..13bd3b4 100644 --- a/src/freesas/test/test_resources.py +++ b/src/freesas/test/test_resources.py @@ -31,16 +31,19 @@ import logging import os from ..resources import resource_filename + logger = logging.getLogger(__name__) -class TestResources(unittest.TestCase): +class TestResources(unittest.TestCase): def test_filename(self): """ - Test for returning the actual path of an existing model + Test for returning the actual path of an existing model """ - self.assertTrue(os.path.exists(resource_filename("keras_models/Rg+Dmax.keras")), - msg="file exists") + self.assertTrue( + os.path.exists(resource_filename("keras_models/Rg+Dmax.keras")), + msg="file exists", + ) def suite(): @@ -49,6 +52,6 @@ def suite(): return test_suite -if __name__ == '__main__': +if __name__ == "__main__": runner = unittest.TextTestRunner() runner.run(suite()) diff --git a/src/freesas/test/test_sas_argparser.py b/src/freesas/test/test_sas_argparser.py index c9a9b4d..713a074 100644 --- a/src/freesas/test/test_sas_argparser.py +++ b/src/freesas/test/test_sas_argparser.py @@ -5,7 +5,7 @@ __authors__ = ["Martha Brennich"] __license__ = "MIT" -__date__ = "25/03/2021" +__date__ = "06/02/2026" import unittest @@ -38,8 +38,7 @@ def test_minimal_guinier_parser_requires_file_argument(self): msg="GuinierParser provides usage if no file provided", ) self.assertTrue( - "the following arguments are required: FILE" - in output_catcher.getvalue(), + "the following arguments are required: FILE" in output_catcher.getvalue(), msg="GuinierParser states that the FILE argument is missing if no file provided", ) @@ -172,7 +171,9 @@ def test_minimal_parser_default_verbosity_level_is_0(self): """ Test that the parser sets the verbosity to 0 if no args are provided """ - basic_parser = SASParser("program", "description", "epilog", exit_on_error=False) + basic_parser = SASParser( + "program", "description", "epilog", exit_on_error=False + ) parsed_arguments = basic_parser.parse_args([]) self.assertEqual( parsed_arguments.verbose, @@ -434,9 +435,11 @@ def test_SASParser_q_unit_argument_does_not_allow_not_predefined_units( _ = basic_parser.parse_args(["-u", "m"]) except SystemExit: pass + output_catcher.seek(0) + output = output_catcher.read() + self.assertTrue( - "argument -u/--unit: invalid choice: 'm' (choose from 'nm', 'Å', 'A')" - in output_catcher.getvalue(), + "argument -u/--unit: invalid choice: 'm'" in output, msg="SASParser does not accept '-u m' argument", ) @@ -500,9 +503,12 @@ def test_GuinierParser_q_unit_argument_does_not_allow_not_predefined_units( _ = basic_parser.parse_args(["afile", "-u", "m"]) except SystemExit: pass + output_catcher.seek(0) + output = output_catcher.read() + # print(output) + valid = "argument -u/--unit: invalid choice: 'm'" in output self.assertTrue( - "argument -u/--unit: invalid choice: 'm' (choose from 'nm', 'Å', 'A')" - in output_catcher.getvalue(), + valid, msg="SASParser does not accept '-u m' argument", ) diff --git a/src/freesas/test/test_sasio.py b/src/freesas/test/test_sasio.py index c897b21..cc641fc 100644 --- a/src/freesas/test/test_sasio.py +++ b/src/freesas/test/test_sasio.py @@ -31,27 +31,34 @@ import logging import io from numpy import array, allclose -from ..sasio import parse_ascii_data, load_scattering_data, \ - convert_inverse_angstrom_to_nanometer +from ..sasio import ( + parse_ascii_data, + load_scattering_data, + convert_inverse_angstrom_to_nanometer, +) + logger = logging.getLogger(__name__) -class TestSasIO(unittest.TestCase): +class TestSasIO(unittest.TestCase): def test_parse_3_ok(self): """ Test for successful parsing of file with some invalid lines """ - file_content = ["Test data for", - "file parsing", - "1 1 1", - "2 a 2", - "3 3 3", - "some stuff at the end", - ] + file_content = [ + "Test data for", + "file parsing", + "1 1 1", + "2 a 2", + "3 3 3", + "some stuff at the end", + ] expected_result = array([[1.0, 1.0, 1.0], [3.0, 3.0, 3.0]]) data = parse_ascii_data(file_content, number_of_columns=3) - self.assertTrue(allclose(data, expected_result, 1e-7), - msg="3 column parse returns expected result") + self.assertTrue( + allclose(data, expected_result, 1e-7), + msg="3 column parse returns expected result", + ) def test_parse_no_data(self): """ @@ -66,71 +73,63 @@ def test_parse_no_valid_data(self): Test that an input list with no valid data raises a ValueError """ file_content = ["a a a", "2 4", "3 4 5 6", "# 3 4 6"] - with self.assertRaises(ValueError, - msg="File with no float float float data" - " cannot be parsed"): + with self.assertRaises( + ValueError, msg="File with no float float float data cannot be parsed" + ): parse_ascii_data(file_content, number_of_columns=3) def test_load_clean_data(self): """ Test that clean float float float data is loaded correctly. """ - file_content = ["# Test data for" - "# file parsing", - "1 1 1", - "2.0 2.0 1.0", - "3 3 3", - "#REMARK some stuff at the end", - ] - expected_result = array([[1.0, 1.0, 1.0], - [2.0, 2.0, 1.0], - [3.0, 3.0, 3.0]]) + file_content = [ + "# Test data for# file parsing", + "1 1 1", + "2.0 2.0 1.0", + "3 3 3", + "#REMARK some stuff at the end", + ] + expected_result = array([[1.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]]) file_data = "\n".join(file_content) mocked_file = io.StringIO(file_data) data = load_scattering_data(mocked_file) - self.assertTrue(allclose(data, expected_result, 1e-7), - msg="Sunny data loaded correctly") + self.assertTrue( + allclose(data, expected_result, 1e-7), msg="Sunny data loaded correctly" + ) def test_load_data_with_unescaped_header(self): """ Test that an unescaped header does not hinder loading. """ - file_content = ["Test data for", - "file parsing", - "1 1 1", - "2.0 2.0 1.0", - "3 3 3", - ] - expected_result = array([[1.0, 1.0, 1.0], - [2.0, 2.0, 1.0], - [3.0, 3.0, 3.0]]) + file_content = [ + "Test data for", + "file parsing", + "1 1 1", + "2.0 2.0 1.0", + "3 3 3", + ] + expected_result = array([[1.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]]) file_data = "\n".join(file_content) mocked_file = io.StringIO(file_data) data = load_scattering_data(mocked_file) - self.assertTrue(allclose(data, expected_result, 1e-7), - msg="Sunny data loaded correctly") + self.assertTrue( + allclose(data, expected_result, 1e-7), msg="Sunny data loaded correctly" + ) def test_load_data_with_unescaped_footer(self): """ Test that an unescaped footer does not hinder loading. """ - file_content = [ - "1 1 1", - "2.0 2.0 1.0", - "3 3 3", - "REMARK some stuff at the end" - ] - expected_result = array([[1.0, 1.0, 1.0], - [2.0, 2.0, 1.0], - [3.0, 3.0, 3.0]]) + file_content = ["1 1 1", "2.0 2.0 1.0", "3 3 3", "REMARK some stuff at the end"] + expected_result = array([[1.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]]) file_data = "\n".join(file_content) mocked_file = io.StringIO(file_data) data = load_scattering_data(mocked_file) - self.assertTrue(allclose(data, expected_result, 1e-7), - msg="Sunny data loaded correctly") - + self.assertTrue( + allclose(data, expected_result, 1e-7), msg="Sunny data loaded correctly" + ) def test_load_invalid_data(self): """ @@ -139,38 +138,33 @@ def test_load_invalid_data(self): file_content = ["a a a", "2 4", "3 4 5 6", "# 3 4 6"] file_data = "\n".join(file_content) mocked_file = io.StringIO(file_data) - with self.assertRaises(ValueError, - msg="File with no float float float " - "data cannot be loaded"): + with self.assertRaises( + ValueError, msg="File with no float float float data cannot be loaded" + ): load_scattering_data(mocked_file) def test_convert_inverse_angstrom_to_nanometer(self): """ Test conversion of data with q in 1/Å to 1/nm """ - input_data = array([[1.0, 1.0, 1.0], - [2.0, 2.0, 1.0], - [3.0, 3.0, 3.0]]) - expected_result = array([[10, 1.0, 1.0], - [20, 2.0, 1.0], - [30, 3.0, 3.0]]) + input_data = array([[1.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]]) + expected_result = array([[10, 1.0, 1.0], [20, 2.0, 1.0], [30, 3.0, 3.0]]) result = convert_inverse_angstrom_to_nanometer(input_data) - self.assertTrue(allclose(result, expected_result, 1e-7), - msg="Converted to 1/nm from 1/Å") + self.assertTrue( + allclose(result, expected_result, 1e-7), msg="Converted to 1/nm from 1/Å" + ) def test_unit_conversion_creates_new_array(self): """ Test conversion of data does not change original data """ - input_data = array([[1.0, 1.0, 1.0], - [2.0, 2.0, 1.0], - [3.0, 3.0, 3.0]]) - expected_data = array([[1.0, 1.0, 1.0], - [2.0, 2.0, 1.0], - [3.0, 3.0, 3.0]]) + input_data = array([[1.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]]) + expected_data = array([[1.0, 1.0, 1.0], [2.0, 2.0, 1.0], [3.0, 3.0, 3.0]]) _ = convert_inverse_angstrom_to_nanometer(input_data) - self.assertTrue(allclose(input_data, expected_data, 1e-7), - msg="Conversion function does not change its input") + self.assertTrue( + allclose(input_data, expected_data, 1e-7), + msg="Conversion function does not change its input", + ) def suite(): @@ -187,6 +181,6 @@ def suite(): return test_suite -if __name__ == '__main__': +if __name__ == "__main__": runner = unittest.TextTestRunner() runner.run(suite()) diff --git a/src/freesas/test/utilstest.py b/src/freesas/test/utilstest.py index 96e0506..df8efce 100644 --- a/src/freesas/test/utilstest.py +++ b/src/freesas/test/utilstest.py @@ -3,14 +3,17 @@ __author__ = "Jérôme Kieffer" __license__ = "MIT" -__date__ = "28/11/2023" -__copyright__ = "2015-2021, ESRF" +__date__ = "06/02/2026" +__copyright__ = "2015-2026, ESRF" import os import logging -logger = logging.getLogger("utilstest") from silx.resources import ExternalResources -downloader = ExternalResources("freesas", "http://www.silx.org/pub/freesas/testdata", "FREESAS_TESTDATA") + +logger = logging.getLogger("utilstest") +downloader = ExternalResources( + "freesas", "http://www.silx.org/pub/freesas/testdata", "FREESAS_TESTDATA" +) def get_datafile(name): @@ -24,13 +27,16 @@ def get_datafile(name): fullpath = downloader.getfile(name) return fullpath + def clean(): pass + class TestOptions: """ Class providing useful stuff for preparing tests. """ + def __init__(self): self.TEST_LOW_MEM = False """Skip tests using too much memory""" @@ -41,8 +47,8 @@ def __init__(self): self.options = None def __repr__(self): - return f"TEST_LOW_MEM={self.TEST_LOW_MEM} "\ - f"TEST_RANDOM={self.TEST_RANDOM} " + return f"TEST_LOW_MEM={self.TEST_LOW_MEM} TEST_RANDOM={self.TEST_RANDOM} " + @property def low_mem(self): """For compatibility""" @@ -55,26 +61,35 @@ def configure(self, parsed_options=None): if parsed_options is not None and parsed_options.low_mem: self.TEST_LOW_MEM = True - elif os.environ.get('FREESAS_LOW_MEM', 'True') == 'False': + elif os.environ.get("FREESAS_LOW_MEM", "True") == "False": self.TEST_LOW_MEM = True if parsed_options is not None and parsed_options.random: self.TEST_RANDOM = True - if os.environ.get('FREESAS_RANDOM', 'False').lower() in ("1", "true", "on"): + if os.environ.get("FREESAS_RANDOM", "False").lower() in ("1", "true", "on"): self.TEST_RANDOM = True - def add_parser_argument(self, parser): """Add extrat arguments to the test argument parser :param ArgumentParser parser: An argument parser """ - parser.add_argument("-l", "--low-mem", dest="low_mem", default=False, - action="store_true", - help="Disable test with large memory consumption (>100Mbyte") - parser.add_argument("-r", "--random", dest="random", default=False, - action="store_true", - help="Enable actual random number to be generated. By default, stable seed ensures reproducibility of tests") - - -test_options = TestOptions() #singleton + parser.add_argument( + "-l", + "--low-mem", + dest="low_mem", + default=False, + action="store_true", + help="Disable test with large memory consumption (>100Mbyte", + ) + parser.add_argument( + "-r", + "--random", + dest="random", + default=False, + action="store_true", + help="Enable actual random number to be generated. By default, stable seed ensures reproducibility of tests", + ) + + +test_options = TestOptions() # singleton diff --git a/src/freesas/transformations.py b/src/freesas/transformations.py index f474c3f..dfd5b58 100644 --- a/src/freesas/transformations.py +++ b/src/freesas/transformations.py @@ -197,8 +197,8 @@ import numpy -__version__ = '2015.03.19' -__docformat__ = 'restructuredtext en' +__version__ = "2015.03.19" +__docformat__ = "restructuredtext en" __all__ = () @@ -329,9 +329,13 @@ def rotation_matrix(angle, direction, point=None): R = numpy.diag([cosa, cosa, cosa]) R += numpy.outer(direction, direction) * (1.0 - cosa) direction *= sina - R += numpy.array([[ 0.0, -direction[2], direction[1]], - [ direction[2], 0.0, -direction[0]], - [-direction[1], direction[0], 0.0]]) + R += numpy.array( + [ + [0.0, -direction[2], direction[1]], + [direction[2], 0.0, -direction[0]], + [-direction[1], direction[0], 0.0], + ] + ) M = numpy.identity(4) M[:3, :3] = R if point is not None: @@ -372,11 +376,11 @@ def rotation_from_matrix(matrix): # rotation angle depending on direction cosa = (numpy.trace(R33) - 1.0) / 2.0 if abs(direction[2]) > 1e-8: - sina = (R[1, 0] + (cosa-1.0)*direction[0]*direction[1]) / direction[2] + sina = (R[1, 0] + (cosa - 1.0) * direction[0] * direction[1]) / direction[2] elif abs(direction[1]) > 1e-8: - sina = (R[0, 2] + (cosa-1.0)*direction[0]*direction[2]) / direction[1] + sina = (R[0, 2] + (cosa - 1.0) * direction[0] * direction[2]) / direction[1] else: - sina = (R[2, 1] + (cosa-1.0)*direction[1]*direction[2]) / direction[0] + sina = (R[2, 1] + (cosa - 1.0) * direction[1] * direction[2]) / direction[0] angle = math.atan2(sina, cosa) return angle, direction, point @@ -456,8 +460,7 @@ def scale_from_matrix(matrix): return factor, origin, direction -def projection_matrix(point, normal, direction=None, - perspective=None, pseudo=False): +def projection_matrix(point, normal, direction=None, perspective=None, pseudo=False): """Return matrix to project onto plane defined by point and normal. Using either perspective point, projection direction, or none of both. @@ -493,14 +496,13 @@ def projection_matrix(point, normal, direction=None, normal = unit_vector(normal[:3]) if perspective is not None: # perspective projection - perspective = numpy.array(perspective[:3], dtype=numpy.float64, - copy=False) - M[0, 0] = M[1, 1] = M[2, 2] = numpy.dot(perspective-point, normal) + perspective = numpy.array(perspective[:3], dtype=numpy.float64, copy=False) + M[0, 0] = M[1, 1] = M[2, 2] = numpy.dot(perspective - point, normal) M[:3, :3] -= numpy.outer(perspective, normal) if pseudo: # preserve relative depth M[:3, :3] -= numpy.outer(normal, normal) - M[:3, 3] = numpy.dot(point, normal) * (perspective+normal) + M[:3, 3] = numpy.dot(point, normal) * (perspective + normal) else: M[:3, 3] = numpy.dot(point, normal) * perspective M[3, :3] = -normal @@ -580,11 +582,10 @@ def projection_from_matrix(matrix, pseudo=False): # perspective projection i = numpy.where(abs(numpy.real(w)) > 1e-8)[0] if not len(i): - raise ValueError( - "no eigenvector not corresponding to eigenvalue 0") + raise ValueError("no eigenvector not corresponding to eigenvalue 0") point = numpy.real(V[:, i[-1]]).squeeze() point /= point[3] - normal = - M[3, :3] + normal = -M[3, :3] perspective = M[:3, 3] / numpy.dot(point[:3], normal) if pseudo: perspective -= normal @@ -631,15 +632,19 @@ def clip_matrix(left, right, bottom, top, near, far, perspective=False): if near <= _EPS: raise ValueError("invalid frustum: near <= 0") t = 2.0 * near - M = [[t/(left-right), 0.0, (right+left)/(right-left), 0.0], - [0.0, t/(bottom-top), (top+bottom)/(top-bottom), 0.0], - [0.0, 0.0, (far+near)/(near-far), t*far/(far-near)], - [0.0, 0.0, -1.0, 0.0]] + M = [ + [t / (left - right), 0.0, (right + left) / (right - left), 0.0], + [0.0, t / (bottom - top), (top + bottom) / (top - bottom), 0.0], + [0.0, 0.0, (far + near) / (near - far), t * far / (far - near)], + [0.0, 0.0, -1.0, 0.0], + ] else: - M = [[2.0/(right-left), 0.0, 0.0, (right+left)/(left-right)], - [0.0, 2.0/(top-bottom), 0.0, (top+bottom)/(bottom-top)], - [0.0, 0.0, 2.0/(far-near), (far+near)/(near-far)], - [0.0, 0.0, 0.0, 1.0]] + M = [ + [2.0 / (right - left), 0.0, 0.0, (right + left) / (left - right)], + [0.0, 2.0 / (top - bottom), 0.0, (top + bottom) / (bottom - top)], + [0.0, 0.0, 2.0 / (far - near), (far + near) / (near - far)], + [0.0, 0.0, 0.0, 1.0], + ] return numpy.array(M) @@ -759,7 +764,7 @@ def decompose_matrix(matrix): if not numpy.linalg.det(P): raise ValueError("matrix is singular") - scale = numpy.zeros((3, )) + scale = numpy.zeros((3,)) shear = [0.0, 0.0, 0.0] angles = [0.0, 0.0, 0.0] @@ -797,15 +802,16 @@ def decompose_matrix(matrix): angles[0] = math.atan2(row[1, 2], row[2, 2]) angles[2] = math.atan2(row[0, 1], row[0, 0]) else: - #angles[0] = math.atan2(row[1, 0], row[1, 1]) + # angles[0] = math.atan2(row[1, 0], row[1, 1]) angles[0] = math.atan2(-row[2, 1], row[1, 1]) angles[2] = 0.0 return scale, shear, angles, translate, perspective -def compose_matrix(scale=None, shear=None, angles=None, translate=None, - perspective=None): +def compose_matrix( + scale=None, shear=None, angles=None, translate=None, perspective=None +): """Return transformation matrix from sequence of transformations. This is the inverse of the decompose_matrix function. @@ -839,7 +845,7 @@ def compose_matrix(scale=None, shear=None, angles=None, translate=None, T[:3, 3] = translate[:3] M = numpy.dot(M, T) if angles is not None: - R = euler_matrix(angles[0], angles[1], angles[2], 'sxyz') + R = euler_matrix(angles[0], angles[1], angles[2], "sxyz") M = numpy.dot(M, R) if shear is not None: Z = numpy.identity(4) @@ -877,11 +883,14 @@ def orthogonalization_matrix(lengths, angles): sina, sinb, _ = numpy.sin(angles) cosa, cosb, cosg = numpy.cos(angles) co = (cosa * cosb - cosg) / (sina * sinb) - return numpy.array([ - [ a*sinb*math.sqrt(1.0-co*co), 0.0, 0.0, 0.0], - [-a*sinb*co, b*sina, 0.0, 0.0], - [ a*cosb, b*cosa, c, 0.0], - [ 0.0, 0.0, 0.0, 1.0]]) + return numpy.array( + [ + [a * sinb * math.sqrt(1.0 - co * co), 0.0, 0.0, 0.0], + [-a * sinb * co, b * sina, 0.0, 0.0], + [a * cosb, b * cosa, c, 0.0], + [0.0, 0.0, 0.0, 1.0], + ] + ) def affine_matrix_from_points(v0, v1, shear=True, scale=True, usesvd=True): @@ -934,11 +943,11 @@ def affine_matrix_from_points(v0, v1, shear=True, scale=True, usesvd=True): # move centroids to origin t0 = -numpy.mean(v0, axis=1) - M0 = numpy.identity(ndims+1) + M0 = numpy.identity(ndims + 1) M0[:ndims, ndims] = t0 v0 += t0.reshape(ndims, 1) t1 = -numpy.mean(v1, axis=1) - M1 = numpy.identity(ndims+1) + M1 = numpy.identity(ndims + 1) M1[:ndims, ndims] = t1 v1 += t1.reshape(ndims, 1) @@ -948,10 +957,10 @@ def affine_matrix_from_points(v0, v1, shear=True, scale=True, usesvd=True): u, s, vh = numpy.linalg.svd(A.T) vh = vh[:ndims].T B = vh[:ndims] - C = vh[ndims:2*ndims] + C = vh[ndims : 2 * ndims] t = numpy.dot(C, numpy.linalg.pinv(B)) t = numpy.concatenate((t, numpy.zeros((ndims, 1))), axis=1) - M = numpy.vstack((t, ((0.0,)*ndims) + (1.0,))) + M = numpy.vstack((t, ((0.0,) * ndims) + (1.0,))) elif usesvd or ndims != 3: # Rigid transformation via SVD of covariance matrix u, s, vh = numpy.linalg.svd(numpy.dot(v1, v0.T)) @@ -959,10 +968,10 @@ def affine_matrix_from_points(v0, v1, shear=True, scale=True, usesvd=True): R = numpy.dot(u, vh) if numpy.linalg.det(R) < 0.0: # R does not constitute right handed system - R -= numpy.outer(u[:, ndims-1], vh[ndims-1, :]*2.0) + R -= numpy.outer(u[:, ndims - 1], vh[ndims - 1, :] * 2.0) s[-1] *= -1.0 # homogeneous transformation matrix - M = numpy.identity(ndims+1) + M = numpy.identity(ndims + 1) M[:ndims, :ndims] = R else: # Rigid transformation matrix via quaternion @@ -970,10 +979,12 @@ def affine_matrix_from_points(v0, v1, shear=True, scale=True, usesvd=True): xx, yy, zz = numpy.sum(v0 * v1, axis=1) xy, yz, zx = numpy.sum(v0 * numpy.roll(v1, -1, axis=0), axis=1) xz, yx, zy = numpy.sum(v0 * numpy.roll(v1, -2, axis=0), axis=1) - N = [[xx+yy+zz, 0.0, 0.0, 0.0], - [yz-zy, xx-yy-zz, 0.0, 0.0], - [zx-xz, xy+yx, yy-xx-zz, 0.0], - [xy-yx, zx+xz, yz+zy, zz-xx-yy]] + N = [ + [xx + yy + zz, 0.0, 0.0, 0.0], + [yz - zy, xx - yy - zz, 0.0, 0.0], + [zx - xz, xy + yx, yy - xx - zz, 0.0], + [xy - yx, zx + xz, yz + zy, zz - xx - yy], + ] # quaternion: eigenvector corresponding to most positive eigenvalue w, V = numpy.linalg.eigh(N) q = V[:, numpy.argmax(w)] @@ -1040,11 +1051,10 @@ def superimposition_matrix(v0, v1, scale=False, usesvd=True): """ v0 = numpy.array(v0, dtype=numpy.float64, copy=False)[:3] v1 = numpy.array(v1, dtype=numpy.float64, copy=False)[:3] - return affine_matrix_from_points(v0, v1, shear=False, - scale=scale, usesvd=usesvd) + return affine_matrix_from_points(v0, v1, shear=False, scale=scale, usesvd=usesvd) -def euler_matrix(ai, aj, ak, axes='sxyz'): +def euler_matrix(ai, aj, ak, axes="sxyz"): """Return homogeneous rotation matrix from Euler angles and axis sequence. ai, aj, ak : Euler's roll, pitch and yaw angles @@ -1070,8 +1080,8 @@ def euler_matrix(ai, aj, ak, axes='sxyz'): firstaxis, parity, repetition, frame = axes i = firstaxis - j = _NEXT_AXIS[i+parity] - k = _NEXT_AXIS[i-parity+1] + j = _NEXT_AXIS[i + parity] + k = _NEXT_AXIS[i - parity + 1] if frame: ai, ak = ak, ai @@ -1080,34 +1090,34 @@ def euler_matrix(ai, aj, ak, axes='sxyz'): si, sj, sk = math.sin(ai), math.sin(aj), math.sin(ak) ci, cj, ck = math.cos(ai), math.cos(aj), math.cos(ak) - cc, cs = ci*ck, ci*sk - sc, ss = si*ck, si*sk + cc, cs = ci * ck, ci * sk + sc, ss = si * ck, si * sk M = numpy.identity(4) if repetition: M[i, i] = cj - M[i, j] = sj*si - M[i, k] = sj*ci - M[j, i] = sj*sk - M[j, j] = -cj*ss+cc - M[j, k] = -cj*cs-sc - M[k, i] = -sj*ck - M[k, j] = cj*sc+cs - M[k, k] = cj*cc-ss + M[i, j] = sj * si + M[i, k] = sj * ci + M[j, i] = sj * sk + M[j, j] = -cj * ss + cc + M[j, k] = -cj * cs - sc + M[k, i] = -sj * ck + M[k, j] = cj * sc + cs + M[k, k] = cj * cc - ss else: - M[i, i] = cj*ck - M[i, j] = sj*sc-cs - M[i, k] = sj*cc+ss - M[j, i] = cj*sk - M[j, j] = sj*ss+cc - M[j, k] = sj*cs-sc + M[i, i] = cj * ck + M[i, j] = sj * sc - cs + M[i, k] = sj * cc + ss + M[j, i] = cj * sk + M[j, j] = sj * ss + cc + M[j, k] = sj * cs - sc M[k, i] = -sj - M[k, j] = cj*si - M[k, k] = cj*ci + M[k, j] = cj * si + M[k, k] = cj * ci return M -def euler_from_matrix(matrix, axes='sxyz'): +def euler_from_matrix(matrix, axes="sxyz"): """Return Euler angles from rotation matrix for specified axis sequence. axes : One of 24 axis sequences as string or encoded tuple @@ -1133,29 +1143,29 @@ def euler_from_matrix(matrix, axes='sxyz'): firstaxis, parity, repetition, frame = axes i = firstaxis - j = _NEXT_AXIS[i+parity] - k = _NEXT_AXIS[i-parity+1] + j = _NEXT_AXIS[i + parity] + k = _NEXT_AXIS[i - parity + 1] M = numpy.array(matrix, dtype=numpy.float64, copy=False)[:3, :3] if repetition: - sy = math.sqrt(M[i, j]*M[i, j] + M[i, k]*M[i, k]) + sy = math.sqrt(M[i, j] * M[i, j] + M[i, k] * M[i, k]) if sy > _EPS: - ax = math.atan2( M[i, j], M[i, k]) - ay = math.atan2( sy, M[i, i]) - az = math.atan2( M[j, i], -M[k, i]) + ax = math.atan2(M[i, j], M[i, k]) + ay = math.atan2(sy, M[i, i]) + az = math.atan2(M[j, i], -M[k, i]) else: - ax = math.atan2(-M[j, k], M[j, j]) - ay = math.atan2( sy, M[i, i]) + ax = math.atan2(-M[j, k], M[j, j]) + ay = math.atan2(sy, M[i, i]) az = 0.0 else: - cy = math.sqrt(M[i, i]*M[i, i] + M[j, i]*M[j, i]) + cy = math.sqrt(M[i, i] * M[i, i] + M[j, i] * M[j, i]) if cy > _EPS: - ax = math.atan2( M[k, j], M[k, k]) - ay = math.atan2(-M[k, i], cy) - az = math.atan2( M[j, i], M[i, i]) + ax = math.atan2(M[k, j], M[k, k]) + ay = math.atan2(-M[k, i], cy) + az = math.atan2(M[j, i], M[i, i]) else: - ax = math.atan2(-M[j, k], M[j, j]) - ay = math.atan2(-M[k, i], cy) + ax = math.atan2(-M[j, k], M[j, j]) + ay = math.atan2(-M[k, i], cy) az = 0.0 if parity: @@ -1165,7 +1175,7 @@ def euler_from_matrix(matrix, axes='sxyz'): return ax, ay, az -def euler_from_quaternion(quaternion, axes='sxyz'): +def euler_from_quaternion(quaternion, axes="sxyz"): """Return Euler angles from quaternion for specified axis sequence. >>> angles = euler_from_quaternion([0.99810947, 0.06146124, 0, 0]) @@ -1176,7 +1186,7 @@ def euler_from_quaternion(quaternion, axes='sxyz'): return euler_from_matrix(quaternion_matrix(quaternion), axes) -def quaternion_from_euler(ai, aj, ak, axes='sxyz'): +def quaternion_from_euler(ai, aj, ak, axes="sxyz"): """Return quaternion from Euler angles and axis sequence. ai, aj, ak : Euler's roll, pitch and yaw angles @@ -1194,8 +1204,8 @@ def quaternion_from_euler(ai, aj, ak, axes='sxyz'): firstaxis, parity, repetition, frame = axes i = firstaxis + 1 - j = _NEXT_AXIS[i+parity-1] + 1 - k = _NEXT_AXIS[i-parity] + 1 + j = _NEXT_AXIS[i + parity - 1] + 1 + k = _NEXT_AXIS[i - parity] + 1 if frame: ai, ak = ak, ai @@ -1211,22 +1221,22 @@ def quaternion_from_euler(ai, aj, ak, axes='sxyz'): sj = math.sin(aj) ck = math.cos(ak) sk = math.sin(ak) - cc = ci*ck - cs = ci*sk - sc = si*ck - ss = si*sk + cc = ci * ck + cs = ci * sk + sc = si * ck + ss = si * sk - q = numpy.empty((4, )) + q = numpy.empty((4,)) if repetition: - q[0] = cj*(cc - ss) - q[i] = cj*(cs + sc) - q[j] = sj*(cc + ss) - q[k] = sj*(cs - sc) + q[0] = cj * (cc - ss) + q[i] = cj * (cs + sc) + q[j] = sj * (cc + ss) + q[k] = sj * (cs - sc) else: - q[0] = cj*cc + sj*ss - q[i] = cj*sc - sj*cs - q[j] = cj*ss + sj*cc - q[k] = cj*cs - sj*sc + q[0] = cj * cc + sj * ss + q[i] = cj * sc - sj * cs + q[j] = cj * ss + sj * cc + q[k] = cj * cs - sj * sc if parity: q[j] *= -1.0 @@ -1244,8 +1254,8 @@ def quaternion_about_axis(angle, axis): q = numpy.array([0.0, axis[0], axis[1], axis[2]]) qlen = vector_norm(q) if qlen > _EPS: - q *= math.sin(angle/2.0) / qlen - q[0] = math.cos(angle/2.0) + q *= math.sin(angle / 2.0) / qlen + q[0] = math.cos(angle / 2.0) return q @@ -1269,11 +1279,14 @@ def quaternion_matrix(quaternion): return numpy.identity(4) q *= math.sqrt(2.0 / n) q = numpy.outer(q, q) - return numpy.array([ - [1.0-q[2, 2]-q[3, 3], q[1, 2]-q[3, 0], q[1, 3]+q[2, 0], 0.0], - [ q[1, 2]+q[3, 0], 1.0-q[1, 1]-q[3, 3], q[2, 3]-q[1, 0], 0.0], - [ q[1, 3]-q[2, 0], q[2, 3]+q[1, 0], 1.0-q[1, 1]-q[2, 2], 0.0], - [ 0.0, 0.0, 0.0, 1.0]]) + return numpy.array( + [ + [1.0 - q[2, 2] - q[3, 3], q[1, 2] - q[3, 0], q[1, 3] + q[2, 0], 0.0], + [q[1, 2] + q[3, 0], 1.0 - q[1, 1] - q[3, 3], q[2, 3] - q[1, 0], 0.0], + [q[1, 3] - q[2, 0], q[2, 3] + q[1, 0], 1.0 - q[1, 1] - q[2, 2], 0.0], + [0.0, 0.0, 0.0, 1.0], + ] + ) def quaternion_from_matrix(matrix, isprecise=False): @@ -1310,7 +1323,7 @@ def quaternion_from_matrix(matrix, isprecise=False): """ M = numpy.array(matrix, dtype=numpy.float64, copy=False)[:4, :4] if isprecise: - q = numpy.empty((4, )) + q = numpy.empty((4,)) t = numpy.trace(M) if t > M[3, 3]: q[0] = t @@ -1340,10 +1353,14 @@ def quaternion_from_matrix(matrix, isprecise=False): m21 = M[2, 1] m22 = M[2, 2] # symmetric matrix K - K = numpy.array([[m00-m11-m22, 0.0, 0.0, 0.0], - [m01+m10, m11-m00-m22, 0.0, 0.0], - [m02+m20, m12+m21, m22-m00-m11, 0.0], - [m21-m12, m02-m20, m10-m01, m00+m11+m22]]) + K = numpy.array( + [ + [m00 - m11 - m22, 0.0, 0.0, 0.0], + [m01 + m10, m11 - m00 - m22, 0.0, 0.0], + [m02 + m20, m12 + m21, m22 - m00 - m11, 0.0], + [m21 - m12, m02 - m20, m10 - m01, m00 + m11 + m22], + ] + ) K /= 3.0 # quaternion is eigenvector of K that corresponds to largest eigenvalue w, V = numpy.linalg.eigh(K) @@ -1363,10 +1380,15 @@ def quaternion_multiply(quaternion1, quaternion0): """ w0, x0, y0, z0 = quaternion0 w1, x1, y1, z1 = quaternion1 - return numpy.array([-x1*x0 - y1*y0 - z1*z0 + w1*w0, - x1*w0 + y1*z0 - z1*y0 + w1*x0, - -x1*z0 + y1*w0 + z1*x0 + w1*y0, - x1*y0 - y1*x0 + z1*w0 + w1*z0], dtype=numpy.float64) + return numpy.array( + [ + -x1 * x0 - y1 * y0 - z1 * z0 + w1 * w0, + x1 * w0 + y1 * z0 - z1 * y0 + w1 * x0, + -x1 * z0 + y1 * w0 + z1 * x0 + w1 * y0, + x1 * y0 - y1 * x0 + z1 * w0 + w1 * z0, + ], + dtype=numpy.float64, + ) def quaternion_conjugate(quaternion): @@ -1482,8 +1504,9 @@ def random_quaternion(rand=None): pi2 = math.pi * 2.0 t1 = pi2 * rand[1] t2 = pi2 * rand[2] - return numpy.array([numpy.cos(t2)*r2, numpy.sin(t1)*r1, - numpy.cos(t1)*r1, numpy.sin(t2)*r2]) + return numpy.array( + [numpy.cos(t2) * r2, numpy.sin(t1) * r1, numpy.cos(t1) * r1, numpy.sin(t2) * r2] + ) def random_rotation_matrix(rand=None): @@ -1524,6 +1547,7 @@ class Arcball: >>> ball.next() """ + def __init__(self, initial=None): """Initialize virtual trackball control. @@ -1542,7 +1566,7 @@ def __init__(self, initial=None): initial = numpy.array(initial, dtype=numpy.float64) if initial.shape == (4, 4): self._qdown = quaternion_from_matrix(initial) - elif initial.shape == (4, ): + elif initial.shape == (4,): initial /= vector_norm(initial) self._qdown = initial else: @@ -1604,7 +1628,7 @@ def drag(self, point): def next(self, acceleration=0.0): """Continue rotation in direction of last drag.""" - q = quaternion_slerp(self._qpre, self._qnow, 2.0+acceleration, False) + q = quaternion_slerp(self._qpre, self._qnow, 2.0 + acceleration, False) self._qpre, self._qnow = self._qnow, q def matrix(self): @@ -1616,11 +1640,11 @@ def arcball_map_to_sphere(point, center, radius): """Return unit sphere coordinates from window coordinates.""" v0 = (point[0] - center[0]) / radius v1 = (center[1] - point[1]) / radius - n = v0*v0 + v1*v1 + n = v0 * v0 + v1 * v1 if n > 1.0: # position outside of sphere n = math.sqrt(n) - return numpy.array([v0/n, v1/n, 0.0]) + return numpy.array([v0 / n, v1 / n, 0.0]) else: return numpy.array([v0, v1, math.sqrt(1.0 - n)]) @@ -1662,14 +1686,31 @@ def arcball_nearest_axis(point, axes): # map axes strings to/from tuples of inner axis, parity, repetition, frame _AXES2TUPLE = { - 'sxyz': (0, 0, 0, 0), 'sxyx': (0, 0, 1, 0), 'sxzy': (0, 1, 0, 0), - 'sxzx': (0, 1, 1, 0), 'syzx': (1, 0, 0, 0), 'syzy': (1, 0, 1, 0), - 'syxz': (1, 1, 0, 0), 'syxy': (1, 1, 1, 0), 'szxy': (2, 0, 0, 0), - 'szxz': (2, 0, 1, 0), 'szyx': (2, 1, 0, 0), 'szyz': (2, 1, 1, 0), - 'rzyx': (0, 0, 0, 1), 'rxyx': (0, 0, 1, 1), 'ryzx': (0, 1, 0, 1), - 'rxzx': (0, 1, 1, 1), 'rxzy': (1, 0, 0, 1), 'ryzy': (1, 0, 1, 1), - 'rzxy': (1, 1, 0, 1), 'ryxy': (1, 1, 1, 1), 'ryxz': (2, 0, 0, 1), - 'rzxz': (2, 0, 1, 1), 'rxyz': (2, 1, 0, 1), 'rzyz': (2, 1, 1, 1)} + "sxyz": (0, 0, 0, 0), + "sxyx": (0, 0, 1, 0), + "sxzy": (0, 1, 0, 0), + "sxzx": (0, 1, 1, 0), + "syzx": (1, 0, 0, 0), + "syzy": (1, 0, 1, 0), + "syxz": (1, 1, 0, 0), + "syxy": (1, 1, 1, 0), + "szxy": (2, 0, 0, 0), + "szxz": (2, 0, 1, 0), + "szyx": (2, 1, 0, 0), + "szyz": (2, 1, 1, 0), + "rzyx": (0, 0, 0, 1), + "rxyx": (0, 0, 1, 1), + "ryzx": (0, 1, 0, 1), + "rxzx": (0, 1, 1, 1), + "rxzy": (1, 0, 0, 1), + "ryzy": (1, 0, 1, 1), + "rzxy": (1, 1, 0, 1), + "ryxy": (1, 1, 1, 1), + "ryxz": (2, 0, 0, 1), + "rzxz": (2, 0, 1, 1), + "rxyz": (2, 1, 0, 1), + "rzyz": (2, 1, 1, 1), +} _TUPLE2AXES = dict((v, k) for k, v in _AXES2TUPLE.items()) @@ -1748,7 +1789,7 @@ def unit_vector(data, axis=None, out=None): if out is not data: out[:] = numpy.array(data, copy=False) data = out - length = numpy.atleast_1d(numpy.sum(data*data, axis)) + length = numpy.atleast_1d(numpy.sum(data * data, axis)) numpy.sqrt(length, length) if axis is not None: length = numpy.expand_dims(length, axis) @@ -1872,7 +1913,7 @@ def is_same_transform(matrix0, matrix1): return numpy.allclose(matrix0, matrix1) -def _import_module(name, package=None, warn=True, prefix='_py_', ignore='_'): +def _import_module(name, package=None, warn=True, prefix="_py_", ignore="_"): """Try import all public attributes from module into global namespace. Existing attributes with name clashes are renamed with prefix. @@ -1883,11 +1924,12 @@ def _import_module(name, package=None, warn=True, prefix='_py_', ignore='_'): """ import warnings from importlib import import_module + try: if not package: module = import_module(name) else: - module = import_module('.' + name, package=package) + module = import_module("." + name, package=package) except ImportError: if warn: warnings.warn("failed to import module %s" % name) @@ -1904,13 +1946,10 @@ def _import_module(name, package=None, warn=True, prefix='_py_', ignore='_'): return True - -#_import_module('_transformations') +# _import_module('_transformations') if __name__ == "__main__": import doctest - import random # used in doctests + numpy.set_printoptions(suppress=True, precision=5) doctest.testmod() - - diff --git a/version.py b/version.py index ed95523..c3e6735 100755 --- a/version.py +++ b/version.py @@ -1,4 +1,4 @@ -#!/usr/bin/env python3 +#!/usr/bin/env python3 # coding: utf-8 # /*########################################################################## # @@ -58,6 +58,10 @@ __all__ = ["date", "version_info", "strictversion", "hexversion", "debianversion", "calc_hexversion", "citation"] + +from collections import namedtuple + + RELEASE_LEVEL_VALUE = {"dev": 0, "alpha": 10, "beta": 11, @@ -68,16 +72,14 @@ "alpha": "a", "beta": "b", "candidate": "rc"} -MAJOR = 2025 -MINOR = 4 +MAJOR = 2026 +MINOR = 2 MICRO = 0 -RELEV = "dev" # <16 +RELEV = "final" # <16 SERIAL = 0 # <16 date = __date__ -from collections import namedtuple - _version_info = namedtuple("version_info", ["major", "minor", "micro", "releaselevel", "serial"]) version_info = _version_info(MAJOR, MINOR, MICRO, RELEV, SERIAL)