
//fix
// 3 modes:
//  read first model
//  overlap
//  many models

#include "Pdb_file.h"

#include "Atom_kit.h"

//fix remove
#include "Create.h"
#include "Bond.h"
#include "Model.h"
#include "Model_kit.h"

#include "File_input.h"
#include "Log.h"
#include "User.h"
#include "Pdb_attributes.h"
#include "All_atoms.h"
#include "Arbitrary_atom_group.h"
#include "Bond_detector.h"
#include "Residue_numerator.h"
#include "Molecular_numerator.h"
#include "Aminoacid_name.h"
#include "Fixed_joint.h"
#include "Float_joint.h"
#include "C_str.h"

#include "Point_3D_impl.h"//fix remove

#include <stdio.h>
#include <string>
#include <fstream>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <map>
#include <ctype.h>
//#include <utility>

#include <stdlib.h>

using namespace std;

namespace MM
{

/*Pdb_file::Pdb_file( Model_accessor_protocol & accessor_protocol) :
   Model_file( accessor_protocol ) {}//*/
Pdb_file Pdb_file::prototype_(it_is_a_prototype);


Pdb_file::
Pdb_file (Prototype const &)
:
    Model_file (Prototype())
{
    add_extension ("pdb");
    add_extension ("ent");
}

Pdb_file::
Pdb_file (Model & model)
:
    Model_file (model)
{
    add_extension ("pdb");
    add_extension ("ent");
}

Model_file * Pdb_file::
clone (Model & model) const
{
    Pdb_file *result = new Pdb_file (model);
    
    //result->extension_ = extension_;   //fix Array += should work
    //for (int i=0;  i<extension_.size();  ++i)
    //    result->extension_.push_back (extension_[i]);
    
    return result;
}

void Pdb_file::check_file_name()
{
    std::string filename = name().c_str();
    if( filename.find( ".pdb" ) == std::string::npos &&
        filename.find( ".ent" ) == std::string::npos)

        FLAW (Text(name()) + 
              " is a wrong file name.\n"
              "Check file name extension.");
}

void Pdb_file::
//save (bool to_comment)
save (bool)
{
    //         1         2         3         4         5         6         7         8
    //12345678901234567890123456789012345678901234567890123456789012345678901234567890
    //ATOM      8  ND2 ASN A   1     -11.806   3.406  -3.543  0.00  0.00           N  
    //ATOM      9 1H   ASN A   1      -8.330   3.957   0.261  0.00  0.00           H  
    //ATOM   1541  HG  SER  1309      18.594 -15.349  26.263  1.00  0.00      1FNF1706
    //HETATM 1361 MG    MG   168       4.750  34.060  19.120  1.00 20.00      521P1458


    
    // set_substr (char *dst, const char *source, int pos, int len)
     
    int i, size;
    Text source;

    check_file_name();
    std::ofstream out (name().c_str());

    Residue_numerator   residue_numerator  (model());
    Molecular_numerator numerate_molecules (model());

    size = model().atom_count();
    for (i=0;  i<size;  ++i)
    {
        Atom const & current_atom = model().atom(i);
        Atom_kit const & kit      = current_atom.kit();
        Text const   current_element = current_atom.element ().c_str();
        char buf[] = "                                                                                ";
        char    format_buffer[13];
        int     length;
        double  value;
        Pdb_attributes          tmp_attributes;
        Pdb_attributes const *  attribute;

        if (current_atom.kit().has_pdb_attributes())
            attribute = &current_atom.kit().pdb();
        else
            attribute = &tmp_attributes;

        if (kit.has_fine_type())
        {
            char const* name = Fine_type::to (kit.fine_type_name(), "PDB", 0);
            if (name != 0)
            {
                attribute = &tmp_attributes;
                tmp_attributes.set_name (name);
            }
        }

        bool not_std_name = true;

        if (&kit.in_residue() != &Residue::none())
        {
            tmp_attributes.set_residue_number 
                (kit.in_residue().local_id().to_int());

            if (kit.in_residue().has_alias())
            {
                Text const& alias = kit.in_residue().alias();
                
                if (Aminoacid_name::is_valid (alias) 
                    || alias == "A" || alias == "C" || alias == "G" 
                    || alias == "I" || alias == "T" || alias == "U" )
                {
                    tmp_attributes.set_residue_name (alias.to_upper());
                    not_std_name = false;
                }
            }
        }

        if (attribute->name() == "" || not_std_name)
        {
            //tmp_attributes.set_name (current_atom.element ().c_str());
            source = "HETATM";
            set_substr (buf, source.c_str(), 0, 6);
        }
        else
        {
            source = "ATOM";
            set_substr (buf, source.c_str(), 0, 4);
        }

        // Atom serial number 
        source = i+1;
        set_substr (buf, source.c_str(), 11-source.length(), source.length());

        // Atom name
        source = attribute->name();
        if (source == "")
        {
            source = current_element;
            source.to_upper();
            set_substr (buf, source.c_str(), 14-source.length(), source.length());
        }
        else
        {
		    if (isdigit (source[0]) || source[0] == '*' || source[0] == '\'') 
                set_substr (buf, source.c_str(), 13-1, source.length());
            else if (current_element.length() == 2)
                set_substr (buf, source.c_str(), 13-1, source.length());
            else
                set_substr (buf, source.c_str(), 14-1, source.length());
        }
        
        // Residue name
        source = attribute->residue_name();
        set_substr (buf, source.c_str(), 18-1, source.length());

        // Residue sequence number
        if (attribute->residue_number() >= 0)
        {
            source = attribute->residue_number();
            set_substr (buf, source.c_str(), 26-source.length(), source.length());
        }

        // xyz
        value = current_atom.x();
        sprintf (format_buffer, "%.3f", value);
        length = strlen (format_buffer);
        set_substr (buf, format_buffer, 38-length, length);

        value = current_atom.y();
        sprintf (format_buffer, "%.3f", value);
        length = strlen (format_buffer);
        set_substr (buf, format_buffer, 46-length, length);

        value = current_atom.z();
        sprintf (format_buffer, "%.3f", value);
        length = strlen (format_buffer);
        set_substr (buf, format_buffer, 54-length, length);

        // Charge to -> Temperature factor
        value = current_atom.charge();
        if (value != 0.)
        {
            sprintf (format_buffer, "%.3f", value);
            length = strlen (format_buffer);
            set_substr (buf, format_buffer, 66-length, length);
        }

        // Element symbol, right-justified
        source = current_element;
        length = source.length();
        source.to_upper();
        set_substr (buf, source.c_str(), 78-length, length);

        out << buf << endl;
    }

    /*size = model().bond_count();
    for (i=0;  i<size;  ++i)
    {
    }//*/
}

/*void Pdb_file::save()
{ 

    check_file_name();

    Pdb_atoms atoms (model());

    std::ofstream out (name().c_str());

    for( int count=0;  count<atoms.atom_count();  ++count )
    {
        Atom & current_atom = atoms[count];

		std::string Temp;
		std::ostringstream ost (Temp);

        ost<< "ATOM  ";

		PutVar (count+1,ost,5,5,'r');

		ost<< " ";

		Text name;
        
        if (current_atom.kit().has_pdb_attributes())
		    name = current_atom.kit().pdb().name();
        else
            name = current_atom.element ().c_str();

		if ( ! isdigit(name[0]) ) {
			ost << " ";
			PutVar(name.c_str(),		ost,	3,3,'l');
		}
		else
			PutVar(name.c_str(),		ost,	4,4,'l');


    	ost << " ";
   		PutVar (current_atom.kit().pdb().residue_name().c_str(),		ost,	3,3,'r');

		//ost << " ";
		ost << "  ";
  		PutVar (current_atom.kit().pdb().residue_number(),		ost,	4,3,'r');
		ost << "    ";


		PutVar (current_atom.x(),			ost,	8,3,'r');
		PutVar (current_atom.y(),			ost,	8,3,'r');
		PutVar (current_atom.z(),			ost,	8,3,'r');

		out << ost.str() << endl;
	}
}//*/

void Pdb_file::
add (Text const & file_name)
{
    /*using namespace std;
    ifstream    in (name().c_str());
    string      _line;
    getline (in, _line);
    while (getline (in, _line))
    {
        _line = "";
    }//*/
    //fix check_file_name();

//     std::ifstream in( name().c_str() );
    //File_input input (name().c_str());
    File_input input (file_name.c_str());
//   if ( ! in ) FLAW (Text("Can't open file ") + name().c_str() );


    char line[81];
    int  i;
/*
     while( std::getline( in, Lbuff, '\n' ) )
     {
        std::string ATOM_Check, HETATM_Check;
        switch( Lbuff[ 0 ] )
        {
            case( 'A' ):
                ATOM_Check.assign( Lbuff.begin(), Lbuff.begin() + 4 );
                if( "ATOM" != ATOM_Check )
                    continue;
                break;

            case( 'H' ):
                HETATM_Check.assign( Lbuff.begin(), Lbuff.begin() + 6 );
                if( "HETATM" != HETATM_Check )
                    continue;
                break;
            default:
                continue;
        }


        //		string  _AtSeNu,_AtNa,_ReNa;
        

        
        char dst[81];
        substr( Lbuff.c_str(), dst, 6, 5);
        int atom_number = atoi(dst);

        substr( Lbuff.c_str(), dst, 12, 4);
        char buff[81];
        sscanf (dst,"%s",buff);
        std::string atom_name = buff;

        substr( Lbuff.c_str(), dst, 17, 3	);
        std::string residue_name = dst;

        substr( Lbuff.c_str(), dst, 21, 1);
        int chain_id = atoi(dst);

        substr( Lbuff.c_str(), dst, 22, 4);
        int residue_number = atoi(dst);
        
        substr( Lbuff.c_str(), dst, 30, 8);
        double x = atof (dst);
        substr( Lbuff.c_str(), dst, 38, 8);
        double y = atof (dst);
        substr( Lbuff.c_str(), dst, 48, 8);
        double z = atof (dst);

        
        
        int element = getElement ( Text(atom_name.c_str())) ;
        Atom * current_atom =  factory().make_atom(Element(element)) ;
        model().add_atom(current_atom);
        Pdb_atom * current_pdb_atom = factory().make_pdb_atom( *current_atom );
        current_atom->aspect().set_pdb_atom(current_pdb_atom);
        
        current_pdb_atom->set_residue_name(residue_name.c_str() );
        current_pdb_atom->set_name(atom_name.c_str() );
        current_pdb_atom->set_residue_number (residue_number);
        current_pdb_atom->set_chain_id(chain_id);
        current_pdb_atom->set_number (atom_number);

        current_pdb_atom->set_x(x);
        current_pdb_atom->set_y(y);
        current_pdb_atom->set_z(z);

   }
//   */

    // map 
    
//     while( std::gets( in, line, '\n' ) )
 /*   char word[81];
    while( input.next_line(line, 80) )
    {
        sscanf (line,"%s", first_word);

        if ( !strcmp (first_word,"ATOM") )
        {
            read_model (first_word,line,input);
        }
        else if ( !strcmp (first_word,"HETATM") )
        {
            read_model (first_word,line,input);
        }
        else if ( !strcmp (first_word,"ENDMDL") )
        {
            FLAW ("ENDMDL must be handled by read_model()");
        }
        else if ( !strcmp (first_word,"TER") )
        {
            FLAW ("TER must be handled by read_model()");
        }
        else if ( !strcmp (first_word,"SEQRES") )
        {
            // fix
        }


    
    }
//*/
    char first_word [81];
    Result result = end_of_model;

    while( result == end_of_model)
    {
    //////////////////////////////////////////////////////////////////////////////
    // fix         Intel ?!! 
    //////////////////////////////////////////////////////////////////////////////
        bool model_was_found = false;
        //bool model_was_found = find_model_line(first_word, line, input);
                //do
                while (input.next_line (line, 80))
                {
                    char dst[81];
                    substr (line, dst, 0, 6);
                    sscanf (dst, "%s", first_word);
                    //sscanf (line,"%s", first_word);
 
                    if ( !strcmp (first_word,"ATOM") 
                      || !strcmp (first_word,"HETATM") )
                    {
                        model_was_found = true;
                        break;
                    } 
                    else if (!strcmp (first_word,"REMARK")) //fix to read remarks
                    {
                        istringstream line_stream(line);
                        string word_2, word_3, word_4;
                        line_stream >> first_word >> word_2 >> word_3 >> word_4;

                        if (word_2=="@mm" && word_3=="joint")
                        {
                            if (word_4 == "fixed" || word_4 == "float")
                            {
                                int first, second, third, atom;
                                line_stream >> first >> second >> third; //fix check errors
                                
                                joint_numbers_.push_back 
                                    (Joint_numbers (word_4.c_str(), first, second, third));

                                string group;
                                line_stream >> group;
                                if (group == "of")
                                {
                                    while (line_stream >> atom)
                                        joint_numbers_.back().
                                            atoms_.push_back (atom);
                                }
                            }
                        }
                    }
                }    
                //while (input.next_line (line, 80));
    //////////////////////////////////////////////////////////////////////////////
        
        if( model_was_found )
            result = read_model (first_word,line,input);
        else
            return;
    }

    //fix to setup joints
    for (i=0;  i<joint_numbers_.size();  ++i)
    {
        if (joint_numbers_[i].type_ == "fixed")
        {
            Atom &first_atom  = model().atom (joint_numbers_[i].first_  -1);
            Atom &second_atom = model().atom (joint_numbers_[i].second_ -1);
            Atom &third_atom  = model().atom (joint_numbers_[i].third_  -1);

            Point_3D_impl                       //fix to Point_3D
                *first  = new Point_3D_impl (first_atom. position()), 
                *second = new Point_3D_impl (second_atom.position()),
                *third  = new Point_3D_impl (third_atom. position()); 

            if (joint_numbers_[i].atoms_.size() == 0)
            {
                model().kit().rigid_body().add_joint (new Fixed_joint 
                    (first, second, third, 
                     first_atom, second_atom, third_atom, 
                     new All_atoms (model())));
            }   
            else
            {
                Arbitrary_atom_group *group = new Arbitrary_atom_group;

                for (int j=0;  j<joint_numbers_[i].atoms_.size();  ++j)
                    group->add (model().atom (joint_numbers_[i].atoms_[j] -1));

                model().kit().rigid_body().add_joint (new Fixed_joint 
                    (first, second, third, 
                     first_atom, second_atom, third_atom, group));
            }
        }
        else if (joint_numbers_[i].type_ == "float")
        {
            Atom &first_atom  = model().atom (joint_numbers_[i].first_  -1);
            Atom &second_atom = model().atom (joint_numbers_[i].second_ -1);
            Atom &third_atom  = model().atom (joint_numbers_[i].third_  -1);

            if (joint_numbers_[i].atoms_.size() == 0)
            {
                model().kit().rigid_body().add_joint (new Float_joint 
                    (first_atom, second_atom, third_atom, 
                     new All_atoms (model())));
            }   
            else
            {
                Arbitrary_atom_group *group = new Arbitrary_atom_group;

                for (int j=0;  j<joint_numbers_[i].atoms_.size();  ++j)
                    group->add (model().atom (joint_numbers_[i].atoms_[j] -1));

                model().kit().rigid_body().add_joint (new Float_joint 
                    (first_atom, second_atom, third_atom, group));
            }
        }
        else
        {
            FLAW ("Wrong joint type.");
        }
    }
    joint_numbers_.clear();

    //Bond_detector bond_detector (model());
    //bond_detector.execute();

    #ifndef NDEBUG
    for (i=0;  i<model().atom_count();  ++i)
        model().atom(i).kit().set_number (-1);
    #endif
}

int Pdb_file::
get_element_N (char const word[])
{
    char char_1 = word[0]; 
    char char_2 = word[1]; 

    if (isspace (char_1))
    {
        switch (char_2)
        {
        case 'C':   return 6;
        case 'O':   return 8;
        case 'N':   return 7;
        case 'H':   return 1;
        default:
            {
                Text text (char_2);
                if (Element::is (text.c_str()))
                    return Element (text.c_str()).number();
                else
                    return Element ("Unknown").number();
            }
        }
    }
    else if (isalpha (char_1))
    {
        Text text (char_1);
        char_2 = tolower (char_2);
        text += char_2;
        if (Element::is (text.c_str()))
            return Element (text.c_str()).number();
        else
            return Element ("Unknown").number();
    }
    else 
    {
        switch (char_1)
        {
        case '1': 
        case '2': 
        case '3': 
        case '*': 
        case '\'':  
            if (char_2 == 'H')
                return 1;
        default:    
            return Element("Unknown").number();
        }
    }
}

Element Pdb_file::
get_element (char const word[])
{
    char char_1 = word[0]; 
    char char_2 = word[1]; 

    if (isspace (char_1))
    {
        switch (char_2)
        {
        case 'C':   return Element(6);
        case 'O':   return Element(8);
        case 'N':   return Element(7);
        case 'H':   return Element(1);
        default:
            {
                Text text (char_2);
                if (Element::is (text.c_str()))
                    return Element (text.c_str());
                else
                    return Element ("Unknown");
            }
        }
    }
    else if (isalpha (char_1))
    {
        Text text (char_1);
        char_2 = tolower (char_2);
        text += char_2;
        if (Element::is (text.c_str()))
            return Element (text.c_str());
        else
            return Element ("Unknown");
    }
    else
    {
        switch (char_1)
        {
        case '1':   
        case '2':   
        case '3':   
        case '*': 
        case '\'':  
            if (char_2 == 'H')
                return Element(1);
        default:    
            return Element("Unknown");
        }
    }
}

Pdb_file::Result Pdb_file::
read_model (char first_word[], char line[], Input & input)
{
/*
    char word[81];
    do 
    {
        sscanf (line,"%s", word);
        if ( !strcmp (word,"ENDMDL") )  
        {
            break;
        }
        else if ( !strcmp (word,"ATOM") 
               || !strcmp (word,"HETATM") )
        {
            read_chain(line, input);
        }
        else if ( !strcmp (word,"SIGATM") 
               || !strcmp (word,"SIGUIJ") 
               || !strcmp (word,"TER") ) 
        {
           // continue;
            FLAW ("SIGATM or SIGUIJ or TER must be handled by read_chain()");
        }
            

    } while( input.next_line(line, 80) );
//*/
//    do  find_model_line(first_word, line, input);
//    while ( ( Result result = read_model( first_word, line, input ) ) == end_of_model) ;
    model();

    Result result = end_of_chain;

    while( result == end_of_chain)
    {
        bool residue_was_found = find_chain_line(first_word, line, input);

        if( residue_was_found )
            result = read_chain( first_word, line, input );
        else
            return  end_of_file;   
    }
    
    return result;
}

Pdb_file::Result Pdb_file::
read_chain(  char first_word[], char line[], Input & input )
{
    Result result = end_of_residue;

    while( result == end_of_residue )
    {
        bool residue_was_found = find_residue_line(first_word, line, input);

        if( residue_was_found )
            result = read_residue( first_word, line, input );
        else
            return  end_of_file;   
    }

//    do  find_residue_line(first_word, line, input);
//    while ( ( Result result = read_residue( first_word, line, input ) ) == end_of_residue) ;

    return result;
}

Pdb_file::Result Pdb_file::
read_residue(  char first_word[], char line[], Input & input )
{
    int reached_residue_serial_number;
    read_atom   ( line, reached_residue_serial_number) ;

    while( input.next_line(line, 80) )
    {
        char dst[81];
        substr (line, dst, 0, 6);
        int result = sscanf (dst, "%s", first_word);
        //int result = sscanf (line,"%s", first_word);
        if (result != 1) 
            continue;
 
        if ( !strcmp (first_word,"TER") )
        {
            return end_of_chain;
        }
        else if (!strcmp (first_word,"ENDMDL") )//fix
        {
            //return end_of_model;
            to_user().warning("Only first model was read.");
            return end_of_file;
        }
        else if ( !strcmp (first_word,"ATOM") 
               || !strcmp (first_word,"HETATM") )
        {
            int current_residue_serial_number;
            read_atom   (line, current_residue_serial_number) ;
            //fix if ( reached_residue_serial_number != current_residue_serial_number )
            //fix     return end_of_residue;
        }
    } 
    return end_of_file;
}


void Pdb_file::
//read_atom   (char rude_line[], int & resudue_serial_number) 
read_atom   (char line[], int & resudue_serial_number) 
{
        char dst[81], buff[81];

        /*//remove wrong symbols
        char *line;
        for (int i=0;  i<100;  ++i)
            if (rude_line[i] == 'A' || rude_line[i] == 'H')
            {
                line = rude_line + i;
                break;
            }//*/

        // Atom serial number 
        substr (line, dst, 7-1, 5);
        int atom_number = atoi(dst);

        // Atom name
        substr (line, dst, 13-1, 4);
        sscanf (dst,"%s",buff);
        Text atom_name = buff;
        int element_N = get_element_N (dst);

        // Residue name
        substr (line, dst, 18-1, 3);
        sscanf (dst,"%s",buff);
        Text residue_name = buff;

        // Chain identifier
        substr (line, dst, 22-1, 1);
        int chain_id = atoi(dst);

        // Residue sequence number
        substr (line, dst, 23-1, 4);
        int residue_number = atoi(dst);
        resudue_serial_number = residue_number;
        
        // xyz
        substr (line, dst, 31-1, 8);
        double x = atof (dst);
        substr (line, dst, 39-1, 8);
        double y = atof (dst);
        substr (line, dst, 47-1, 8);
        double z = atof (dst);//*/

        // Temperature factor to -> Charge
        substr (line, dst, 61-1, 7);
        //string str (dst);
        istringstream charge_str (dst);
        double charge;
        charge_str >> charge;

        // Element symbol, right-justified
        substr( line, dst, 77-1, 2);
        if (Element:: is (dst))
            element_N = Element (dst).number();
        //int element_N_2 = get_element_N (dst);
        //if (element_N_2 != Element ("Unknown").number())
        //    element_N = element_N_2;

//        int element = getElement ( Text(atom_name.c_str())) ;
//        int element = getElement (atom_name) ;
///        Atom * current_atom = create_atom( element );
        Atom & current_atom = new_atom (element_N);
        //MM_atom * current_atom = new MM_atom( element );
//        Atom *    tmp_atom      = create <Atom> (element);
//        MM_atom * current_atom  = static_cast <MM_atom*> (tmp_atom);

        CHECK("", !current_atom.kit().has_pdb_attributes() );
        current_atom.kit().add_pdb_attributes( new Pdb_attributes );
        //model().add_atom( current_atom );
        current_atom.kit().pdb().set_residue_name(residue_name.c_str() );
        current_atom.kit().pdb().set_name(atom_name.c_str() );
        current_atom.kit().pdb().set_residue_number (residue_number);
        current_atom.kit().pdb().set_chain_id(chain_id);
        current_atom.kit().pdb().set_number (atom_number);

        current_atom.set_x(x);
        current_atom.set_y(y);
        current_atom.set_z(z);

        if (charge_str.good())
            current_atom.set_charge (charge);
}

bool   Pdb_file::
find_residue_line(char first_word[], char line[], Input & input) 
{
    do
    {
        char dst[81];
        substr (line, dst, 0, 6);
        sscanf (dst, "%s", first_word);
        //sscanf (line,"%s", first_word);
 
        if ( !strcmp (first_word,"ATOM") 
          || !strcmp (first_word,"HETATM") )
        {
            return true;
        }

    }    
    while( input.next_line(line, 80) );

    return false;       
}  

bool   Pdb_file::
find_chain_line  (char first_word[], char line[], Input & input)
{
    return find_residue_line (first_word, line, input);
} 
 
bool   Pdb_file::
find_model_line  (char first_word[], char line[], Input & input)
{
    return find_residue_line (first_word, line, input);
}  

}//MM



