// FeatureRangeAlphabet.cc
// Thu Aug 13 16:56:28 PDT 2009 Kevin Karplus

// This file contains the implementation of the FeatureRangeAlphabet class;
// see FeatureRangeAlphabet.h for a description of this class.
//
//	Note: the actual assignment of letters based on Conformations is
//	deliberately not part of this class, but is handled by
//	    friend Base FeatureRangeClassify(const Conformation* conf,
//			    const ChainsResiduesAndAtoms *cra,
//			    int res_idx,
//			    const FeatureRangeAlphabet &al);



#include <sstream>

#include <math.h>	// for isfinite() and floor()

#include "EqualStrings/EqualStrings.h"
#include "FeatureRangeAlphabet.h"

using namespace std;

// the command table
NameToPtr *FeatureRangeAlphabet::CommandTable = NULL;

// required stuff for things derived from NamedClass
static NamedClass *create_FeatureRangeAlphabet(void)
{   
    return new FeatureRangeAlphabet;
}

IdObject FeatureRangeAlphabet::ID("FeatureRangeAlphabet",
				 create_FeatureRangeAlphabet,
				 FeatureRangeAlphabet::init_is_a);


// initialize the command table with all the commands from the base class
// plus some new ones
void FeatureRangeAlphabet::init_command_table(void)
{   
    assert(! CommandTable);
    CommandTable = new NameToPtr(30);
    CommandTable->ignore_case();

    // commands from the base class
    CommandTable->AddName(new AlphabetInputCommand("Comment",
						   Alphabet::ReadComment));
    CommandTable->AddName(new AlphabetInputCommand("ClassName",
						   Alphabet::VerifyClassName));
    CommandTable->AddName(new AlphabetInputCommand("EndClassName",
						   Alphabet::VerifyClassName));
    CommandTable->AddName(new AlphabetInputCommand("Name",
						   Alphabet::ReadName));
    CommandTable->AddName(new AlphabetInputCommand("NormalChars",
						   Alphabet::ReadNormalChars));
    CommandTable->AddName(new AlphabetInputCommand("Alias",
						   Alphabet::ReadAlias));
    CommandTable->AddName(new AlphabetInputCommand("Wildcard",
						   Alphabet::ReadWildcard));
    CommandTable->AddName(new AlphabetInputCommand("AllMatch",
						   Alphabet::ReadAllMatch));
    CommandTable->AddName(new AlphabetInputCommand("CharName",
						   Alphabet::ReadCharName));
    CommandTable->AddName(new AlphabetInputCommand("IsNucleic",
						   Alphabet::ReadIsNucleic));
    CommandTable->AddName(new AlphabetInputCommand("IsRNucleic",
						   Alphabet::ReadIsRNucleic));
    CommandTable->AddName(new AlphabetInputCommand("IsAmino",
						   Alphabet::ReadIsAmino));


    // new commands for FeatureRangeAlphabet
    CommandTable->AddName(new AlphabetInputCommand("Feature",
				       FeatureRangeAlphabet::ReadFeature));
    CommandTable->AddName(new AlphabetInputCommand("Modulus",
				       FeatureRangeAlphabet::ReadModulus));
    CommandTable->AddName(new AlphabetInputCommand("Offset",
				       FeatureRangeAlphabet::ReadOffset));
    CommandTable->AddName(new AlphabetInputCommand("Range",
				       FeatureRangeAlphabet::ReadRange));

}


// required command for reading the alphabet
int FeatureRangeAlphabet::read_knowing_type(istream &in)
{
    if (!CommandTable)
        init_command_table();

    char word[500];
    int non_end_class = 1;

    while (in.good() && non_end_class)
    {   
	get_word(in,word, '=');	
	AlphabetInputCommand *comm = dynamic_cast<AlphabetInputCommand *>
			(CommandTable->FindOldName(word, ZeroIfNew));

        if (comm != 0)
        {   non_end_class = comm->execute(in, this);
	}
	else
        {   cerr << "# error: unknown FeatureRangeAlphabet command: "
                 << word << endl;
        }
    }

    if (non_end_class)
    {
        cerr << "Warning: FeatureRangeAlphabet "
		<< name() << " failed command: "
	     << word << " EndClassName command\n";
    }

    if (! ranges_ok())
    {   cerr << "ERROR: FeatureRangeAlphabet "
    		<<name() << " has invalid ranges ";
	print_modulus(cerr);
	print_ranges(cerr);
	return 0;
    }
    
    if (name() == NULL)
    {
        cerr << "Warning: FeatureRangeAlphabet::read_knowing_type(): alphabet "
	     << "has no name, not registered in lookup table.\n";
    }

    return !non_end_class;
}

