mapSearch.cpp

Go to the documentation of this file.

//============================================ ProSHADE
#include "../proshade/ProSHADE.h"
#include "../proshade/ProSHADE_internal.h"

//============================================ Main
int main ( int argc, char **argv )
{
    //======================================== Initialise variables
    ProSHADE::ProSHADE_settings* setUp        = new ProSHADE::ProSHADE_settings ( );
    setUp->verbose                            = -1;
    setUp->extraSpace                         = 0.0;
    setUp->volumeTolerance                    = 10.0;
    setUp->useCOM                             = false;
    
    std::string mapName                       = "../xx_searchMap_data/02/emd_7780.map";
    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainA.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainB.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainC.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainD.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainE.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainG.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainH.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainI.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainJ.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainK.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainL.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainM.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainN.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainO.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainP.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainQ.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainS.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainT.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainU.pdb" );
//    setUp->structFiles.emplace_back           ( "../xx_searchMap_data/02/ChainV.pdb" );
    
    std::vector<std::array<double,2> > strRadii;
    std::vector<double> shSpacing             ( setUp->structFiles.size(), setUp->shellSpacing );
    std::vector<unsigned int> bandVec         ( setUp->structFiles.size(), setUp->bandwidth );
    std::vector<unsigned int> thetaVec        ( setUp->structFiles.size(), setUp->theta );
    std::vector<unsigned int> phiVec          ( setUp->structFiles.size(), setUp->phi );
    std::vector<unsigned int> glIntegVec      ( setUp->structFiles.size(), setUp->glIntegOrder );
    std::vector<double> extraVec              ( setUp->structFiles.size(), setUp->extraSpace );
    std::array<double,2> hlpArr;
    
    double maxRange                           = 0.0;
    double secRange                           = 0.0;
    
    //======================================== For each input structure
    for ( unsigned int strIt = 0; strIt < static_cast<unsigned int> ( setUp->structFiles.size ( ) ); strIt++ )
    {
        //==================================== Initialise structure reading
        ProSHADE_internal::ProSHADE_data *str = new ProSHADE_internal::ProSHADE_data ( );
        
        //==================================== Read in the structures
        if ( ProSHADE_internal::checkFileType ( setUp->structFiles.at ( strIt ) ) == 1 )
        {
            str->getDensityMapFromPDB         ( setUp->structFiles.at    ( strIt ),
                                               &shSpacing.at             ( strIt ),
                                                setUp->mapResolution,
                                               &bandVec.at               ( strIt ),
                                               &thetaVec.at              ( strIt ),
                                               &phiVec.at                ( strIt ),
                                               &glIntegVec.at            ( strIt ),
                                               &extraVec.at              ( strIt ),
                                                setUp->mapResDefault,
                                                setUp,
                                                setUp->bFactorValue,
                                                false,
                                                setUp->overlayDefaults );
        }
        else
        {
            str->getDensityMapFromMAP         ( setUp->structFiles.at    ( strIt ),
                                               &shSpacing.at             ( strIt ),
                                                setUp->mapResolution,
                                               &bandVec.at               ( strIt ),
                                               &thetaVec.at              ( strIt ),
                                               &phiVec.at                ( strIt ),
                                               &glIntegVec.at            ( strIt ),
                                               &extraVec.at              ( strIt ),
                                                setUp->mapResDefault,
                                                setUp->rotChangeDefault,
                                                setUp,
                                                setUp->overlayDefaults );
        }
        
        //==================================== Find input structure radius
        maxRange                              = std::max ( str->getMapXRange ( ), std::max ( str->getMapYRange ( ), str->getMapZRange ( ) ) );
        if ( maxRange == str->getMapXRange ( ) )
        {
            secRange                          = std::max ( str->getMapYRange ( ), str->getMapZRange ( ) );
        }
        else if ( maxRange == str->getMapYRange ( ) )
        {
            secRange                          = std::max ( str->getMapXRange ( ), str->getMapZRange ( ) );
        }
        else
        {
            secRange                          = std::max ( str->getMapXRange ( ), str->getMapYRange ( ) );
        }
        
        hlpArr[0]                             = strIt;
        hlpArr[1]                             = sqrt ( pow ( maxRange, 2.0 ) + pow ( secRange, 2.0 ) );
        strRadii.emplace_back                 ( hlpArr );

        //==================================== Free memory
        delete str;
    }
    
    //======================================== Sort by radii
    std::sort ( strRadii.begin (), strRadii.end(), [](const std::array<double,2>& a, const std::array<double,2>& b) { return a[1] > b[1]; } );
    
    //======================================== Find settings max's and min's
    double minSpacing                         = std::numeric_limits<double>::infinity ( );
    unsigned int maxBand                      = 0;
    unsigned int maxTheta                     = 0;
    unsigned int maxPhi                       = 0;
    unsigned int maxGlInt                     = 0;
    double maxExtraSpace                      = 0.0;
    for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( setUp->structFiles.size ( ) ); iter++ )
    {
        if ( shSpacing.at ( iter )  < minSpacing )    { minSpacing = shSpacing.at ( iter );   }
        if ( bandVec.at ( iter )    > maxBand    )    { maxBand = bandVec.at ( iter );        }
        if ( thetaVec.at ( iter )   > maxTheta   )    { maxTheta = thetaVec.at ( iter );      }
        if ( phiVec.at ( iter )     > maxPhi     )    { maxPhi = phiVec.at ( iter );          }
        if ( glIntegVec.at ( iter ) > maxGlInt   )    { maxGlInt = glIntegVec.at ( iter );    }
        if ( extraVec.at ( iter ) > maxExtraSpace)    { maxExtraSpace = extraVec.at ( iter ); }
    }
   
    double minSpacingCopy                     = minSpacing;
    unsigned int maxBandCopy                  = maxBand;
    unsigned int maxThetaCopy                 = maxTheta;
    unsigned int maxPhiCopy                   = maxPhi;
    unsigned int maxGlIntCopy                 = maxGlInt;
    double maxExtraSpaceCopy                  = maxExtraSpace;
    
    //======================================== Build the database
    setUp->taskToPerform                      = ProSHADE::BuildDB;
    setUp->databaseName                       = "testDB.bin";

    setUp->shellSpacing                       = minSpacing;
    setUp->bandwidth                          = maxBand;
    setUp->theta                              = maxTheta;
    setUp->phi                                = maxPhi;
    setUp->glIntegOrder                       = maxGlInt;
    setUp->extraSpace                         = maxExtraSpace;
    
    setUp->clearMapData                       = false;
    setUp->saveWithAndWithout                 = true;
    
    ProSHADE::ProSHADE *run                   = new ProSHADE::ProSHADE ( setUp );
    delete run;
    
    //======================================== Read in the map
    ProSHADE_internal::ProSHADE_data *inMap   = new ProSHADE_internal::ProSHADE_data ( );
    inMap->getDensityMapFromMAP               ( mapName,
                                               &minSpacingCopy,
                                                setUp->mapResolution,
                                               &maxBandCopy,
                                               &maxThetaCopy,
                                               &maxPhiCopy,
                                               &maxGlIntCopy,
                                               &maxExtraSpaceCopy,
                                                setUp->mapResDefault,
                                                setUp->rotChangeDefault,
                                                setUp,
                                                setUp->overlayDefaults );
    
    //======================================== Map the query to the shells
    inMap->keepPhaseInMap                     ( setUp->alpha,
                                                setUp->bFactorChange,
                                               &maxBandCopy,
                                               &maxThetaCopy,
                                               &maxPhiCopy,
                                               &maxGlIntCopy,
                                               setUp,
                                               setUp->useCOM,
                                               setUp->noIQRsFromMap,
                                               setUp->verbose,
                                               setUp->clearMapData,
                                               setUp->rotChangeDefault,
                                               setUp->overlayDefaults,
                                               setUp->maskBlurFactor,
                                               setUp->maskBlurFactorGiven );

    inMap->mapPhaselessToSphere               ( setUp,
                                                maxThetaCopy,
                                                maxPhiCopy,
                                                minSpacingCopy,
                                                setUp->manualShells,
                                                true,
                                                false );
    
    //======================================== Database comparison settings
    ProSHADE::ProSHADE_settings* dbCompSetUp  = new ProSHADE::ProSHADE_settings ( );
    
    // ... Settings regarding which distances to compute
    dbCompSetUp->energyLevelDist              = true;
    dbCompSetUp->traceSigmaDist               = true;
    dbCompSetUp->fullRotFnDist                = true;
    
    // ... Settings regarding hierarchical distance computation
    dbCompSetUp->enLevelsThreshold            = -999.9;
    dbCompSetUp->trSigmaThreshold             = -999.9;
    
    // ... Settings regarding the task
    dbCompSetUp->taskToPerform                = ProSHADE::Distances;
    
    // ... Settings regarding the database
    dbCompSetUp->databaseName                 = setUp->databaseName;
    dbCompSetUp->volumeTolerance              = setUp->volumeTolerance;
    dbCompSetUp->clearMapData                 = false;
    dbCompSetUp->extraSpace                   = setUp->extraSpace;
    
    // ... Settings regarding loudness
    dbCompSetUp->verbose                      = setUp->verbose;
    dbCompSetUp->htmlReport                   = setUp->htmlReport;

    //======================================== For each input structure, fragment the map
    for ( unsigned int strIt = 0; strIt < static_cast<unsigned int> ( strRadii.size ( ) ); strIt++ )
    {
        //==================================== Fragmentation set up
        setUp->mapFragBoxSize                 = strRadii.at(strIt)[1] * 2.0;
        
        //==================================== Do the fragmentation
        std::vector<ProSHADE_internal::ProSHADE_data*> frags = inMap->fragmentMap ( setUp, false );
        
        //==================================== For each passing fragment, check against the input db
        for ( unsigned int frIt = 0; frIt < static_cast<unsigned int> ( frags.size() ); frIt++ )
        {
            //================================ Initialise fragment comparison settings
            std::stringstream frgName;
            frgName << "proshade_tmp_frag" << strIt << "_" << frIt << ".map";
            dbCompSetUp->structFiles.clear    ( );
            dbCompSetUp->structFiles.emplace_back ( frgName.str() );
            
            //================================ Print the fragment - this could be made fasted by using the internal representation, but it would also require more work. Will get back to this in due time.
            frags.at(frIt)->writeMap          ( frgName.str(), frags.at(frIt)->getMap ( ) );
            
            //================================ Now, compare the fragment against the database
            ProSHADE::ProSHADE *compareFrag   = new ProSHADE::ProSHADE ( dbCompSetUp );
            
            //================================ And get distances
            std::vector<double> crossCorrDists = compareFrag->getCrossCorrDists ( );
            std::vector<double> traceSigmaDists = compareFrag->getTraceSigmaDists ( );
            std::vector<double> rotFunDists   = compareFrag->getRotFunctionDists ( );
            
            //================================ Release the computation object
            delete compareFrag;
            
            for ( int i = 0; i < static_cast<int> ( crossCorrDists.size() ); i++ )
            {
//                if ( rotFunDists.at(i) > 0.1 )
                {
                    printf ( "%i\t%i\t||\t%+.3f\t\t%+.3f\t\t%+.3f\t\t%s\n", strIt, frIt, crossCorrDists.at(i), traceSigmaDists.at(i), rotFunDists.at(i), dbCompSetUp->structFiles.at(i+1).c_str() );
                }
            }
        }
        
        std::cout << " ... " << std::endl; exit (0);
        
        //==================================== Clean up
        for ( unsigned int iter = 0; iter < static_cast<unsigned int> ( frags.size ( ) ); iter++ )
        {
            delete frags.at(iter);
        }
        frags.clear                           ( );
    }
    
    //======================================== Free database comparison memory
    delete dbCompSetUp;

    //======================================== Free memory
    delete inMap;
    delete setUp;
    
    //======================================== DONE
    return 0;
}