#include <list>
#include <algorithm>
#include <math.h>
#include "../QueryAtom.h"
#include "../QueryBond.h"
#include "../SmilesParse/SmilesWrite.h"
#include "../SmilesParse/SmartsWrite.h"
//#include "../SmilesParse/SmilesParse.h"
//#include "../Substruct/SubstructMatch.h"
#include "SubstructMatchCustom.h"
#include "MaximumCommonSubgraph.h"

#ifdef VERBOSE_STATISTICS_ON
    ExecStatistics stat;
#endif
namespace RDKit
{
 namespace FMCS
 {
    bool ConsoleOutputEnabled = false;
}}

namespace RDKit
{
 namespace FMCS
 {

    struct LabelDefinition
    {
        unsigned ItemIndex;   // item with this label value
        unsigned Value;
        LabelDefinition() : ItemIndex(-1), Value(-1) {}
        LabelDefinition(unsigned i, unsigned value) : ItemIndex(i), Value(value) {}
    };
        
MaximumCommonSubgraph::MaximumCommonSubgraph(const MCSParameters* params)
{
    Parameters = ( 0 != params ? *params : MCSParameters());
    if (Parameters.ProgressCallback == MCSProgressCallbackTimeout)
        Parameters.ProgressCallbackUserData = &To;
        To = time(0);
}

static
bool molPtr_NumBondLess (const ROMol* l, const ROMol* r)    // need for sorting the source molecules by size
{
    return l->getNumBonds() < r->getNumBonds();
}

void MaximumCommonSubgraph::init()
{
    QueryMolecule = Molecules.front();

    void* userData = Parameters.CompareFunctionsUserData;
    if(Parameters.BondCompareParameters.CompleteRingsOnly || Parameters.BondCompareParameters.RingMatchesRingOnly)
    {
#ifdef FAST_SUBSTRUCT_CACHE
        RingMatchTables.init(QueryMolecule);
        Parameters.CompareFunctionsUserData = &RingMatchTables;
#endif
    }
#ifdef FAST_SUBSTRUCT_CACHE
#ifdef PRECOMPUTED_TABLES_MATCH
        // fill out RingMatchTables for check cache Hash collision by checking match a part of Query to Query
        if(!userData // predefined functor - compute RingMatchTable for all targets
        && (Parameters.BondCompareParameters.CompleteRingsOnly || Parameters.BondCompareParameters.RingMatchesRingOnly))
            RingMatchTables.computeRingMatchTable(QueryMolecule, QueryMolecule, Parameters);

        // fill out match tables
        size_t nq = QueryMolecule->getNumAtoms();
        QueryAtomMatchTable.resize(nq, nq);
        for(size_t aj = 0; aj < nq; aj++)
            for(size_t ai = 0; ai < nq; ai++)
                QueryAtomMatchTable.set(ai, aj, Parameters.AtomTyper(Parameters.AtomCompareParameters, 
                    *QueryMolecule, ai, *QueryMolecule, aj, Parameters.CompareFunctionsUserData));

        nq = QueryMolecule->getNumBonds();
        QueryBondMatchTable.resize(nq, nq);
        for(size_t aj = 0; aj < nq; aj++)
            for(size_t ai = 0; ai < nq; ai++)
                QueryBondMatchTable.set(ai, aj, Parameters.BondTyper(Parameters.BondCompareParameters,
                    *QueryMolecule, ai, *QueryMolecule, aj, Parameters.CompareFunctionsUserData));
#endif // PRECOMPUTED_TABLES_MATCH
        // Compute label values based on current functor and parameters for code Morgan correct computation.
        unsigned currentLabelValue = 1;
        std::vector<LabelDefinition> labels;
        nq = QueryMolecule->getNumAtoms();
        QueryAtomLabels.resize(nq);
        for(size_t ai = 0; ai < nq; ai++)
        {
            if(MCSAtomCompareAny == Parameters.AtomTyper) // predefined functor without atom compare parameters
                QueryAtomLabels[ai] = 1;
            else
            {
                const Atom* atom = QueryMolecule->getAtomWithIdx(ai);
                if(MCSAtomCompareElements == Parameters.AtomTyper) // predefined functor without atom compare parameters
                    QueryAtomLabels[ai] = atom->getAtomicNum()|(Parameters.AtomCompareParameters.MatchValences ? (atom->getTotalValence()>>8):0);
                else if(MCSAtomCompareIsotopes == Parameters.AtomTyper) // predefined functor without atom compare parameters
                    QueryAtomLabels[ai] = atom->getAtomicNum()|(atom->getIsotope()>>8)|(Parameters.AtomCompareParameters.MatchValences ? (atom->getTotalValence()>>16):0);
                else // custom user defined functor
                {
                    QueryAtomLabels[ai] = -1;
                    for(size_t i = 0; i < labels.size(); i++)
                     if(Parameters.AtomTyper(Parameters.AtomCompareParameters, 
                         *QueryMolecule, labels[i].ItemIndex, *QueryMolecule, ai, userData)) // equal itoms
                    {
                        QueryAtomLabels[ai] = labels[i].Value;
                        break;
                    }
                    if(-1 == QueryAtomLabels[ai]) // not found -> create new label
                    {
                        QueryAtomLabels[ai] = ++currentLabelValue;
                        labels.push_back(LabelDefinition(ai, currentLabelValue));
                    }
                }
            }
        }
        labels.clear();
        currentLabelValue = 1;
        nq = QueryMolecule->getNumBonds();
        QueryBondLabels.resize(nq);
        for(size_t aj = 0; aj < nq; aj++)
         {
            const Bond* bond = QueryMolecule->getBondWithIdx(aj);
            unsigned ring = 0;
            if(Parameters.BondCompareParameters.CompleteRingsOnly || Parameters.BondCompareParameters.RingMatchesRingOnly)
            {
                ring = RingMatchTables.isQueryBondInRing(aj) ? 0 : 1;  // is bond in ring
            }
            if(MCSBondCompareAny == Parameters.BondTyper) // predefined functor without atom compare parameters
                QueryBondLabels[aj] = 1 | (ring>>8);
            else if(MCSBondCompareOrderExact == Parameters.BondTyper) // predefined functor without compare parameters
                QueryBondLabels[aj] = (bond->getBondType() + 1) | (ring>>8);
            else if(MCSBondCompareOrder == Parameters.BondTyper)     // predefined functor, ignore Aromatization
            {
                unsigned order = bond->getBondType();
                if(Bond::AROMATIC == order || Bond::ONEANDAHALF == order) // ignore Aromatization
                    order = Bond::SINGLE;
                else if(Bond::TWOANDAHALF == order)
                    order = Bond::DOUBLE;
                else if(Bond::THREEANDAHALF == order)
                    order = Bond::TRIPLE;
                else if(Bond::FOURANDAHALF == order)
                    order = Bond::QUADRUPLE;
                else if(Bond::FIVEANDAHALF == order)
                    order = Bond::QUINTUPLE;
                QueryBondLabels[aj] = (order + 1) | (ring>>8);
            }
            else // custom user defined functor
            {
                QueryBondLabels[aj] = -1;
                for(size_t i = 0; i < labels.size(); i++)
                 if(Parameters.BondTyper(Parameters.BondCompareParameters, 
                                        *QueryMolecule, labels[i].ItemIndex, *QueryMolecule, aj, userData)) // equal bonds + ring ...
                {
                    QueryBondLabels[aj] = labels[i].Value;
                    break;
                }
                if(-1 == QueryAtomLabels[aj]) // not found -> create new label
                {
                    QueryBondLabels[aj] = ++currentLabelValue;
                    labels.push_back( LabelDefinition(aj, currentLabelValue));
                }
            }
        }
#endif
    Targets.resize(Molecules.size()-1);
    size_t i=0;
    for(std::vector<const ROMol*>::iterator it = Molecules.begin()+1; it != Molecules.end(); it++, i++)
    {
        Targets[i].Molecule = *it;
        // build Target Topology ADD ATOMs
        size_t j=0;    // current item
        for(ROMol::ConstAtomIterator a = Targets[i].Molecule->beginAtoms(); a != Targets[i].Molecule->endAtoms(); a++, j++)
        {
            Targets[i].Topology.addAtom((*a)->getIdx());
#ifdef FAST_INCREMENTAL_MATCH
            Targets[i].AtomAdjacency.resize((*it)->getNumAtoms());
            ROMol::OEDGE_ITER beg,end;
            for(boost::tie(beg,end) = (*it)->getAtomBonds(*a); beg!=end; beg++)  // all bonds from the atom
            {
                AtomAdjacency aa;
                aa.Bond = (*(*it))[*beg].get();
                aa.BondIdx = aa.Bond->getIdx();
                aa.ConnectedAtomIdx = aa.Bond->getEndAtomIdx();
                if((*a)->getIdx() == aa.ConnectedAtomIdx)
                    aa.ConnectedAtomIdx = aa.Bond->getBeginAtomIdx();
                Targets[i].AtomAdjacency[j].push_back(aa);
            }
#endif // FAST_INCREMENTAL_MATCH
        }
        // build Target Topology ADD BONDs
        for(ROMol::ConstBondIterator b = Targets[i].Molecule->beginBonds(); b != Targets[i].Molecule->endBonds(); b++)
        {
            const Bond* bond = *b;
            unsigned ii = bond->getBeginAtomIdx();
            unsigned jj = bond->getEndAtomIdx();
            Targets[i].Topology.addBond((*b)->getIdx(), ii, jj);
        }

        // fill out RingMatchTables
        if(!userData // predefined functor - compute RingMatchTable for all targets
        && (Parameters.BondCompareParameters.CompleteRingsOnly || Parameters.BondCompareParameters.RingMatchesRingOnly))
        {
#ifdef FAST_SUBSTRUCT_CACHE
            RingMatchTables.addTargetBondRingsIndeces(Targets[i].Molecule);
            RingMatchTables.computeRingMatchTable(QueryMolecule, Targets[i].Molecule, Parameters);
#endif
        }

#ifdef PRECOMPUTED_TABLES_MATCH
        // fill out match tables
        size_t nq = QueryMolecule->getNumAtoms();
        size_t nt = (*it)->getNumAtoms();
        Targets[i].AtomMatchTable.resize(nq, nt);

        for(size_t aj = 0; aj < nt; aj++)
            for(size_t ai = 0; ai < nq; ai++)
                Targets[i].AtomMatchTable.set(ai, aj, Parameters.AtomTyper(Parameters.AtomCompareParameters, 
                    *QueryMolecule, ai, *Targets[i].Molecule, aj, Parameters.CompareFunctionsUserData));

        nq = QueryMolecule->getNumBonds();
        nt = (*it)->getNumBonds();
        Targets[i].BondMatchTable.resize(nq, nt);
        for(size_t aj = 0; aj < nt; aj++)
            for(size_t ai = 0; ai < nq; ai++)
                Targets[i].BondMatchTable.set(ai, aj, Parameters.BondTyper(Parameters.BondCompareParameters,
                    *QueryMolecule, ai, *Targets[i].Molecule, aj, Parameters.CompareFunctionsUserData));
#endif
    }

    Parameters.CompareFunctionsUserData = userData; // restore
}

void MaximumCommonSubgraph::makeInitialSeeds()
{
    // build a set of initial seeds as "all" single bonds from query molecule
    std::vector<bool> excludedBonds(QueryMolecule->getNumBonds());
    for(size_t i = 0; i < excludedBonds.size(); i++)
        excludedBonds[i] = false;

    Seeds.clear();

    //R1 additional performance OPTIMISATION
    //if(Parameters.BondCompareParameters.CompleteRingsOnly)
    // disable all mismatched rings, and do not generate initial seeds from such disabled bonds
    //  for(  rings .....) for(i......) 
    //   if(mismatched) excludedBonds[i.......] = true;

    for(RWMol::ConstBondIterator bi = QueryMolecule->beginBonds(); bi != QueryMolecule->endBonds(); bi++)
    {
        //R1 additional performance OPTIMISATION
        //if(excludedBonds[(*bi)->getIdx()])
        //    continue;
        Seed seed;
#ifdef FAST_INCREMENTAL_MATCH
        seed.MatchResult.resize(Targets.size());
#endif

#ifdef VERBOSE_STATISTICS_ON
        ++stat.Seed;
        ++stat.InitialSeed;
#endif
        seed.addAtom((*bi)->getBeginAtom());
        seed.addAtom((*bi)->getEndAtom());
        seed.ExcludedBonds = excludedBonds; // all bonds from first to current
        seed.addBond (*bi);
        excludedBonds[(*bi)->getIdx()] = true;

        seed.computeRemainingSize(*QueryMolecule);

        if( ! checkIfMatchAndAppend(seed))
        {
#ifdef VERBOSE_STATISTICS_ON
            ++stat.MismatchedInitialSeed;
#endif
            // optionally remove all such bonds from all targets TOPOLOGY where it exists.
            //..........

            // disable (mark as already processed) mismatched bond in all seeds
            for(SeedSet::iterator si = Seeds.begin(); si != Seeds.end(); si++)
                si->ExcludedBonds[(*bi)->getIdx()] = true;

        }
    }
    Seeds.MaxBonds = 1;
    Seeds.MaxAtoms = 2;
}

void MaximumCommonSubgraph::growSeeds(MolFragment& mcsIdx, MCSResult& res)
{
    unsigned steps = 99999; // steps from last progress callback call. call it immediately in the begining

    // Find MCS -- SDF Seed growing OPTIMISATION (it works in 3 times faster)
    while(!Seeds.empty())
    {
        ++steps;
        stat.TotalSteps++;
        SeedSet::iterator si = Seeds.begin();

        si->grow(*this, *QueryMolecule);

        const Seed& fs = Seeds.front();
        // bigger substructure found
        if((!Parameters.MaximizeBonds && (fs.getNumAtoms() > res.NumAtoms || (fs.getNumAtoms() == res.NumAtoms && fs.getNumBonds() > res.NumBonds)))
         ||( Parameters.MaximizeBonds && (fs.getNumBonds() > res.NumBonds || (fs.getNumBonds() == res.NumBonds && fs.getNumAtoms() > res.NumAtoms)))
         )
        {
            stat.MCSFoundStep = stat.TotalSteps;
            res.NumAtoms    = fs.getNumAtoms();
            res.NumBonds    = fs.getNumBonds();
            mcsIdx.Atoms    = fs.MoleculeFragment.Atoms;
            mcsIdx.Bonds    = fs.MoleculeFragment.Bonds;
            mcsIdx.AtomsIdx = fs.MoleculeFragment.AtomsIdx;
            mcsIdx.BondsIdx = fs.MoleculeFragment.BondsIdx;
//TMP DEBUG
/*    std::cout<<"MCS atoms=(";
    for(unsigned seedAtomIdx = 0; seedAtomIdx < mcs.getNumAtoms(); seedAtomIdx++)
    {
        const Atom* atom = mcs.MoleculeFragment.Atoms[seedAtomIdx];
        std::cout<<atom->getIdx()<<", ";
    }
    std::cout<<") new atoms=[";
    for(unsigned seedAtomIdx = mcs.LastAddedAtomsBeginIdx; seedAtomIdx < mcs.getNumAtoms(); seedAtomIdx++)
    {
        const Atom* atom = mcs.MoleculeFragment.Atoms[seedAtomIdx];
        std::cout<<atom->getIdx()<<", ";
    }
    std::cout<<"] bonds=(";
    for(int i=0; i<mcsIdx.BondsIdx.size(); i++)
        std::cout<<mcsIdx.BondsIdx[i]<<", ";
    std::cout<<") Size="<< mcs.getNumAtoms() <<", "<< mcs.getNumBonds() <<" Remain=" << mcs.RemainingAtoms <<", "<<mcs.RemainingBonds<<" = ";
    std::cout<< MolFragmentToSmiles(*QueryMolecule, *(const std::vector<int>*) &mcs.MoleculeFragment.AtomsIdx, (const std::vector<int>*) &mcs.MoleculeFragment.BondsIdx) <<"\n";  // unsigned
*/
/*
if(0==si->MoleculeFragment.BondsIdx[0])
{
    std::cout<<"\n"
        <<" -----------------"<<(void*)&*si 
        <<"Remaining "<<si->RemainingBonds<< ",  " << si->RemainingAtoms
        <<": LastAddedAtomsBeginIdx = "<<si->LastAddedAtomsBeginIdx<<", LastAddedBondsBeginIdx = "<<si->LastAddedBondsBeginIdx<<"\n";
    for(size_t i = 0; i < si->MoleculeFragment.Bonds.size(); i++)
        std::cout << i << " "<<si->MoleculeFragment.Bonds[i]->getIdx()<<" : "
                        <<" "<<si->MoleculeFragment.Bonds[i]->getBeginAtom()->getIdx()
                        <<" "<<si->MoleculeFragment.Bonds[i]->getEndAtom()->getIdx()
                        <<"\n";
}
*/            
/*
#ifdef xxxWIN32 // TEMP DEBUG !!!!!!!!!!!!!
  if(30 == res.NumAtoms)
  {
    std::cout<<"-------------\n";
      for(int i=0; i<mcs.getNumBonds(); i++)
        std::cout<<i<<" bIdx=" << mcs.MoleculeFragment.BondsIdx[i]//<<"\n";
            <<" =("<< QueryMolecule->getBondWithIdx(mcs.MoleculeFragment.BondsIdx[i])->getBeginAtomIdx()
            <<" , "<< QueryMolecule->getBondWithIdx(mcs.MoleculeFragment.BondsIdx[i])->getEndAtomIdx()<<")\n";
    std::cout<<"-------------\n";
    for(int i=0; i<30; i++)
    {
        const Atom* a = QueryMolecule->getAtomWithIdx(mcs.MoleculeFragment.AtomsIdx[i]);
        std::cout<<i<<" aIdx=" << mcs.MoleculeFragment.AtomsIdx[i]<<" : "<<a->getAtomicNum()<<"\n";
        if(a->getAtomicNum()==7)   // N
        {
            ROMol::OEDGE_ITER beg,end;
            for(boost::tie(beg,end) = QueryMolecule->getAtomBonds(a); beg!=end; beg++)
            {
                const Bond* bond = &*((*QueryMolecule)[*beg]);
                std::cout<<"  bond="<<bond->getIdx()<<" =("<< bond->getBeginAtomIdx()
                    <<" , "<< bond->getEndAtomIdx()<<")\n";
            }
        }
    }
    std::cout<<"-------------\n";
  }
#endif
*/
        }

        if(-1 == si->GrowingStage) //finished
            Seeds.erase(si);

        if(Parameters.ProgressCallback && (steps > 777))// || res.NumAtoms > Stat.NumAtoms))
        {
            steps = 0;
            Stat.NumAtoms = res.NumAtoms;
            Stat.NumBonds = res.NumBonds;
#ifdef VERBOSE_STATISTICS_ON
            Stat.SeedProcessed = stat.Seed;
#endif
            if(!Parameters.ProgressCallback(Stat, Parameters, Parameters.ProgressCallbackUserData))
            {
                res.Canceled = true;
                break;
            }
        }
    }

}

