main page namespaces classes files
CPSP home

src/cpsp/gecode/GC_ThreadingSpace.cc

Go to the documentation of this file.

/*
 *  Main authors:
 *     Martin Mann http://www.bioinf.uni-freiburg.de/~mmann/
 *
 *  Contributing authors:
 *     Sebastian Will http://www.bioinf.uni-freiburg.de/~will/
 *
 *  Copyright:
 *     Martin Mann, 2007
 *
 *  This file is part of the CPSP-tools package:
 *     http://www.bioinf.uni-freiburg.de/sw/cpsp/
 *
 *  See the file "LICENSE" for information on usage and
 *  redistribution of this file, and for a
 *     DISCLAIMER OF ALL WARRANTIES.
 *
 */



//  #define TMPOUT

#include "cpsp/gecode/GC_ThreadingSpace.hh"
#include "cpsp/gecode/GC_Search.hh"
#include "cpsp/gecode/GC_StlSetRangeIterator.hh"
#include "cpsp/gecode/GC_LatticeNeighbored2.hh"
#include "cpsp/gecode/GC_ValNearCenter.hh"
#include "cpsp/gecode/GC_RankViewSel.hh"
#include "cpsp/gecode/GC_FirstSolBranching.hh"

#include <gecode/kernel.hh>
#include <gecode/int.hh>
#include <gecode/search.hh>
#include <gecode/minimodel.hh>
#include <gecode/int/branch.hh>

#include <biu/LatticeFrame.hh>

#ifdef TMPOUT
    #include <iostream>
#endif

namespace cpsp
{
  namespace gecode
  {

        // just a helper function ... TODO : delete ;)
    void print(biu::IndexSet& x, const biu::LatticeFrame* latFrame) {
        std::cout<<"( ";
        for (biu::IndexSet::const_iterator it = x.begin(); it != x.end(); it++)
            std::cout <<latFrame->getPoint(*it) <<"| ";
        std::cout <<")";
        std::cout.flush();
    }

        // just a helper function ... TODO : delete ;)
    void GC_ThreadingSpace::print() const {
      std::cerr<<"\n";
      for (int i=0; i< domains.size(); i++)  {
        const Gecode::Int::IntView vd(domains[i]);
        Gecode::Int::ViewRanges< Gecode::Int::IntView > v(vd);
        std::cerr <<" " <<i <<"(";
        while(v()) {
            if (v.min() != v.max())
                std::cerr <<v.min() <<".." <<v.max()<<"," ;
            else
                std::cerr <<v.min() <<",";
            ++v;
        }
        std::cerr <<")\n";
      }
      std::cerr<<"\n";
        }

