#! /usr/bin/env starta
/*--*- Java -*---------------------------------------------------------------*\
**$Author: saulius $
**$Date: 2025-07-02 17:59:59 +0000 (Wed, 02 Jul 2025) $ 
**$Revision: 12285 $
**$URL: file:///home/saulius/svn-repositories/paskaitos/VU/software/bioinformatics-III-assignments/trunk/bin/pdbhull $
\*---------------------------------------------------------------------------*/
//*
// Starta programa, kuri kiekvienam PDB(x) failui suskaičiuoja jo
// atomų iškiląjį apvalkalėlį. Algoritmas: "gift wrapping".
//
// Programos iškvietimas:
// pdbhull < 1zyx.pdb
// pdbhull 1zyx.pdb 2zyx*.cif
//**
    
pragma append  "$P/lib/starta";
pragma prepend "$P/lib";
    
use std;
use Strings;
use BiomoleculeAtom;
use PDBline;
use PDBxline;
use Det3_4;
use Vectors;
use EdgeSet;

var null_Atom = new Atom;

var Id = "$Id: pdbhull 12285 2025-07-02 17:59:59Z saulius $";

. "#", Id[1:-2];

// XAtom is an atom that can be sorted by x coordinates. If x
// coordinates are equal, then y coordinates are taken into account,
// and finally z.

class XAtom : Atom {
    constructor() {}
    constructor fromAtom( Atom a )
    {
        self := a;
    }
    operator ">" ( XAtom a; XAtom b ): bool
    {
        return a.x > b.x ? true :
            (a.x < b.x ? false :
             (a.y > b.y ? true : (a.y < b.y ? false : a.z > b.z)
              )
             );
    }
    operator "<" ( XAtom a; XAtom b ): bool
    {
        return a.x < b.x ? true :
            (a.x > b.x ? false :
             (a.y < b.y ? true : (a.y > b.y ? false : a.z < b.z)
              )
             );
    }
};

var null_XAtom = new XAtom;

// Flag which file format we are currently reading:
type format = enum byte( PDB, PDBx );
var current_format = PDB format;

exception FORMAT_ERROR;

var atom = new XAtom();
var molecule: array of XAtom = new XAtom[0];

var long n = 1;

var bool in_text_field = false;

while(<>) {
    if( !strstart( ";", $_ ) && in_text_field ) {
        continue;
    }
    if( strstart( ";", $_ )) {
        if( in_text_field ) {
            in_text_field = false;
        } else {
            in_text_field = true;
        }
        continue;
    }
    if( strstart( "data_", $_ )) {
        current_format = PDBx format;
        continue;
    }
    if( strstart( "HEADER", $_ ) || strstart( "CRYST1", $_ ) ) {
        current_format = PDB format;
        continue;
    }
    if( current_format == PDB format && strstart( "ATOM  ", $_ ) ||
        current_format == PDBx format && strstart( "ATOM", $_ ) ||
        strstart( "HETATM", $_ )) {
        
        if( current_format == PDB format ) then
            atom = new XAtom.fromAtom( PDBline::split( $_ ));
        elsif( current_format == PDBx format ) then
            atom = new XAtom.fromAtom( PDBxline::split( $_ ));
        else
            raise FORMAT_ERROR( "unsupported format %d" %% current_format@int );
        endif;
        molecule = push( molecule, atom );
        atom.atomNumber = n;
        n++;
    }
}

function cmpx( a, b: XAtom ): int
{
    return a.x < b.x ? -1 : 1;
}

function cmpy( a, b: XAtom ): int
{
    return a.y < b.y ? -1 : 1;
}

function cmpz( a, b: XAtom ): int
{
    return a.z < b.z ? -1 : 1;
}

function minmax( a: array of XAtom; cmp: function(XAtom a; XAtom b)->(int) ): XAtom, XAtom, int, int
{
    var min = a[0];
    var max = a[0];
    var imin, imax: int;

    for var i = 1 to last(a) {
            if( cmp(min, a[i]) > 0 ) {
                min = a[i];
                imin = i;
            }
            if( cmp(max, a[i]) < 0 ) {
                max = a[i];
                imax = i;
            }
        }
    return min, max, imin, imax
}

var top, bottom, left, right, front, back: XAtom =
    null_XAtom, null_XAtom, null_XAtom, null_XAtom, null_XAtom, null_XAtom;

var itop, ibottom, ileft, iright, ifront, iback: int;

bottom, top, ibottom, itop = minmax( molecule, cmpx );
left, right, ileft, iright = minmax( molecule, cmpy );
front, back, ifront, iback = minmax( molecule, cmpz );

// . ">>> bot:", bottom->asString();
// . ">>> top:", top->asString();
// . ">>> left:", left->asString();
// . ">>> rght:", right->asString();
// . ">>> frnt:", front->asString();
// . ">>> back:", back->asString();
// . "";

// Eliminate all atoms inside the top-bottom/left-right/front-back
// "diamond":

var outsiders = new XAtom[0];

var seen_atoms = new bool[length(molecule)];

for( var i in [ileft, iright, itop, ibottom, ifront, iback] ) {
    if( !seen_atoms[i] ) {
        outsiders = push( outsiders, molecule[i] );
        seen_atoms[i] = true;
    }
}

function abs( x : double ): double
{
    return x >= 0.0D ? x : -x;
}