// read the Feature name
int FeatureRangeAlphabet::ReadFeature(istream &in,
				       Alphabet *change,
				       AlphabetInputCommand *self) 
{   assert(change && self);

    // make an input stream for the current line
    string line;
    getline(in, line);
    istringstream linestream(line);

    // make sure this command is only called for FeatureRangeAlphabet
    FeatureRangeAlphabet *ara = dynamic_cast<FeatureRangeAlphabet*>(change);
    assert(ara);

    linestream >> ara->FeatName >> ws;
    if ( ara->FeatName.length()==0 
	|| linestream.fail() || ! linestream.eof())
    {   
	cerr << "ERROR: bad value or format for " << self->name()
		<< " in FeatureRangeAlphabet " << ara->name()     << endl;
        return 0;
    }

    return 1;
}

// read the Modulus
int FeatureRangeAlphabet::ReadModulus(istream &in,
					Alphabet *change,
					AlphabetInputCommand *self) 
{   
    assert(change && self);

    // make sure this is command is only called for FeatureRangeAlphabet
    FeatureRangeAlphabet *ara = dynamic_cast<FeatureRangeAlphabet*>(change);
    assert(ara);

    in >> ara->Modulus;    
    return 1;
}

// read the Offset
int FeatureRangeAlphabet::ReadOffset(istream &in,
					Alphabet *change,
					AlphabetInputCommand *self) 
{   
    assert(change && self);

    // make sure this is command is only called for FeatureRangeAlphabet
    FeatureRangeAlphabet *ara = dynamic_cast<FeatureRangeAlphabet*>(change);
    assert(ara);

    in >> ara->Offset;    
    return 1;
}

// read the Ranges for the various letters
int FeatureRangeAlphabet::ReadRange(istream &in,
					Alphabet *change,
					AlphabetInputCommand *self) 
{   
    assert(change && self);

    // make an input stream for the current line
    string line;
    getline(in, line);
    istringstream linestream(line);

    // make sure this is command is only called for FeatureRangeAlphabet
    FeatureRangeAlphabet *ara = dynamic_cast<FeatureRangeAlphabet*>(change);
    assert(ara);

    // make sure the match all wildcard (and indirectly normal chars) are read
    if ((ara->first_all_match_wildcard()).is_null())
    {
	cerr << "error: wildcard not present before " << self->name()
		<< " in FeatureRangeAlphabet " << ara->name()     << endl;
	return 0;
    }

    // read the letter to which this commands pertains
    char inputbase;
    if (! (linestream >> inputbase))
    {   cerr << "ERROR: format error"
		<< " in FeatureRangeAlphabet " << ara->name()
		<< " for "  << self->name()     << endl;
	
	return 0;
    }
    
    int alpha_idx = ara->index_from_char(inputbase);
    if (alpha_idx < 0 || alpha_idx > ara->norm_length())
    {   cerr << "ERROR: character " << inputbase << " is not a normal character"
		<< " in FeatureRangeAlphabet " << ara->name()
		<< " for "  << self->name()     << endl;
	return 0;
    }

    unsigned int subscript = static_cast<unsigned int>(alpha_idx);
    
    // Make sure vectors are big enough
    if (subscript>= ara->LowEnd.size() || subscript>=ara->HighEnd.size())
    {   ara->LowEnd.resize(ara->norm_length());
    	ara->HighEnd.resize(ara->norm_length());
    }
    
    string word;
    
    linestream >> word >> ws;
    if (EqualStrings(word.c_str(), "-inf", 1))
    {	ara->LowEnd[subscript] = 0-INFINITY;
    }
    else
    {	stringstream tonum(word);
        tonum>>ara->LowEnd[subscript];
    }

    linestream >> word >> ws;
    if (EqualStrings(word.c_str(), "+inf", 1))
    {	ara->HighEnd[subscript] = INFINITY;
    }
    else
    {	stringstream tonum(word);
        tonum>>ara->HighEnd[subscript];
    }

    return 1;
}



