CoalescentSimTypes.cc

/*

Copyright (C) 2003-2009 Kevin Thornton, krthornt[]@[]uci.edu

Remove the brackets to email me.

This file is part of libsequence.

libsequence is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

libsequence 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 General Public License for more details.

You should have received a copy of the GNU General Public License
long with libsequence.  If not, see <http://www.gnu.org/licenses/>.

*/

#include <Sequence/Coalescent/SimTypes.hpp>
#include <cassert>
#include <iostream>
#include <cstdlib>

namespace Sequence
{
  segment::segment() : beg(0),end(0),desc(0)
  {
  }

  segment::segment(const int & b,
                   const int & e,
                   const int & d = 0) : beg(b),end(e),desc(d)
  {
    assert ( end >= beg );
  }

  chromosome::chromosome() : segs(NULL),
                             pop(0),
                             nsegs(0)
  {
  }

  chromosome::chromosome( const std::vector<segment> & initial_segs,
                          const int & population ) : 
    segs((segment *)malloc(initial_segs.size()*
                           sizeof(segment))), 
    pop(population),
    nsegs(int(initial_segs.size()))
  {
    std::copy(initial_segs.begin(),initial_segs.end(),segs);
  }

  chromosome::~chromosome()
  {
    if(nsegs>0) free(segs);
  }

  chromosome::chromosome( const chromosome & ch ) :
    segs((segment *)malloc(ch.nsegs*sizeof(segment))),
    pop(ch.pop),
    nsegs(ch.nsegs)
  {
    std::copy(ch.segs,ch.segs+ch.nsegs,segs);
  }

  chromosome & chromosome::operator=(const chromosome & ch)
  {
    if(this == &ch) return *this;
    if(ch.nsegs>this->nsegs)
      {
        this->segs = (segment *)realloc(this->segs,ch.nsegs*sizeof(segment));
      }
    std::copy(ch.segs,ch.segs+ch.nsegs,segs);
    this->nsegs = ch.nsegs;
    this->pop = ch.pop;
    return *this;
  }


  void chromosome::swap_with( chromosome & ch )
  {
    std::swap(this->segs,ch.segs);
    std::swap(this->nsegs,ch.nsegs);
    std::swap(this->pop,ch.pop);
  }

  void chromosome::assign_allocated_segs( segment * newsegs,
                                          const int & new_nsegs )
  {
    free(segs);
    segs=newsegs;
    nsegs=new_nsegs;
  }

  chromosome::iterator chromosome::begin()
  {
    return segs;
  }

  chromosome::iterator chromosome::end()
  {
    return segs+nsegs;
  }

  chromosome::const_iterator chromosome::begin() const
  {
    return segs;
  }

  chromosome::const_iterator chromosome::end() const
  {
    return segs+nsegs;
  }

  int chromosome::links() const
  {
    return( (segs+nsegs-1)->end - segs->beg );
  }

  node::node(const double & t,
             const int & a) : time(t),abv(a)
  {
  }

  marginal::marginal(const int & b, const int & ns,
                     const int & nn,
                     const std::vector<node> & t) :
    beg(b),nsam(ns),nnodes(nn),tree(t)
  {
  }

  marginal::iterator marginal::begin()
  {
    return tree.begin();
  }

  marginal::iterator marginal::end()
  {
    return tree.end();
  }

  marginal::const_iterator marginal::begin() const
  {
    return tree.begin();
  }

  marginal::const_iterator marginal::end() const
  {
    return tree.end();
  }

  marginal::reference marginal::operator[](const std::vector<node>::size_type &i)
  {
    assert( int(i) <= 2*nsam-2 );
    return tree[i];
  }

  marginal::const_reference marginal::operator[](const std::vector<node>::size_type &i) const
  {
    assert( int(i) <= 2*nsam-2 );
    return tree[i];
  }

  bool marginal::operator<(const marginal & rhs) const
  {
    return this->beg < rhs.beg;
  }

  std::ostream & operator<<(std::ostream & s,const  chromosome & c)
  {
    s << '(';
    chromosome::const_iterator seg = c.begin();
    for( ; seg < c.end() ; 
         ++seg)
      {
        s << seg->beg << ':' << seg->end << ',';
      }
    s << ' ' << c.nsegs << ' ' << c.links() << ')';
    return s;
  }

