package SmilesScripts::Chirality;

use strict;
use warnings;

use parent Exporter::;
our @EXPORT = qw(
    chiral_volume
    order_points_around_center
    square_planar
    trigonal_bipyramidal
);

use List::Util qw( all sum );
use Math::MatrixReal;
use Scalar::Util qw( blessed );

my $PI = 4 * atan2( 1, 1 );

sub chiral_volume
{
    my @coords = @_[0..3];
    my $normalise_vectors = $_[4];

    if(      all { blessed $_ && $_->isa( 'Chemistry::Atom' ) } @coords ) {
        @coords = map { $_->coords } @coords;
    } elsif( all { blessed $_ && $_->isa( 'Math::VectorReal' ) } @coords ) {
        # No need to do anything
    } else {
        die "chiral_volume() called with parameters of unknown type\n";
    }

    my( $a, $b, $c ) = map { $_ - $coords[0] } @coords[1..3];
    ( $a, $b, $c ) = map { $_->norm } ( $a, $b, $c ) if $normalise_vectors;
    return $a . ( $b x $c );
}

# Project the given 3D points to a minimal distance plane.
# Return the order of points in which they are "encountered" around their center.
sub order_points_around_center
{
    my( @points3D ) = @_;

    my( $A, $B, $C, $D ) = plane_equation( @points3D );
    my @points2D = project_3D_to_plane( $A, $B, $C, $D, @points3D );
    my $reference = shift @points2D; # first point is the reference

    my @angles = map { $_ >= 0 ? $_ : $_ + 2*$PI }
                 map { atan2 $_->[1] - $reference->[1], $_->[0] - $reference->[0] }
                     @points2D;

    return sort { $angles[$a] <=> $angles[$b] } 0..$#angles;
}

sub plane_equation
{
    my( @points ) = @_;

    my( $A, $B, $C, $D );
    if( @points == 3 ) {
        # From https://math.stackexchange.com/questions/2686606/equation-of-a-plane-passing-through-3-points
        $A = $points[0]->[1] * $points[1]->[2] - $points[1]->[2] * $points[2]->[1];
        $B = $points[0]->[2] * $points[2]->[0] - $points[0]->[0] * $points[2]->[2];
        $C = $points[1]->[0] * $points[2]->[1] - $points[0]->[1] * $points[1]->[0];
        $D = $points[0]->[0] * $points[1]->[1] * $points[2]->[2] -
             $points[0]->[2] * $points[1]->[1] * $points[2]->[0];
    } else {
        # Calculate the centroid of points
        my @centroid = ( sum( map { $_->[0] } @points ) / @points,
                         sum( map { $_->[1] } @points ) / @points,
                         sum( map { $_->[2] } @points ) / @points );

        # Calculate the covariance matrix
        my $xx = sum map { ($_->[0] - $centroid[0])**2 } @points;
        my $yy = sum map { ($_->[1] - $centroid[1])**2 } @points;
        my $zz = sum map { ($_->[2] - $centroid[2])**2 } @points;
        my $xy = sum map { ($_->[0] - $centroid[0]) * ($_->[1] - $centroid[1]) } @points;
        my $xz = sum map { ($_->[0] - $centroid[0]) * ($_->[2] - $centroid[2]) } @points;
        my $yz = sum map { ($_->[1] - $centroid[1]) * ($_->[2] - $centroid[2]) } @points;

        my @matrix = ( [ $xx, $xy, $xz ],
                       [ $xy, $yy, $yz ],
                       [ $xz, $yz, $zz ] );
        my $matrix = Math::MatrixReal->new_from_rows( \@matrix );

        # Find eigenvalues and eigenvectors, select smallest eigenvalue
        my( $l, $V ) = $matrix->sym_diagonalize;
        my $min = 1;
        for (2..3) {
            $min = $_ if $l->element( $min, 1 ) > $l->element( $_, 1 );
        }

        # Calculate the plane coefficients
        ( $A, $B, $C ) = @{$V->row( $min )->[0][0]};
        $D = $A * $centroid[0] + $B * $centroid[1] + $C * $centroid[2];
    }

    return ( $A, $B, $C, $D );
}

sub project_3D_to_plane
{
    my $A = shift;
    my $B = shift;
    my $C = shift;
    my $D = shift;

    my @points = @_;

    return map { [ $_->[0] - $A * ( $A * $_->[0] + $B * $_->[1] + $C * $_->[2] + $D ) / ( $A**2 + $B**2 + $C**2 ),
                   $_->[1] - $B * ( $A * $_->[0] + $B * $_->[1] + $C * $_->[2] + $D ) / ( $A**2 + $B**2 + $C**2 ) ] }
               @points;
}

sub square_planar
{
    my @points3D = @_;

    my @order = order_points_around_center( @points3D );
    while( $order[0] != 0 ) {
        push @order, shift @order;
    }

    return '@SP1' if $order[2] == 2;
    return '@SP2' if $order[2] == 1;
    return '@SP3' if $order[2] == 3;
}

sub trigonal_bipyramidal
{
    my @points3D = @_;

    my @order = order_points_around_center( @points3D );
    while( $order[0] != 0 ) {
        push @order, shift @order;
    }

    return $order[1] == 1 ? '@TB1' : '@TB2';
}

1;