src/cpsp/SideChain/chain_threading.cpp

Go to the documentation of this file.
00001 /*
00002  *  Main authors:
00003  *     Mohamad Rabbath <rabbath@informatik.uni-freiburg.de>
00004  *
00005  *  Contributing authors:
00006  *     Martin Mann <mmann@informatik.uni-freiburg.de>
00007  *     Sebastian Will <will@informatik.uni-freiburg.de>
00008  *
00009  *  This file is part of the CPSP-tools package:
00010  *     http://www.bioinf.uni-freiburg.de/sw/cpsp/
00011  *
00012  *  See the file "LICENSE" for information on usage and
00013  *  redistribution of this file, and for a
00014  *     DISCLAIMER OF ALL WARRANTIES.
00015  *
00016  */
00017 
00018 #include "chain_threading.h"
00019 #include <biu/Point.hh>
00020 #include<iostream>
00021 #include <cpsp/gecode/GC_ThreadingSymmBreaker.hh>
00022 #include <cpsp/gecode/GC_RankViewSel.hh>
00023 #include "options/chain_option.h"
00024 #include <biu/LatticeDescriptorCUB.hh>
00025 namespace cpsp{
00026 //side chain threading
00027 namespace scth{
00028 
00029 SideChainThreading::SideChainThreading(const std::string *sequence,
00030         const biu::LatticeFrame* latFrame, 
00031         const biu::IndexVec* neighVecs,cpsp::HCore* hCore,
00032         const biu::IPointVec* shiftVec,
00033         bool withBreakingSymmetry){
00034     //Set all the domains
00035     assert(sequence!=NULL);
00036 
00037     //Settting the domain of the Hs-Array
00038     setHDomain(latFrame,hCore);
00039     //Setting the domain of the BHs-Array
00040     setBHDomain(latFrame,hCore);
00041     //Setting the domain of the Ps-Array
00042     setOutCoreDomain(sequence,latFrame,hCore);
00043     //Now applying all constrains
00044     //All H's are in the HCore
00045     inCore(hCore);
00046 
00047     //All the back bones of H's are 1 distance unit far from the H-Core
00048     attachedToCore(hCore );
00049 
00050     //All none H are out of the HCore
00051     outCore(sequence,hCore );
00052 
00053     //Just for test
00054     /*biu::IndexSet test= hCore->getIndexedHull(latFrame, 0);
00055        for (biu::IndexSet::const_iterator t=test.begin(); t!= test.end();t++){
00056           std::cout<<(*t)<<std::endl;
00057           biu::IntPoint tpoint=latFrame->getPoint(*t);
00058           std::cout<<" " <<tpoint.getX()<<" "<<tpoint.getY()<<" "<<tpoint.getZ()<<std::endl;
00059        }
00060       biu::IntPoint ph1(7,8,9);
00061       int phi1=latFrame->getIndex(ph1);
00062       Hs[3]=IntVar(this,phi1,phi1);
00063       biu::IntPoint ph2(8,8,8);
00064       int phi2=latFrame->getIndex(ph2);
00065       Hs[0]=IntVar(this,phi2,phi2);
00066       biu::IntPoint ph3(8,8,9);
00067       int phi3=latFrame->getIndex(ph3);
00068       Hs[2]=IntVar(this,phi3,phi3);
00069       biu::IntPoint ph4(9,8,9);
00070       int phi4=latFrame->getIndex(ph4);
00071       Hs[1]=IntVar(this,phi4,phi4);*/
00072     /*biu::IntPoint p1(7,9,9);
00073       int pi1=latFrame->getIndex(p1);
00074       BHs[3]=IntVar(this,pi1,pi1);
00075       biu::IntPoint p2(8,9,8);
00076       int pi2=latFrame->getIndex(p2);
00077       BHs[0]=IntVar(this,pi2,pi2);
00078       biu::IntPoint p3(8,9,9);
00079       int pi3=latFrame->getIndex(p3);
00080       BHs[2]=IntVar(this,pi3,pi3);
00081       biu::IntPoint p4(9,9,9);
00082       int pi4=latFrame->getIndex(p4);
00083       BHs[1]=IntVar(this,pi4,pi4); 
00084       biu::IntPoint p5(8,8,7);
00085       int pi5=latFrame->getIndex(p5);
00086       Ps[0]=IntVar(this,pi5,pi5);
00087       biu::IntPoint p6(9,8,8);
00088       int pi6=latFrame->getIndex(p6);
00089       Ps[1]=IntVar(this,pi6,pi6);
00090       biu::IntPoint p7(8,9,7);
00091       int pi7=latFrame->getIndex(p7);
00092       BPs[0]=IntVar(this,pi7,pi7);
00093       biu::IntPoint p8(9,9,8);
00094       int pi8=latFrame->getIndex(p8);
00095       BPs[1]=IntVar(this,pi8,pi8);*/
00096     //End test
00097 
00098     //Neighbouring constrains for the back bones
00099     BsNeighbours(sequence,neighVecs);
00100 
00101     //Neighbourig constrains between back bones and side chains P
00102     B_PsNeighbours(neighVecs);
00103 
00104     //Neighbourig constrains between back bones and side chains P
00105     B_HsNeighbours(neighVecs);
00106 
00107     //Self avoiding
00108     selfAvoid();
00109     //Branching
00110     ViewArray<Int::IntView> allView=ViewArray<Int::IntView>(this,Hs.size()+BHs.size()+BPs.size()+Ps.size());
00111     //Converting to the index type
00112     typedef cpsp::gecode::IdxIntView IndexType;
00113     Gecode::ViewArray<IndexType> viewIndex(this,allView.size());
00114     int qindx=0;
00115 
00116     for(int indx=0;indx<Hs.size();indx++){
00117         allView[qindx]=Hs[indx];
00118         viewIndex[qindx]=cpsp::gecode::IdxIntView(allView[qindx], qindx);
00119         qindx++;
00120     }
00121 
00122     for(int indx=0;indx<BHs.size();indx++){
00123         allView[qindx]=BHs[indx];
00124         viewIndex[qindx]=cpsp::gecode::IdxIntView(allView[qindx], qindx);
00125         qindx++;
00126     }
00127 
00128     for(int indx=0;indx<BPs.size();indx++){
00129         allView[qindx]=BPs[indx];
00130         viewIndex[qindx]=cpsp::gecode::IdxIntView(allView[qindx], qindx);
00131         qindx++;
00132     }
00133 
00134     for(int indx=0;indx<Ps.size();indx++){
00135         allView[qindx]=Ps[indx];
00136         viewIndex[qindx]=cpsp::gecode::IdxIntView(allView[qindx], qindx);
00137         qindx++;
00138     }
00139 
00140     cpsp::gecode::GC_ValNearCenter::latFrame = latFrame;
00141     if(withBreakingSymmetry){
00142         Gecode::BoolVarArray cpInit = cpsp::gecode::GC_ThreadingSymmBreaker::initCpByCore(this, hCore, latFrame, shiftVec);
00143         cpsp::gecode::GC_ThreadingSymmBreaker symmBreaker(this,latFrame, shiftVec, cpInit);
00144         new (this)
00145         Gecode::Search::SBViewValBranching
00146         //                  <Gecode::Int::IntView,
00147         <IndexType,
00148         int,
00149         //                  Gecode::Int::Branch::BySizeMin,
00150         cpsp::gecode::GC_RankViewSel<Gecode::Int::Branch::ByDegreeMax>,
00151         cpsp::gecode::GC_ValNearCenter,
00152         cpsp::gecode::GC_ThreadingSymmBreaker>(this, viewIndex,symmBreaker);
00153     }
00154     //This is branching without breaking the symmetry
00155     new (this) Gecode::ViewValBranching< Int::IntView,int,Int::Branch::ByDegreeMax,
00156     cpsp::gecode::GC_ValNearCenter>
00157     (this, allView);
00158 }
00159 
00160 SideChainThreading::SideChainThreading(bool share, SideChainThreading& old):
00161     Gecode::Space(share, old){
00162     Hs.update(this, share, old.Hs);
00163     Ps.update(this, share, old.Ps);
00164     BHs.update(this, share, old.BHs);
00165     BPs.update(this, share, old.BPs);
00166 }
00167 
00168 Gecode::Space* SideChainThreading::copy(bool share){
00169     return new SideChainThreading(share, *this);
00170 }
00171 
00172 
00173 SideChainThreading::~SideChainThreading(){}
00174 
00175 void SideChainThreading::inCore(cpsp::HCore* hCore){
00176     Hs=IntVarArray(this,hCore->getSize(),HDomain);
00177     //Hs=IntVarArray(this,hCore->getSize(),1000,3000);  
00178 }
00179 
00180 void SideChainThreading::attachedToCore(cpsp::HCore* hCore){
00181     BHs=IntVarArray(this,hCore->getSize(),BHDomain);
00182     //BHs=IntVarArray(this,hCore->getSize(),1000,3000); 
00183 }
00184 //passing the sequence as parameter for the effeciency
00185 void SideChainThreading::printSol(const std::string* seq,const biu::LatticeFrame* latFrame){
00186     int hIndx=0;
00187     int pIndx=0;
00188     //     std::cout<<"#Cubic Protein Solution:"<<std::endl;
00189     for(size_t i=0;i<seq->size();i++){
00190         if( (*seq)[i]=='H' ){
00191             biu::IntPoint point=latFrame->getPoint(Hs[hIndx].val());
00192             std::cout<<"h: " <<point.getX()<<" "<<point.getY()<<" "<<point.getZ()<<std::endl;
00193 
00194             point=latFrame->getPoint(BHs[hIndx].val());
00195             std::cout<<"bh: " <<point.getX()<<" "<<point.getY()<<" "<<point.getZ()<<std::endl;
00196             hIndx++;
00197         }
00198         else{
00199             biu::IntPoint point=latFrame->getPoint(Ps[pIndx].val());
00200             std::cout<<"p: " <<point.getX()<<" "<<point.getY()<<" "<<point.getZ()<<std::endl;
00201 
00202             point= latFrame->getPoint(BPs[pIndx].val());
00203             std::cout<<"bp: " <<point.getX()<<" "<<point.getY()<<" "<<point.getZ()<<std::endl;
00204             pIndx++;
00205         }
00206     }
00207 }
00208 
00209 std::string SideChainThreading::printAsMoves(const std::string* seq,const biu::LatticeFrame* latFrame){
00210     std::string latName=latFrame->getDescriptor()->getName();
00211     const biu::LatticeDescriptorCUB* cubicDescInstance=
00212         dynamic_cast<const biu::LatticeDescriptorCUB*>(latFrame->getDescriptor());
00213         if(cubicDescInstance){
00214             int hIndx=0;
00215             int pIndx=0;
00216             std::string s="";
00217 
00218             //      std::cout<<"#"<< latName<<" Protein Solution:"<<std::endl;
00219             for(size_t i=0;i<seq->size();i++){
00220                 biu::IntPoint point;
00221                 biu::IntPoint bpoint;
00222                 if( (*seq)[i]=='H' ){
00223                     point=latFrame->getPoint(Hs[hIndx].val());
00224                     bpoint=latFrame->getPoint(BHs[hIndx].val());
00225                     hIndx++;
00226                 }
00227                 else{
00228                     point=latFrame->getPoint(Ps[pIndx].val());
00229                     bpoint= latFrame->getPoint(BPs[pIndx].val());
00230                     pIndx++;
00231                 }
00232                 if( bpoint.getX()!=point.getX() ){
00233                     if( bpoint.getX() > point.getX() )
00234                         s+="(L)";
00235                     else s+="(R)";
00236                 }
00237                 else if( bpoint.getY()!=point.getY() ){
00238                     if( bpoint.getY() > point.getY() )
00239                         s+="(B)";
00240                     else s+="(F)";
00241 
00242                 }
00243                 else{
00244                     if( bpoint.getZ() > point.getZ() )
00245                         s+="(D)";
00246                     else s+="(U)";
00247                 }
00248                 //This is for the move to next backbone
00249                 if(i<seq->size()-1){
00250                     biu::IntPoint nextbpoint;
00251                     if( (*seq)[i+1]=='H' ){
00252                         nextbpoint=latFrame->getPoint(BHs[hIndx].val());
00253                     }
00254                     else{
00255                         nextbpoint=latFrame->getPoint(BPs[pIndx].val());
00256                     }
00257                     if( bpoint.getX()!=nextbpoint.getX() ){
00258                         if( bpoint.getX() > nextbpoint.getX() )
00259                             s+="L";
00260                         else s+="R";
00261                     }
00262                     else if( bpoint.getY()!=nextbpoint.getY() ){
00263                         if( bpoint.getY() > nextbpoint.getY() )
00264                             s+="B";
00265                         else s+="F";
00266 
00267                     }
00268                     else{
00269                         if( bpoint.getZ() > nextbpoint.getZ() )
00270                             s+="D";
00271                         else s+="U";
00272                     }
00273 
00274                 }//End if
00275             }//End fore
00276             std::cout<<s<<std::endl;
00277             return s;
00278         }
00279 
00280         else{
00281             int hIndx=0;
00282             int pIndx=0;
00283             std::string s="";
00284             //      std::cout<<"#"<< latName<<" Protein Solution:"<<std::endl;
00285             for(size_t i=0;i<seq->size();i++){
00286                 biu::IntPoint point;
00287                 biu::IntPoint bpoint;
00288                 if( seq->at(i)=='H' ){
00289                     point=latFrame->getPoint(Hs[hIndx].val());
00290                     bpoint=latFrame->getPoint(BHs[hIndx].val());
00291                     hIndx++;
00292                 }
00293                 else{
00294                     point=latFrame->getPoint(Ps[pIndx].val());
00295                     bpoint= latFrame->getPoint(BPs[pIndx].val());
00296                     pIndx++;
00297                 }
00298                 biu::IntPoint absVec = point - bpoint;
00299                 //Get the moves for the FCC
00300                 biu::Move absMove =
00301                     latFrame->getDescriptor()->getNeighborhood().getElement(absVec).getMove();
00302 
00303                 std::string moveStr =
00304                     latFrame->getDescriptor()->getAlphabet()->getString( absMove );
00305                 //
00306                 s+="("+moveStr+")";
00307                 //This is for the move to next backbone
00308                 if(i<seq->size()-1){
00309                     biu::IntPoint nextbpoint;
00310                     if( (*seq)[i+1]=='H' ){
00311                         nextbpoint=latFrame->getPoint(BHs[hIndx].val());
00312                     }
00313                     else{
00314                         nextbpoint=latFrame->getPoint(BPs[pIndx].val());
00315                     }
00316                     biu::IntPoint nextAbsVec = nextbpoint - bpoint;
00317                     //Get the moves for the FCC
00318                     biu::Move nextAbsMove =
00319                         latFrame->getDescriptor()->getNeighborhood().getElement(nextAbsVec).getMove();
00320 
00321                     std::string nextMoveStr =
00322                         latFrame->getDescriptor()->getAlphabet()->getString( nextAbsMove );
00323                     //
00324                     s+=nextMoveStr;
00325 
00326                 }//End if
00327 
00328             }
00329             std::cout<<s<<std::endl;
00330             return s;
00331         }
00332 }
00333 
00334 //Printing for drawing
00335 void SideChainThreading::printAsDrowMoves(const std::string* seq,const biu::LatticeFrame* latFrame){
00336     int hIndx=0;
00337     int pIndx=0;
00338     std::string s="";
00339     //     std::cout<<"#"<< latFrame->getDescriptor()->getName()<<" Protein Solution:"<<std::endl;
00340     for(size_t i=0;i<seq->size();i++){
00341         biu::IntPoint point;
00342         biu::IntPoint bpoint;
00343         if( (*seq)[i]=='H'){
00344             point=latFrame->getPoint(Hs[hIndx].val());
00345             bpoint=latFrame->getPoint(BHs[hIndx].val());
00346             hIndx++;
00347         }
00348         else{
00349             point=latFrame->getPoint(Ps[pIndx].val());
00350             bpoint= latFrame->getPoint(BPs[pIndx].val());
00351             pIndx++;
00352         }
00353         if( bpoint.getX()!=point.getX() ){
00354             if( bpoint.getX() > point.getX() )
00355                 s+="LR";
00356             else s+="RL";
00357         }
00358         else if( bpoint.getY()!=point.getY() ){
00359             if( bpoint.getY() > point.getY() )
00360                 s+="BF";
00361             else s+="FB";
00362 
00363         }
00364         else{
00365             if( bpoint.getZ() > point.getZ() )
00366                 s+="DU";
00367             else s+="UD";
00368         }
00369         //This is for the move to next backbone
00370         if(i<seq->size()-1){
00371             biu::IntPoint nextbpoint;
00372             if( (*seq)[i+1]=='H' ){
00373                 nextbpoint=latFrame->getPoint(BHs[hIndx].val());
00374             }
00375             else{
00376                 nextbpoint=latFrame->getPoint(BPs[pIndx].val());
00377             }
00378             if( bpoint.getX()!=nextbpoint.getX() ){
00379                 if( bpoint.getX() > nextbpoint.getX() )
00380                     s+="L";
00381                 else s+="R";
00382             }
00383             else if( bpoint.getY()!=nextbpoint.getY() ){
00384                 if( bpoint.getY() > nextbpoint.getY() )
00385                     s+="B";
00386                 else s+="F";
00387 
00388             }
00389             else{
00390                 if( bpoint.getZ() > nextbpoint.getZ() )
00391                     s+="D";
00392                 else s+="U";
00393 
00394             }
00395 
00396         }//End if
00397 
00398     }//End fore
00399     std::cout<<s<<std::endl;
00400 }
00401 
00402 std::string SideChainThreading::print(const std::string* seq,const biu::LatticeFrame* latFrame,int option){
00403     assert(latFrame!=NULL);
00404     assert(seq!=NULL);
00405     //Option 0 means print the coordination
00406     if(option==cpsp::scth::NORMAL){
00407         printSol(seq,latFrame);
00408         return "";
00409     }
00410     else if(option==cpsp::scth::MOVES){
00411         return printAsMoves(seq,latFrame);
00412     }
00413     else{
00414         printAsDrowMoves(seq,latFrame);
00415         return "";
00416     }
00417 }
00418 //All none H are out of the HCore
00419 void  SideChainThreading::outCore(const std::string *sequence,cpsp::HCore* hCore){
00420     assert(sequence!=NULL);
00421     Ps=IntVarArray(this,sequence->size()-hCore->getSize(),OutCoreDomain);
00422     BPs=IntVarArray(this,sequence->size()-hCore->getSize(),OutCoreDomain);
00423     //Ps=IntVarArray(this,sequence->size()-hCore->getSize(),1000,3000);
00424     //BPs=IntVarArray(this,sequence->size()-hCore->getSize(),1000,3000);
00425 }
00426 
00427 void SideChainThreading::setHDomain(const biu::LatticeFrame* latFrame,cpsp::HCore* hCore){
00428     biu::IndexSet hCoreDomain= hCore->getIndexedHull(latFrame, 0);
00429     biu::LatticeFrame::index_type HDomainArr[hCoreDomain.size()];
00430     int i=0;
00431     for (biu::IndexSet::const_iterator it=hCoreDomain.begin(); it!= hCoreDomain.end();it++){
00432         HDomainArr[i]=*it;
00433         i++;
00434     }
00435     HDomain=IntSet(HDomainArr,hCoreDomain.size());
00436 }
00437 
00438 void SideChainThreading::setBHDomain(const biu::LatticeFrame* latFrame,cpsp::HCore* hCore){
00439     biu::IndexSet bhCoreDomain= hCore->getIndexedHull(latFrame, 1);
00440     biu::LatticeFrame::index_type BHDomainArr[bhCoreDomain.size()];
00441     int i=0;
00442     for (biu::IndexSet::const_iterator it=bhCoreDomain.begin(); it!= bhCoreDomain.end();it++){
00443         BHDomainArr[i]=*it;
00444         i++;
00445     }
00446     BHDomain=IntSet(BHDomainArr,bhCoreDomain.size());
00447 }
00448 
00449 //Setting the out-core domain
00450 void SideChainThreading::setOutCoreDomain(const std::string *sequence,const biu::LatticeFrame* latFrame,cpsp::HCore* hCore){
00451     int maxHull=getMaxHull(sequence);
00452     biu::IndexSet outCoreDomain= hCore->getIndexedHull(latFrame, maxHull);
00453     biu::LatticeFrame::index_type outDomainArr[outCoreDomain.size()];
00454     int i=0;
00455     for (biu::IndexSet::const_iterator it=outCoreDomain.begin(); it!= outCoreDomain.end();it++){
00456         outDomainArr[i]=*it;
00457         i++;
00458     }
00459     OutCoreDomain=IntSet(outDomainArr,outCoreDomain.size());
00460 }
00461 
00462 //Neighbourig constrains between back bones and side chains P
00463 void SideChainThreading::B_PsNeighbours(const biu::IndexVec* neighVecs){
00464     assert(Ps.size()==BPs.size());
00465     for(int i=0; i<Ps.size(); i++) {
00466         GECODE_ES_FAIL(this,cpsp::gecode::GC_LatticeNeighbored2::post(this,Ps[i],BPs[i],neighVecs,Gecode::ICL_DOM,150));
00467     }
00468 }
00469 
00470 //Neighbourig constrains between back bones and side chains H
00471 void SideChainThreading::B_HsNeighbours(const biu::IndexVec* neighVecs){
00472     assert(Hs.size()==BHs.size());
00473     for(int i=0; i<Hs.size(); i++) {
00474         GECODE_ES_FAIL(this,cpsp::gecode::GC_LatticeNeighbored2::post(this,Hs[i],BHs[i],neighVecs,Gecode::ICL_DOM,150));
00475     }
00476 }
00477 
00478 //Neighbouring constrains for the back bones
00479 void SideChainThreading::BsNeighbours(const std::string *sequence,const biu::IndexVec* neighVecs){
00480     assert(sequence!=NULL);
00481     //The view of all back bones
00482     ViewArray<Int::IntView> backBones=ViewArray<Int::IntView>(this,sequence->size());
00483     int hIndex=0;
00484     int pIndex=0;
00485     for(size_t i=0;i<sequence->size();i++){
00486         if((*sequence)[i]=='H'){
00487             backBones[i]=BHs[hIndex];
00488             hIndex++;
00489         }
00490         else{
00491             backBones[i]=BPs[pIndex];
00492             pIndex++;
00493         }
00494     }
00495 
00496     for(int i=1; i<backBones.size(); i++) {
00497         GECODE_ES_FAIL(this,cpsp::gecode::GC_LatticeNeighbored2::post(this,backBones[i],backBones[i-1],neighVecs,Gecode::ICL_DOM,150));
00498     }
00499 }
00500 
00501 //Self avoiding
00502 void SideChainThreading::selfAvoid(){
00503     distinct(this,Hs, Gecode::ICL_DOM);
00504     //Everything except the HCore
00505     IntVarArgs B_PsVars(Ps.size()+BPs.size()+BHs.size());
00506     int k=0;
00507     for(int i=0;i<Ps.size();i++){
00508         B_PsVars[k]=Ps[i];
00509         k++;
00510     }
00511     for(int i=0;i<BHs.size();i++){
00512         B_PsVars[k]=BHs[i];
00513         k++;
00514     }
00515 
00516     for(int i=0;i<BPs.size();i++){
00517         B_PsVars[k]=BPs[i];
00518         k++;
00519     }
00520     distinct(this,B_PsVars,Gecode::ICL_DOM);
00521 }
00522 
00523 //Get the maximum hull possible hull
00524 int SideChainThreading::getMaxHull(const std::string *sequence){
00525     assert(sequence!=NULL);
00526     int maxHull=0;
00527     bool first=true;
00528     for(size_t i=0;i<sequence->size();){
00529         //Escaping the hs
00530         while((*sequence)[i]=='H' ){
00531 
00532             i++;
00533             if(i>=sequence->size())
00534                 break;
00535         }
00536 
00537         //Break when finishing the sequence
00538         if(i>=sequence->size())
00539             break;
00540 
00541         //Count the ps
00542         int pNumber=0;
00543         int newDistance=0;
00544         while((*sequence)[i]=='P' ){
00545             pNumber++;
00546             i++;
00547             if(i>=sequence->size())
00548                 break;
00549         }//End while
00550 
00551         //If this is the first sequence of ps then take the distance according to the most left p
00552         if(first){
00553             first=false;
00554             newDistance=pNumber;
00555         }
00556         //If this is the last sequence of ps then take the distance according to the most right p
00557         else if(i>=sequence->size()){
00558             newDistance=pNumber;
00559         }
00560         //Take the distance according to the middle p
00561         else{
00562             if(pNumber % 2==0)
00563                 newDistance=pNumber/2;
00564             else newDistance=(pNumber+1)/2;
00565         }
00566         maxHull=std::max(maxHull,newDistance);
00567         if(i>=sequence->size())
00568             break;
00569 
00570     }//end for
00571     //Because of the side chain
00572     maxHull=maxHull+2;
00573     return maxHull;
00574 }
00575 
00576 }
00577 
00578 }
00579 
00580 /*int main(int arg, char** args){
00581 
00582    return 0;
00583 }*/