// print the FeatureName
void FeatureRangeAlphabet::print_feature(std::ostream &out) const
{   
    out << "Feature = " << FeatName << endl;
}

// print the Modulus
void FeatureRangeAlphabet::print_modulus(std::ostream &out) const
{   
    out << "Modulus = " << Modulus << endl;
}

// print the Offset
void FeatureRangeAlphabet::print_offset(std::ostream &out) const
{   
    out << "Offset = " << Offset << endl;
}

// print the ranges
void FeatureRangeAlphabet::print_ranges(std::ostream &out) const
{   
    for (unsigned int i=0; i<LowEnd.size(); i++)
    {    out << "Range = " << to_char(unindex(i)) << " ";
         if (isfinite(LowEnd[i]))
	 {    out << LowEnd[i];
	 }
	 else 
	 {    out << "-inf";
	 }
	 out << " ";
         if (isfinite(HighEnd[i]))
	 {    out << HighEnd[i];
	 }
	 else 
	 {    out << "+inf";
	 }
	 out << endl;
    }
}


// required function to write angle vector alphabets
void FeatureRangeAlphabet::write_knowing_type(ostream &out) const
{
    out << "\nName = " << name() << "\n";

    Alphabet::print_is_values(out);
    Alphabet::print_normal_chars(out);
    Alphabet::print_wildcards(out);
    Alphabet::print_alias_chars(out);
    Alphabet::print_char_names(out);

    out << endl;
    FeatureRangeAlphabet::print_feature(out);
    FeatureRangeAlphabet::print_offset(out);
    FeatureRangeAlphabet::print_modulus(out);
    
    out << endl;
    FeatureRangeAlphabet::print_ranges(out);
    
    out << endl;
}

bool FeatureRangeAlphabet::ranges_ok(void) const
{
    
    // check that sizes of LowEnd and HighEnd are correct
    unsigned int num_norm = norm_length();
    if (num_norm==0) return false;	// must have some characters
    if (num_norm != LowEnd.size()) return false;
    if (num_norm != HighEnd.size()) return false;
    
    // Note: we constrain alphabets so that 
    //	LowEnd[i+1] = HighEnd[i]  (mod modulus) 
    // to  allow for sanity checks (like no missing classification).

    if (Modulus < 0) return false;	// Modulus must be >=0
    
    if (Modulus==0)
    {	// If Modulus==0, we need to have LowEnd[i]<HighEnd[i];
    	// LowEnd[0] and HighEnd[num_norm-1] should be effectively -inf and inf,
	//	but this is not checked, as for some ResidueFeatures finite
	//	values may make good sense.
	for (unsigned int i=0; i< num_norm; i++)
	{   if (LowEnd[i] >= HighEnd[i]) return false;
	}
	
	// Make sure that range is covered without holes
	for (unsigned int i=1; i< num_norm; i++)
	{    if (LowEnd[i] != HighEnd[i-1]) return false;
	}
	
	return true;	// no more checks
    }
    
    // With positive Modulus, we need to check
    // that ranges match modulo the Modulus.
    // x mod y is x-y*floor(x/y).

    // Make sure that range is covered without holes (including wraparound)
    for (unsigned int i=0; i< num_norm; i++)
    {   double diff= (HighEnd[i] -LowEnd[(i+1)%num_norm]);
        if (diff -Modulus*floor(diff/Modulus) > 1.e-07*Modulus) return false;
    }
	
    if (num_norm <=1) return true;
    
    // Make sure that adjacent ranges don't overlap
    for (unsigned int i=0; i< num_norm; i++)
    {   if (in_range(LowEnd[(i+2) % num_norm], i)) return false;
    }
    
    return true;

}

// CHANGE LOG:
//  Fri Aug 14 17:03:19 PDT 2009 Kevin Karplus
//	Changed sanity test to use new in_range() function