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
{
  namespace coalsim {
    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)))), 
      segs( static_cast<segment*>(malloc(initial_segs.size()*sizeof(segment)))), 
      pop(population),
      nsegs(unsigned(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))),
      segs(static_cast<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));
          this->segs = static_cast<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 unsigned & 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;
    }

    class newick_stream_marginal_tree_impl
    {
    public:
      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(std::vector<int>::size_type(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(std::vector<int>::size_type(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(std::vector<int>::size_type(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(std::vector<int>::size_type(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
    {
      if( left[std::vector<int>::size_type(noden)] == -1 ) 
        {
          assert( (mi+((mi+noden)->abv))->time >= 0. );
          o << noden+1 << ':' <<(mi+((mi+noden)->abv))->time;
        }
      else
        {
          o << '(';
          parens( left[std::vector<int>::size_type(noden)], o ) ;
          o << ',';
          parens( right[std::vector<int>::size_type(noden)], o ) ;
          if( (mi+noden)->abv == -1 ) 
            {
              o << ");";
            }
          else 
            {
              double 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[ std::vector<int>::size_type( (mi+i)->abv ) ] == -1 ) 
            {
              left[std::vector<int>::size_type( (mi+i)->abv ) ] = i;
            }
          else 
            {
              right[ std::vector<int>::size_type( (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( )
    {
    }

    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);
    }
  }
}