molecule.h

/****************************************************************************************
                                          molecule.h 
                                        --------------
    copyright            : (C) 2006 Jean-Luc Perret - Pierre Mahé
    email                : jean-luc.perret@unine.ch - pierre.mahe@ensmp.fr
 ***************************************************************************************/

/****************************************************************************************
 *                                                                                      *
 *      This program is free software; you can redistribute it and/or                   *
 *      modify it under the terms of the GNU Lesser General Public                      *
 *      License as published by the Free Software Foundation; either                    *
 *      version 2.1 of the License, or (at your option) any later version.              *
 *                                                                                      *
 *      This program is distributed in the hope that it will be useful,                 *
 *      but WITHOUT ANY WARRANTY; without even the implied warranty of                  *
 *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU               *
 *      Lesser General Public License for more details.                                 *
 *                                                                                      *
 *      You should have received a copy of the GNU Lesser General Public                *
 *      License along with this library; if not, write to the Free Software             *
 *      Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA      *
 *                                                                                      *
 ****************************************************************************************/



#ifndef MOLECULE_H
#define MOLECULE_H

#include <vector>
#include <map>
#include <sstream>
#include <locale>
#include <string>
#include <algorithm>


#include <datacontainer.h>
#include <atom.h>
#include <elements.h>
#include <constant.h>

//#include <kashimathreadstruct.h>

extern Elements elements;


class Molecule : public DataContainer  {

public:



        Molecule();

        virtual ~Molecule();

        Molecule( Molecule& aMolecule , bool bool_resetMorganIndex = true );

        Molecule(
                                Molecule& molecule1,
                                Molecule& molecule2,
                                double (*pt2AtomKernel)( Atom*, Atom* ),
                                double (*pt2BondKernel)( Bond*, Bond* )
                                );

        Molecule(
                        Molecule& m1,
                        Molecule& m2,
                        double (*pt2AtomKernel)( Atom*, Atom*),
                        double (*pt2BondKernel)( float, float, float ),
                        float edgeKernelParameter
                        );


        Molecule& operator=( const Molecule& aMolecule );

        virtual Atom* addAtom( string aSymbol, bool resetSSSR = true ) throw( CError );

        Atom* addAtom( Atom*, bool resetSSSR = true, bool resetMorganIndex = true ) throw( CError);


        Bond* linkAtoms( Atom* aSource, Atom* aTarget, int aBondLabel, int aBondStereo = 0,  int aBondNotUsed = 0, int aBondTopology = 0, int aBondReactionCenter = 0, bool resetSSSR = true );

        Bond* linkAtoms( string aSource, string aTarget, int aBondLabel, int aBondStereo = 0,  int aBondNotUsed = 0, int aBondTopology = 0, int aBondReactionCenter = 0, bool resetSSSR = true ) throw( CError );

        Bond* linkAtoms( int firstAtom, int secondAtom, int aBondLabel, int aBondStereo = 0,  int aBondNotUsed = 0, int aBondTopology = 0, int aBondReactionCenter = 0, bool resetSSSR = true ) throw( CError );


        void readMOL( string aFilename, bool genericAtomType = false ) throw(CError);

        void readMOLOld( string aFilename ) throw(CError);

        void erase();
        
        void eraseHiddenAtoms();

        void eraseRings();

        void eraseAdjacency();
        
        void eraseWalks();






        int getId() const { return( id ); }

        string getIdString();

        string getName(){ return( getStringDescriptor( "name" )->getValue() ); }

        void setName( string aName ) { setStringDescriptor( "name", aName, "", "", true, true ); }

        void setActivity( float aNumber ){ activity = aNumber; flagActivity = true; }

        void setActivity( string aNumber ){ activity = atof( aNumber.c_str() ); flagActivity = true; };

        float getActivity( bool silentMode = false ) throw( CError );
        
        float hasActivity(){ return( flagActivity ); }

        void unsetActivity(){ flagActivity = false; }


        vector<Atom*>& getAtoms(){ return( atoms ); }

        Atom* getAtom(int anId) throw( CError );

        Atom* getAtomByIndex(int ind);

        bool atomExists(Atom* anAtom);

        int numAtoms(){ return( atoms.size() ); }

        int numHiddenAtoms(){ return( hiddenAtoms.size() ); }


        int numBonds();
        int numHiddenBonds();

        long bondSum();

        void unsetBondFlags();
        void unsetBondFlagsOriginal();

        bool isChiral(){ return( chiral ); }

        void select(){ selectedFlag = true; }
        void unSelect(){ selectedFlag = false; }
        bool isSelected(){ return( selectedFlag ); }

        void setSortDescriptor( string aName, int aType );


