import numbers
import traceback
from typing import List, Tuple, Union

import numpy as np

from geogst.core.geometries.shapes.points import Point
from geogst.core.utils.inspections.errors import Error
from geogst.core.utils.inspections.functions import caller_name
from geogst.core.mathematics.arrays import singular_value_decomposition
from geogst.core.mathematics.vectors3d import Vect3D
from geogst.core.geology.orientations import Plane, Direct


def calculate_best_fit_plane_from_coordinates(
    coordinates: List[Tuple[numbers.Real, numbers.Real, numbers.Real]],
) -> Tuple[Union[None, Tuple[Plane, Point, np.ndarray]], Error]:

    try:

        xyz_array = np.array(
            coordinates,
            dtype=np.float64
        )

        xyz_mean = np.mean(
            xyz_array,
            axis=0
        )

        (_, sorted_eigenvalues, sorted_eigenvectors), err = singular_value_decomposition(
            xyz_array=xyz_array - xyz_mean
        )

        if err:
            return None, err

        lowest_eigenvector = sorted_eigenvectors[-1, : ]  # Solution is last row
        normal = lowest_eigenvector[: 3 ] / np.linalg.norm(lowest_eigenvector[: 3 ])
        normal_vector = Vect3D(normal[0], normal[1], normal[2])
        direction_normal_to_best_fit_plane, err = Direct.from_vector(normal_vector)

        if err:
            return None, err

        best_fit_geological_plane = direction_normal_to_best_fit_plane.normal_plane()
        central_point = Point(xyz_mean[0], xyz_mean[1], xyz_mean[2])

        return (best_fit_geological_plane, central_point, sorted_eigenvalues), Error()

    except Exception as e:

        return None, Error(
            True,
            caller_name(),
            e,
            traceback.format_exc()
        )
