Source code for TEMPy.protein.scoring_functions

# =============================================================================
#     This file is part of TEMPy.
#
#     TEMPy is a software designed to help the user in the manipulation
#     and analyses of macromolecular assemblies using 3D electron microscopy
#     maps.
#
#     Copyright  2015 Birkbeck College University of London.
#
#     Authors: Maya Topf, Daven Vasishtan, Arun Prasad Pandurangan,
#     Irene Farabella, Agnel-Praveen Joseph, Harpal Sahota
#
#     This software is made available under GPL V3 license
#     http://www.gnu.org/licenses/gpl-3.0.html
#
#
#     Please cite your use of TEMPy in published work:
#
#     Farabella, I., Vasishtan, D., Joseph, A.P., Pandurangan, A.P., Sahota, H.
#     & Topf, M. (2015). J. Appl. Cryst. 48.
#
# =============================================================================

from collections import OrderedDict

import itertools
import math
import sys
import gc
import warnings

import numpy as np
from scipy.fftpack import (
    ifftn,
    fftshift,
    ifftshift,
)
from scipy.spatial import KDTree
from scipy import stats
import pyfftw

from TEMPy.protein.structure_blurrer import StructureBlurrer
from TEMPy.protein.rigid_body_parser import RBParser
from TEMPy.math import vector
from TEMPy.maps.em_map import Map
from TEMPy.math import fourier


def _MI_C_helper(
        arr1,
        arr2,
        layers1=20,
        layers2=20,
        N=0,
        lc1=0.0,
        lc2=0.0):
    ly1 = int(layers1)  # noqa:F841
    ly2 = int(layers2)  # noqa:F841
    # input 3D arrays
    nz = int(arr1.shape[0])  # noqa:F841
    ny = int(arr1.shape[1])  # noqa:F841
    nx = int(arr1.shape[2])  # noqa:F841
    # min and max to set left and right bound
    ma1 = np.ma.masked_less_equal(arr1, lc1, copy=False)
    min1 = float(ma1.min())
    max1 = float(ma1.max())

    ma2 = np.ma.masked_less_equal(arr2, lc2, copy=False)
    min2 = float(ma2.min())
    max2 = float(ma2.max())

    min1 = float(min1 - ((max1 - min1) / layers1) * 0.0001)
    min2 = float(min2 - ((max2 - min2) / layers2) * 0.0001)
    # bin width
    step1 = (max1-min1) / float(layers1)  # noqa:F841
    step2 = (max2-min2) / float(layers2)  # noqa:F841
    # histogram freq in bins
    freq1 = np.zeros(layers1, dtype=float)  # noqa:F841
    freq2 = np.zeros(layers2, dtype=float)  # noqa:F841
    comb_freq = np.zeros((layers1, layers2), dtype=float)  # noqa:F841

    s1 = 0
    s2 = 0
    p1 = 0.0
    p2 = 0.0
    pcomb = 0.0
    Hxy = 0.0
    Hy = 0.0
    Hx = 0.0
    # /*long index = 0;
    # long indexend = nz * ny * nx;
    # while (index < indexend){
    #     va1  = arr1[index];
    #     va2 = arr2[index];*/
    # /* use 3d array loop */
    for z in range(nz):
        for y in range(ny):
            for x in range(nx):
                va1 = arr1[z, y, x] ** 2
                va2 = arr2[z, y, x] ** 2
                i = 0
                for i in range(ly1):
                    if (
                            (va1 > (min1 + i * step1)) and
                            (va1 <= (min1 + (i + 1) * step1))
                    ):
                        freq1[i] += 1.0
                        s1 += 1
                        break
                if i == ly1:
                    i = i-1
                for j in range(ly2):
                    if (
                            (va2 > (min2+j * step2)) and
                            (va2 <= (min2+(j+1) * step2))
                    ):
                        freq2[j] += 1.0
                        s2 += 1
                        comb_freq[i, j] += 1.0
                        break
    for i in range(ly1):
        p1 = freq1[i] / s1
        # /*std::cout << s1 << ' ' << s2 << std::endl;*/
        for j in range(ly2):
            p2 = freq2[j] / s2
            pcomb = comb_freq[i, j] / s1
            if (pcomb != 0.0):
                Hxy += (-pcomb*np.log2(pcomb))
            if ((i == 0) and (p2 != 0.0)):
                Hy += (-p2*np.log2(p2))
        if (p1 != 0.0):
            Hx += (-p1*np.log2(p1))
    P, e1, e2 = np.histogram2d(ma1.ravel(), ma2.ravel(), bins=20)
    P /= P.sum()
    p1 = P.sum(axis=0)
    p2 = P.sum(axis=1)
    p1 = p1[p1 > 0]
    p2 = p2[p1 > 0]
    P = P[P > 0]
    Hx_ = (-p1 * np.log2(p1)).sum()
    Hy_ = (-p2 * np.log2(p2)).sum()
    Hxy_ = (-P * np.log2(P)).sum()
    t2 = Hx_ + Hy_ - Hxy_  # noqa: F841
    import pdb
    pdb.set_trace()
    # /*std::cout << Hxy << ' ' << Hx << ' ' << Hy << ' ' << std::endl;*/
    if (N == 1):
        if (Hxy != 0.0):
            return (Hx+Hy)/Hxy
        else:
            return 0.0
    else:
        return Hx+Hy-Hxy


def _MI_C_helper_new(
        arr1,
        arr2,
        layers1=20,
        layers2=20,
        N=0,
        lc1=0.0,
        lc2=0.0):
    ma1 = np.ma.masked_less_equal(arr1, lc1, copy=False)
    ma2 = np.ma.masked_less_equal(arr2, lc2, copy=False)

    P, xe, ye = np.histogram2d(
        ma1.ravel(),
        ma2.ravel(),
        bins=(layers1, layers2)
    )
    P /= P.sum()
    p1 = P.sum(axis=0)
    p2 = P.sum(axis=1)
    p1 = p1[p1 > 0]
    p2 = p2[p2 > 0]
    P = P[P > 0]
    Hx_ = (-p1 * np.log2(p1)).sum()
    Hy_ = (-p2 * np.log2(p2)).sum()
    Hxy_ = (-P * np.log2(P)).sum()
    return Hx_ + Hy_ - Hxy_


