From 63d3e22dd717e57a2c426bf2ae6732b49bd50a01 Mon Sep 17 00:00:00 2001 From: Doeke Hekstra Date: Wed, 13 Dec 2023 08:21:20 -0500 Subject: [PATCH 01/10] clip I for Sigma calc --- reciprocalspaceship/algorithms/scale_merged_intensities.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/reciprocalspaceship/algorithms/scale_merged_intensities.py b/reciprocalspaceship/algorithms/scale_merged_intensities.py index d4e57349..fcb55ffd 100644 --- a/reciprocalspaceship/algorithms/scale_merged_intensities.py +++ b/reciprocalspaceship/algorithms/scale_merged_intensities.py @@ -95,6 +95,7 @@ def _french_wilson_posterior_quad(Iobs, SigIobs, Sigma, centric, npoints=100): return mean, np.sqrt(variance), mean_F, np.sqrt(variance_F) + def mean_intensity_by_miller_index(I, H, bandwidth): """ Use a gaussian kernel smoother to compute mean intensities as a function of miller index. @@ -115,6 +116,7 @@ def mean_intensity_by_miller_index(I, H, bandwidth): """ H = np.array(H, dtype=np.float32) I = np.array(I, dtype=np.float32) + I = np.clip(I, a_min=0,a_max=1e20) bandwidth = np.float32(bandwidth) ** 2.0 n = len(I) From c43c76d92a79bd6fdc7e59d27752b86ebfb2c68b Mon Sep 17 00:00:00 2001 From: Doeke Hekstra Date: Sun, 17 Dec 2023 00:12:53 -0500 Subject: [PATCH 02/10] added flag --- .../algorithms/scale_merged_intensities.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/reciprocalspaceship/algorithms/scale_merged_intensities.py b/reciprocalspaceship/algorithms/scale_merged_intensities.py index fcb55ffd..8608508a 100644 --- a/reciprocalspaceship/algorithms/scale_merged_intensities.py +++ b/reciprocalspaceship/algorithms/scale_merged_intensities.py @@ -96,7 +96,7 @@ def _french_wilson_posterior_quad(Iobs, SigIobs, Sigma, centric, npoints=100): -def mean_intensity_by_miller_index(I, H, bandwidth): +def mean_intensity_by_miller_index(I, H, bandwidth, clip_neg_Iobs): """ Use a gaussian kernel smoother to compute mean intensities as a function of miller index. @@ -108,6 +108,8 @@ def mean_intensity_by_miller_index(I, H, bandwidth): Nx3 array of miller indices bandwidth : float(optional) Kernel bandwidth in miller units + clip_neg_Iobs : bool(optional) + If true, clips Iobs to be at least 0 before calculating Sigma Returns ------- @@ -116,7 +118,8 @@ def mean_intensity_by_miller_index(I, H, bandwidth): """ H = np.array(H, dtype=np.float32) I = np.array(I, dtype=np.float32) - I = np.clip(I, a_min=0,a_max=1e20) + if clip_neg_Iobs: + I = np.clip(I, a_min=0,a_max=1e20) bandwidth = np.float32(bandwidth) ** 2.0 n = len(I) @@ -187,6 +190,7 @@ def scale_merged_intensities( mean_intensity_method="isotropic", bins=100, bw=2.0, + clip_neg_Iobs=False, ): """ Scales merged intensities using Bayesian statistics in order to @@ -242,6 +246,11 @@ def scale_merged_intensities( parameter controls the distance that each reflection impacts in reciprocal space. Only affects output if mean_intensity_method is \"anisotropic\". + clip_neg_Iobs : bool + Will set negative Iobs to 0 for the purpose of calculating Sigma. + Addresses rare cases where local average Iobs is negative. + Default: False. + Returns ------- @@ -288,7 +297,7 @@ def scale_merged_intensities( ) elif mean_intensity_method == "anisotropic": Sigma = ( - mean_intensity_by_miller_index(I / multiplicity, ds.get_hkls(), bw) + mean_intensity_by_miller_index(I / multiplicity, ds.get_hkls(), bw, clip_neg_Iobs) * multiplicity ) From 95642de4b95a420d223edd2bb82bc317161fdbd7 Mon Sep 17 00:00:00 2001 From: Doeke Hekstra Date: Sun, 17 Dec 2023 12:27:09 -0500 Subject: [PATCH 03/10] revert to c43c76d92a79bd6fdc7e59d27752b86ebfb2c68b --- .../algorithms/scale_merged_intensities.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/reciprocalspaceship/algorithms/scale_merged_intensities.py b/reciprocalspaceship/algorithms/scale_merged_intensities.py index 8608508a..10c5549c 100644 --- a/reciprocalspaceship/algorithms/scale_merged_intensities.py +++ b/reciprocalspaceship/algorithms/scale_merged_intensities.py @@ -95,8 +95,12 @@ def _french_wilson_posterior_quad(Iobs, SigIobs, Sigma, centric, npoints=100): return mean, np.sqrt(variance), mean_F, np.sqrt(variance_F) +<<<<<<< HEAD def mean_intensity_by_miller_index(I, H, bandwidth, clip_neg_Iobs): +======= +def mean_intensity_by_miller_index(I, H, bandwidth): +>>>>>>> d25d7da19615acf5c6694af20766c9a7cbeb80b2 """ Use a gaussian kernel smoother to compute mean intensities as a function of miller index. @@ -118,8 +122,12 @@ def mean_intensity_by_miller_index(I, H, bandwidth, clip_neg_Iobs): """ H = np.array(H, dtype=np.float32) I = np.array(I, dtype=np.float32) +<<<<<<< HEAD if clip_neg_Iobs: I = np.clip(I, a_min=0,a_max=1e20) +======= + I = np.clip(I, a_min=0, a_max=1e20) +>>>>>>> d25d7da19615acf5c6694af20766c9a7cbeb80b2 bandwidth = np.float32(bandwidth) ** 2.0 n = len(I) From 1ae48c302991bea138385a139a1b8499ecb14c22 Mon Sep 17 00:00:00 2001 From: Doeke Hekstra Date: Sun, 17 Dec 2023 12:39:44 -0500 Subject: [PATCH 04/10] manually resolved merge conflict --- .../algorithms/scale_merged_intensities.py | 9 --------- 1 file changed, 9 deletions(-) diff --git a/reciprocalspaceship/algorithms/scale_merged_intensities.py b/reciprocalspaceship/algorithms/scale_merged_intensities.py index 10c5549c..32726c69 100644 --- a/reciprocalspaceship/algorithms/scale_merged_intensities.py +++ b/reciprocalspaceship/algorithms/scale_merged_intensities.py @@ -95,12 +95,7 @@ def _french_wilson_posterior_quad(Iobs, SigIobs, Sigma, centric, npoints=100): return mean, np.sqrt(variance), mean_F, np.sqrt(variance_F) -<<<<<<< HEAD - def mean_intensity_by_miller_index(I, H, bandwidth, clip_neg_Iobs): -======= -def mean_intensity_by_miller_index(I, H, bandwidth): ->>>>>>> d25d7da19615acf5c6694af20766c9a7cbeb80b2 """ Use a gaussian kernel smoother to compute mean intensities as a function of miller index. @@ -122,12 +117,8 @@ def mean_intensity_by_miller_index(I, H, bandwidth): """ H = np.array(H, dtype=np.float32) I = np.array(I, dtype=np.float32) -<<<<<<< HEAD if clip_neg_Iobs: I = np.clip(I, a_min=0,a_max=1e20) -======= - I = np.clip(I, a_min=0, a_max=1e20) ->>>>>>> d25d7da19615acf5c6694af20766c9a7cbeb80b2 bandwidth = np.float32(bandwidth) ** 2.0 n = len(I) From 1041521c7b398018d1982d5f7def8133a36a6ccd Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 17 Dec 2023 17:41:49 +0000 Subject: [PATCH 05/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../algorithms/scale_merged_intensities.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/reciprocalspaceship/algorithms/scale_merged_intensities.py b/reciprocalspaceship/algorithms/scale_merged_intensities.py index 32726c69..7b8267a1 100644 --- a/reciprocalspaceship/algorithms/scale_merged_intensities.py +++ b/reciprocalspaceship/algorithms/scale_merged_intensities.py @@ -118,7 +118,7 @@ def mean_intensity_by_miller_index(I, H, bandwidth, clip_neg_Iobs): H = np.array(H, dtype=np.float32) I = np.array(I, dtype=np.float32) if clip_neg_Iobs: - I = np.clip(I, a_min=0,a_max=1e20) + I = np.clip(I, a_min=0, a_max=1e20) bandwidth = np.float32(bandwidth) ** 2.0 n = len(I) @@ -246,7 +246,7 @@ def scale_merged_intensities( reciprocal space. Only affects output if mean_intensity_method is \"anisotropic\". clip_neg_Iobs : bool - Will set negative Iobs to 0 for the purpose of calculating Sigma. + Will set negative Iobs to 0 for the purpose of calculating Sigma. Addresses rare cases where local average Iobs is negative. Default: False. @@ -296,7 +296,9 @@ def scale_merged_intensities( ) elif mean_intensity_method == "anisotropic": Sigma = ( - mean_intensity_by_miller_index(I / multiplicity, ds.get_hkls(), bw, clip_neg_Iobs) + mean_intensity_by_miller_index( + I / multiplicity, ds.get_hkls(), bw, clip_neg_Iobs + ) * multiplicity ) From 3b1f155a0a9a52c464ceea954e6d2d402526cc55 Mon Sep 17 00:00:00 2001 From: Doeke Hekstra Date: Tue, 19 Dec 2023 17:21:52 -0500 Subject: [PATCH 06/10] switch to clipping Sigma --- .../algorithms/scale_merged_intensities.py | 27 +++++++++++-------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/reciprocalspaceship/algorithms/scale_merged_intensities.py b/reciprocalspaceship/algorithms/scale_merged_intensities.py index 7b8267a1..f76cb0f3 100644 --- a/reciprocalspaceship/algorithms/scale_merged_intensities.py +++ b/reciprocalspaceship/algorithms/scale_merged_intensities.py @@ -40,6 +40,7 @@ def _french_wilson_posterior_quad(Iobs, SigIobs, Sigma, centric, npoints=100): Iobs = np.array(Iobs, dtype=np.float64) SigIobs = np.array(SigIobs, dtype=np.float64) Sigma = np.array(Sigma, dtype=np.float64) + print(np.amin(Sigma)) # Integration window based on the normal, likelihood distribution window_size = 20.0 # In standard devs @@ -95,7 +96,7 @@ def _french_wilson_posterior_quad(Iobs, SigIobs, Sigma, centric, npoints=100): return mean, np.sqrt(variance), mean_F, np.sqrt(variance_F) -def mean_intensity_by_miller_index(I, H, bandwidth, clip_neg_Iobs): +def mean_intensity_by_miller_index(I, H, bandwidth, clip_neg_Sigma=False, eps=0.001): """ Use a gaussian kernel smoother to compute mean intensities as a function of miller index. @@ -107,8 +108,9 @@ def mean_intensity_by_miller_index(I, H, bandwidth, clip_neg_Iobs): Nx3 array of miller indices bandwidth : float(optional) Kernel bandwidth in miller units - clip_neg_Iobs : bool(optional) - If true, clips Iobs to be at least 0 before calculating Sigma + clip_neg_Sigma : bool(optional, default False) + If true, clips Sigma to be at least 0.0001 before calculating Sigma + eps (float) : small positive number to clip Sigma to (if clip_neg_Sigma=True). Default=0.001. Returns ------- @@ -117,8 +119,6 @@ def mean_intensity_by_miller_index(I, H, bandwidth, clip_neg_Iobs): """ H = np.array(H, dtype=np.float32) I = np.array(I, dtype=np.float32) - if clip_neg_Iobs: - I = np.clip(I, a_min=0, a_max=1e20) bandwidth = np.float32(bandwidth) ** 2.0 n = len(I) @@ -126,7 +126,10 @@ def mean_intensity_by_miller_index(I, H, bandwidth, clip_neg_Iobs): for i in range(n): K = np.exp(-0.5 * ((H - H[i]) * (H - H[i])).sum(1) / bandwidth) S[i] = (I * K).sum() / K.sum() - + + if clip_neg_Sigma: + S = np.clip(S, a_min=eps, a_max=1e20) + return S @@ -189,7 +192,8 @@ def scale_merged_intensities( mean_intensity_method="isotropic", bins=100, bw=2.0, - clip_neg_Iobs=False, + clip_neg_Sigma=False, + eps=0.001, ): """ Scales merged intensities using Bayesian statistics in order to @@ -245,10 +249,11 @@ def scale_merged_intensities( parameter controls the distance that each reflection impacts in reciprocal space. Only affects output if mean_intensity_method is \"anisotropic\". - clip_neg_Iobs : bool - Will set negative Iobs to 0 for the purpose of calculating Sigma. + clip_neg_Sigma : bool + Will set negative Sigma to 0 for the purpose of calculating Sigma. Addresses rare cases where local average Iobs is negative. - Default: False. + Default: False. *Used by valdo*. + eps : float Floor imposed on Sigma if clipping. Returns @@ -297,7 +302,7 @@ def scale_merged_intensities( elif mean_intensity_method == "anisotropic": Sigma = ( mean_intensity_by_miller_index( - I / multiplicity, ds.get_hkls(), bw, clip_neg_Iobs + I / multiplicity, ds.get_hkls(), bw, clip_neg_Sigma, eps ) * multiplicity ) From 74e2fbbc85b5eebe5e3f930f95dd57cf12d08035 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 19 Dec 2023 22:22:37 +0000 Subject: [PATCH 07/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- reciprocalspaceship/algorithms/scale_merged_intensities.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/reciprocalspaceship/algorithms/scale_merged_intensities.py b/reciprocalspaceship/algorithms/scale_merged_intensities.py index f76cb0f3..37761b99 100644 --- a/reciprocalspaceship/algorithms/scale_merged_intensities.py +++ b/reciprocalspaceship/algorithms/scale_merged_intensities.py @@ -126,10 +126,10 @@ def mean_intensity_by_miller_index(I, H, bandwidth, clip_neg_Sigma=False, eps=0. for i in range(n): K = np.exp(-0.5 * ((H - H[i]) * (H - H[i])).sum(1) / bandwidth) S[i] = (I * K).sum() / K.sum() - + if clip_neg_Sigma: S = np.clip(S, a_min=eps, a_max=1e20) - + return S From 42469888a61065d89941bbb4fb23292aabce77f1 Mon Sep 17 00:00:00 2001 From: Doeke Hekstra Date: Tue, 19 Dec 2023 17:26:12 -0500 Subject: [PATCH 08/10] removed print statement --- reciprocalspaceship/algorithms/scale_merged_intensities.py | 1 - 1 file changed, 1 deletion(-) diff --git a/reciprocalspaceship/algorithms/scale_merged_intensities.py b/reciprocalspaceship/algorithms/scale_merged_intensities.py index f76cb0f3..3616f704 100644 --- a/reciprocalspaceship/algorithms/scale_merged_intensities.py +++ b/reciprocalspaceship/algorithms/scale_merged_intensities.py @@ -40,7 +40,6 @@ def _french_wilson_posterior_quad(Iobs, SigIobs, Sigma, centric, npoints=100): Iobs = np.array(Iobs, dtype=np.float64) SigIobs = np.array(SigIobs, dtype=np.float64) Sigma = np.array(Sigma, dtype=np.float64) - print(np.amin(Sigma)) # Integration window based on the normal, likelihood distribution window_size = 20.0 # In standard devs From 9339461a0ffa000416a9a8bc1f7d88e169d4832a Mon Sep 17 00:00:00 2001 From: Doeke Hekstra Date: Thu, 21 Dec 2023 15:37:49 -0500 Subject: [PATCH 09/10] switch to minimum_sigma --- .../algorithms/scale_merged_intensities.py | 25 ++++++------------- 1 file changed, 8 insertions(+), 17 deletions(-) diff --git a/reciprocalspaceship/algorithms/scale_merged_intensities.py b/reciprocalspaceship/algorithms/scale_merged_intensities.py index cc3cad1b..17c6f6ca 100644 --- a/reciprocalspaceship/algorithms/scale_merged_intensities.py +++ b/reciprocalspaceship/algorithms/scale_merged_intensities.py @@ -95,7 +95,7 @@ def _french_wilson_posterior_quad(Iobs, SigIobs, Sigma, centric, npoints=100): return mean, np.sqrt(variance), mean_F, np.sqrt(variance_F) -def mean_intensity_by_miller_index(I, H, bandwidth, clip_neg_Sigma=False, eps=0.001): +def mean_intensity_by_miller_index(I, H, bandwidth): """ Use a gaussian kernel smoother to compute mean intensities as a function of miller index. @@ -107,9 +107,6 @@ def mean_intensity_by_miller_index(I, H, bandwidth, clip_neg_Sigma=False, eps=0. Nx3 array of miller indices bandwidth : float(optional) Kernel bandwidth in miller units - clip_neg_Sigma : bool(optional, default False) - If true, clips Sigma to be at least 0.0001 before calculating Sigma - eps (float) : small positive number to clip Sigma to (if clip_neg_Sigma=True). Default=0.001. Returns ------- @@ -126,9 +123,6 @@ def mean_intensity_by_miller_index(I, H, bandwidth, clip_neg_Sigma=False, eps=0. K = np.exp(-0.5 * ((H - H[i]) * (H - H[i])).sum(1) / bandwidth) S[i] = (I * K).sum() / K.sum() - if clip_neg_Sigma: - S = np.clip(S, a_min=eps, a_max=1e20) - return S @@ -191,8 +185,7 @@ def scale_merged_intensities( mean_intensity_method="isotropic", bins=100, bw=2.0, - clip_neg_Sigma=False, - eps=0.001, + minimum_sigma=-np.inf, ): """ Scales merged intensities using Bayesian statistics in order to @@ -248,11 +241,8 @@ def scale_merged_intensities( parameter controls the distance that each reflection impacts in reciprocal space. Only affects output if mean_intensity_method is \"anisotropic\". - clip_neg_Sigma : bool - Will set negative Sigma to 0 for the purpose of calculating Sigma. - Addresses rare cases where local average Iobs is negative. - Default: False. *Used by valdo*. - eps : float Floor imposed on Sigma if clipping. + minimum_sigma : float + Minimum value imposed on Sigma (default: -np.inf, that is: no minimum). Returns @@ -296,15 +286,16 @@ def scale_merged_intensities( if mean_intensity_method == "isotropic": dHKL = ds["dHKL"].to_numpy(dtype=np.float64) Sigma = ( - mean_intensity_by_resolution(I / multiplicity, dHKL, bins) * multiplicity + mean_intensity_by_resolution(I / multiplicity, dHKL, bins) ) elif mean_intensity_method == "anisotropic": Sigma = ( mean_intensity_by_miller_index( - I / multiplicity, ds.get_hkls(), bw, clip_neg_Sigma, eps + I / multiplicity, ds.get_hkls(), bw ) - * multiplicity ) + Sigma = np.clip(Sigma, a_min=minimum_sigma, a_max=np.inf) + Sigma = Sigma * multiplicity # Initialize outputs ds[outputI] = 0.0 From ca49796202330e61dacc7be0fc638eb5d16c8c27 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 21 Dec 2023 20:38:29 +0000 Subject: [PATCH 10/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../algorithms/scale_merged_intensities.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/reciprocalspaceship/algorithms/scale_merged_intensities.py b/reciprocalspaceship/algorithms/scale_merged_intensities.py index 17c6f6ca..1f9f697c 100644 --- a/reciprocalspaceship/algorithms/scale_merged_intensities.py +++ b/reciprocalspaceship/algorithms/scale_merged_intensities.py @@ -241,7 +241,7 @@ def scale_merged_intensities( parameter controls the distance that each reflection impacts in reciprocal space. Only affects output if mean_intensity_method is \"anisotropic\". - minimum_sigma : float + minimum_sigma : float Minimum value imposed on Sigma (default: -np.inf, that is: no minimum). @@ -285,15 +285,9 @@ def scale_merged_intensities( I, Sig = ds[intensity_key].to_numpy(), ds[sigma_key].to_numpy() if mean_intensity_method == "isotropic": dHKL = ds["dHKL"].to_numpy(dtype=np.float64) - Sigma = ( - mean_intensity_by_resolution(I / multiplicity, dHKL, bins) - ) + Sigma = mean_intensity_by_resolution(I / multiplicity, dHKL, bins) elif mean_intensity_method == "anisotropic": - Sigma = ( - mean_intensity_by_miller_index( - I / multiplicity, ds.get_hkls(), bw - ) - ) + Sigma = mean_intensity_by_miller_index(I / multiplicity, ds.get_hkls(), bw) Sigma = np.clip(Sigma, a_min=minimum_sigma, a_max=np.inf) Sigma = Sigma * multiplicity