# Path: pyem/geom/quat_numba.py
# @numba.jit(cache=False, nopython=True, nogil=True)
# def distq(q1, q2):
#     pi_half = np.pi / 2
#     v = np.abs(np.sum(q1 * q2))
#     v = np.clip(0, 1.0, v)
#     v = np.arccos(v)
#     v *= 2
#     # msk = v > pi_half
#     # v[msk] = np.pi - v[msk]
#     return v
# 
# Path: pyem/geom/quat_numba.py
# @numba.jit(cache=True, nopython=True, nogil=True)
# def qslerp(q1, q2, t, longest=False):
#     cos_half_theta = np.dot(q1, q2)
#     if cos_half_theta >= 1.0:
#         return q1.copy()
#     if longest:
#         if cos_half_theta > 0:
#             cos_half_theta = -cos_half_theta
#             q1 = -q1
#     elif cos_half_theta < 0:
#         cos_half_theta = -cos_half_theta
#         q1 = -q1
#     half_theta = np.arccos(cos_half_theta)
#     sin_half_theta = np.sqrt(1 - cos_half_theta * cos_half_theta)
#     if np.abs(sin_half_theta) < 1E-12:
#         return (q1 + q2) / 2
#     a = np.sin((1 - t) * half_theta)
#     b = np.sin(t * half_theta)
#     return (q1 * a + q2 * b) / sin_half_theta
# 
# Path: pyem/geom/quat_numba.py
# @numba.guvectorize(["void(float64[:], float64[:], float64[:])"],
#                    "(m),(m)->(m)", nopython=True, cache=False)
# def qtimes(q1, q2, q3):
#     _qtimes(q1, q2, q3)
# 
# Path: pyem/geom/quat.py
# def meanq(q, w=None):
#     if w is None:
#         return np.linalg.eigh(np.einsum('ij,ik->...jk', q, q))[1][:, -1]
#     else:
#         return np.linalg.eigh(np.einsum('ij,ik,i->...jk', q, q, w))[1][:, -1]


import numpy as np
from .quat_numba import distq
from .quat_numba import qslerp
from .quat_numba import qtimes
from .quat import meanq
#
# Library of miscellaneous utility functions.
# See help text and README file for more information.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.


def double_center(arr, reference=None, inplace=False):
    assert arr.ndim == 2
    if reference is None:
        mu0 = np.mean(arr, axis=0, keepdims=True)
        mu = np.mean(arr)
    else:
        mu0 = np.mean(reference, axis=0, keepdims=True)
        mu = np.mean(reference)
    mu1 = np.mean(arr, axis=1, keepdims=True)
    if inplace:
        arr -= mu0
        arr -= mu1

# Path: pyem/geom/quat_numba.py
# @numba.jit(cache=False, nopython=True, nogil=True)
# def distq(q1, q2):
#     pi_half = np.pi / 2
#     v = np.abs(np.sum(q1 * q2))
#     v = np.clip(0, 1.0, v)
#     v = np.arccos(v)
#     v *= 2
#     # msk = v > pi_half
#     # v[msk] = np.pi - v[msk]
#     return v
# 
# Path: pyem/geom/quat_numba.py
# @numba.jit(cache=True, nopython=True, nogil=True)
# def qslerp(q1, q2, t, longest=False):
#     cos_half_theta = np.dot(q1, q2)
#     if cos_half_theta >= 1.0:
#         return q1.copy()
#     if longest:
#         if cos_half_theta > 0:
#             cos_half_theta = -cos_half_theta
#             q1 = -q1
#     elif cos_half_theta < 0:
#         cos_half_theta = -cos_half_theta
#         q1 = -q1
#     half_theta = np.arccos(cos_half_theta)
#     sin_half_theta = np.sqrt(1 - cos_half_theta * cos_half_theta)
#     if np.abs(sin_half_theta) < 1E-12:
#         return (q1 + q2) / 2
#     a = np.sin((1 - t) * half_theta)
#     b = np.sin(t * half_theta)
#     return (q1 * a + q2 * b) / sin_half_theta
# 
# Path: pyem/geom/quat_numba.py
# @numba.guvectorize(["void(float64[:], float64[:], float64[:])"],
#                    "(m),(m)->(m)", nopython=True, cache=False)
# def qtimes(q1, q2, q3):
#     _qtimes(q1, q2, q3)
# 
# Path: pyem/geom/quat.py
# def meanq(q, w=None):
#     if w is None:
#         return np.linalg.eigh(np.einsum('ij,ik->...jk', q, q))[1][:, -1]
#     else:
#         return np.linalg.eigh(np.einsum('ij,ik,i->...jk', q, q, w))[1][:, -1]


# Path: pyem/geom/geom.py
import numpy as np
from .quat_numba import distq
from .quat_numba import qslerp
from .quat_numba import qtimes
from .quat import meanq

#
# Library of miscellaneous utility functions.
# See help text and README file for more information.
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.


def double_center(arr, reference=None, inplace=False):
    assert arr.ndim == 2
    if reference is None:
        mu0 = np.mean(arr, axis=0, keepdims=True)
        mu = np.mean(arr)
    else:
        mu0 = np.mean(reference, axis=0, keepdims=True)
        mu = np.mean(reference)
    mu1 = np.mean(arr, axis=1, keepdims=True)
    if inplace:
        arr -= mu0
        arr -= mu1

        arr += mu