class ScoringFunctions:
    """Class implementing most TEMPy scoring functions.

    Args:
        None.

    Returns:
        A ScoringFunction object. See
    """
    def __init__(self):
        warnings.warn("ScoringFunctions class is deprecated, call scoring "
            "functions directly using TEMPy.protein.scores module.",
            category=DeprecationWarning, stacklevel=1)
        pass

    def _overlap_map_samebox(self, map1, map2):
        """ Is this function used somewhere - candidate for deletion?

        volume overlap within 2 maps with same box size

        Return:
           % of overlap
        """
        binmap1 = map1.fullMap > 0.0
        binmap2 = map2.fullMap > 0.0
        mask_array = (binmap1 * binmap2) > 0.0
        return[
            np.count_nonzero(binmap1),
            np.count_nonzero(binmap2),
            np.count_nonzero(mask_array),
            mask_array.size
        ]

    def calculate_overlap_scores(
                                self,
                                map_target,
                                map_probe,
                                map_target_threshold,
                                map_probe_threshold
                                ):
        """
            mask maps with 2 cut-off map_target_threshold and
            map_probe_threshold (vol thr.)

            return:
            fraction of overlap with respect to target, with respect to
            probe and with respect to the total (union)
        """
        if not self.mapComparison(map_target, map_probe):
            raise ValueError(
                f"Could not calculate scores for input maps with pixel sizes {map_target.apix}, {map_probe.apix},"
                f" shapes {map_target.box_size()}, {map_probe.box_size()} and "
                f"origins {map_target.origin}, {map_probe.origin}."
            )
        binmap1 = map_target.fullMap > float(map_target_threshold)
        binmap2 = map_probe.fullMap > float(map_probe_threshold)
        mask_array = (binmap1 * binmap2) > 0

        size1 = np.sum(binmap1)
        size2 = np.sum(binmap2)
        return (
            float(np.sum(mask_array)) / size1,
            float(np.sum(mask_array)) / size2,
            float(np.sum(mask_array)) / (size1 + size2),
        )

    def _overlap_map_array(
            self,
            map_target,
            map_target_threshold,
            map_probe,
            map_probe_threshold,
    ):
        """Helper function for mutual information calculation.

        Masks maps using unique cutoff for each (map_target_threshold and
        map_probe_threshold) and returns a final mask which is True where
        both maps are above their respective thresholds.

        Args:
            map_target, map_probe: Input Map instance
            map_target_threshold, map_probe_threshold: Threshold used to mask
            respective map instance

        Returns:
            mask_array: Mask which is true where both maps are above their
            respective thresholds.
        """
        binmap1 = map_target.fullMap > float(map_target_threshold)
        binmap2 = map_probe.fullMap > float(map_probe_threshold)
        mask_array = (binmap1 * binmap2) > 0
        return mask_array

    def calculate_map_threshold(self, map_target):
        """Calculates map threshold as mean value + (2 * map standard deviation)

        Args:
            map_target: Map instance

        Returns:
            float: Threshold value
        """
        # NOTE: MOVE TO TEMPy.maps.em_map.Map??
        try:
            peak, ave, sigma = map_target._peak_density()
            vol_threshold = float(ave) + (2.0 * float(sigma))
        except Exception:
            if len(map_target.header) == 0:
                amean = map_target.mean()
                rms = map_target.std()
                vol_threshold = float(amean) + (1.5 * float(rms))
            else:
                amean = map_target.mean()
                rms = map_target.std()
                vol_threshold = float(amean) + (1.5 * float(rms))
        return vol_threshold

    def mapComparison(self, map_target, map_probe):
        """
        Checks if properties (sampling rate, box size and origin) of two maps
        are equal.

        Args:
            map_target, map_probe
                Map instance to compare.

        Returns:
            True if the map properties are the same between two maps, False
            otherwise.
        """
        if (
                np.isclose(map_target.apix,
                           map_probe.apix,
                           atol=1E-6).all() and
                map_target.box_size() == map_probe.box_size()
        ):
            if (
                    round(map_target.origin[0], 2) == round(map_probe.origin[0], 2) and  # noqa:E501
                    round(map_target.origin[1], 2) == round(map_probe.origin[1], 2) and  # noqa:E501
                    round(map_target.origin[2], 2) == round(map_probe.origin[2], 2)  # noqa:E501
            ):
                return True
            else:
                return False
        else:
            return False

    def _failed_match(self):
        # should implement TEMPy errors for more elegant error handling
        print(
            'Warning: can\'t match the map at the moment,'
            'use map with same box size.'
        )
        sys.exit()

    def scale_median(self, arr1, arr2):
        """Scale one list/array of scores with respect to another based on
        distribution around median.

        Arguments:
            arr1, arr2: Array/list of scores

        Returns:
            Scaled arr1, based on the values in arr2
        """
        nparr1 = np.array(arr1)
        nparr2 = np.array(arr2)
        med_dev_1 = np.median(np.absolute(nparr1 - np.median(nparr1)))
        med_dev_2 = np.median(np.absolute(nparr2 - np.median(nparr2)))
        if med_dev_1 == 0.0:
            scale_factor = 0.0
        else:
            scale_factor = med_dev_2 / med_dev_1
        shift_factor = np.median(nparr2) - (scale_factor * np.median(nparr1))

        # TODO: find a better way to avoid outliers in general
        if (max(nparr1) - min(nparr1)) > 0.1:
            scaled_arr = ((scale_factor * nparr1 + shift_factor) + nparr2) / 2.
        else:
            scaled_arr = nparr2
        return scaled_arr

    def _CCC_calc(self, m1, m2):
        arr1 = m1.view(float)
        arr2 = m2.view(float)  # noqa:F841
        nd = len(arr1.shape)
        if nd == 2 and len(arr1.shape)[1] == 0:
            nd = 1
        l = 1  # noqa:E741
        dim = np.zeros(3, dtype=int)
        for i in range(nd):
            l *= arr1.shape[i]  # noqa:E741
            dim[i] = arr1.shape[i]
        corr = 0.0
        numer = 0.0
        var1 = 0.0
        var2 = 0.0
        if nd == 1:
            for z in range(dim[0]):
                numer += arr1[z]*arr2[z]
                var1 += arr1[z]**2
                var2 += arr2[z]**2
        elif nd == 3:
            for z in range(dim[0]):
                for y in range(dim[1]):
                    for x in range(dim[2]):
                        numer += arr1[z, y, x]*arr2[z, y, x]
                        var1 += arr1[z, y, x]**2
                        var2 += arr2[z, y, x]**2
        corr = numer/math.sqrt(var1*var2)
        corr = min(1.0, corr)
        corr = max(-1.0, corr)
        return corr

    def CCC_map(
            self,
            map_target,
            map_probe,
            map_target_threshold=0.0,
            map_probe_threshold=0.0,
            mode=1,
            meanDist=False,
            cmode=True,
    ):
        """Calculate cross-correlation between two Map instances using one of
        3 methods.

        Args:
            map_target, map_probe: EMMap instance to compare.

            map_target_threshold, map_probe_threshold:
                Map threshold, used to select pixels when calculating CCC with
                method 2. If not specified, will be calculated using
                :meth:`calculate_map_threshold <TEMPy.protein.scoring_functions.ScoringFunctions.calculate_map_threshold>`  # noqa:E501

            mode:
                1. CCC calculation on all pixels in map
                2. CCC only on pixels with density value above respective thresholds
                3. CCC only on pixels within the mask

            meanDist: True if the deviation from mean needs to be calculated

            cmode: If True, use pure python implementation to calculate CCC.
                If False, use numpy implementation.

        Returns:

            CCC, percentage_overlap: The cross correlation score and the
            percentage overlap between the two input maps

        Raises:

            Warning: Will return -1, 0 if the two input maps do not match
                shape, pixel size and origin.
        """
        if self.mapComparison(map_target, map_probe):
            if not mode == 1:
                if map_target_threshold == 0 and map_probe_threshold == 0:
                    map_target_threshold = self.calculate_map_threshold(
                        map_target
                    )
                    map_probe_threshold = self.calculate_map_threshold(
                        map_probe
                    )
                bin_map1 = map_target.fullMap > float(map_target_threshold)
                bin_map2 = map_probe.fullMap > float(map_probe_threshold)
                minim = np.sum(bin_map1)
                minim2 = np.sum(bin_map2)
                if minim2 < minim:
                    minim = minim2
                mask_array = (bin_map1 * bin_map2) > 0
                if not minim == 0.0:
                    perc_ovr = float(np.sum(mask_array)) / minim
                else:
                    print(
                        'No map overlap (Cross correlation score), exiting '
                        'score calculation..'
                    )
                    return -1.0, 0.0
                if perc_ovr < 0.02:
                    return -1.0, 0.0
            else:
                perc_ovr = 1.0

            # calculate CCC within volume of overlap
            if mode == 3:
                if np.sum(mask_array) == 0:
                    print(
                        'No map overlap (Cross correlation score), '
                        'exiting score calculation..'
                    )
                    return -1.0, 0.0
                map1_mask = map_target.fullMap[mask_array]
                map2_mask = map_probe.fullMap[mask_array]
                if meanDist:
                    map1_mask = map1_mask - np.mean(map1_mask)
                    map2_mask = map2_mask - np.mean(map2_mask)
                if cmode:
                    corr = self._CCC_calc(
                        map1_mask.flatten(),
                        map2_mask.flatten(),
                    )
                else:
                    corr = None
                if corr is None:
                    return (
                        (
                            np.sum(map1_mask * map2_mask) /
                            np.sqrt(
                                np.sum(np.square(map1_mask)) *
                                np.sum(np.square(map2_mask))
                            )
                        ),
                        perc_ovr
                    )
                else:
                    return corr, perc_ovr

            # calculate CCC for contoured maps based on threshold
            elif mode == 2:
                map1_mask = map_target.fullMap * bin_map1
                map2_mask = map_probe.fullMap * bin_map2
                if meanDist:
                    map1_mask = (
                            map1_mask -
                            np.mean(map_target.fullMap[bin_map1])
                    )
                    map2_mask = (
                            map2_mask -
                            np.mean(map_probe.fullMap[bin_map2])
                    )
                    map1_mask = map1_mask * bin_map1
                    map2_mask = map2_mask * bin_map2
                else:
                    map1_mask = map_target.fullMap * bin_map1
                    map2_mask = map_probe.fullMap * bin_map2
                if cmode:
                    corr = self._CCC_calc(map1_mask, map2_mask)
                else:
                    corr = None
                if corr is None:
                    return (
                        (
                            np.sum(map1_mask * map2_mask) /
                            np.sqrt(
                                np.sum(np.square(map1_mask)) *
                                np.sum(np.square(map2_mask))
                            )
                        ),
                        perc_ovr
                    )
                else:
                    return corr, perc_ovr
            # calculate on the complete map
            if meanDist:
                if cmode:
                    corr = self._CCC_calc(
                        map_target.fullMap - np.mean(map_target.fullMap),
                        map_probe.fullMap - np.mean(map_probe.fullMap)
                    )
                else:
                    corr = None
                if corr is None:
                    return (
                        (
                            np.sum(
                                (
                                    map_target.fullMap -
                                    np.mean(map_target.fullMap)
                                ) *
                                (
                                    map_probe.fullMap -
                                    np.mean(map_probe.fullMap)
                                )
                            ) /
                            (
                                np.sqrt(
                                    np.sum(
                                        np.square(
                                            map_target.fullMap -
                                            np.mean(map_target.fullMap)
                                        )
                                    ) *
                                    np.sum(
                                        np.square(
                                            map_probe.fullMap -
                                            np.mean(map_probe.fullMap)
                                        )
                                    )
                                )
                            )
                        ),
                        perc_ovr
                    )
                else:
                    return corr, perc_ovr
            if cmode:
                corr = self._CCC_calc(map_target.fullMap, map_probe.fullMap)
            else:
                corr = None
            if corr is None:
                return (
                    np.sum(map_target.fullMap * map_probe.fullMap) /
                    np.sqrt(
                        np.sum(np.square(map_target.fullMap)) *
                        np.sum(np.square(map_probe.fullMap))
                    ),
                    perc_ovr
                )
            else:
                return corr, perc_ovr
        else:
            print('@@@ Maps could not be matched')
            return -1., 0.

    def CCC(self, map_target, map_probe):
        """Calculate cross-correlation between two Map instances.

        Args:
            map_target, map_probe
                Map instance to compare.
        Returns:
            CCC score
        Raises:
            sys.exit: Input maps do not match in shape, pixel size or origin
        """
        if self.mapComparison(map_target, map_probe):
            return (
                map_target.normalise().getMap()*map_probe.normalise().getMap()
            ).mean()
        else:
            self._failed_match()

    def LSF(self, map_target, map_probe):
        """Calculate least-squares between two Map instances.

        Arguments:
            map_target, map_probe: Map instance to compare.
        Return:
            Least-squares value
        Raises:
            sys.exit: Input maps do not match in shape, pixel size or origin
        """
        if self.mapComparison(map_target, map_probe):
            map_target, map_probe = map_target, map_probe
        else:
            self._failed_match()
        return ((map_target.getMap() - map_probe.getMap())**2).mean()

    def laplace_CCC(
            self,
            map_target,
            map_probe,
            prefil=(False, False)
    ):
        """Calculate CCC between two laplacian filtered Map instances.

        The laplacian filter is an edge detection filter, which computes the
        second derviatives of an image/volume. Using Laplacian filters in
        combination with CCC scoring of cryo EM models  was shown to be more
        accurate than CCC scoring alone for low resolution maps, in
        `(Chacon and Wriggers, 2002) <https://pubmed.ncbi.nlm.nih.gov/11922671/>`_

        Arguments:
            map_target, map_probe: Map instances to compare.

            prefil: 2-tuple of boolean values, one for each map respectively.
                True if Map instance is already Laplacian-filtered. False
                otherwise.
        Return:
            Laplacian cross-correlation score
        Raises:
            sys.exit: Input maps do not match in shape, pixel size or origin
        """

        if self.mapComparison(map_target, map_probe):
            # I don't understand this - can we replace it with pass?
            _, _ = map_target, map_probe  # noqa:F841
        else:
            self._failed_match()

        if not prefil[0]:
            map_target = map_target.laplace_filtered()
        if not prefil[1]:
            map_probe = map_probe.laplace_filtered()
        map_target = map_target.normalise()
        map_probe = map_probe.normalise()
        return self.CCC(map_target, map_probe)

    def normal_vector_score(
            self,
            map_target,
            map_probe,
            primary_boundary,
            secondary_boundary=0.0,
            Filter=None,
    ):
        """Calculate the Normal Vector Score between two Map surfaces.

        Based on 3SOM algorithm (Ceulemans and Russell, 2004).

        Arguments:
            map_target, map_probe:
                Map instances to compare.
            primary_boundary, secondary_boundary:
                If a filter is selected, just input a contour level as primary
                threshold. Otherwise, need to run get_primary_boundary and
                get_second_boundary based on map target.
            Filter:
                Filter to be applied to inputs maps before calculating score.
                    1. Sobel Filter (Filter=='Sobel')
                    2. Laplace Filter (Filter=='Laplace')
                    3. Minimum Filter (Filter=='Minimum')
                    4. Mean Filter (Filter=='Mean')
                    5. No filtering (Filter==None)
        Return:
            Normal vector score.
        Raises:
            sys.exit: Input maps do not match in shape, pixel size or origin
        """
        if Filter not in ['Sobel', 'Laplace', 'Mean', 'Minimum', None]:
            print('Incorrect name of filter: %s" % Filter')
            print(
                'Select one of the following Filters if applicable: %s\n' %
                ', '.join(['Sobel', 'Laplace'])
            )
            sys.exit()

        scores = []
        if not self.mapComparison(map_target, map_probe):
            self._failed_match()
        assert isinstance(primary_boundary, float)
        assert isinstance(secondary_boundary, float)
        if primary_boundary > secondary_boundary:
            temp_thr = secondary_boundary
            secondary_boundary = primary_boundary
            primary_boundary = temp_thr

        points = np.argwhere(
            (map_target.fullMap > primary_boundary) &
            (map_target.fullMap < secondary_boundary)
        )

        if Filter == 'Sobel':
            map1_surface = map_target._sobel_filter_contour(primary_boundary)
            points = np.argwhere(
                map1_surface.fullMap > (map1_surface.max() / 2.0)
            )
        elif Filter == 'Laplace':
            map1_surface = map_target._laplace_filtered_contour(
                primary_boundary
            )
            points = np.argwhere(
                map1_surface.fullMap > (map1_surface.max() / 2.0)
            )
        elif Filter == 'Minimum':
            map1_surface = map_target._surface_minimum_filter(
                float(primary_boundary)
            )
            points = np.argwhere(map1_surface == 1)
        elif Filter == 'Mean':
            map1_filter = map_target._surface_features(float(primary_boundary))
            bin_test = [0.0001]
            for ii in range(1, 41):
                bin_test.append(0.025 * ii)
            freq_test = np.histogram(map1_filter.fullMap, bin_test)[0]
            sum_freq = 0.0
            for fr in range(len(freq_test)):
                sum_freq += float(freq_test[fr])
                if (
                        sum_freq / np.sum(freq_test) > 0.05 and
                        bin_test[fr+1] >= 0.3
                ):
                    t1 = bin_test[fr+1]
                    break
                if (
                        sum_freq / np.numsum(freq_test) > 0.10 or
                        sum_freq > 100000
                ):
                    t1 = bin_test[fr+1]
                    break
            points = np.argwhere(
                (map1_filter.fullMap > 0.0) & (map1_filter.fullMap < t1)
            )
        # C++ calculation
        flagc = 1
        try:
            vecnorm_target = map_target._get_normal_vector(points)
            vecnorm_probe = map_probe._get_normal_vector(points)
        except Exception:
            flagc = 0
        if vecnorm_target is None or vecnorm_probe is None:
            flagc = 0
        ct = 0

        if flagc == 1:
            for i in range(len(vecnorm_target)):
                ct += 1
                nvec = vecnorm_target[i]
                ovec = vecnorm_probe[i]
                # add max value for regions of null variation
                if (
                        nvec[0] == 0. and
                        nvec[1] == 0. and
                        nvec[2] == 0.
                ):
                    if (
                            ovec[0] == 0. and
                            ovec[1] == 0. and
                            ovec[2] == 0.0
                    ):
                        continue
                    else:
                        scores.append(3.14)
                        continue
                else:
                    if (
                            ovec[0] == 0. and
                            ovec[1] == 0. and
                            ovec[2] == 0.
                    ):
                        scores.append(3.14)
                        continue
                try:
                    dotprod = (
                            ovec[0] *
                            nvec[0] +
                            ovec[1] *
                            nvec[1] +
                            ovec[2] *
                            nvec[2]
                    )
                    den = (
                            np.sqrt(nvec[0]**2 + nvec[1]**2 + nvec[2]**2) *
                            np.sqrt(ovec[0]**2 + ovec[1]**2 + ovec[2]**2)
                    )
                    if abs(dotprod - den) < 0.00001:
                        ang = 0.0
                    else:
                        ang = math.acos(min(max(dotprod / den, -1.0), 1.0))
                    if den == 0.0:
                        print(dotprod, den, nvec, ovec)
                    scores.append(abs(ang))
                except ValueError:
                    print(
                        'Error: Angle could not be calculated: ',
                        nvec,
                        ' ',
                        ovec
                    )
            if len(scores) == 0:
                print(
                    'There are no points to be scored! The threshold values or'
                    ' the number of points to be considered needs to be'
                    ' changed.'
                )
                return None
            else:
                if sum(scores) == 0:
                    return 0.0
                else:
                    return 1 - (sum(scores) / (len(points) * 3.14))
        scores = []
        ct1 = 0
        if flagc == 0:
            for v in points:
                n_vec = map_target.get_normal_vector(v[2], v[1], v[0])
                o_vec = map_probe.get_normal_vector(v[2], v[1], v[0])
                ct1 += 1
                if (
                    n_vec.x == -9 and
                    n_vec.y == -9 and
                    n_vec.z == -9
                ):
                    if (
                        o_vec.x == -9 and
                        o_vec.y == -9 and
                        o_vec.z == -9
                    ):
                        continue
                    else:
                        scores.append(3.14)
                        continue
                else:
                    if (
                        o_vec.x == -9 and
                        o_vec.y == -9 and
                        o_vec.z == -9
                    ):
                        scores.append(3.14)
                        continue
                try:
                    scores.append(abs(n_vec.arg(o_vec)))
                except ValueError:
                    print(
                        'Error: Angle between ' +
                        str(n_vec) +
                        ', ' +
                        str(o_vec) +
                        ' for point %d, %d, %d cannot be calculated.' %
                        (v.x, v.y, v.z)
                    )
        if len(scores) == 0:
            print(
                'There are no points to be scored! The threshold values or '
                'the number of points to be considered needs to be changed.'
            )
        else:
            if sum(scores) == 0:
                return 0
            else:
                return 1 - (sum(scores) / (len(points) * 3.14))

    def get_partial_DLSF(
            self,
            num_of_points,
            map_target,
            map_probe
    ):
        """Calculate the DLSF score between two Map instances.

        The DLSF is similar to the LSF;
        whereas the LSF compares absolute density values,
        the DLSF compares the difference between pairs of values.

        Arguments:
            map_target, map_probe: Map instances to compare.
            num_of_points: Number of significant points.
        Return:
            DLSF score, or string "can't match the map" if Input maps do not
            match in shape, pixel size or origin
        """

        if not self.mapComparison(map_target, map_probe):
            return "can't Match the map"

        map_target_sig_pairs = map_target._get_random_significant_pairs(
            int(num_of_points)
        )
        otherMap = map_probe
        score = 0.0
        for p in map_target_sig_pairs:
            z1 = p[0]
            y1 = p[1]
            x1 = p[2]
            z2 = p[3]
            y2 = p[4]
            x2 = p[5]
            dens = p[6]
            prot_dens = (
                otherMap.fullMap[z1][y1][x1] - otherMap.fullMap[z2][y2][x2]
            )
            score += (dens - prot_dens)**2
        return score / map_target.fullMap.size

    def _MI(self, map_target, map_probe, layers=20):
        """Calculate the mutual information score between two Maps.

        Old version of the score?

        Arguments:
            map_target, map_probe: Map instances to compare.

            layers: Number of layers used to bin the map. Default is 20 as in
                Shatsky et al., 2008.
        Return: MI score
        Raises:
            sys.exit: Input maps do not match in shape, pixel size or origin
        """

        if self.mapComparison(map_target, map_probe):
            m1, m2 = map_target, map_probe
        else:
            self._failed_match()
        score = 0
        m1_levels = (m1.max() - m1.min()) / layers
        m2_levels = (m2.max() - m2.min()) / layers
        for x in range(layers):
            for y in range(layers):
                m1_level_map = (
                        (m1.getMap() >= m1.min() + (x * m1_levels)) *
                        (m1.getMap() <= m1.min() + ((x + 1) * m1_levels))
                )
                m2_level_map = (
                        (m2.getMap() >= m2.min() + (y * m2_levels)) *
                        (m2.getMap() <= m2.min() + ((y + 1) * m2_levels))
                )
                comb_level_map = m1_level_map*m2_level_map
                p_m1 = float(m1_level_map.sum())/m1_level_map.size
                p_m2 = float(m2_level_map.sum())/m2_level_map.size
                p_comb = float(comb_level_map.sum())/comb_level_map.size
                if p_comb == 0:
                    mi_score = 0.0
                else:
                    mi_score = p_comb * math.log(p_comb / (p_m1 * p_m2), 2)
                score += mi_score
        return score

    def _MI_C(
            self,
            m1,
            m2,
            layers1=20,
            layers2=20,
            N=0,
            lc1=0.0,
            lc2=0.0,
    ):
        arr1 = (m1).view(float)
        arr2 = (m2).view(float)
        return _MI_C_helper_new(arr1, arr2, layers1, layers2, N, lc1, lc2)

    def MI(
            self,
            map_target,
            map_probe,
            map_target_threshold=0.0,
            map_probe_threshold=0.0,
            mode=1,
            layers1=None,
            layers2=None,
            weight=False,
            cmode=True,
    ):
        """Calculate the mutual information score between two Map instances.

        Arguments:
            map_target, map_probe: Map instance to compare.
            map_target_threshold, map_probe_threshold: Threshold used for
                contouring
            mode:
                1. use complete map for calculation
                3. use overlap region for calculation
            layers1: Number of layers used to bin map_target. If None, value
                is calculated based on Sturges, 1926 method.
            layers2: Number of layers used to bin map_probe. If None, value
                is calculated based on Sturges, 1926 method.
            weight: If True: normalised MI (Studholme et al.) is used to
                account for overlap of 'contours'
            cmode:
        Returns:
            MI score
        Raises:
            sys.exit: Input maps do not match in shape, pixel size or origin
        """
        if not self.mapComparison(map_target, map_probe):
            self._failed_match()
        if map_target_threshold == 0.0:
            map_target_threshold = self.calculate_map_threshold(map_target)
        if map_probe_threshold == 0.0:
            map_probe_threshold = self.calculate_map_threshold(map_probe)
        # calculation on the complete map
        if mode == 1:
            if weight:
                wt = 1
            else:
                wt = 0
            if layers1 is None:
                layers1 = 20
            if layers2 is None:
                layers2 = 20
            min1 = (
                    np.amin(map_target.fullMap) -
                    0.00001 *
                    (np.amax(map_target.fullMap) - np.amin(map_target.fullMap))
            )
            min2 = (
                    np.amin(map_probe.fullMap) -
                    0.00001 *
                    (np.amax(map_probe.fullMap) - np.amin(map_probe.fullMap))
            )
            if cmode:
                mic = self._MI_C(
                    map_target.fullMap,
                    map_probe.fullMap,
                    layers1,
                    layers2,
                    wt,
                    min1,
                    min2
                )
            else:
                mic = None
            if mic is not None:
                return mic
            # digitize whole map based on layers
            map1_bin = map_target._map_digitize(
                map_target.min(),
                layers1,
                True,
            )
            map2_bin = map_probe._map_digitize(map_probe.min(), layers2, True)
            bins1 = []
            for i in range(layers1+2):
                bins1.append(i)
            bins2 = []
            for i in range(layers2+2):
                bins2.append(i)
            # calculate frequency of bins
            map1_freq = np.histogram(map1_bin.fullMap, bins1)[0][1:]
            map2_freq = np.histogram(map2_bin.fullMap, bins2)[0][1:]
        elif mode == 3:
            mask_array = self._overlap_map_array(
                map_target,
                map_target_threshold,
                map_probe,
                map_probe_threshold
            )
            if np.sum(mask_array) == 0:
                print(
                    'No map overlap (Mutual information score), exiting score '
                    'calculation..'
                )
                return 0.0
            # sturges rule provides a way of calculating number of bins :
            # 1 + math.log(number of points)
            if layers1 is None:
                try:
                    layers1 = int(1 + math.log(np.sum(mask_array), 2))
                except ValueError:
                    print(
                        'No map overlap (Mutual information score), '
                        'exiting score calculation..'
                    )
                    return 0.0
            if layers2 is None:
                try:
                    layers2 = int(1 + math.log(np.sum(mask_array), 2))
                except ValueError:
                    print(
                        'No map overlap (Mutual information score), exiting '
                        'score calculation..'
                    )
                    return 0.0
            layers1 = max(layers1, 15)
            layers2 = max(layers2, 15)
            if weight:
                wt = 1
            else:
                wt = 0
            if cmode:
                mic = self._MI_C(
                    np.array(map_target.fullMap * mask_array),
                    np.array(map_probe.fullMap * mask_array),
                    layers1,
                    layers2,
                    wt,
                )
            else:
                mic = None
            if mic is not None:
                return mic
            # digitize masked map based on layers
            map1_bin = map_target.copy()
            map2_bin = map_probe.copy()
            map1_bin.fullMap = map1_bin.fullMap * mask_array
            map2_bin.fullMap = map2_bin.fullMap * mask_array
            map1_bin = map1_bin._map_digitize(
                map_target.fullMap[mask_array].min(),
                layers1,
                True
            )
            map2_bin = map2_bin._map_digitize(
                map_probe.fullMap[mask_array].min(),
                layers2,
                True
            )
            # make sure the outside region is filled with zeros
            map1_bin.fullMap = map1_bin.fullMap * mask_array
            map2_bin.fullMap = map2_bin.fullMap * mask_array
            # background frequencies from the whole map
            bins1 = []
            for i in range(layers1 + 2):
                bins1.append(i)
            bins2 = []
            for i in range(layers2 + 2):
                bins2.append(i)
            # calculate frequency of bins
            map1_freq = np.histogram(map1_bin.fullMap, bins1)[0][1:]
            map2_freq = np.histogram(map2_bin.fullMap, bins2)[0][1:]

        score = 0.0
        total = 0

        if np.sum(map1_freq) == 0:
            print(
                'No map overlap (Mutual information score), exiting score '
                'calculation..'
            )
            return 0.0
        if np.sum(map2_freq) == 0:
            print(
                'No map overlap (Mutual information score), exiting score '
                'calculation..'
            )
            return 0.0
        list_overlaps = []
        for x in range(layers1):
            mask_array = map1_bin.fullMap == float(x+1)
            overlap_freq = np.histogram(
                map2_bin.fullMap[mask_array], bins2
            )[0][1:]
            total += float(np.sum(overlap_freq))
            list_overlaps.append(overlap_freq)

        if total == 0:
            print(
                'No map overlap (Mutual information score), exiting score '
                'calculation..'
            )
            return 0.0
        enter = 0
        Hxy = 0.0
        Hx = 0.0
        Hy = 0.0
        mi_score = 0.0
        p_comb = 0.0
        for x in range(layers1):
            p_m1 = map1_freq[x]/float(np.sum(map1_freq))
            for y in range(layers2):
                enter = 1
                p_comb = list_overlaps[x][y]/total
                p_m2 = map2_freq[y] / float(np.sum(map2_freq))
                if p_comb == 0:
                    mi_score = 0.0

                else:
                    Hxy += -p_comb * math.log(p_comb, 2)
                score += mi_score
                if x == 0 and not p_m2 == 0.0:
                    Hy += (-p_m2 * math.log(p_m2, 2))
            if not p_m1 == 0.0:
                Hx += (-p_m1 * math.log(p_m1, 2))
        if enter == 1:
            if weight:
                if Hxy == 0.0:
                    return 0.0
                return (Hx + Hy) / Hxy
            return Hx + Hy - Hxy
        else:
            return None

    def MI_new(self, exp_map, sim_map, bins1=20, bins2=20, normalise=False):
        """
        Taken from https://stackoverflow.com/questions/20491028/optimal-way-to-compute-pairwise-mutual-information-using-numpy
        """

        c_xy, x, y = np.histogram2d(exp_map.fullMap.ravel(),
                            sim_map.fullMap.ravel(),
                            bins=[bins1, bins2])

        g, _, _, _ = stats.chi2_contingency(c_xy, lambda_='log-likelihood')

        mi_nl = 0.5 * g / c_xy.sum()

        return np.log2(np.exp(mi_nl))

    def _hausdorff_list(
            self,
            primary_boundary,
            secondary_boundary,
            kdtree,
            map_probe,
    ):
        """
        This is for the chamdef distance def chamfer_distance, min max density
        value that define the surface of the protein

        Arguments:
            *kdtree*
                (there are 2 of them in numpy one Cbased on py-based, the
                latter is better, ctrl) this have to be one of the input.
                kdtree from map_target
            *primary_boundary, secondary_boundary*
                need to run get_primary_boundary and get_second_boundary for
                map_probe

            NOTE: if you keep the kdtree as parametre out os less time
            consuming as building it takes time.
        """
        points = map_probe.get_pos(primary_boundary, secondary_boundary)
        return kdtree.query(points)[0]

    def chamfer_distance(
            self,
            map_target,
            map_probe,
            primary_boundary,
            secondary_boundary,
            kdtree=None,
    ):
        """Calculate the chamfer distance Score between two Map instances.

        Currently unstable and not recommended for use.

        Arguments:
            map_target, map_probe: Map instances to compare.
            primary_boundary: Threshold that defined the EM map surface, based
                on the volume of the protein model. Can be calculated using
                :meth:`Map.makeKDTree <TEMPy.maps.em_map.Map.get_primary_boundary>`
            secondary_boundary: The second bound density point.
                Can be calculated using
                :meth:`Map.makeKDTree <TEMPy.maps.em_map.Map.get_second_boundary>`
            kdtree: Precalculated kdtree. If None, KDTree is constructed using
                :meth:`Map.makeKDTree <TEMPy.maps.em_map.Map.makeKDTree>`
        """
        if self.mapComparison(map_target, map_probe):
            m1, m2 = map_target, map_probe
        else:
            self._failed_match()
        if kdtree:
            return self._hausdorff_list(
                primary_boundary, secondary_boundary, kdtree, m2
            ).mean()
        else:
            print(m1, primary_boundary, secondary_boundary)
            kdtree = m1.makeKDTree(primary_boundary, secondary_boundary)
            if kdtree is None:
                print('Error. No points selected, change boundary parameters.')
                sys.exit()
            else:
                return self._hausdorff_list(
                    primary_boundary,
                    secondary_boundary,
                    kdtree,
                    m2
                ).mean()

    def _surface_distance_score(
            self,
            map_target,
            map_probe,
            map_target_threshold1=0.0,
            map_probe_threshold=0.0,
            Filter=None,
            map_target_threshold2=0.0,
            weight=False
    ):
        """
        Calculate the chamfer distance Score between two Map instances.

        Arguments:
            *map_target, map_probe*
                EMMap instances to compare.
            *map_target_threshold1*
                contour threshold of the target map.
                This value is used the primary boundary if
                map_target_threshold2 is given.
            *map_probe_threshold*
                contour threshold for the probe map.
            *Filter*
                definition of the surface:
                    1) None : surface defined by known boundaries -
                              map_target_threshold1 & map_target_threshold2
                    If the boundaries are not known and target&probe map
                    contour levels are known:
                        2) Std : to define the boundaries, contour level +- 5%
                                 sigma is calculated. 5%sigma is used to limit
                                 the number of points picked as surface. For
                                 small maps, higher values (eg: 10%sigma)
                                 can be used.
                        3) Mean: a mean filter is applied on the binary contour
                                 mask over a long window. The resulting mask
                                 has values between 0 and 1. Points with values
                                 less than 0.3 is used to represent surface.
                                 As the average is calculated on a long window,
                                 highly exposed surface points have very low
                                 values and partially exposed surfaces/grooves
                                 have relatively higher values.
                                 This definition is useful especially when the
                                 map surface has many features/projections.
                        4) Minimum: a minimum filter is applied on a binary
                                    contour mask to locate surface points.
                                    Voxels surrounded by points outside the
                                    contour (zeroes) are detected as surface.
                                    Voxels surrounded by points outside the
                                    contour (zeroes) are detected as surface.
                        5) Sobel: sobel filter is applied on the map to detect
                                  high density gradients. Before applying the
                                  sobel filter, it is important to reduce the
                                  noise density and large variations
                                  (gradients) in the noise region.
            *weight*
                If set true, the distances between the surface points is
                normalized in a way similar to GDT (Zemla 2007)
                calculation for atomic co-ordinate alignments.
        """
        # check if both maps are on the same grid
        if not self.mapComparison(map_target, map_probe):
            print('@@@ Maps could not be matched')
            return -999.
        # if the boundaries are known, calculate the kdtree
        if Filter is None:
            kdtree = map_target.makeKDTree(
                map_target_threshold1,
                map_target_threshold2
            )
            probe_points = map_probe.get_pos(
                map_target_threshold1,
                map_target_threshold2
            )
        elif Filter == 'Std':
            target_points = np.argwhere(
                (map_target.fullMap > (
                        float(map_target_threshold1)-(map_target.std()*0.10))
                 ) &
                (map_target.fullMap < (
                        float(map_target_threshold1)+(map_target.std()*0.10))
                 )
            )
            probe_points = np.argwhere(
                (map_probe.fullMap > (
                        float(map_probe_threshold)-(map_probe.std()*0.10))
                 ) &
                (map_probe.fullMap < (
                        float(map_probe_threshold)+(map_probe.std()*0.10))
                 )
            )
            if len(target_points) < len(probe_points):
                probe_points1 = np.copy(target_points)
                target_points = np.copy(probe_points)
                probe_points = np.copy(probe_points1)
            if len(target_points) == 0 or len(probe_points) == 0:
                print('Surface detection failed (Std filter), exiting..')
                return None
            try:
                from scipy.spatial import cKDTree
                try:
                    kdtree = cKDTree(target_points)
                except RuntimeError:
                    return None
            except ImportError:
                try:
                    kdtree = KDTree(target_points)
                except RuntimeError:
                    return None
        elif Filter == 'Mean':
            map1_filter = map_target._surface_features(
                float(map_target_threshold1)
            )
            map2_filter = map_probe._surface_features(
                float(map_probe_threshold)
            )
            """
            define surface based on the filtered mask values.
            points with values less than 0.3 are usually preferred. But in
            some cases like viruses, most surface points are highly exposed
            and a large number of points are returned and the calculation
            becomes slow.
            Hence an additional filter is added: the maximum allowed points
            is 10% of box size.
            The minimum number of points is kept as 7%. This mode is less
            sensitive to the number of surface points chosen
            as the extent of exposure is used for defining surface. Hence
            thick surface is not usually required.
            calculate frequencies in bins for filtered mask.
            The smaller the bins, more precise will be the calculation of
            points allowed based on percent of points chosen.
            As this is just an additional filter and doesn't affect the
            calculations drastically, 40 bins are used to calculate
            frequencies.
            """
            bin_test = [0.0001]
            for ii in range(1, 41):
                bin_test.append(0.025 * ii)
            freq_test = np.histogram(map1_filter.fullMap, bin_test)[0]
            map1_filled = np.sum(map1_filter.fullMap > 0)

            sum_freq = 0.0
            for fr in range(len(freq_test)):
                sum_freq += float(freq_test[fr])
                if sum_freq / map1_filled > 0.05 and bin_test[fr+1] >= 0.3:
                    t1 = bin_test[fr+1]
                    break
                if sum_freq / map1_filled > 0.10 or sum_freq > 200000:
                    t1 = bin_test[fr+1]
                    break
            sum_freq = 0.0
            freq_test = np.histogram(map2_filter.fullMap, bin_test)[0]
            map2_filled = np.sum(map2_filter.fullMap > 0)
            for fr in range(len(freq_test)):
                sum_freq += float(freq_test[fr])
                if sum_freq/map2_filled > 0.05 and bin_test[fr+1] >= 0.3:
                    t2 = bin_test[fr+1]
                    break
                if sum_freq/map2_filled > 0.10 or sum_freq > 200000:
                    t2 = bin_test[fr+1]
                    break
            # t1 and t2 are the selected levels based on filtered values and
            # percent of points
            target_points = np.argwhere(
                (map1_filter.fullMap > 0.0) & (map1_filter.fullMap <= t1)
            )
            probe_points = np.argwhere(
                (map2_filter.fullMap > 0.0) & (map2_filter.fullMap <= t2)
            )
            if len(target_points) == 0 or len(probe_points) == 0:
                print('Surface detection failed (Mean filter), exiting..')
                return None
            if len(target_points) < len(probe_points):
                probe_points1 = np.copy(target_points)
                target_points = np.copy(probe_points)
                probe_points = np.copy(probe_points1)
            try:
                from scipy.spatial import cKDTree
                try:
                    kdtree = cKDTree(target_points)
                except RuntimeError:
                    return None
            except ImportError:
                try:
                    kdtree = KDTree(target_points)
                except RuntimeError:
                    return None
        elif Filter == 'Minimum':
            map1_surface = map_target._surface_minimum_filter(
                float(map_target_threshold1)
            )
            map2_surface = map_probe._surface_minimum_filter(
                float(map_probe_threshold)
            )
            # select the surface points represented by the mask
            target_points = np.argwhere(map1_surface == 1)
            probe_points = np.argwhere(map2_surface == 1)
            if len(target_points) == 0 or len(probe_points) == 0:
                print('Surface detection failed (Minimum filter), exiting..')
                return None
            if len(target_points) + len(probe_points) > 250000:
                return None
            if len(target_points) < len(probe_points):
                probe_points1 = np.copy(target_points)
                target_points = np.copy(probe_points)
                probe_points = np.copy(probe_points1)
            try:
                from scipy.spatial import cKDTree
                try:
                    kdtree = cKDTree(target_points)
                except RuntimeError:
                    return None
            except ImportError:
                try:
                    kdtree = KDTree(target_points)
                except RuntimeError:
                    return None
        elif Filter == 'Sobel':
            map1_surface = map_target._sobel_filter_contour(
                float(map_target_threshold1)
            )
            map2_surface = map_probe._sobel_filter_contour(
                float(map_probe_threshold)
            )

            target_points = np.argwhere(
                map1_surface.fullMap > map1_surface.max() / float(2)
            )
            probe_points = np.argwhere(
                map2_surface.fullMap > map2_surface.max() / float(2)
            )
            if len(target_points) == 0 or len(probe_points) == 0:
                print('Surface detection failed (Sobel filter), exiting..')
                return None
            if len(target_points) < len(probe_points):
                probe_points1 = np.copy(target_points)
                target_points = np.copy(probe_points)
                probe_points = np.copy(probe_points1)
            try:
                from scipy.spatial import cKDTree
                try:
                    kdtree = cKDTree(target_points)
                except RuntimeError:
                    return None
            except ImportError:
                try:
                    kdtree = KDTree(target_points)
                except RuntimeError:
                    return None
        distances = kdtree.query(probe_points)[0]
        if len(distances) == 0:
            return None
        if not weight:
            if not np.mean(distances) <= 0.05:
                return 1 / np.mean(distances)
            else:
                return 1/0.05

        x = int(30.0 / map_target.apix)
        if np.amin(distances) < x / 2:
            distances = distances - np.amin(distances)
        bins = []
        i = 0
        while i <= float(x):
            bins.append(i * 1.0)
            i += 1
        num_distances = len(distances)
        overlap_freq = np.histogram(distances, bins)[0]
        for fr_i in range(len(overlap_freq)):
            if overlap_freq[fr_i] > np.amax(overlap_freq) / 3.:
                break
        total_ext = fr_i
        bins = bins[fr_i:]
        if cl:  # noqa:F821
            points_cl = probe_points
            if len(points_cl) == 0:
                return None, None
            try:
                kdtree = cKDTree(points_cl)
            except Exception:
                return None, None
            neighbors_num = 20
            distance_lim = 3.0
            neigh = kdtree.query(
                points_cl,
                k=neighbors_num,
                distance_upper_bound=distance_lim,
            )[1]

            cl_weight = (
                    np.numsum(np.sum(neigh < len(neigh), axis=1) > 17) /
                    float(len(probe_points))
            )
            distances_align = distances
            distances_sel = distances_align[
                np.sum(neigh < len(neigh), axis=1) > 17
            ]
            distances = distances_sel[:]

        overlap_freq = np.histogram(distances, bins)[0]
        total = total_ext
        cumul_freq = 0.0
        enter = 0
        sum_sc = 0.0
        for i in range(len(overlap_freq)):
            w = len(overlap_freq) - (i)
            try:
                cumul_freq += overlap_freq[i]
            except IndexError:
                pass
            try:
                perc_equiv = float(cumul_freq) / num_distances
            except ZeroDivisionError:
                print('Distance weighting failed!!. Check surface defined')
                return None, None
            sum_sc = sum_sc + ((w) * perc_equiv)
            total += (w)
            enter = 1
        score = float(sum_sc) / total
        if cl:  # noqa:F821
            if enter == 1:
                if len(distances_sel) == 0.0:
                    return 0.0
                if np.mean(distances_sel) == 0.0:
                    return 0.0
                if cl_weight == 0.0:
                    return 0.0
                return score
            else:
                return None, None
        if enter == 1:
            if np.mean(distances) <= 0.05:
                return 1.0
            if np.mean(distances) == 0.0:
                return 1.0
            return score
        else:
            return None, None

    def envelope_score(
            self,
            map_target,
            primary_boundary,
            structure_instance,
            norm=True,
            ):
        """
        Calculate the envelope score between a target Map and a Structure
        Instances.

        Arguments:
            map_target: Target Map Instance.
            primary_boundary: Value specified is calculated with
                primary_boundary of the map object.
            structure_instance: Model structure instance.
            norm: Normalise the score between 0 and 1 if True.

        Returns:
            Envelope score, normalised between 0 and 1 if norm is True

        """
        binMap = map_target.make_bin_map(primary_boundary)
        max_score = float(-2 * np.sum(binMap.fullMap))
        min_score = float(np.sum(binMap.fullMap)-2 * np.sum(binMap.fullMap+1))

        blurrer = StructureBlurrer()
        struct_binMap = blurrer.make_atom_overlay_map1(
            map_target,
            structure_instance
        )
        grid = struct_binMap.get_pos(0.9, 1.1)
        for x, y, z in grid:
            g = binMap[z][y][x]
            if g == -1:
                binMap[z][y][x] = 2
            elif g == 0:
                binMap[z][y][x] = -2
        score = float(np.sum(binMap.fullMap))
        if norm:
            norm_score = float((score-min_score) / (max_score-min_score))
            return norm_score
        else:
            return score

    def envelope_score_map(
            self,
            map_target,
            map_probe,
            map_target_threshold=0,
            map_probe_threshold=0,
            norm=True,
    ):
        """Calculate the envelope score between two maps

        Arguments:
            map_target, map_probe: Map instance to compare.
            map_target_threshold, map_probe_threshold:
                Map threshold that defines the envelope, or surface, of the
                map. Will be automatically calculated using
                :meth:`self.calculate_map_threshold <TEMPy.protein.scoring_functions.ScoringFunctions.calculate_map_threshold>`
                if not specified.
            norm: Normalise the score between 0 and 1 if True.

        Returns:
            Envelope score, normalised between 0 and 1 if norm is True
        """
        if self.mapComparison(map_target, map_probe):
            if map_target_threshold == 0:
                map_target_threshold = self.calculate_map_threshold(map_target)
            if map_probe_threshold == 0:
                map_probe_threshold = self.calculate_map_threshold(map_probe)

        binMap = map_target.make_bin_map(map_target_threshold)
        max_score = float(-2 * np.sum(binMap.fullMap))
        min_score = float(
            np.sum(binMap.fullMap)-2 * np.sum(binMap.fullMap + 1)
        )
        struct_binMap = map_probe.make_bin_map(map_probe_threshold)
        newMap = binMap.fullMap + 2 * struct_binMap.fullMap
        hist_array = np.histogram(newMap, 4)
        score = (
                2 *
                hist_array[0][0] -
                (2 * (hist_array[0][1])) -
                (hist_array[0][2])
        )
        if norm:
            norm_score = float((score - min_score)) / (max_score - min_score)
            return norm_score
        else:
            return score

    def _percent_overlap(
            self,
            map_target,
            map_probe,
            map_target_threshold,
            map_probe_threshold,
            flagsize=0,
    ):
        """Calculate the fraction of overlap between two map grids.

        Arguments:
            map_target, map_probe: Map instance to compare.
            map_target_threshold, map_probe_threshold:
                map contour thresholds for map_target and map_probe.
        Return:
            Percent overlap with respect to smaller grid
        """
        if self.mapComparison(map_target, map_probe):
            binmap1 = map_target.fullMap > float(map_target_threshold)
            binmap2 = map_probe.fullMap > float(map_probe_threshold)
            minim = len(map_target.fullMap[binmap1])
            if len(map_probe.fullMap[binmap2]) < minim:
                minim = len(map_probe.fullMap[binmap2])
            maskmap = (binmap1 * binmap2) > 0
            if flagsize == 1:
                return np.sum(maskmap), np.sum(binmap1), np.sum(binmap2)
            if not minim == 0.0:
                return float(len(map_target.fullMap[maskmap]))/minim
            else:
                print('Check map contour!!')
                return 0.0
        else:
            print('@@@ Maps could not be matched')
            return -1.0

    def SCCC(
            self,
            map_target,
            resolution_densMap,
            sigma_map,
            structure_instance,
            rigid_body_structure,
            write=False,
            c_mode=False,
    ):
        """Calculate local cross-correlation for a segment of a protein model

        Segments are defined using a rigid-body structure instance, which can
        be generated using `RIBFIND <https://ribfind.ismb.lon.ac.uk/>`_.
        The CCC is calculated using
        :meth:`CCC_map <TEMPy.protein.scoring_functions.ScoringFunctions.CCC_map>`.

        Implementation discussed in more detail at the following reference:
            * `Pandurangan et al., J Struct Biol. 2013 Dec 12 <https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3988922/>`_

        Args:
            map_target: Map Instance.
            resolution_densMap: Resolution of the input map.
            sigma_map: Parameter need for Structure Blurrer. Full explanation
                :ref:`here<Note on real space blurring:>`
            structure_instance: Structure instance to compare
            rigid_body_structure: Rigid-body Structure instance that defines
                the segment.
            write: If True, a string is returned with format :code:`f'SCCC for segment {SCCC_score}'`
            c_mode: No longer in use
        Returns:
            float if :code:`write=False`, else string:
                SCCC Score
        """
        blurrer = StructureBlurrer()
        scorer = ScoringFunctions()
        outline = ""
        resolution_densMap = float(resolution_densMap)
        whole_fit_map = blurrer.gaussian_blur(
            structure_instance,
            resolution_densMap,
            densMap=map_target,
            sigma_coeff=sigma_map,
            normalise=True,
        )
        sim_map = blurrer.gaussian_blur(
            rigid_body_structure,
            resolution_densMap,
            densMap=map_target,
            sigma_coeff=sigma_map,
            normalise=True,
        )
        minDens = sim_map.std()
        sim_mask_array = sim_map._get_maskArray(minDens)
        mask_emMap = map_target._get_maskMap(sim_mask_array)
        mask_simMap = whole_fit_map._get_maskMap(sim_mask_array)
        sse_lccf, ov = scorer.CCC_map(mask_emMap, mask_simMap, cmode=c_mode)
        if write is True:
            outline += 'SCCC for segment %f\n' % (sse_lccf)
            return outline
        return sse_lccf

    def SCCC_LAP(
            self,
            map_target,
            resolution_densMap,
            sigma_map,
            structure_instance,
            rigid_body_structure,
            write=False,
    ):
        """Calculate the Laplacian filtered cross correlation for a segment of
        a protein model.

        Segments are defined using a rigid-body structure instance, which can
        be generated using `RIBFIND <https://ribfind.ismb.lon.ac.uk/>`_.
        Score calculated using
        :meth:`laplace_CCC <TEMPy.protein.scoring_functions.ScoringFunctions.laplace_CCC>`

        Arguments:
            map_target: Map Instance.
            resolution_densMap: Resolution of the map instance.
            sigma_map: Parameter need for Structure Blurrer. Full explanation
                :ref:`here<Note on real space blurring:>`
            structure_instance: Structure instance to compare
            rigid_body_structure: Rigid-body Structure instance.
            write: If True, a string is returned with format
                :code:`f'SCCC_LAP for segment {SCCC_LAP_score}'`
        Return:
            float if :code:`write=False`, else string:
                SCCC_LAP score
        """
        blurrer = StructureBlurrer()
        scorer = ScoringFunctions()
        outline = ""
        resolution_densMap = float(resolution_densMap)
        whole_fit_map = blurrer.gaussian_blur(
            structure_instance,
            resolution_densMap,
            densMap=map_target,
            sigma_coeff=sigma_map,
            normalise=True,
        )
        sim_map = blurrer.gaussian_blur(
            rigid_body_structure,
            resolution_densMap,
            densMap=map_target,
            sigma_coeff=sigma_map,
            normalise=True,
        )
        minDens = sim_map.std()
        sim_mask_array = sim_map._get_maskArray(minDens)
        mask_emMap = map_target._get_maskMap(sim_mask_array)
        mask_simMap = whole_fit_map._get_maskMap(sim_mask_array)
        sse_lccf = scorer.laplace_CCC(mask_emMap, mask_simMap)
        if write is True:
            outline += 'SCCC for segment %f\n' % (sse_lccf)
            return outline
        return sse_lccf

    def SCCC_MI(
            self,
            map_target,
            resolution_densMap,
            sigma_map,
            structure_instance,
            rigid_body_structure,
            write=False,
    ):
        """Calculate the mutual information score for a segment of a protein
        model.

        Segments are defined using a rigid-body structure instance, which can
        be generated using `RIBFIND <https://ribfind.ismb.lon.ac.uk/>`_.
        The score for the segment is calculated using
        :meth:`MI <TEMPy.protein.scoring_functions.ScoringFunctions.MI>`

        Arguments:
            map_target: Map Instance.
            resolution_densMap: Resolution of the map instance.
            sigma_map: Parameter need for Structure Blurrer. Full explanation
                :ref:`here<Note on real space blurring:>`
            structure_instance: Structure instance to compare
            rigid_body_structure: Rigid-body Structure instance.
            write: If True, a string is returned with format
                :code:`f'SCCC_LAP for segment {SCCC_LAP_score}'`
        Return:
            float if :code:`write=False`, else string:
                SCCC_MI score
        """
        blurrer = StructureBlurrer()
        scorer = ScoringFunctions()
        outline = ""
        resolution_densMap = float(resolution_densMap)
        whole_fit_map = blurrer.gaussian_blur(
            structure_instance,
            resolution_densMap,
            densMap=map_target,
            sigma_coeff=sigma_map,
            normalise=True
        )
        sim_map = blurrer.gaussian_blur(
            rigid_body_structure,
            resolution_densMap,
            densMap=map_target,
            sigma_coeff=sigma_map,
            normalise=True
        )
        minDens = sim_map.std()
        sim_mask_array = sim_map._get_maskArray(minDens)
        mask_emMap = map_target._get_maskMap(sim_mask_array)
        mask_simMap = whole_fit_map._get_maskMap(sim_mask_array)
        sse_lccf = scorer.MI(mask_emMap, mask_simMap)
        if write is True:
            outline += 'SCCC for segment %f\n' % (sse_lccf)
            return outline
        return sse_lccf

    def get_sccc(
            self,
            map_target,
            structure_instance,
            rigid_body_file,
            resolution,
            sigma_map=0.187
    ):
        """Read rigid body file and calculate sccc scores for each rigid body

        Args:
            map_target: Map instance
            structure_instance: Model structure instance
            rigid_body_file: Path to rigid body file
            resolution: Resolution of map instance
            sigma_map: Parameter need for Structure Blurrer. Full explanation
                :ref:`here<Note on real space blurring:>`

        Returns:
            Two dictionaries; dict_chain_scores, dict_rigid_scores
                * dict_chain_scores is a dictionary with keys equal to chain
                  labels (e.g. :code:`['A', 'B', 'C']`). Each value is another
                  dictionary, with keys as residue number for each residue in a
                  rigid body, and values the CCC score for that residue.

                * dict_rigid_scores is a dictionary with keys equal to each
                  rigid body (string?? not sure) and values the CCC score for
                  that body.

        """
        dict_rbf = RBParser.get_rigid_body_str_instances(
            rigid_body_file,
            structure_instance,
        )
        dict_chain_details = structure_instance.get_dict_str()
        dict_chain_scores = {}
        dict_rigid_scores = {}
        minscore = 1.0
        maxscore = -1.0
        for ch in dict_chain_details:
            dict_res_scores = {}
            for res in dict_chain_details[ch]:
                dict_res_scores[res] = -2.0
            if ch not in dict_rbf:
                continue
            for rb in dict_rbf[ch]:
                res_list = []
                for seg in dict_rbf[ch][rb][0]:
                    try:
                        for r in range(int(seg[0]), int(seg[1]) + 1):
                            res_list.append(r)
                    except TypeError:
                        continue
                score_SCCC = self.SCCC(
                    map_target,
                    resolution,
                    sigma_map,
                    structure_instance,
                    dict_rbf[ch][rb][1]
                )
                if rb not in dict_rigid_scores:
                    dict_rigid_scores[rb] = score_SCCC
                    if score_SCCC < minscore:
                        minscore = score_SCCC
                    if score_SCCC > maxscore:
                        maxscore = score_SCCC
                for res in res_list:
                    dict_res_scores[res] = score_SCCC

            dict_chain_scores[ch] = dict_res_scores

        minscore = max(minscore - (maxscore - minscore) / 5., -1.0)
        for ch in dict_chain_scores:
            for res in dict_chain_scores[ch]:
                if dict_chain_scores[ch][res] == -2.0:
                    dict_chain_scores[ch][res] = minscore

        return dict_chain_scores, dict_rigid_scores

    def calc_moc(
            self,
            indices,
            map_probe,
            map_target,
    ):
        map_target_mask = map_target.fullMap[indices]
        map_probe_mask = map_probe.fullMap[indices]
        num = np.sum(map_target_mask * map_probe_mask)
        den = np.sqrt(
            np.sum(np.square(map_target_mask)) *
            np.sum(np.square(map_probe_mask))
        )
        if den == 0.0:
            return -1.0
        return num / den

    def calc_ccc(
            self,
            indices,
            map_probe,
            map_target
    ):
        map_target_mask = map_target.fullMap[indices]
        map_target_mask = map_target_mask - float(
            map_target_mask.sum()/len(map_target_mask)
        )
        map_probe_mask = map_probe.fullMap[indices]
        map_probe_mask = map_probe_mask - float(
            map_probe_mask.sum()/len(map_probe_mask)
        )
        num = np.sum(map_target_mask * map_probe_mask)
        den = np.sqrt(
            np.sum(np.square(map_target_mask)) *
            np.sum(np.square(map_probe_mask))
        )
        if den == 0.0:
            return -1.0
        return num / den

    def indices_smoc(
            self,
            score_indices,
            grid_indices
    ):
        tmplist = score_indices[:]
        setlist = set(tmplist)
        indices = list(setlist)
        sc_indices = []
        for ii in indices:
            sc_indices.append(grid_indices[ii])
        array_indices = np.array(sc_indices)
        ind_arrxyz = np.transpose(array_indices)
        ind_arrzyx = (ind_arrxyz[2], ind_arrxyz[1], ind_arrxyz[0])
        return ind_arrzyx

    def set_score_as_bfactor(
            self,
            structure_instance,
            dict_scores,
            default=0.0,
            chid=None,
            outfile='score.pdb',
    ):
        # include scores as b-factor records
        for x in structure_instance.atomList:
            cur_chain = x.chain
            cur_res = x.get_res_no()
            if cur_chain in dict_scores:
                try:
                    x.temp_fac = "{0:.2f}".format(
                        float(dict_scores[cur_chain][cur_res])
                    )
                except (KeyError, IndexError):
                    print('Residue missing: ', cur_res, chid)
                    x.temp_fac = default
            else:
                x.temp_fac = default
        structure_instance.write_to_PDB(outfile)

    def smoc_score_rigid_bodies(
            self,
            dict_res_indices,
            rigid_body_file,
            ch,
            indi,
            map_target,
            sim_map
    ):
        dict_res_scores = {}
        dict_rbf = RBParser.read_rigid_body_file(rigid_body_file)
        chainid_set = ch
        if len(ch) == 0:
            chainid_set = '-'
        if chainid_set in dict_rbf:
            for rb in dict_rbf[chainid_set]:
                # each rigid body
                indices = []
                res_list = []
                for seg in dict_rbf[chainid_set][rb]:
                    try:
                        for r in range(int(seg[0]), int(seg[1]) + 1):
                            indices.extend(dict_res_indices[r])
                            res_list.append(r)
                    except TypeError:
                        continue
                # minimum value set as 0 here (assuming only positive values!)
                if len(indices) == 0:
                    for res in res_list:
                        dict_res_scores[res] = 0.0  # -0.99
                    continue
                try:
                    ind_arrzyx = self.indices_smoc(indices, indi)
                    smoc = self.calc_moc(ind_arrzyx, sim_map, map_target)
                except IndexError:
                    smoc = 0.0
                # save scores
                for res in res_list:
                    dict_res_scores[res] = smoc
        return dict_res_scores

    def smoc_get_indices_fragment_window(self, dict_res_indices,
                                         dict_chain_res, ch, res,
                                         win):
        indices = dict_res_indices[res][:]
        # consider residues on both sides. NOTE: wont work for insertion codes!
        # need to rewite res numbers to avoid insertion codes
        for ii in range(1, int(round((win + 1) / 2))):
            try:
                res_idx = dict_chain_res[ch].index(res) - ii
                if res_idx < 0:
                    continue
                # get prev residue indices
                indices.extend(
                    dict_res_indices[
                        dict_chain_res[ch]
                        [
                            res_idx
                        ]
                    ]
                )
            except (KeyError, IndexError):
                pass
        for ii in range(1, int(round((win + 1) / 2))):
            try:
                indices.extend(
                    dict_res_indices[
                        dict_chain_res[ch]
                        [
                            dict_chain_res[ch].index(res) + ii
                        ]
                    ]
                )
            except (IndexError, KeyError):
                pass
        return indices

    def smoc_get_indices_sphere(self, gridtree, coord, dist=5.0):
        list_points = gridtree.query_ball_point(
            [coord[0], coord[1], coord[2]],
            dist
        )
        return list_points