    GC_ThreadingSpace::GC_ThreadingSpace(const std::string *sequence,
            const biu::LatticeFrame* latFrame,
            const biu::IndexVec* neighVecs,
            SeqFeatureMap* seqFeatureMap,
            const HullLevel *hullLvl,
            HCore* hCore,
            const GC_ThreadingSymmBreaker::GlobalShiftVec* shiftVec,
            int branchingType,
            NeighPropLvl neighPropLvl ) :
        domains(    this, sequence->size(),
                    0, Gecode::Limits::Int::int_max), rank(sequence->size(), 999)
    {
        assert(hCore != NULL);
            // temp variables
        biu::IndexSet data;
        GC_StlSetRangeIterator dataIt;
        Gecode::IntSet dataSet;
        Gecode::IntVarArgs hVars(hCore->getSize()), pVars(domains.size() - hCore->getSize());
        int hVarI = 0, pVarI = 0;

            // domain init

        int i=0;

        for (i=0; i<domains.size(); i++) {
                // get domain data
            data = hCore->getIndexedHull(latFrame, hullLvl->at(i));
                // create iterator 
            dataIt.init(&data);
                // shrink domains to given data in dataIt via IntView
            Gecode::Int::IntView actDom(domains[i]);
                // use macro for security
            GECODE_ME_FAIL(this, actDom.inter(this,dataIt));
                // init distinct viewvecs
            if (hullLvl->at(i) == 0) {
                hVars[hVarI++] = domains[i];
            } else {
                pVars[pVarI++] = domains[i];
            }
        }

            // p-singlet init
        std::vector<unsigned int>::const_iterator it;
            // get domain data
        data = hCore->getIndexedPSinglets(latFrame);
            // create iterator 
        dataIt.init(&data);
            // shrink domains 
        for (   it = (*seqFeatureMap)[P_SINGLET].begin();
                it != (*seqFeatureMap)[P_SINGLET].end(); it++)
        {
                // view to singlet domain
            Gecode::Int::IntView actDom(domains[*it]);
                // use macro for security
            GECODE_ME_FAIL(this, actDom.inter(this,dataIt));
                // set rank 
            rank[*it] = 0;
        }

            // h-singlet init
            // get domain data
        data = hCore->getIndexedHSinglets(latFrame);
            // create iterator 
        dataIt.init(&data);
            // shrink domains 
        for (   it = (*seqFeatureMap)[H_SINGLET].begin();
                it != (*seqFeatureMap)[H_SINGLET].end(); it++)
        {
                // view to singlet domain
            Gecode::Int::IntView actDom(domains[*it]);
                // use macro for security
            GECODE_ME_FAIL(this, actDom.inter(this,dataIt));
                // set rank 
            rank[*it] = 0;
        }

            // init remaining rank stuff
        {
            std::string seqPattern = "HHPP";
            std::string::size_type pos = 0, shift = 1;
            while( (pos = sequence->find(seqPattern,pos)) != std::string::npos) {
                pos += shift;
                rank [pos] = 1;
            }
            seqPattern = "PPHH";
            pos = 0, shift = 2;
            while( (pos = sequence->find(seqPattern,pos)) != std::string::npos) {
                pos += shift;
                rank [pos] = 1;
            }
            if (domains.size() > 1) {
                if ((*sequence)[0] != (*sequence)[1])
                    rank[0] = 1;
                if ((*sequence)[sequence->size()-1] != (*sequence)[sequence->size()-2])
                    rank[sequence->size()-1] = 1;
            }
        }

            // init constraints 

            // ALL DIFFERENT
            // Supports value (ICL_VAL, default), bounds (ICL_BND),
            // and domain-consistency (ICL_DOM)
        Gecode::distinct(this, hVars, Gecode::ICL_DOM);
        Gecode::distinct(this, pVars, Gecode::ICL_DOM);

            // binary lattice neighbor constraints along the sequence
        switch(neighPropLvl) {
        
        case CUSTOM_2 : {   
                        for(i=1; i<domains.size(); i++) {
                            GECODE_ES_FAIL(this,
                            GC_LatticeNeighbored2::post(    this,
                                                    domains[i-1], domains[i],
                                                    neighVecs,
                                                    Gecode::ICL_DOM,
                                                    150
                                                ));
                        }
                        break; }
        case NO_PROP : 
        default :   //  NO NEIGHBOR PROPAGATION INSERTED
                    break;
        }

            // GC_ValNearCenter lattice frame init
        GC_ValNearCenter::latFrame = latFrame;

        // branching
        
        if ( branchingType != BR_NONE ) {

            assertbiu( !(branchingType & BR_DDS) , "DDS branching currently not supported"); 

                // viewarray for branching
            Gecode::ViewArray<MyViewType> xv(this,domains.size());
            for (i=0; i< domains.size(); i++)
                xv[i] = IdxIntView(domains[i], i);
                
                // symmetry breaking
            if ( branchingType & BR_SYM ) {
//              Gecode::ViewArray<Gecode::Int::IntView> xv(this,(Gecode::IntVarArgs)domains);
    
                Gecode::BoolVarArray cpInit = GC_ThreadingSymmBreaker::initCpByCore(this, hCore, latFrame, shiftVec);
                GC_ThreadingSymmBreaker symmBreaker(this,latFrame, shiftVec, cpInit);
    
                new (this)
                Gecode::Search::SBViewValBranching
//                  <Gecode::Int::IntView,
                    <MyViewType,
                    int,
//                  Gecode::Int::Branch::BySizeMin,
                    GC_RankViewSel<Gecode::Int::Branch::ByDegreeMax>,
                    GC_ValNearCenter,
                    GC_ThreadingSymmBreaker>(this, xv,symmBreaker);
            }
    
                // DFS branching
            if ( branchingType & BR_DFS ) {
                typedef Gecode::ViewValBranching<   MyViewType,
                                                    int,
                                                    GC_RankViewSel<Gecode::Int::Branch::ByDegreeMax>,
                                                    GC_ValNearCenter
                                                > DFS;
                                                
                // the branching
                new (this) DFS(this, xv);
            }
        }
    }