        int getSortDescriptorType(){ return( sortDescriptorType ); }

        string getSortDescriptorName(){ return( sortDescriptorName ); }

        Ring* getRingWithID( int anID, bool createIfMissing) throw( CError );

        bool hasRing( Ring* newRing, bool detectingRing = false ) throw( CError );

        bool hasRing() throw( CError );

        int numRings() throw( CError ){
                if( !hasSSSRDetected() ){
                        //detectSSSR();
                        CError e( SSSRNOTDETECTED, "Smallest Set of Smallest Rings was not detected before calling Molecule::numRings()" );
                        //e.describe();
                        throw(e);
                }
                return( sssr.size() );
        }


        bool hasSSSRDetected(){ return( flagHasSSSRDetected ); }

        float getMW( bool silentError = false ) throw( CError );

        int getNumAtoms(){
                return( atoms.size() );
        }

        string getLocation(){ return( location ); }

        int getOriginalFormat(){ return( originalFormat ); }
        
        void setOriginalFormat( int a ){ originalFormat = a; }



        





        int hideAtomsByIntDescriptor( string aDescriptorName, int aValue, bool refreshBonds = true );

        void hideAtom( vector< Atom* >::iterator anAtomI );

        void hideAtomAndToFromBonds( vector< Atom* >::iterator anAtomI );

        void hideAtomAndToFromBonds( Atom* anAtom );

        int hideHydrogens();

        bool isHiddenAtom( Atom* anAtom );

        int restoreHiddenAtoms( bool flagRestoreBonds = true );

        int restoreHiddenBonds();

        int refreshBonds();

        void eraseAtom( Atom* anAtom ) throw( CError );

        void deleteBonds();

        void deleteHiddenAtoms();

        void setMorganLabels( int anOrder );

        void setPerretLabels();

        int setUniqueMorganIndices();

        void resetMorganIndex();

        int getNumberOfDistinctMorganIndices( int anOrder );

        int getMaxMorganIteration();


        int getNumCarbonsOfComponent( string aDescriptorName, int aValue );

        int getNumNitrogensOfComponent( string aDescriptorName, int aValue );

        int numAtomsNonCSkeleton();     

        float atomicDistance(Atom* atom1, Atom* atom2);

        void compute();

        void binClassifyFromDescriptor( string descriptorName, float value, bool smallerOrEqual = true );

        
        void atomsLabelsListing( vector<string>* );

        void atomsSymbolsListing( vector<string>* );

        void bondsListing( vector<int>* );

        void noTottersTransform();

        void threeDtransform(int nBins, double distMin, double distMax);

        void readPartialCharges(string charges);
        
        void setMorganChargesLabels(double threshold);

        int hideSalts( stringstream* out );

        int markFragments();

        void unmarkFragments();

        void hideAllFragmentsBut( int aFragmentNumber );

        void writeFragments( ofstream* outStream );

        int DFS( Atom* startAtom, string intDescriptorName, int markValue );

        int detectSSSR();

        void resetSSSR(){ flagHasSSSRDetected = false; }

        void setHasSSSR(){ flagHasSSSRDetected = true; }

        Bond* checkEdges( Atom* );









        vector<Atom*>::iterator beginAtom(){ return( atoms.begin() ); }
        vector<Atom*>::iterator endAtom(){ return( atoms.end() ); }

        map<Atom*, Bond*>::iterator beginBond( Atom* anAtom ){ return( anAtom->beginBond() ); }
        map<Atom*, Bond*>::iterator endBond( Atom* anAtom ){ return( anAtom->endBond() ); }

        map<Atom*, Bond*>::iterator beginBond( int anId ){ return( getAtom( anId )->beginBond() ); }
        map<Atom*, Bond*>::iterator endBond( int anId ){ return( getAtom( anId )->endBond() ); }

        map<int, int>::iterator beginComponentSizes(){
                return( componentSizes.begin() );
        }

        map<int, int>::iterator endComponentSizes(){
                return( componentSizes.end() );
        }






        void setKashimaKernelProb( double aPq, bool skipSkeleton = false );

        double sumPT();

        double sumProbabilities();
        
        double sumProbabilitiesFast();

        double sumPQPS();

        double sumPQPSFast();


        void raisePowerFast();

        double computeKernel(
                                Molecule* anotherMolecule,
                                double (*pt2GraphKernel)(
                                        Molecule* mol1,
                                        Molecule* mol2,
                                        double (*pt2AtomKernel) ( Atom*, Atom* ),
                                        double (*pt2BondKernel)( Bond*, Bond* ),
                                        int parameter1, int parameter2
                                ),
                                double (*pt2AtomKernel)( Atom*, Atom* ),
                                double (*pt2BondKernel)( Bond*, Bond* ),
                                int parameter1, int parameter2 = 1 );