    struct AtomMatch // for each seed atom (matched)
    {
        unsigned QueryAtomIdx;
        unsigned TargetAtomIdx;
        AtomMatch() : QueryAtomIdx(-1), TargetAtomIdx(-1) {}
    };
    typedef std::vector<AtomMatch> AtomMatchSet; 

std::string MaximumCommonSubgraph::generateSMARTS(const MolFragment& mcsIdx)
{
    // match the result MCS with all targets to check if it is exact match or template
    Seed seed; // result MCS
    seed.ExcludedBonds.resize(QueryMolecule->getNumBonds());
    for(size_t i = 0; i < seed.ExcludedBonds.size(); i++)
        seed.ExcludedBonds[i] = false;
    std::vector<AtomMatchSet> atomMatchResult(Targets.size());
    std::vector<unsigned> atomIdxMap(QueryMolecule->getNumAtoms());
    std::vector<std::map<unsigned, const Bond*> > bondMatchSet (mcsIdx.Bonds.size()); //key is unique BondType
    std::vector<std::map<unsigned, const Atom*> > atomMatchSet (mcsIdx.Atoms.size()); //key is unique atomic number
//isotope, mass, charge, Hs    std::vector<std::map<unsigned, const Atom*> > atomMatchSet (mcsIdx.Atoms.size()); //key is unique atomic number

    for(std::vector<const Atom*>::const_iterator atom = mcsIdx.Atoms.begin(); atom != mcsIdx.Atoms.end(); atom++)
    {
        atomIdxMap[(*atom)->getIdx()] = seed.getNumAtoms();
        seed.addAtom((*atom));
    }
    for(std::vector<const Bond*>::const_iterator bond = mcsIdx.Bonds.begin(); bond != mcsIdx.Bonds.end(); bond++)
        seed.addBond((*bond));

    unsigned itarget = 0;
    for(std::vector<Target>::const_iterator tag = Targets.begin(); tag != Targets.end(); tag++, itarget++)
    {
        match_V_t match;    // THERE IS NO Bonds match INFO !!!! 
        bool target_matched = 
#ifdef PRECOMPUTED_TABLES_MATCH
            SubstructMatchCustomTable(tag->Topology, seed.Topology, tag->AtomMatchTable, tag->BondMatchTable, &match);
#else //noticable slowly:
            SubstructMatchCustom(
            tag->Topology.GraphTopology
            , *tag->Molecule
            , seed.GraphTopology
            , *QueryMolecule
            , Parameters.AtomTyper, Parameters.BondTyper
            , Parameters.AtomCompareParameters, Parameters.BondCompareParameters, Parameters.CompareFunctionsUserData
            );
#endif
        atomMatchResult[itarget].resize(seed.getNumAtoms());
        for(match_V_t::const_iterator mit = match.begin(); target_matched && mit != match.end(); mit++)
        {
            unsigned ai = mit->first;  // SeedAtomIdx
            atomMatchResult[itarget][ai].QueryAtomIdx  = seed.Topology[mit->first];
            atomMatchResult[itarget][ai].TargetAtomIdx = tag->Topology[mit->second];
            const Atom* ta = tag->Molecule->getAtomWithIdx(tag->Topology[mit->second]);
            if(ta && ta->getAtomicNum() != seed.MoleculeFragment.Atoms[ai]->getAtomicNum())
                atomMatchSet[ai][ta->getAtomicNum()] = ta; // add
        }
        // AND BUILD BOND MATCH INFO
        unsigned bi=0;
        for(std::vector<const Bond*>::const_iterator bond = mcsIdx.Bonds.begin(); target_matched && bond != mcsIdx.Bonds.end(); bond++, bi++)
        {
            unsigned i = atomIdxMap[(*bond)->getBeginAtomIdx()];
            unsigned j = atomIdxMap[(*bond)->getEndAtomIdx()];
            unsigned ti= atomMatchResult[itarget][i].TargetAtomIdx;
            unsigned tj= atomMatchResult[itarget][j].TargetAtomIdx;
            const Bond* tBond = tag->Molecule->getBondBetweenAtoms(ti, tj);
            if(tBond && (*bond)->getBondType() != tBond->getBondType())
                bondMatchSet[bi] [tBond->getBondType()] = tBond; // add

        }
    }

    // Generate result's SMARTS

    RWMol mol;  // create molecule from MCS for MolToSmarts()
    unsigned ai = 0;  // SeedAtomIdx
    for(std::vector<const Atom*>::const_iterator atom = mcsIdx.Atoms.begin(); atom != mcsIdx.Atoms.end(); atom++, ai++)
    {
        QueryAtom* a = new QueryAtom(*(*atom));  // generate [#6] instead of C or c !
        //for all atomMatchSet[ai] items add atom query to template
        for(std::map<unsigned, const Atom*>::const_iterator am = atomMatchSet[ai].begin(); am != atomMatchSet[ai].end(); am++)
        {
            ATOM_OR_QUERY* a2 = new ATOM_OR_QUERY();
            a2->addChild(QueryAtom::QUERYATOM_QUERY::CHILD_TYPE(makeAtomNumQuery((*atom)->getAtomicNum())));
            a2->addChild(QueryAtom::QUERYATOM_QUERY::CHILD_TYPE(makeAtomNumQuery(am->second->getAtomicNum())));      
            a2->setDescription("AtomOr");
//ATOM_EQUALS_QUERY *makeAtomIsotopeQuery(int what)
//....
            a->setQuery(a2);
        }
        mol.addAtom(a, true);
    }
    unsigned bi = 0;  // Seed Idx
    for(std::vector<const Bond*>::const_iterator bond = mcsIdx.Bonds.begin(); bond != mcsIdx.Bonds.end(); bond++, bi++)
    {
        unsigned i = atomIdxMap[(*bond)->getBeginAtomIdx()];
        unsigned j = atomIdxMap[(*bond)->getEndAtomIdx()];
        QueryBond* b = new QueryBond(*(*bond));
        // add OR template if need
        for(std::map<unsigned, const Bond*>::const_iterator bm = bondMatchSet[bi].begin(); bm != bondMatchSet[bi].end(); bm++)
        {
            BOND_OR_QUERY* b2 = new BOND_OR_QUERY();
            b2->addChild(QueryBond::QUERYBOND_QUERY::CHILD_TYPE(makeBondOrderEqualsQuery((*bond)->getBondType())));
            b2->addChild(QueryBond::QUERYBOND_QUERY::CHILD_TYPE(makeBondOrderEqualsQuery(bm->second->getBondType())));      
            b2->setDescription("BondOr");
            b->setQuery(b2);
        }
        b->setBeginAtomIdx(i);
        b->setEndAtomIdx(j);
        mol.addBond(b, true);
    }
    return MolToSmarts(mol, true);
}

MCSResult MaximumCommonSubgraph::find(const std::vector<ROMOL_SPTR>& src_mols)
{
    clear();
    MCSResult res;

    if(src_mols.size() < 2)
        throw std::runtime_error("FMCS. Invalid argument. mols.size() must be at least 2");
    if (Parameters.Threshold > 1.0)
        throw std::runtime_error("FMCS. Invalid argument. Parameter Threshold must be 1.0 or less.");

    ThresholdCount = (unsigned) ceil((src_mols.size()-1) * Parameters.Threshold);   // minimal required number of matched targets
    if (ThresholdCount < 1) // at least one target
        ThresholdCount = 1;
    if (ThresholdCount > src_mols.size()-1) // max all targets
        ThresholdCount = src_mols.size()-1;
    
    // Selecting CompleteRingsOnly option also enables --ring-matches-ring-only. ring--ring and chain bonds only match chain bonds.
    if(Parameters.BondCompareParameters.CompleteRingsOnly)
        Parameters.BondCompareParameters.RingMatchesRingOnly = true;

    for(std::vector<ROMOL_SPTR>::const_iterator it = src_mols.begin(); it != src_mols.end(); it++)
    {
        Molecules.push_back((*it).get());
        if(!Molecules.back()->getRingInfo()->isInitialized())
            Molecules.back()->getRingInfo()->initialize();  // but do not fill out !!!
    }

    // sort source set of molecules by their 'size' and assume the smallest molecule as a query
    std::stable_sort(Molecules.begin(), Molecules.end(), molPtr_NumBondLess);

    init();
    makeInitialSeeds();

/* // TEMP DEBUG !!!!!!!!!!!!!
  {
    std::cout<<"---Bonds----------\n";
      for(int i=0; i<QueryMolecule->getNumBonds(); i++)
        std::cout<<i<<" "
            <<" =("<< QueryMolecule->getBondWithIdx(i)->getBeginAtomIdx()
            <<" , "<< QueryMolecule->getBondWithIdx(i)->getEndAtomIdx()<<")\n";
    std::cout<<"---Atoms----------\n";
    for(int i=0; i<QueryMolecule->getNumAtoms(); i++)
    {
        const Atom* a = QueryMolecule->getAtomWithIdx(i);
        std::cout<<i<<" : elem="<<a->getAtomicNum()<<" bonds: ";
        ROMol::OEDGE_ITER beg,end;
        for(boost::tie(beg,end) = QueryMolecule->getAtomBonds(a); beg!=end; beg++)
        {
            const Bond* bond = &*((*QueryMolecule)[*beg]);
            std::cout<<bond->getIdx()<<"=("<< bond->getBeginAtomIdx()<<", "<<bond->getEndAtomIdx()<<") ";
        }
        std::cout<<"\n";
    }
    std::cout<<"-------------\n";
  }
// !!!!!!!!!!!!!!!!!!!!!!
*/
    MolFragment mcsIdx; // current MCS
    growSeeds(mcsIdx, res);

    res.SmartsString = generateSMARTS(mcsIdx);

/* // TMP DEBUG 
    std::cout<<"---MCS Bonds------\n";
    for(int i=0; i<mcsIdx.BondsIdx.size(); i++)
        std::cout<<i<<" bIdx=" << mcsIdx.BondsIdx[i]
            <<" =("<< QueryMolecule->getBondWithIdx(mcsIdx.BondsIdx[i])->getBeginAtomIdx()
            <<" , "<< QueryMolecule->getBondWithIdx(mcsIdx.BondsIdx[i])->getEndAtomIdx()<<")\n";
// TMP DEBUG 
//---------------------------
std::cout<<"Query      "<< MolToSmiles(*QueryMolecule) <<"\n";
std::cout<<"MCS Smiles "<< MolFragmentToSmiles(*QueryMolecule, *(const std::vector<int>*) &mcsIdx.AtomsIdx, (const std::vector<int>*) &mcsIdx.BondsIdx) <<"\n";  // unsigned
*/
/*
    unsigned itarget = 0;
    for(std::vector<Target>::const_iterator tag = Targets.begin(); tag != Targets.end(); tag++, itarget++)
    {
        MatchVectType match;
        bool target_matched = SubstructMatch(*tag->Molecule, *QueryMolecule, match);
        std::cout<<"Target "<< itarget+1 << (target_matched ? " matched" : " MISMATCHED") <<"\n";
    }
*/
//---------------------------

#ifdef VERBOSE_STATISTICS_ON
if(ConsoleOutputEnabled)
{
    std::cout << "STATISTICS:\n";
    std::cout << "Total Growing Steps  = " << stat.TotalSteps<<",  MCS found on "<<stat.MCSFoundStep<<" step\n";
    std::cout << "Initial   Seeds      = " << stat.InitialSeed << ",  Mismatched " << stat.MismatchedInitialSeed<<"\n";
    std::cout << "Inspected Seeds      = " << stat.Seed<<"\n";
    std::cout << "Rejected by BestSize = " << stat.RemainingSizeRejected << "\n";
#ifdef EXCLUDE_WRONG_COMPOSITION   
    std::cout << "Rejected by WrongComposition = " << stat.WrongCompositionRejected
                        << " [ " << stat.WrongCompositionDetected << " Detected ]\n";
#endif
    std::cout << "MatchCheck Seeds     = " << stat.SeedCheck      <<"\n";
    std::cout //<< "\n"
              << "     MatchCalls = " << stat.MatchCall      <<"\n"
              << "     MatchFound = " << stat.MatchCallTrue  <<"\n";
#ifdef FAST_INCREMENTAL_MATCH
    std::cout << " fastMatchCalls = " << stat.FastMatchCall <<"\n"
              << " fastMatchFound = " << stat.FastMatchCallTrue <<"\n";
    std::cout << " slowMatchCalls = " << stat.MatchCall     - stat.FastMatchCallTrue <<"\n"
              << " slowMatchFound = " << stat.MatchCallTrue - stat.FastMatchCallTrue <<"\n";
#endif
#ifdef PRECOMPUTED_TABLES_MATCH
    std::cout << "--- USED PreComputed Match TABLES ! ---\n";
#endif
#ifdef VERBOSE_STATISTICS_FASTCALLS_ON
    std::cout << "AtomFunctorCalls = " << stat.AtomFunctorCalls << "\n";
    std::cout << "BondCompareCalls = " << stat.BondCompareCalls << "\n";
#endif
//    std::cout << "\n";
    std::cout << "  DupCacheFound = " << stat.DupCacheFound 
        <<"   "<< stat.DupCacheFoundMatch<<" matched, "
        <<stat.DupCacheFound - stat.DupCacheFoundMatch <<" mismatched\n";
#ifdef FAST_SUBSTRUCT_CACHE
    std::cout << "HashCache size  = " << HashCache.keyssize() << " keys\n";
    std::cout << "HashCache size  = " << HashCache.fullsize() << " entries\n";
    std::cout << "FindHashInCache = " << stat.FindHashInCache << "\n";
    std::cout << "HashFoundInCache= " << stat.HashKeyFoundInCache << "\n";
    std::cout << "ExactMatchCalls = " << stat.ExactMatchCall  <<"\n"
              << "ExactMatchFound = " << stat.ExactMatchCallTrue <<"\n";
#endif
}
#endif
//---------------------
    clear();
    return res;
}

bool MaximumCommonSubgraph::checkIfMatchAndAppend(Seed& seed)
{
#ifdef TRACE_ON
        TRACE() << "CHECK ";    // print out time
        for(std::vector<const Bond*>::const_iterator bi = seed.MoleculeFragment.Bonds.begin(); bi != seed.MoleculeFragment.Bonds.end(); bi++)
            TRACE(0) << (*bi)->getIdx() << " ";
        TRACE(0) << "\n";
#endif
#ifdef VERBOSE_STATISTICS_ON
    ++stat.SeedCheck;
#endif

    bool foundInCache = false;
    bool foundInDupCache = false;

#ifdef DUP_SUBSTRUCT_CACHE
    if(DuplicateCache.find(seed.DupCacheKey, foundInCache))
    {
    // duplicate found. skip match() but store both seeds, because they will grow by different paths !!!
        #ifdef VERBOSE_STATISTICS_ON
            stat.DupCacheFound++;
            stat.DupCacheFoundMatch += foundInCache ? 1 : 0;
        #endif
        if(!foundInCache) // mismatched !!!
            return false;
    }
    foundInDupCache = foundInCache;
#endif
#ifdef FAST_SUBSTRUCT_CACHE
    SubstructureCache::HashKey      cacheKey;
    SubstructureCache::TIndexEntry* cacheEntry = 0;
    bool cacheEntryIsValid = false;
    if(!foundInCache)
    {
    #ifdef VERBOSE_STATISTICS_ON
        ++stat.FindHashInCache;
    #endif
        cacheEntry = HashCache.find(seed, QueryAtomLabels, QueryBondLabels, cacheKey);
        cacheEntryIsValid = true;
        if(cacheEntry)   // possibly found. check for hash collision
        {
            #ifdef VERBOSE_STATISTICS_ON
                ++stat.HashKeyFoundInCache;
            #endif
            // check hash collisions (time +3%):
            for(SubstructureCache::TIndexEntry::const_iterator g = cacheEntry->begin(); !foundInCache && g != cacheEntry->end(); g++)
            {
                if(g->m_vertices.size() != seed.getNumAtoms() || g->m_edges.size() != seed.getNumBonds())
                    continue;
            #ifdef VERBOSE_STATISTICS_ON
                ++stat.ExactMatchCall;
            #endif
            // EXACT MATCH
            #ifdef PRECOMPUTED_TABLES_MATCH
                foundInCache = SubstructMatchCustomTable((*g), seed.Topology, QueryAtomMatchTable, QueryBondMatchTable);
            #else //..................
            #endif
            #ifdef VERBOSE_STATISTICS_ON
                if(foundInCache)
                    ++stat.ExactMatchCallTrue;
            #endif
            }
        }
    }
    #endif
    bool found = foundInCache;

    if(!found) 
    {
        found = match(seed);
    }

    if(found)  // Store new generated seed, if found in cache or in all(- threshold) targets
    {
        Seed& new_seed = Seeds.add(seed);

    #ifdef DUP_SUBSTRUCT_CACHE
        if(!foundInDupCache && seed.getNumBonds() >= 3)  // only seed with a ring can be duplicated - do not store very small seed in cache
            DuplicateCache.add(seed.DupCacheKey, true);
    #endif
    #ifdef FAST_SUBSTRUCT_CACHE
        if(!foundInCache)
            HashCache.add(seed, cacheKey, cacheEntry);
    #endif
    }
    else
    {
    #ifdef DUP_SUBSTRUCT_CACHE
        if(seed.getNumBonds() > 3)
            DuplicateCache.add(seed.DupCacheKey, false);   //opt. cache mismatched duplicates too
    #endif
    }
    return found;  // new matched seed has been actualy added
}

bool MaximumCommonSubgraph::match(Seed& seed)
{
    unsigned max_miss = Targets.size() - ThresholdCount;
    unsigned missing  = 0;
    unsigned passed   = 0;
    unsigned itarget  = 0;

    for(std::vector<Target>::const_iterator tag = Targets.begin(); tag != Targets.end(); tag++, itarget++)
    {
#ifdef VERBOSE_STATISTICS_ON
    ++stat.MatchCall;
#endif
        bool target_matched = false;
#ifdef FAST_INCREMENTAL_MATCH
        if(!seed.MatchResult[itarget].empty())
            target_matched = matchIncrementalFast(seed, itarget);
#endif
        if(!target_matched) // slow full match
        {
            match_V_t match;    // THERE IS NO Bonds match INFO !!!! 
            target_matched = 
#ifdef PRECOMPUTED_TABLES_MATCH
            SubstructMatchCustomTable(tag->Topology, seed.Topology, tag->AtomMatchTable, tag->BondMatchTable, &match);
#else //noticable slowly:
            SubstructMatchCustom(
            tag->Topology
            , *tag->Molecule
            , seed.GraphTopology
            , *QueryMolecule
            , Parameters.AtomTyper, Parameters.BondTyper
            , Parameters.AtomCompareParameters, Parameters.BondCompareParameters, Parameters.CompareFunctionsUserData
            );
#endif
#ifdef FAST_INCREMENTAL_MATCH   // save current match info
            seed.MatchResult[itarget].clear();//resize(target_matched ? match.size() : 0);
            for(match_V_t::const_iterator mit = match.begin(); target_matched && mit != match.end(); mit++)
            {
                BondMatch m;
                m.SeedAtomIdx   = mit->first;
                m.QueryAtomIdx  = seed.GraphTopology[mit->first];
                m.TargetAtomIdx = mit->second; //==target.Topology[it->second]);
                seed.MatchResult[itarget].push_back(m); // add for each matched bond
            }
            //AND RESTORE BOND MATCH INFO !!!!!!!!!!!!!!!!!!!!!
#endif
        }

        if(target_matched)
        {
            if(++passed >= ThresholdCount) // it's enought
                break;
        }
        else // mismatched
        {
            if(++missing > max_miss)
                break;
        }
    }
    if(missing <= max_miss)
    {
        #ifdef VERBOSE_STATISTICS_ON
            ++stat.MatchCallTrue;
        #endif
        return true;
    }
    return false;
}

#ifdef FAST_INCREMENTAL_MATCH
// call it by target, if fail perform full match check
bool MaximumCommonSubgraph::matchIncrementalFast(Seed& seed, unsigned itarget)
{
    // use and update results of previous match stored in the seed
#ifdef VERBOSE_STATISTICS_ON
    ++stat.FastMatchCall;
#endif
    Target& t = Targets[itarget];
    BondMatchSet& match = seed.MatchResult[itarget];
    size_t previousLen = match.size();
    for(unsigned newBondSeedIdx = seed.LastAddedBondsBeginIdx; newBondSeedIdx < seed.getNumBonds(); newBondSeedIdx++)
    {
        bool found=false;
        const Bond* bond = seed.MoleculeFragment.Bonds[newBondSeedIdx];
        unsigned newBondQueryIdx   = seed.MoleculeFragment.BondsIdx[newBondSeedIdx];
        unsigned i = seed.MoleculeFragment.SeedAtomIdxMap[bond->getBeginAtomIdx()];
        if(i >= seed.LastAddedAtomsBeginIdx) // old atom in the seed
            i = seed.MoleculeFragment.SeedAtomIdxMap[bond->getEndAtomIdx()];
        unsigned newBondSourceAtomSeedIdx   = i; // seed's index of atom from which new bond was added
        const BondMatch& ma = seed.MatchResult[itarget][i];
        unsigned newBondSourceAtomTargetIdx = ma.TargetAtomIdx; // corresponding atom in the target
unsigned newBondSourceAtomQueryIdx  = ma.QueryAtomIdx;
//const Atom* newBondSourceAtomTarget = t.Molecule->getAtomWithIdx(newBondSourceAtomTargetIdx);
        // check all unvisited bonds from newBondSourceAtomTargetIdx
        const AtomAdjacencyList& taAdj = t.AtomAdjacency[newBondSourceAtomTargetIdx];
        for(size_t itaAdji = 0; itaAdji < taAdj.size(); itaAdji++)
        {
            unsigned tbi = taAdj[itaAdji].BondIdx;
            if(t.BondMatchTable.at(newBondQueryIdx, tbi))
            {
/*
                // and check ending atom of the bond too 
                //..........
                bool visited = false;
                //match.find
                //MULTIPLE various bonds for one atom !!!
                for(size_t im = 0; !visited && im < match.size(); im++)
                    if(seed.MatchResult[itarget][im]. == )
                        visited = true;
                if(visited) // visited == targer bons already presents in the match of this seed
                    continue;
                // append to match result
                seed.MatchResult[itarget][]

                found = true;
                break;
*/
            }
        }
    }
/*/-----------------------
//Seed:
unsigned LastAddedBondBeginIdx;
//history for each tag
seedBond, tagBond,
both (seedAtom, tagAtom) of

        Graph q = seed.GraphTopology;
        for(unsigned ib = seed.LastAddedBondBeginIdx; ib < seed.getNumBonds(); ib++)   // all added bonds

            ib
            endAtom // of ib
            unsigned tagBeginAtomIdx = ...; // already verified. matched with seed.
            for() all outgoing bonds except already matched to seed
                check match of bond to ib and endAtom to ia if ia is not already verified (ring back to seed)
        //        for(unsigned srcAtomIdx = seed.LastAddedAtomsBeginIdx; srcAtomIdx < seed.getNumAtoms(); srcAtomIdx++)   // all added atoms

//-----------------------*/
    if(false // match_ok
        )
    {
    #ifdef VERBOSE_STATISTICS_ON
        ++stat.FastMatchCallTrue;
    #endif
        return true;
    }
    match.resize(previousLen);  // clean up - remove new items only
    return false;
}
#endif // FAST_INCREMENTAL_MATCH

}}   // namespace RDKit