    GC_ThreadingSpace::GC_ThreadingSpace(bool share, GC_ThreadingSpace& toCopy):
        SuperSpace(share, toCopy), domains(), rank(toCopy.rank)
    {
        domains.update(this, share, toCopy.domains);

    }


    GC_ThreadingSpace::~GC_ThreadingSpace()
    {
    }


    SuperSpace*
    GC_ThreadingSpace::copy(bool share) {
        return new GC_ThreadingSpace(share, *this);
    }

    biu::IndexVec
    GC_ThreadingSpace::getSolution() const {
        biu::IndexVec retVec;
        for (int i=0; i<domains.size();i++) {
            if (domains[i].assigned()) {
                retVec.push_back((biu::LatticeFrame::index_type)domains[i].val());
            } else {
                std::cerr <<"Error: Variable "<<i<<" is not assigned !!! size = "<<domains[i].size()<<"\n"; 
                std::cerr <<"       Are you using DDS ?\n";
            }
        }
        return retVec;
    }


    int
    GC_ThreadingSpace::getRank(int index) const {
        assert(index >= 0 && index < (int)rank.size());
        return rank[index];
    }



// GC_ThreadingSpaceShapes



    GC_ThreadingSpaceShapes::GC_ThreadingSpaceShapes(
        const std::string *sequence,
        const biu::LatticeFrame* latFrame,
        const biu::IndexVec* neighVecs,
        SeqFeatureMap* seqFeatureMap,
        const HullLevel *hullLvl,
        HCore* hCore,
        const GC_ThreadingSymmBreaker::GlobalShiftVec* shiftVec,
        int branchingType)
      : GC_ThreadingSpace( sequence, latFrame, neighVecs, seqFeatureMap,
            hullLvl, hCore, shiftVec, BR_NONE) 
    {
        // posting branchings
        if ( branchingType != BR_NONE ) {

            assertbiu( !(branchingType & BR_DDS) , "DDS branching currently not supported"); 
        
                // viewarrays for branching
            Gecode::ViewArray<MyViewType> xvH(this,hCore->getSize());
            Gecode::ViewArray<MyViewType> xvP(this,domains.size() - hCore->getSize());
                // fill arrays
            unsigned int iH=0, iP=0;
            for (int i=0; i< domains.size(); i++) {
                if ((*sequence)[i] == 'H') {
                    xvH[iH++] = IdxIntView(domains[i], i);
                } else {
                    xvP[iP++] = IdxIntView(domains[i], i);
                }
            }
            if ( iH != hCore->getSize() || iP != domains.size()-hCore->getSize() )
             std::cerr <<" init failed: iH=" <<iH <<"/" <<hCore->getSize() <<" iP=" <<iP <<"/" <<(domains.size()-hCore->getSize())<<std::endl;
                
                // symmetry breaking on H-monomers
            if ( branchingType & BR_SYM ) {
//              Gecode::ViewArray<Gecode::Int::IntView> xv(this,(Gecode::IntVarArgs)domains);
    
                Gecode::BoolVarArray cpInit = GC_ThreadingSymmBreaker::initCpByCore(this, hCore, latFrame, shiftVec);
                GC_ThreadingSymmBreaker symmBreaker(this,latFrame, shiftVec, cpInit);
    
                new (this)
                Gecode::Search::SBViewValBranching
//                  <Gecode::Int::IntView,
                    <MyViewType,
                    int,
//                  Gecode::Int::Branch::BySizeMin,
                    GC_RankViewSel<Gecode::Int::Branch::ByDegreeMax>,
                    GC_ValNearCenter,
                    GC_ThreadingSymmBreaker>(this, xvH,symmBreaker);
            }
    
                // DFS branching on remaining H-monomers
            if ( branchingType & BR_DFS ) {
                new (this) Gecode::ViewValBranching<   MyViewType,
                                                    int,
                                                    GC_RankViewSel<Gecode::Int::Branch::ByDegreeMax>,
                                                    GC_ValNearCenter
                                                >(this, xvH);
            }
                // find only one valid solution for the remaining P-monomers
            new (this) FirstSolBranching<   MyViewType,
                                            int,
                                            GC_RankViewSel<Gecode::Int::Branch::ByDegreeMax>,
                                            GC_ValNearCenter
                                        >( this, xvP );
        }
    }




  } // namespace gecode
}  // namespace cpsp