        double getSelfKernel( double(*pt2GraphKernel)( Molecule* mol1,
                Molecule* mol2, double(*pt2AtomKernel)(Atom*, Atom*),
                double(*pt2BondKernel)( Bond*, Bond* ), int, int ),
                double(*pt2AtomKernel)( Atom*, Atom* ),
                double(*pt2BondKernel)( Bond*, Bond* ) ,
                int parameter1, int parameter2 = 1
                ){
                        #ifdef DEBUG
                                cout << "Molecule::getSelfKernel(...) " << endl;
                        #endif

                        if(selfKernelCalculated == false){
                                #ifdef DEBUG
                                        cout << "  sk not computed, computing" << endl;
                                #endif
                                calculateSelfKernel( pt2GraphKernel, pt2AtomKernel, pt2BondKernel, parameter1, parameter2 );
                        }
                        #ifdef DEBUG
                                cout << "  returning " << selfKernel << endl;
                        #endif
                        return( selfKernel );
                }


        double calculateSelfKernel
        (
                                double(*pt2GraphKernel)
                                (
                                                Molecule* mol1,
                                                Molecule* mol2,
                                                double(*pt2AtomKernel)(Atom*, Atom*),
                                                double(*pt2BondKernel)(Bond*, Bond*),
                                                int, int
                                ),
                                double(*pt2AtomKernel)( Atom*, Atom* ),
                                double(*pt2BondKernel)( Bond*, Bond* ),
                                int paramter1, int parameter2
        );


        double getSelfKernel() throw( CError );

        void setSelfKernel(double value);

        void addToSelfKernel( double );
        
        void substractToSelfKernel( double );

        void resetSelfKernel(){ selfKernelCalculated = false; }






        void setAdjacency(int i, int j, double value);

        double getAdjacency(int i, int j);

        void setWalks(int i, int j, double value);

        double getWalks(int i, int j);

        void raisePowerAdjacency();
        
        double traceWalks();

        double traceDiagWalks();







        string toString();

        string toStringShort();

        string toStringLong();


        void describe();

        void describeShort();

        void describeLong();

        void describeEachAtom();

        void writeMOL( string aFileName );

        void writeDOT( string aFilename, bool perretLabels = false );

        void writeSD( string aFileName );



// ******************************* //
// **** DEPRECATED FUNCTIONS **** //
// ****************************** //


        //void saveAllBonds();





        //virtual float moleculeKernel( Molecule* anotherMolecule, float (*pt2AtomKernel)( Atom*, Atom* ), float (*pt2BondKernel)( Bond*, Bond* ), int convergenceCondition = 1000 );


        //virtual float powerKernelUntilN( Molecule* anotherMolecule, float (*pt2AtomKernel)( Atom*, Atom* ), int maxPower = 4 );

        //virtual float powerKernelConverge( Molecule* anotherMolecule, float (*pt2AtomKernel)( Atom*, Atom* ), int converg = 1000 );

        //virtual float powerKernelOrderN( Molecule* anotherMolecule, float (*pt2AtomKernel)( Atom*, Atom* ), int length );




        //  void* kashimaKernelThread(void* arg);




        //void raisePower();





        //void exportFragments( string sdFileName, int minAtoms = 1 );



protected:

        void moleculeChanged( bool resetSSSR = true, bool resetMorganIndex = true );


        vector<Atom*> atoms;

        vector<Atom*> hiddenAtoms;


        vector<Ring*> sssr;

        bool flagHasSSSRDetected;

        static int counter;

        int id;

        bool selectedFlag;

        double selfKernel;

        bool selfKernelCalculated;

        bool chiral;

        bool flagActivity;

        int sortDescriptorType;

        string sortDescriptorName;

        int maxMorganIteration;

        

private: // Private methods

        //void hideAllBonds();

        Bond* linkAtomsNoReturn( Atom* aSource, Atom* aTarget, int aBondLabel, int aBondStereo = 0,  int aBondNotUsed = 0, int aBondTopology = 0, int aBondReactionCenter = 0 );

        void linkAtomsNoReturn( int firstAtom, int secondAtom, int aBondLabel ) throw( CError );


        map<Atom*, map<Atom*, double>* >* fastPT;
        map<Atom*, map<Atom*, double>* >* fastPTNext;
        map<Atom*, map<Atom*, double>* >* fastPTSave;   
        map<Atom*, double> fastPQ;
        map<Atom*, double> fastPS;


        map< int, int > componentSizes;

        float activity;

        string location;        
        
        int originalFormat;

        vector< vector<double> >* adjacency;
        
        vector< vector<double> >* walks;

};

#endif

---