[docs] def SMOC(self, map_target, resolution_densMap, structure_instance=None, win=11, rigid_body_file=None, sim_map=None, sigma_map=0.225, write=False, c_mode=True, sigma_thr=2.5, fragment_score=True, dist=5.0, atom_centre='CA', calc_metric='smoc', get_coord=True ): """Calculate the local Mander's Overlap for each residue in a protein model SMOC can be calculated in two ways: * **SMOCf**: The SMOC score is calculated for residues adjacent in sequence, using a sliding and overlapping window. * **SMOCd**: Pixels are scored in a spherical radius around each residue Arguments: map_target: Target Map Instance. resolution_densMap: Resolution of the target map. structure_instance: Model structure instance. win: Overlapping Window length to calculate the score rigid_body_file: Path to rigid-body file. sim_map: Precomputed simulated map. If :code:`None`, sim_map is calculated automatically. sigma_map: Parameter need for Structure Blurrer. Full explanation :ref:`here<Note on real space blurring:>` write: Deprecated, not used. c_mode: Deprecated, not used. sigma_thr: Parameter used to label pixels near each atom in the model. Explained in further detail :meth:`here <TEMPy.protein.structure_blurrer.StructureBlurrer.get_indices>` fragment_score: If True, use SMOCf method. If False, use SMOCd method. dist: Deprecated, not used. atom_centre: Which atom type should be considered the centre of residues. calc_metric: Which method to use for calculating the local score at each residue. Can be :code:`"smoc"` or :code:`"sccc"` get_coord: Return additional dict_chain_CA dictionary. Returns: Dictionary: The dictionaries dict_chain_scores, dict_chain_res. dict_chain_CA is returned additionally if get_coord = True * dict_chain_scores: 2D dictionary containing the SMOC scores for each residue. Keys are chain_label for first dictionary and residue_number for second dictionary, e.g.: :code:`dict_chain_scores[chain_label][residue_number] = SMOC_score` * dict_chain_res: Dictionary with chain labels (e.g. :code:`"A"`) as keys and a list of all residue numbers for a given chain as values. * dict_chain_CA: 2D dictionary containing a list for each residue, containing the amino acid type of each residue and its 3D coordinates. Keys are chain_label for first dictionary and residue_number for second dictionary, e.g.: :code:`dict_chain_CA[chain_label][residue_number] = residue_info` where, :code:`residue_info = [residue_type, x, y, z]` """ blurrer = StructureBlurrer() if sim_map is None: sim_map = blurrer.gaussian_blur_real_space( structure_instance, resolution_densMap, densMap=map_target, sigma_coeff=sigma_map, normalise=True) # mask out background try: peak, ave, sigma = sim_map.peak_density() if peak >= ave - sigma and peak < (ave + 5 * sigma): sim_map.fullMap = sim_map.fullMap * (sim_map.fullMap > peak) except: # noqa: E722 pass dict_chain_indices, dict_chain_res, dict_chain_CA, gridtree = \ blurrer.get_indices( structure_instance, map_target, resolution_densMap, sim_sigma_coeff=sigma_map, sigma_thr=sigma_thr, atom_centre=atom_centre) del structure_instance gc.collect() # get details of map nz, ny, nx = map_target.fullMap.shape zg, yg, xg = np.mgrid[0:nz, 0:ny, 0:nx] # x y,z coordinates of the grid indi = list(zip(xg.ravel(), yg.ravel(), zg.ravel())) list_sccc = [] # save scores for each chain and res dict_chain_scores = OrderedDict() for ch in dict_chain_indices: dict_res_scores = OrderedDict() dict_res_indices = dict_chain_indices[ch] # multi-chain rigid body parser if rigid_body_file is not None: dict_res_scores = self.smoc_score_rigid_bodies( dict_res_indices, rigid_body_file, ch, indi, map_target, sim_map) # for residues not in rigid bodies for res in dict_res_indices: # if not dict_res_scores.has_key(res): if res not in dict_res_scores: if fragment_score: try: indices = self.smoc_get_indices_fragment_window( dict_res_indices, dict_chain_res, ch, res, win) except TypeError: indices = [] # unusual cases if len(indices) > 0 and len(indices) < 10: try: dict_res_scores[res] = \ dict_res_scores[dict_chain_res[ch][dict_chain_res[ch].index(res)-1]] # noqa: E501 try: dict_res_scores[res] = \ (dict_res_scores[res] + dict_res_scores[dict_chain_res[ch][dict_chain_res[ch].index(res)+1]])/2.0 # noqa: E501 except (IndexError, KeyError): pass except (IndexError, KeyError): try: dict_res_scores[res] = \ dict_res_scores[dict_chain_res[ch][dict_chain_res[ch].index(res)+1]] # noqa: E501 except (IndexError, KeyError): dict_res_scores[res] = 0.0 continue else: indices = dict_res_indices[res][:] # for lower resolutions calculate score in a sphere if resolution_densMap > 3.5: dist = 1.0*resolution_densMap try: indices_local = \ self.smoc_get_indices_sphere( gridtree, dict_chain_CA[ch][res][1:], dist=dist) indices.extend(indices_local) except (TypeError, KeyError): pass if len(indices) == 0: print('Score calculation failed for: ', ch, res) smoc = -1.0 else: try: ind_arrzyx = self.indices_smoc(indices, indi) if calc_metric == 'smoc': smoc = self.calc_moc( ind_arrzyx, sim_map, map_target) elif calc_metric == 'sccc': smoc = self.calc_ccc( ind_arrzyx, sim_map, map_target) except IndexError: smoc = 0.0 dict_res_scores[res] = smoc # error in score calculation if smoc == -1.0: dict_res_scores[res] = 0.0 try: dict_res_scores[res] = \ dict_res_scores[dict_chain_res[ch][dict_chain_res[ch].index(res)-1]] # noqa: E501 try: dict_res_scores[res] = \ (dict_res_scores[res] + dict_res_scores[dict_chain_res[ch][dict_chain_res[ch].index(res)+1]])/2.0 # noqa: E501 except (IndexError, KeyError): pass except (IndexError, KeyError): try: dict_res_scores[res] = \ dict_res_scores[dict_chain_res[ch][dict_chain_res[ch].index(res)+1]] # noqa: E501 except (IndexError, KeyError): dict_res_scores[res] = 0.0 continue list_sccc.append(smoc) dict_chain_scores[ch] = OrderedDict() dict_chain_scores[ch] = dict_res_scores del gridtree del dict_chain_indices gc.collect() if get_coord: return dict_chain_scores, dict_chain_res, dict_chain_CA else: del dict_chain_CA return dict_chain_scores, dict_chain_res
def _SMOC1( # can we delete this? self, map_target, resolution_densMap, structure_instance, win=11, rigid_body_file=None, sigma_map=0.225, write=False, ): """ Calculate Local cross correlation (Mander's Overlap) It is a local Overlap Coefficient calculated on atoms in sliding residue windows along the chain. Arguments: *map_target* Target Map Instance. *resolution_densMap* Parameter need for Structure Blurrer. Resolution of the target map. *structure_instance* Model structure instance. *win* Overlapping Window length to calculate the score *rigid_body_file* Rigid-body file. sigma_map: Parameter need for Structure Blurrer. Full explanation :ref:`here<Note on real space blurring:>` Return: Dictionary of smoc scores for residues in the chain """ blurrer = StructureBlurrer() sim_map = blurrer.gaussian_blur_real_space( structure_instance, resolution_densMap, densMap=map_target, sigma_coeff=sigma_map, normalise=True, ) peak, ave, sigma = sim_map._peak_density() # NOTE: filter background _, dict_res_indices, dict_res_dist = blurrer.get_indices( structure_instance, map_target, resolution_densMap ) # get details of map # origin = map_target.origin # apix = map_target.apix # box_size = map_target.box_size() nz, ny, nx = map_target.fullMap.shape zg, yg, xg = np.mgrid[0:nz, 0:ny, 0:nx] indi = zip(xg.ravel(), yg.ravel(), zg.ravel()) # save rigid body details dict_rf_res = {} dict_rf_sc = {} res_list = [] rb_list = [] list_sccc = [] # save scores for each res dict_res_scores = {} r_ct = 0 if rigid_body_file is not None: inp = open(rigid_body_file, 'r') for ln in inp: if ln[0] != '#': score_indices = [] lrb = ln.split() if len(lrb) == 0: continue r_ct += 1 res_list = [] rb_pairs = [] # get scores for each res and each rigid body for i in range(max((len(lrb) / 2) - 1, 1)): rb_pairs.append([int(lrb[2 * i]), int(lrb[2 * i + 1])]) # NOTE: wont work for insertion codes for r in range(int(lrb[2 * i]), int(lrb[2*i + 1]) + 1): score_indices.extend(dict_res_indices[r]) res_list.append(r) rb_list.append(lrb) dict_rf_res[r_ct] = rb_pairs if len(score_indices) == 0: dict_rf_sc[r_ct] = 0.0 for res in res_list: dict_res_scores[res] = 0.0 continue tmplist = score_indices[:] setlist = set(tmplist) score_indices = list(setlist) sc_indices = [] for ii in score_indices: sc_indices.append(indi[ii]) array_indices = np.array(sc_indices) ind_arrxyz = np.transpose(array_indices) ind_arrzyx = (ind_arrxyz[2], ind_arrxyz[1], ind_arrxyz[0]) sccc = self.calc_moc(ind_arrzyx, sim_map, map_target) dict_rf_sc[r_ct] = sccc # save scores for res in res_list: dict_res_scores[res] = sccc list_sccc.append(sccc) inp.close() # for residues not in rigid bodies: consider pentapeptides for res in dict_res_indices: if res in dict_res_scores: indices = dict_res_indices[res][:] # consider residues on both sides. NOTE: wont work for # insertion codes! # need to rewite res numbers to avoid insertion codes for ii in range(1, int(round((win + 1) / 2))): try: indices.extend(dict_res_indices[res - ii]) except Exception: pass for ii in range(1, int(round((win + 1) / 2))): try: indices.extend(dict_res_indices[res + ii]) except Exception: pass tmplist = indices[:] setlist = set(tmplist) indices = list(setlist) sc_indices = [] for ii in indices: sc_indices.append(indi[ii]) if len(indices) == 0: dict_res_scores[res] = 0.0 continue array_indices = np.array(sc_indices) ind_arrxyz = np.transpose(array_indices) ind_arrzyx = (ind_arrxyz[2], ind_arrxyz[1], ind_arrxyz[0]) sccc = self.calc_moc(ind_arrzyx, sim_map, map_target) dict_res_scores[res] = sccc list_sccc.append(sccc) return dict_res_scores def _get_shell(self, dist1, maxlevel, step, x): # indices between upper and lower shell bound return np.argwhere((dist1 < min(maxlevel, x + step)) & (dist1 >= x)) # match power spectra for two maps def _amplitude_match( self, map_1, map_2, shellmin, shellmax, step=0.005, c1=0, c2=0, reso=None, lpfiltb=False, lpfilta=False, ref=False, ): """ Scale amplitudes to the average in each resolutions shell Arguments: *step : shell width (1/A) """ # fourier transform: use pyfftw if available pyfftw_flag = 1 try: import pyfftw except ImportError: pyfftw_flag = 0 try: if pyfftw_flag == 0: raise ImportError inputa1 = pyfftw.n_byte_align_empty( map_1.fullMap.shape, 16, 'complex128', ) outputa1 = pyfftw.n_byte_align_empty( map_1.fullMap.shape, 16, 'complex128', ) # fft planning, set planning_timelimit or flags to make it faster fft = pyfftw.FFTW( inputa1, outputa1, direction='FFTW_FORWARD', axes=(0, 1, 2), flags=['FFTW_ESTIMATE'], ) inputa1[:, :, :] = map_1.fullMap[:, :, :] fft() ft1 = Map( fftshift(outputa1), map_1.origin, map_1.apix, map_1.filename, map_1.header[:], ) except Exception: # use numpy fft instead ft1 = map_1.fourier_transform() try: if pyfftw_flag == 0: raise ImportError inputa2 = pyfftw.n_byte_align_empty( map_2.fullMap.shape, 16, 'complex128', ) outputa2 = pyfftw.n_byte_align_empty( map_2.fullMap.shape, 16, 'complex128', ) fft = pyfftw.FFTW( inputa2, outputa2, direction='FFTW_FORWARD', axes=(0, 1, 2), flags=['FFTW_ESTIMATE'], ) inputa2[:, :, :] = map_2.fullMap[:, :, :] fft() ft2 = Map( fftshift(outputa2), map_2.origin, map_2.apix, map_2.filename, map_2.header[:], ) except Exception: ft2 = map_2.fourier_transform() # low pass filter before scaling if reso is not None: cutoff1 = map_1.apix / float(reso) cutoff2 = map_2.apix / float(reso) if lpfiltb and not lpfilta: ft1._tanh_lowpass(cutoff1, fall=0.2, ftmap=True) ft2._tanh_lowpass(cutoff2, fall=0.2, ftmap=True) # size1 = max(map_1.x_size(),map_1.y_size(),map_1.z_size()) dist1 = map_1._make_fourier_shell(1)/map_1.apix[0] # size2 = max(map_2.x_size(), map_2.y_size(), map_2.z_size()) dist2 = map_2._make_fourier_shell(1) / map_2.apix[0] ft1_avg = [] ft2_avg = [] ft1_avg_new = [] lfreq = [] # select max spatial frequency to iterate to. low resolution map maxlevel = 0.5 / np.max((map_1.apix[0], map_2.apix[0]), axis=0) # loop over freq shells, shellwidth=0.005 # for x in arange(0,maxlevel+step,step): nc = 0 x = 0.0 highlevel = x + step while (x < maxlevel): # print x,highlevel, maxlevel # indices between upper and lower shell bound fshells1 = ((dist1 < min(maxlevel, highlevel)) & (dist1 >= x)) # radial average shellvec1 = ft1.fullMap[fshells1] # indices between upper and lower shell bound fshells2 = ((dist2 < min(maxlevel, highlevel)) & (dist2 >= x)) # radial average shellvec2 = ft2.fullMap[fshells2] abs1 = abs(shellvec1) abs2 = abs(shellvec2) ns1 = len(np.nonzero(abs1)[0]) ns2 = len(np.nonzero(abs2)[0]) if ns1 < 10 or ns2 < 10: nc += 1 highlevel = min(maxlevel, x + (nc + 1) * step) x = max(0.0, x - nc * step) continue else: nc = 0 mft1 = np.mean(abs1) mft2 = np.mean(abs2) if mft1 == 0.0 and mft2 == 0.0: continue # sq of radial avg amplitude ft1_avg.append(np.log10(np.mean(np.square(abs1)))) ft2_avg.append(np.log10(np.mean(np.square(abs2)))) # scale to amplitudes of the ref map if ref: if mft1 == 0.0: continue ft1.fullMap[fshells1] = shellvec1 * (mft2 / mft1) else: # replace with avg amplitudes for the two maps ft1.fullMap[fshells1] = shellvec1 * (mft2 + mft1) / (2 * mft1) ft2.fullMap[fshells2] = shellvec2 * (mft2 + mft1) / (2 * mft2) # new radial average (to check) mft1 = np.mean(abs(ft1.fullMap[fshells1])) ft1_avg_new.append( np.log10( np.mean( np.square(abs(ft1.fullMap[fshells1])) ) ) ) lfreq.append(highlevel) sampling_frq = highlevel cutoff_freq = min((1.0 / reso) + 0.25, maxlevel) # scale the rest and break after relevant frequencies if sampling_frq > cutoff_freq: fshells1 = (dist1 >= highlevel) shellvec1 = ft1.fullMap[fshells1] mft1 = np.mean(abs(shellvec1)) fshells2 = (dist2 >= highlevel) shellvec2 = ft2.fullMap[fshells2] mft2 = np.mean(abs(shellvec2)) if mft1 == 0.0 and mft2 == 0.0: break ft1_avg.append(np.log10(np.mean(np.square(abs(shellvec1))))) ft2_avg.append(np.log10(np.mean(np.square(abs(shellvec2))))) if ref: if mft1 == 0.0: break ft1.fullMap[fshells1] = shellvec1*(mft2/mft1) else: ft1.fullMap[fshells1] = shellvec1*(mft2+mft1)/(2*mft1) ft2.fullMap[fshells2] = shellvec2*(mft2+mft1)/(2*mft2) mft1 = np.mean(abs(ft1.fullMap[fshells1])) ft1_avg_new.append( np.log10( np.mean( np.square(abs(ft1.fullMap[fshells1])) ) ) ) lfreq.append((highlevel + step / 2)) break x = highlevel highlevel = x + step # low pass filter after? # low pass filter before scaling if reso is not None: if lpfilta and not lpfiltb: ft1._tanh_lowpass(cutoff1, fall=0.2, ftmap=True) ft2._tanh_lowpass(cutoff2, fall=0.2, ftmap=True) try: if pyfftw_flag == 0: raise ImportError ifft = pyfftw.FFTW( inputa1, outputa1, direction='FFTW_BACKWARD', axes=(0, 1, 2), flags=['FFTW_ESTIMATE'], ) inputa1[:, :, :] = ifftshift(ft1.fullMap)[:, :, :] ifft() map1_filt = Map( outputa1.real.astype('float'), map_1.origin, map_1.apix, map_1.filename, map_1.header[:], ) except Exception: # use numpy ifft instead map1_filt = map_1.copy() map1_filt.fullMap = np.real(ifftn(ifftshift(ft1.fullMap))) try: if pyfftw_flag == 0: raise ImportError ifft = pyfftw.FFTW( inputa2, outputa2, direction='FFTW_BACKWARD', axes=(0, 1, 2), flags=['FFTW_ESTIMATE'], ) inputa2[:, :, :] = ifftshift(ft2.fullMap)[:, :, :] ifft() map2_filt = Map( outputa2.real.astype('float'), map_2.origin, map_2.apix, map_2.filename, map_2.header[:], ) except Exception: map2_filt = map_2.copy() map2_filt.fullMap = np.real(ifftn(ifftshift(ft2.fullMap))) try: import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt from matplotlib import pylab try: plt.style.use('ggplot') except AttributeError: pass plt.rcParams.update({'font.size': 18}) plt.rcParams.update({'legend.fontsize': 18}) plt.plot(lfreq, ft1_avg, 'r--', label='map1') plt.plot(lfreq, ft2_avg, 'bs', label='map2') plt.plot(lfreq, ft1_avg_new, 'g^', label='scaled') leg = plt.legend(loc='upper right') for legobj in leg.legendHandles: legobj.set_linewidth(2.0) pylab.savefig("spectra.png") plt.close() except Exception: pass return map1_filt.fullMap, map2_filt.fullMap def get_clash_map(self, emmap, apix): template_grid = emmap._make_clash_map(apix) return template_grid def get_sm_score( self, struct, ncomp, template_grid, cvol, apix, ): """Calculate the clash score between chains within a protein model The clash score calculates the amount that chains in a protein overlap one another. The structure is mapped into discrete bins in 3D space, defined by the template_grid input. The score is proportional to the number of these bins occupied by more than one of the protein chains. Recommended pixel spacing of the template_grid is 3.5 Å. Args: struct: Structure instance. ncomp: Number of chains in the structure. template_grid: Map object with empty 3D-grid (filled with zeros). This can be calculated using :meth:`StructureBlurrer.protMap <TEMPy.protein.structure_blurrer.StructureBlurrer.protMap>` cvol: List containing the volume of each chain. This can be calculated using :meth:`Map._get_component_volumes <TEMPy.maps.em_map.Map._get_component_volumes>` apix: Pixel size of template_grid, used for scaling the volume calculations. Returns: The clash score for the structure (close to 0 is good, large negative score is bad) """ overlay_maplist = [] overlay_maplist = self.get_overlay_comp_maplist(struct, template_grid) nc = range(ncomp) cpair = list(itertools.combinations(nc, 2)) score = 0.0 n_overlap_voxel = 0 overlap_volume = 0.0 for i in cpair: n_overlap_voxel = ( overlay_maplist[i[0]].fullMap * overlay_maplist[i[1]].fullMap ).sum() overlap_volume = ((apix.prod()) * n_overlap_voxel) * 2 clash_percent = ( float(overlap_volume / (cvol[i[0]] + cvol[i[1]])) ) score = score + clash_percent return -(score) def get_overlay_comp_maplist(self, struct, template_grid): blurrer = StructureBlurrer() overlay_maplist = [] ssplit = struct.split_into_chains() for x in ssplit: overlay_maplist.append( blurrer.make_atom_overlay_map1( template_grid, x, ) ) return overlay_maplist def calculate_alcps(self, ref_chain, probe_chain): """ Calculate the Arc Length Component Placement Score between two models. Generally, this score is calculated between individual chains in a structure. Inputs: ref_chain: One of the input structure instances being compared probe_chain: The other input structure instance being compared Output: alcps: Arc Length Component Placement Score """ diff_vec = ref_chain.CoM - probe_chain.CoM r = diff_vec.mod() ref_matrix = vector.get_coordinate_mass_array(ref_chain) probe_matrix = vector.get_coordinate_mass_array(probe_chain) rotation_matrix, _t = vector.get_transform(ref_matrix, probe_matrix) theta, n1, n2, n3 = vector._rotmat_to_axisangle(rotation_matrix) theta = math.degrees(theta) alcps = 2 * math.pi * r * theta / 360 return alcps def get_log10_alcps(self, ref_chain, probe_chain, lowest_value=-2): """ Calculate the log10 of the ALCPS score. Inputs: ref_chain: Reference structure instance probe_chain: Probe structure instance lowest_value: Value to use if two structures are perfectly aligned (i.e. if alcps = 0) - default value is -2 Outputs: log10(alcps): log10 of the ALCPS lowest_value: returned if ALCPS = 0, default value is -2 """ alcps = self.calculate_alcps(ref_chain, probe_chain) if alcps == 0.0: return lowest_value else: return math.log10(alcps) def _get_max_distance_between_residues(atoms): max_dist = 0 atom_list = atoms[:] for atom in atoms: atom_list.pop(atom_list.index(atom)) for atom2 in atom_list: vec = np.linalg.norm([atom.x, atom.y, atom.z], [atom2.x, atom2.y, atom2.z]) if vec > max_dist: max_dist = vec return max_dist def map_model_fsc_resolution( self, structure, target_map, map_resolution, mask=None, sim_map=None, cutoff_correlation=0.5, use_highest_resolution=True): """Calculate the Map-model "resolution" using Fourier Shell Correlation (FSC). Arguments: structure: A model structure instance. target_map: A Map instance. map_resolution: The resolution of target_map mask: A Map instance of a soft-edged mask. sim_map: A Map instance of a pre-calculated simulated map, if not explicitly supplied it will be calculated automatically from the structure instance. cutoff_correlation: Threshold correlation used to define the model-map "resolution". use_highest_resolution: If True, the highest resolution shell that crosses the cutoff_correlation threshold is taken. Otherwise, the frequency that first crosses the threshold is taken. Returns: Map-model FSC "resolution": The fourier shell frequency (in angstroms) where the correlation crosses the stated threshold. """ if sim_map is None: blurrer = StructureBlurrer(with_vc=True) sim_map = blurrer.gaussian_blur_real_space( structure, map_resolution, densMap=target_map, ) fsc = self.map_model_fsc(structure, target_map, map_resolution, mask=mask, sim_map=sim_map) max_frequency_at_cutoff = fourier.get_fsc_resolution( fsc, fsc_threshold=cutoff_correlation, interpolate_cutoff=False, use_highest_resolution=use_highest_resolution ) return max_frequency_at_cutoff def map_model_fsc(self, structure, target_map, map_resolution, mask=None, sim_map=None): """Calculate the map-model Fourier Shell Correlation (FSC). Arguments: structure: A model structure instance. target_map: A Map instance. map_resolution: The resolution of target_map mask: A Map instance of a soft-edged mask. sim_map: A Map instance of a pre-calculated simulated map, if not explicitly supplied it will be calculated automatically from the structure instance Returns: The map-model FSC as an np.array with shape (2, (N / 2) + 1) for an input map with shape (N, N, N). Output has format: [[shell frequency (angstroms), correlation] ........ [shell frequency (angstroms), correlation]] """ if sim_map is None: blurrer = StructureBlurrer(with_vc=True) sim_map = blurrer.gaussian_blur_real_space( structure, map_resolution, densMap=target_map, ) if mask: fsc = fourier.calculate_fsc_with_mask( target_map.fullMap, sim_map.fullMap, mask, pixel_size=target_map.apix[0], ) else: fsc = fourier.calculate_fsc( target_map.fullMap, sim_map.fullMap, pixel_size=target_map.apix[0], ) return fsc def LoQFit( self, structure, target_map, resolution_map_target, sim_map=None, max_res=18., verbose=True, half_map_1=None, half_map_2=None, ): """Calculate the LoQFit score for each residue in a model. The LoQFit score is essentially a local implementation of the map-model FSC. More information can be found here: https://www.biorxiv.org/content/10.1101/2022.06.08.495259v1 Arguments: structure: Model structure instance target_map: Experimental Map instance resolution_map_target: The resolution of the target map sim_map: Pre-calculated simulated map - will be calculated automatically if not supplied max_res: The maximum allowed LoQFit score for any residue, this parameter prevents verbose: If True, prints a progress bar and warnings half_map_1: A map instance of a half-map from reconstruction of the target map. If supplied (with second half map), this will be used to normalise the LoQFit score half_map_2: A map instance of the second half map from target_map reconstruction. Returns: The LoQFit score for each residue in structure instance. The output is a dictionary with keys: (chain_id, residue_number) where each value is the LoQFit score for a given residue. If half_map_1 and half_map_2 arguments are supplied the LoQFit score will be normalised based on the local resolution. """ print(f"Calculating LoQFit score for model {structure.filename} into " f"map {target_map.filename}") # set up some stuff target_map.apix = np.array(target_map.apix) blurrer = StructureBlurrer(with_vc=True) if sim_map is None: sim_map = blurrer.gaussian_blur_real_space( structure, resolution_map_target, densMap=target_map) fsc_threshold = 0.5 sim_map.normalise_to_1_minus1(in_place=True) target_map.normalise_to_1_minus1(in_place=True) # approximate mask sphere diameter sphere_diameter_angstroms = resolution_map_target * 5 sphere_diameter_pixels = int(sphere_diameter_angstroms / target_map.apix[0]) # ensure that parameters are sensible i.e. not a 1 pixel wide mask if sphere_diameter_pixels < 6: sphere_diameter_pixels = 6 cos_edge = int(round(sphere_diameter_pixels / 3)) if cos_edge < 10: cos_edge = 10 if half_map_1 is not None and half_map_2 is not None: interpolate_cutoff = False normalise_with_half_map = True half_map_1.origin = target_map.origin half_map_2.origin = target_map.origin else: interpolate_cutoff = True normalise_with_half_map = False # get an appropriate box size to crop from main map # will be 2^n, 3^n or 5^n size where n is the number required for # the minimum box size needed to encompass entire mask in box mask_dims = fourier.get_fft_optimised_box_size( sphere_diameter_pixels + (2*cos_edge), target_map.fullMap.shape) if mask_dims == [0, 0, 0]: mask = Map(np.ones(target_map.fullMap.shape), target_map.origin, target_map.apix, target_map.filename) mask_dims = mask.fullMap.shape else: centre_of_mask_map = [mask_dims[0]/2, mask_dims[1]/2, mask_dims[2]/2] mask = blurrer.make_spherical_mask(sphere_diameter_pixels, mask_dims, centre_of_mask_map, cos_edge) # save scores for each chain and res dict_res_scores = OrderedDict() chains = structure.split_into_chains() number_of_residues = structure.number_of_residues() # set up pyfftw for very fast ffts input1 = pyfftw.empty_aligned(mask_dims, dtype='float64') fftw_object1 = fourier.get_pyfftw_object(mask_dims) fftw_object2 = fourier.get_pyfftw_object(mask_dims) # precompute these to speed up calculations qbins, labelled_shells = fourier.get_labelled_shells( input1.shape) warnings = [] res_no = 0 for chain in chains: for atom in chain.atomList: if atom.atom_name == 'CA' or atom.atom_name == "C1'": pass else: continue cropped_map_1 = target_map.get_cropped_box_around_atom( atom, mask_dims, ).fullMap cropped_map_2 = sim_map.get_cropped_box_around_atom( atom, mask_dims, ).fullMap fsc = fourier.calculate_fsc_with_mask( cropped_map_1, cropped_map_2, mask.fullMap, pixel_size=target_map.apix[0], labelled_shells=labelled_shells, qbins=qbins, fftw_object1=fftw_object1, fftw_object2=fftw_object2, ) # find highest frequency shell with 0.5 correlation freq_at_corr_cutoff = fourier.get_fsc_resolution( fsc, fsc_threshold=fsc_threshold, interpolate_cutoff=interpolate_cutoff, use_highest_resolution=True, ) # freq_at_cutoff can = 0 if correlation never drops below 0.5 if freq_at_corr_cutoff < min(target_map.apix) * 2: freq_at_corr_cutoff = min(target_map.apix) * 2 warnings.append(f"Warning: FSC at residue {atom.res_no} in" f" chain {atom.chain} never falls below " f"0.5. LoQFit score for this residue set " f"to {freq_at_corr_cutoff}. This is NOT " f"expected and NOT an indication of a " f"good fit between map and model.") # If there is poor correlation in low resolution shells the # freq_at_cutoff can be ridiculously high (i.e. >100 angstroms) # which ruins the scaling when the score is plotted if freq_at_corr_cutoff > max_res: warnings.append(f"Warning: Residue {atom.res_no} in chain " f"{atom.chain} has local resolution of " f" {freq_at_corr_cutoff}, which is set " f"to {max_res}") freq_at_corr_cutoff = max_res if normalise_with_half_map: crop_hmap_1 = half_map_1.get_cropped_box_around_atom( atom, mask_dims, ) crop_hmap_2 = half_map_2.get_cropped_box_around_atom( atom, mask_dims, ) fsc = fourier.calculate_fsc_with_mask( crop_hmap_1.fullMap, crop_hmap_2.fullMap, mask.fullMap, pixel_size=target_map.apix[0], labelled_shells=labelled_shells, qbins=qbins, fftw_object1=fftw_object1, fftw_object2=fftw_object2, ) local_res = fourier.get_fsc_resolution( fsc, fsc_threshold=0.5, interpolate_cutoff=False, use_highest_resolution=True, ) freq_at_corr_cutoff = freq_at_corr_cutoff / local_res dict_res_scores[(atom.chain, int(atom.res_no))] = \ freq_at_corr_cutoff res_no += 1 if verbose: zfill_len = len(str(number_of_residues)) print(f"\rProcessed residue {str(res_no).zfill(zfill_len)}" f"/{number_of_residues}", sep='', end='', flush=True) if verbose: print("") for warning in warnings: print(warning) return dict_res_scores def map_local_resolution( self, half_map_1, half_map_2, structure=None, protein_mask=None, verbose=True, fsc_threshold=0.5, ): """Calculate the local resolution using the local FSC method. Arguments: half_map_1: A map instance of the first half-map. half_map_2: A map instance of the second half map. structure: Model structure instance. If defined, the local resolution is calculated at the CA atom for each residue in the structure only. protein_mask: A Map instance of a mask. If defined, the local resolution is calculated at each point within the mask only. verbose: If True, prints a progress bar and warnings fsc_threshold: Correlation threshold value that defines how the map resolution is found. Returns: The local resolution at each position in the map, or as defined by the structure instance or mask. The output is a dictionary where each key has the format: (z, y, x) where z, y, x is a 3D-index from the half maps, and the value is the local resolution (in angstroms). """ print(f"Calculating Local Resolution for half map:" f"{half_map_1.filename} and map {half_map_2.filename}") # set up some stuff apix = np.array(half_map_1.apix)[0] map_res = self.map_resolution_FSC(half_map_1, half_map_2) blurrer = StructureBlurrer() # approximate mask sphere diameter sphere_diameter_angstroms = map_res * 7 sphere_diameter_pixels = int(sphere_diameter_angstroms/apix) # ensure that parameters are sensible i.e. not a 1 pixel wide mask if sphere_diameter_pixels < 6: sphere_diameter_pixels = 6 cos_edge = int(round(sphere_diameter_pixels / 3)) if cos_edge < 10: cos_edge = 10 # get an appropriate box size to crop from main map # will be 2^n, 3^n or 5^n size where n is the number required for # the minimum box size needed to encompass entire mask in box mask_dims = fourier.get_fft_optimised_box_size( sphere_diameter_pixels + (cos_edge * 2), half_map_1.fullMap.shape, ) if mask_dims == [0, 0, 0]: mask = Map(np.ones(half_map_1.fullMap.shape), half_map_1.origin, half_map_1.apix, half_map_1.filename) mask_dims = mask.fullMap.shape else: centre_of_mask_map = [mask_dims[0]/2, mask_dims[1]/2, mask_dims[2]/2] mask = blurrer.make_spherical_mask(sphere_diameter_pixels, mask_dims, centre_of_mask_map, cos_edge) # save scores for each chain and res dict_res_scores = OrderedDict() # set up pyfftw for very fast ffts input1 = pyfftw.empty_aligned(mask_dims, dtype='float64') fftw_object1 = fourier.get_pyfftw_object(mask_dims) fftw_object2 = fourier.get_pyfftw_object(mask_dims) qbins, labelled_shells = fourier.get_labelled_shells( input1.shape) # get coordinates where local FSC should be calculated if structure is not None: coordinates = structure.get_CAonly().get_atom_zyx_map_indexes( half_map_1) elif protein_mask is not None: (z, y, x) = np.where(protein_mask.fullMap > 0) coordinates = np.stack((z, y, x), axis=-1) else: (z, y, x) = np.where(half_map_1.fullMap < np.inf) coordinates = np.stack((z, y, x), axis=-1) warnings = [] n = 0 for coordinate in coordinates: cropped_map_1 = half_map_1.get_cropped_box( coordinate, mask_dims, ) cropped_map_2 = half_map_2.get_cropped_box( coordinate, mask_dims, ) fsc = fourier.calculate_fsc_with_mask( cropped_map_1.fullMap, cropped_map_2.fullMap, mask.fullMap, pixel_size=apix, labelled_shells=labelled_shells, qbins=qbins, fftw_object1=fftw_object1, fftw_object2=fftw_object2, ) # find highest frequency shell with 0.5 correlation freq_at_corr_cutoff = fourier.get_fsc_resolution( fsc, fsc_threshold=fsc_threshold, interpolate_cutoff=False, use_highest_resolution=True, ) # FSC can go below nyquist frequency when shells never drop # below 0.5 correlation and freq_at_corr_cutoff still equals 0 if freq_at_corr_cutoff < min(half_map_1.apix) * 2: freq_at_corr_cutoff = min(half_map_1.apix) * 2 warnings.append(f"Warning: FSC at coordinate {coordinate}" f" never falls below 0.5. SLRS " f"score for this residue set to " f"{freq_at_corr_cutoff}. This is NOT " f"expected and NOT an indication of a good " f"fit between map and model.") dict_res_scores[tuple(coordinate)] = freq_at_corr_cutoff n += 1 if verbose: zfill_len = len(str(len(coordinates))) print(f"\rProcessed residue {str(n).zfill(zfill_len)}" f"/{len(coordinates)}", sep='', end='', flush=True) if verbose: print("") for warning in warnings: print(warning) return dict_res_scores def map_resolution_FSC(self, half_map1, half_map2, mask=None, fsc_threshold=0.143, use_highest_resolution=True): """Calculate the map resolution from two half-maps using the Fourier Shell Correlation. Arguments: half_map1: Map instance of first half map half_map2: Map instance of second half map mask: A Map instance of a soft edged mask fsc_threshold: Correlation threshold value that defines how the map resolution is found. use_highest_resolution: If True, takes the highest resolution shell that crosses the fsc_threshold, if there are multiple crossovers. Returns: The FSC resolution (in angstroms). """ if mask is not None: fsc = fourier.calculate_fsc_with_mask( half_map1.fullMap, half_map2.fullMap, mask.fullMap, pixel_size=half_map1.apix[0],) else: fsc = fourier.calculate_fsc(half_map1.fullMap, half_map2.fullMap, pixel_size=half_map1.apix[0],) resolution = fourier.get_fsc_resolution( fsc, fsc_threshold=fsc_threshold, use_highest_resolution=use_highest_resolution, interpolate_cutoff=False) return resolution class SlidingResidueWindow: def __init__(self, accumulator): self.accumulator = accumulator def scores(self, residue_atoms_list, size): self.accumulator.zero() half_window = int((size - 1) / 2) scores = [] for head in range(0, len(residue_atoms_list) + half_window): tail = head - size # While the head is with in the allowed range, use all the voxels # associated with the residues atoms. if head < len(residue_atoms_list): for atom in residue_atoms_list[head]: self.accumulator.add_atom(atom) # Drop the tail if it within the range of residues. if tail >= 0: for atom in residue_atoms_list[tail]: self.accumulator.del_atom(atom) # Extract the SMOC score if head - half_window >= 0: scores.append(self.accumulator.get_val()) return scores class _AbstractSlidingScore: def __init__(self, struct, accumulator): mapping = {} for a in struct.atomList: if mapping.get(a.chain) is None: mapping[a.chain] = {} if mapping[a.chain].get(a.res_no) is None: mapping[a.chain][a.res_no] = [] mapping[a.chain][a.res_no].append(a) self.mapping = mapping self.slider = SlidingResidueWindow(accumulator) def _get_contigs(self, chain): contigs = [[]] for res_num in self.mapping[chain]: last = contigs[-1] if len(last) == 0: last.append(res_num) else: if last[-1] + 1 == res_num: last.append(res_num) else: contigs.append([res_num]) return contigs def score_chain_contig(self, chain, size): """Generate sliding window scores where windows do not span chain breaks." This method computes sliding windows over each contiguous region of the chain. Args: chain (str): The chain name to generate window scores. size (int): The size of the window. Should be an odd number. Returns: A dictionary mapping each residue to its score. """ scores = {} res_mapping = self.mapping[chain] contigs = self._get_contigs(chain) for contig in contigs: residue_atoms_list = [res_mapping[r] for r in contig] contig_scores = self.slider.scores(residue_atoms_list, size) for k, v in zip(contig, contig_scores): scores[k] = v return scores def score_chain_span(self, chain, size): """Generate sliding window scores where windows do span chain breaks. This method ignores any breaks in the chain, so windows may include distant residues. Args: chain (str): The chain name to generate window scores. size (int): The size of the window. Should be an odd number. Returns: A dictionary mapping each residue to its score. """ residue_mapping = self.mapping[chain] residue_atoms_list = list(residue_mapping.values()) numbers = residue_mapping.keys() scores = self.slider.scores(residue_atoms_list, size) return {k: v for k, v in zip(numbers, scores)} class _Accumulator: """An accumulator computes a score over a sliding window of atoms""" def add_atom(self, atom): pass def del_atom(self, atom): pass def get_val(self): pass def zero(self): pass class _MOCAccumulator(_Accumulator): """Placeholder docstring """ def __init__(self, exp_map, sim_map, radius): from math import ceil, log, pow import voxcov as vc max_dim = max(exp_map.box_size()) next_pow2 = int(pow(2, ceil(log(max_dim) / log(2)))) self.radius = radius self.moc = vc.SMOC( exp_map.apix, exp_map.origin, next_pow2, [exp_map.x_size(), exp_map.y_size(), exp_map.z_size()], np.ascontiguousarray(exp_map.fullMap, dtype="float64"), np.ascontiguousarray(sim_map.fullMap, dtype="float64"), ) def add_atom(self, atom): self.moc.add_sphere([atom.x, atom.y, atom.z], self.radius) def del_atom(self, atom): self.moc.del_sphere([atom.x, atom.y, atom.z], self.radius) def get_val(self): return self.moc.get_val() def zero(self): self.moc.zero() class _MIAccumulator(_Accumulator): """Incremental calculation of MI score over voxels covered by atoms.""" def __init__(self, exp_map, sim_map, radius, x_bins, y_bins): from math import ceil, log, pow import voxcov as vc max_dim = max(exp_map.box_size()) next_pow2 = int(pow(2, ceil(log(max_dim) / log(2)))) self.radius = 3 self.exp_map = exp_map self.sim_map = sim_map self.mi = vc.SMI( exp_map.apix, exp_map.origin, next_pow2, [exp_map.x_size(), exp_map.y_size(), exp_map.z_size()], np.ascontiguousarray(sim_map.fullMap, dtype="float64"), np.ascontiguousarray(exp_map.fullMap, dtype="float64"), np.min(exp_map.fullMap), np.max(exp_map.fullMap), np.min(sim_map.fullMap), np.max(sim_map.fullMap), x_bins, y_bins ) def add_atom(self, atom): self.mi.add_sphere( [atom.x, atom.y, atom.z], self.radius, ) def del_atom(self, atom): self.mi.del_sphere( [atom.x, atom.y, atom.z], self.radius, ) # Not strictly essential, but I should add zeroing to DiffMap in voxcov. def zero(self): self.mi.zero() # Return current MI value def get_val(self): return self.mi.get_val()
[docs]class LocalMI(_AbstractSlidingScore): """ """
[docs] def __init__(self, struct, exp_map, sim_map, resolution, radius, no_bins=20): sigma_coeff=0.225 self.exp_map = exp_map self.sim_map = sim_map self.radius = radius self.no_bins = no_bins if sim_map is None: sim_map = StructureBlurrer().gaussian_blur_real_space( struct, resolution, exp_map, sigma_coeff=sigma_coeff ) mi = _MIAccumulator(exp_map, sim_map, radius, no_bins, no_bins) super().__init__(struct, mi)
[docs]class FastSMOC(_AbstractSlidingScore): """A faster version of the sliding window SMOC score. Unlike the original algorithm, this ones is optimized in the following ways: * lazily determines required voxels saving space and time. * The sliding-window calculation is performed in time proportional only to chain size and not window size. This method does not offer all the bells and whistles of old SMOC function such as rigid body parameters. However, two alternative methods are provided for dealing with SMOC scores for unmodelled regions. The function `score_chain_contig`, will score all contigs in a chain where each window only contains residues belonging to the contig. The function `score_chain_span` mimicks the old SMOC behaviour. Windows ignore contigs boundaries (they 'span' them), thus residues in windows may be physically distant. In general, the consider using the `score_chain_contig` method. Args: struct (:class:`BioPy_Structure <TEMPy.protein.prot_rep_biopy.BioPy_Structure>`): Model to score exp_map (:class:`EMMap <TEMPy.maps.em_map.Map>`): Experimental density map resolution (float): The resolution of the experimental map sigma_coeff (float, optional): The sigma coefficient for simulating the map. Defaults to 0.225. sigma_threshold (float, optional): The sigma threshold determines the radius around the atoms to consider. Defaults to 2.5. sim_map (EMMap, optional): A simulated map to use instead of generating one. Returns: FastSMOC object """
[docs] def __init__( self, struct, exp_map, resolution, sigma_coeff=0.225, sigma_threshold=2.5, sim_map=None, locres=None, ): radius = sigma_threshold * max(sigma_coeff * resolution, 1.0) if sim_map is None: sim_map = StructureBlurrer(locres=locres).gaussian_blur_real_space( struct, resolution, exp_map, sigma_coeff=sigma_coeff ) moc = _MOCAccumulator(exp_map, sim_map, radius) super().__init__(struct, moc)
[docs]class FastSCCC:
[docs] def __init__( self, structure, map_target, resolution, sigma_coeff=0.225, ): """Calculate local cross-correlation for a segment of a protein model. Args: structure: The entire atomic structure map_target: An experimental map. resolution: The sesolution of the experimental map. sigma_coeff: The sigma coefficient for blurring. :ref:`here<Note on real space blurring:>` Returns: A SCCC object. """ self.blurrer = StructureBlurrer() self.scorer = ScoringFunctions() self.map_target = map_target self.resolution = resolution self.sigma_coeff = sigma_coeff self.whole_sim_map = self.blurrer.gaussian_blur_real_space( structure, resolution, densMap=map_target, sigma_coeff=sigma_coeff, normalise=True, )
def score_segment(self, segment_structure): """Score a segment of the structure. Args: segment: An atomic structure which should be a component Returns: CCC score: (float) """ segment_sim_map = self.blurrer.gaussian_blur_real_space( segment_structure, self.resolution, self.map_target, self.sigma_coeff, normalise=True, ) minDens = segment_sim_map.std() segment_mask = segment_sim_map._get_maskArray(minDens) target_masked = self.map_target._get_maskMap(segment_mask) sim_masked = self.whole_sim_map._get_maskMap(segment_mask) sse, _ = self.scorer.CCC_map(target_masked, sim_masked, cmode=False) return sse