  std::ostream & operator<<(std::ostream & s, const marginal & m)
  {
    s << "marginal tree begins at: " << m.beg << '\n';
    int seg=0;
    for( ;seg < m.nnodes ; ++seg)
      {
        s << seg << ' ' 
          << (m.begin()+seg)->time << ' '
          << (m.begin()+seg)->abv << '\n';
      }
    s << seg << ' '
      << (m.begin()+seg)->time << ' '
      << (m.begin()+seg)->abv;
    return s;
  }

  struct newick_stream_marginal_tree_impl
  {
    marginal::const_iterator mi;
    const int nsam;
    std::vector<int> left,right;
    std::vector<node> tree;
    void init();
    std::ostream & parens( const int & noden,
                           std::ostream & o) const;
    newick_stream_marginal_tree_impl(const marginal & m);
    newick_stream_marginal_tree_impl(const marginal * m);
    newick_stream_marginal_tree_impl(arg::const_iterator m);
    newick_stream_marginal_tree_impl(arg::iterator m);
  };

  newick_stream_marginal_tree_impl::newick_stream_marginal_tree_impl(const marginal & m) :
    mi(m.begin()),nsam(m.nsam),
    left(2*nsam-1,-1),right(left),
    tree(std::vector<node>())
  {
    init();
  }

  newick_stream_marginal_tree_impl::newick_stream_marginal_tree_impl(const marginal * m) :
    mi(m->begin()),nsam(m->nsam),
    left(2*nsam-1,-1),right(left),
    tree(std::vector<node>())
  {
    init();
  }

  newick_stream_marginal_tree_impl::newick_stream_marginal_tree_impl(arg::const_iterator m) :
    mi(m->begin()),nsam(m->nsam),
    left(2*nsam-1,-1),right(left),
    tree(std::vector<node>())
  {
    init();
  }

  newick_stream_marginal_tree_impl::newick_stream_marginal_tree_impl(arg::iterator m) :
    mi(m->begin()),nsam(m->nsam),
    left(2*nsam-1,-1),right(left),
    tree(std::vector<node>())
  {
    init();
  }

  std::ostream &
  newick_stream_marginal_tree_impl::parens( const int & noden,
                                            std::ostream & o) const
  {
    double time;

    if( left[noden] == -1 ) 
      {
        assert( (mi+((mi+noden)->abv))->time >= 0. );
        o << noden+1 << ':' <<(mi+((mi+noden)->abv))->time;
      }
    else
      {
        o << '(';
        parens( left[noden], o ) ;
        o << ',';
        parens( right[noden], o ) ;
        if( (mi+noden)->abv == -1 ) 
          {
            o << ");";
          }
        else 
          {
            time = (mi + (mi+noden)->abv )->time - (mi+noden)->time ;
            assert(time >= 0.);
            o << "):"<<time;
          }
      }
    return o;
  }

  void newick_stream_marginal_tree_impl::init()
  {
    for(int i=0;i<2*nsam-2;++i)
      {
        if( left[ (mi+i)->abv ] == -1 ) 
          {
            left[(mi+i)->abv] = i;
          }
        else 
          {
            right[ (mi+i)->abv] = i;
          }
      }

  }

  newick_stream_marginal_tree::newick_stream_marginal_tree( const marginal & m ) :
    impl( new newick_stream_marginal_tree_impl(m) )
  {
  }

  newick_stream_marginal_tree::newick_stream_marginal_tree( const marginal * m ) :
    impl( new newick_stream_marginal_tree_impl(m) )
  {
  }

  newick_stream_marginal_tree::newick_stream_marginal_tree( arg::const_iterator m ) :
    impl( new newick_stream_marginal_tree_impl(m) )
  {
  }

  newick_stream_marginal_tree::newick_stream_marginal_tree( arg::iterator m ) :
    impl( new newick_stream_marginal_tree_impl(m) )
  {
  }

  newick_stream_marginal_tree::~newick_stream_marginal_tree( )
  {
    delete impl;
  }

  std::vector<node> newick_stream_marginal_tree::get_tree() const
  {
    return impl->tree;
  }

  std::ostream & newick_stream_marginal_tree::print( std::ostream & o ) const
  {
    impl->parens( 2*(impl->nsam)-2,o);
    return o;
  }

  std::istream & newick_stream_marginal_tree::read( std::istream & i )
  {
    return i;
  }

  std::ostream & operator<<(std::ostream & o, const newick_stream_marginal_tree & n)
  {
    return n.print(o);
  }

  std::istream & operator>>(std::istream & i, newick_stream_marginal_tree & n)
  {
    return n.read(i);
  }
}

---