{
    var vtop    = [ top.x, top.y, top.z, 1.0D ];
    var vbottom = [ bottom.x, bottom.y, bottom.z, 1.0D ];
    var vleft   = [ left.x, left.y, left.z, 1.0D ];
    var vright  = [ right.x, right.y, right.z, 1.0D ];
    var vfront  = [ front.x, front.y, front.z, 1.0D ];
    var vback   = [ back.x, back.y, back.z, 1.0D ];

    for( var a in molecule ) {
        if( a == top || a == bottom ||
            a == left || a == right ||
            a == front || a == back ) {
            continue;
        }
        var va = [ a.x, a.y, a.z, 1.0D ];
        var d1 = det4( [ vfront, vright, vtop, va ] );
        var d2 = det4( [ vright, vback, vtop, va ] );
        var d3 = det4( [ vback, vleft, vtop, va ] );
        var d4 = det4( [ vleft, vfront, vtop, va ] );

        var d5 = det4( [ vfront, vright, vbottom, va ] );
        var d6 = det4( [ vright, vback, vbottom, va ] );
        var d7 = det4( [ vback, vleft, vbottom, va ] );
        var d8 = det4( [ vleft, vfront, vbottom, va ] );

        // . ">>> dets:", d1, d2, d3, d4, d5, d6, d7, d8, a->asString();

        var eps = 1E-4D;
        
        pragma double;
        if( d1 <= 0 || d2 <= 0 || d3 <= 0 || d4 <= 0 ||
            d5 >= 0 || d6 >= 0 || d7 >= 0 || d8 >= 0 ||
            abs(d1) < eps || abs(d2) < eps || abs(d3) < eps || abs(d4) < eps ||
            abs(d5) < eps || abs(d6) < eps || abs(d7) < eps || abs(d8) < eps ) {
            outsiders = push( outsiders, a );
        }
    }
}

// for( var a in outsiders ) {
//     . ">>>", a->asString();
// }

import HeapSort(XAtom) as SortAtomsByX;

var sorted_atoms_by_x = SortAtomsByX::sorted( outsiders );

var a2 = sorted_atoms_by_x[0];
var a2i = 0;

// The first flip axis:
var a3 = new XAtom.fromAtom( a2 );
a3.y += 1.0D;

// Another atom to add to the hull, potentially:
var a1 = sorted_atoms_by_x[1];
var a1i = 1;

var v1 = [a1.x,a1.y,a1.z,1.0D];
var v2 = [a2.x,a2.y,a2.z,1.0D];
var v3 = [a3.x,a3.y,a3.z,1.0D];

for var i = 2 to last(sorted_atoms_by_x) {
        var a4 = sorted_atoms_by_x[i];
        var v4 = [a4.x,a4.y,a4.z,1.0D];

        var d = det4( [v1, v2, v3, v4] );

        if( d < 0.0D ) {
            a1 = a4;
            v1 = v4;
            a1i = i;
        }
    };

var used_edges = new EdgeSet( length(sorted_atoms_by_x) );
var potential_edges = new EdgeSet( length(sorted_atoms_by_x) );

potential_edges->add_edge([a1i,a2i]);
used_edges->add_edge([a1i,a2i]);

var hull = new XAtom[0];
var in_hull = new bool[length(sorted_atoms_by_x)];

in_hull[a1i] = true;
in_hull[a2i] = true;

hull = push( hull, a1 );
hull = push( hull, a2 );

procedure find_edge_in_array( int[][][] edges ): int[]
{
    for( var edge_array in edges ) {
        if( length(edge_array) > 0 ) {
            return edge_array[0];
        }
    }
    return null;
}

procedure find_edge( EdgeSet es ): int[]
{
    return find_edge_in_array( es->get_edge_table() );
}

function lookup_atom_not_in_hull( XAtom[] atoms; bool[] in_hull ): XAtom, int
{
    for( var i = 0; i < length(atoms); i++ ) {
        if( !in_hull[i] ) {
            return atoms[i], i;
        }
    }
    return null_XAtom, -1
}

repeat {
    var b1i: int;
    var b1: XAtom = null_XAtom;
    b1, b1i = lookup_atom_not_in_hull( sorted_atoms_by_x, in_hull );
    if( b1i < 0 ) {
        break;
    }

    var flip_edge = find_edge( potential_edges );

    potential_edges->remove_edge( flip_edge );
    used_edges->add_edge( flip_edge );

    var b2i = flip_edge[0];
    var b3i = flip_edge[1];
    var b2 = sorted_atoms_by_x[b2i];
    var b3 = sorted_atoms_by_x[b3i];

    v1 = [b1.x,b1.y,b1.z,1.0D];
    v2 = [b2.x,b2.y,b2.z,1.0D];
    v3 = [b3.x,b3.y,b3.z,1.0D];
    for var i = 1 to last(sorted_atoms_by_x) {
            if( i == b1i || i == b2i || i == b3i ) {
                continue;
            }
            var b4 = sorted_atoms_by_x[i];
            var v4 = [b4.x,b4.y,b4.z,1.0D];
            var d = det4( [v1, v2, v3, v4] );
            if( d < 0.0D ) {
                b1 = b4;
                b1i = i;
                v1 = v4;
            }
        }
    if( !in_hull[b1i] ) {
        in_hull[b1i] = true;
        hull = push( hull, b1 );
    }
    // Insert ordered edges, with the directions as specified below:
    var e1 = [b2i,b1i];
    var e2 = [b1i,b3i];
    if( !used_edges->has_edge( e1 ) ) {
        potential_edges->add_edge( e1 );
    }
    if( !used_edges->has_edge( e2 ) ) {
        potential_edges->add_edge( e2 );
    }
} while( potential_edges->edge_count() > 0 );

for( var a in hull ) {
    . a->asString();
}