/*------------------------------------------------------
 Maximum likelihood estimation 
 of migration rate  and effectice population size
 using a Metropolis-Hastings Monte Carlo algorithm                            
 -------------------------------------------------------                        
 M C M C 2   R O U T I N E S 
 
 Tree changing routines
 
 
 Peter Beerli 1996, Seattle
 beerli@scs.fsu.edu
 
Copyright 1996-2002 Peter Beerli and Joseph Felsenstein, Seattle WA
Copyright 2003-2004 Peter Beerli, Tallahassee FL
 
 This software is distributed free of charge for non-commercial use
 and is copyrighted. Of course, we do not guarantee that the software
 works and are not responsible for any damage you may cause or have.
 
$Id: mcmc2.c 733 2007-05-31 02:28:34Z beerli $
-------------------------------------------------------*/
/* \file mcmc2.c
Contains low-level routines to manipulate the treechanger
*/

#include "migration.h"
#include "sighandler.h"
#include "tree.h"
#include "uep.h"
#include "tools.h"
#include "migrate_mpi.h"

#ifdef DMALLOC_FUNC_CHECK
#include <dmalloc.h>
#endif

#define comment(a) FPRINTF(stderr,a)


/* prototypes ------------------------------------------- */
/* changes the tree to the proposed tree by rebuilding parts
   of the tree and also inserts new migration events. */
void coalesce1p (proposal_fmt * proposal);
/* and its pretend version */
void pretendcoalesce1p (proposal_fmt * proposal);

/* private functions */
/* family of tree manipulation routines: dependent on the position of
   target and origin in the tree different routines are used.
   oback: there is no change in the tree topology.
   ocousin: target is the cousin-node, the branch is inserted in its
   former sister branch.
   rbcoa: target is the bottommost node, the ripped branch will be 
   bring the new bottommost node.
   stancoa: all other cases. */
void coat_oback (proposal_fmt * proposal);
void coat_ocousin (proposal_fmt * proposal);
void coat_obbcoa (proposal_fmt * proposal);
void coat_rbcoa (proposal_fmt * proposal);
void coat_stancoa (proposal_fmt * proposal);
/* pretend versions of the above */
void target_oback (proposal_fmt * proposal);
void target_obbcoa (proposal_fmt * proposal);
void target_ocousin (proposal_fmt * proposal);
void target_rbcoa (proposal_fmt * proposal);
void target_stancoa (proposal_fmt * proposal);
/* subroutines in target_stancoa */
void t_s_upper (proposal_fmt * proposal, node * connect, node * obb);
node *t_s_obbncon (proposal_fmt * proposal, node * obb);
void t_s_tncon (proposal_fmt * proposal, node * obb);
void t_s_tcon (proposal_fmt * proposal, node * uppern);
long erase_migr_nodes (node * up);
node *findcrossing (node ** ptrl1, node ** ptrl2);
void connectnodes (node * mother, node * brother, node * sister);
void gotoroot (node * origin, node ** ptrlist);
void adjust (node * theNode, MYREAL tyme, long level);
void localevaluate (node * mother);
void copy_x (proposal_fmt * proposal, xarray_fmt xx1, xarray_fmt xx2);
void fix_root_pop (node * p);
void pseudoevaluate (proposal_fmt * proposal, xarray_fmt x, MYREAL *lx,
                     node * mother, node * newdaughter, MYREAL v);
node *crawl_down (node * theNode, MYREAL tyme);


/*=========================================================*/

void
coalesce1p (proposal_fmt * proposal)
{
    node *obb = NULL;
    //traverse_check(crawlback (proposal->root->next));
    erase_migr_nodes (proposal->origin);
    if (proposal->migr_table_counter2 > 0) //|| proposal->mig_removed)
        erase_migr_nodes (proposal->realtarget);
    if (proposal->target == proposal->ocousin)
    {
        coat_ocousin (proposal);
        if (proposal->oback->tyme > showtop (proposal->oback->back)->tyme
                && crawlback (proposal->oback)->type != 'r')
            error ("problem with time encountered"); /* was >= ** */
    }
    else
    {
        if (((proposal->target == proposal->oback)
                || (proposal->target == proposal->osister)
                || (proposal->target == proposal->origin)))
        {
            obb = proposal->oback->back;
            coat_oback (proposal);
            if (proposal->oback->tyme > showtop (proposal->oback->back)->tyme)
            {   /* was >= ** */
                if (proposal->oback->back->type == 'r')
                {
                    adjust (proposal->root, proposal->oback->tyme + 10000, 0);
                    comment ("Root time adjusted: was to small");
                }
                else
                    error ("problem with time encountered");
            }
        }
        else
        {
            if (crawlback (proposal->target)->type == 'r')
            {
                coat_rbcoa (proposal);
                if (proposal->oback->tyme >
                        showtop (proposal->oback->back)->tyme)
                {  /* was >= ** */
                    if (proposal->oback->back->type == 'r')
                    {
                        adjust (proposal->root, proposal->oback->tyme + 10000,
                                0);
                        comment ("Root time adjusted: was to small");
                    }
                    else
                        error ("problem with time encountered");
                }
            }
            else
            {
                obb = showtop (crawlback (proposal->oback));
                if ((obb != NULL) && (obb == proposal->target))
                {
                    coat_obbcoa (proposal);
                    if (proposal->oback->tyme > showtop (proposal->oback->back)->tyme) /* was >= ** */
                        error ("problem with time encountered");
                }
                else
                {
                    coat_stancoa (proposal);
                    if (proposal->oback->tyme > showtop (proposal->oback->back)->tyme) /* was >= ** */
                        error ("problem with time encountered");
                }
            }
        }
    }
    set_pop (proposal->oback, proposal->origin->actualpop,
             proposal->origin->actualpop);
  //  free_migr_node(proposal->origin, crawlback(proposal->origin));
    if (proposal->migr_table_counter > 0)
    {
        set_pop (proposal->oback,
                 proposal->migr_table[proposal->migr_table_counter - 1].from,
                 proposal->migr_table[proposal->migr_table_counter - 1].from);
        insert_migr_node (proposal->world, proposal->origin,
                          crawlback (proposal->origin), proposal->migr_table,
                          &(proposal->migr_table_counter));
    }
    proposal->migr_table_counter = 0;
    if (proposal->migr_table_counter2)
    {
        insert_migr_node (proposal->world, proposal->realtarget,
                          crawlback (proposal->target), proposal->migr_table2,
                          &(proposal->migr_table_counter2));
        proposal->migr_table_counter2 = 0;
    }
    /*    if (proposal->oback == crawlback(proposal->root->next)) */
    fix_root_pop (crawlback (proposal->root->next));
    //traverse_check(crawlback (proposal->root->next));
}



void
pretendcoalesce1p (proposal_fmt * proposal)
{
    if (proposal->target == proposal->ocousin)
    {
        target_ocousin (proposal);
    }
    else
    {
        if (((proposal->target == proposal->oback)
                || (proposal->target == proposal->osister)
                || (proposal->target == proposal->origin)))
        {
            target_oback (proposal);
        }
        else
        {
            if (crawlback (proposal->target)->type == 'r')
            {
                target_rbcoa (proposal);
            }
            else
            {
                if ((showtop (crawlback (proposal->oback)) != NULL)
                        && (showtop (crawlback (proposal->oback)) ==
                            proposal->target))
                {
                    target_obbcoa (proposal);
                }
                else
                {
                    target_stancoa (proposal);
                }
            }
        }
    }
}

/*================================================================*/

void
coat_oback (proposal_fmt * proposal)
{
    node *obb, *obm;
    obb = showtop (crawlback (proposal->oback));
    connectnodes (showtop (proposal->oback->back), proposal->osister, NULL);
    adjust_time (proposal->oback, proposal->time);
    obm = showtop (crawl_down (proposal->osister, proposal->time)->back);
    connectnodes (proposal->oback, proposal->origin, proposal->osister);
    proposal->oback->back = NULL;
    if (obb->type == 'r')
        connectnodes (obb, proposal->oback, NULL);
    else
    {
        if (obm->type == 'm')
            connectnodes (obm, proposal->oback, NULL);
        connectnodes (obb, proposal->oback, proposal->ocousin);
    }
    adjust (obb, obb->tyme, 2);
    set_dirty (proposal->oback);
    localevaluate (proposal->oback);
}

void
coat_ocousin (proposal_fmt * proposal)
{
    node *obb, *rcback;
    obb = showtop (crawlback (proposal->oback));
    connectnodes (showtop (proposal->oback->back), proposal->osister, NULL);
    rcback = showtop (crawl_down (proposal->ocousin, proposal->time)->back);
    adjust_time (proposal->oback, proposal->time);
    proposal->oback->back = NULL;
    connectnodes (rcback, proposal->oback, NULL);
    connectnodes (proposal->oback, proposal->ocousin, proposal->origin);
    connectnodes (obb, proposal->oback, proposal->osister);
    adjust (obb, obb->tyme, 2);
    set_dirty (proposal->oback);
    set_dirty (obb);
    if (crawlback (obb)->type != 'r')
        localevaluate (showtop (crawlback (obb)));
}

void
coat_rbcoa (proposal_fmt * proposal)
{
    node *obb, *root;
    if (proposal->ocousin != NULL)
    {
        obb = showtop (crawlback (proposal->oback));
        root = proposal->root;
        connectnodes (showtop (proposal->oback->back), proposal->osister, NULL);
        connectnodes (obb, proposal->ocousin, proposal->osister);
        adjust (obb, obb->tyme, 2);
        adjust_time (proposal->oback, proposal->time);
        proposal->oback->back = NULL;
        connectnodes (proposal->oback, proposal->origin, proposal->target);
        connectnodes (root, proposal->oback, NULL);
        adjust (proposal->oback, proposal->time, 2);
        set_dirty (obb);
        set_dirty (proposal->oback);
        localevaluate (showtop (crawlback (obb)));
    }
    else
    {
        error ("error in coat_rbcoa\n");
    }
}

void
coat_obbcoa (proposal_fmt * proposal)
{
    node /**obb,*/  * proposalack;
    if (proposal->ocousin != NULL)
    {
        /*   obb = showtop(crawlback(proposal->oback)); */
        connectnodes (showtop (proposal->oback->back), proposal->osister, NULL);
        proposalack = showtop (crawlback (proposal->target));
        adjust_time (proposal->oback, proposal->time);
        connectnodes (proposal->target, proposal->ocousin, proposal->osister);
        adjust (proposal->target, proposal->target->tyme, 1);
        proposal->oback->back = NULL;
        adjust_time (proposal->oback, proposal->time);
        if (showtop (proposal->realtarget->back)->type == 'm')
        {
            connectnodes (showtop (proposal->realtarget->back), proposal->oback,
                          NULL);
        }
        connectnodes (proposal->oback, proposal->origin, proposal->target);
        connectnodes (proposalack, proposal->tsister, proposal->oback);
        adjust (proposalack, proposalack->tyme, 2);
        set_dirty (proposal->target);
        set_dirty (proposal->oback);
        localevaluate (proposalack);
    }
    else
    {
        error ("error in coat_obbcoa\n");
    }
}

void
coat_stancoa (proposal_fmt * proposal)
{
    node *obb, *proposalack;
    if (proposal->ocousin != NULL)
    {
        obb = showtop (crawlback (proposal->oback));
        proposalack = showtop (crawlback (proposal->target));
        connectnodes (showtop (proposal->oback->back), proposal->osister, NULL);
        proposal->oback->back = NULL;
        connectnodes (obb, proposal->osister, proposal->ocousin);
        adjust_time (proposal->oback, proposal->time);
        if (showtop (proposal->realtarget->back)->type == 'm')
            connectnodes (showtop (proposal->realtarget->back), proposal->oback,
                          NULL);
        connectnodes (proposal->oback, proposal->origin, proposal->target);
        connectnodes (proposalack, proposal->oback, proposal->tsister);
        adjust (proposalack, proposalack->tyme, 3);
        adjust (obb, obb->tyme, 3);
        set_dirty (obb);
        set_dirty (proposal->oback);
        localevaluate (proposalack);
        localevaluate (obb);
    }
    else
    {
        adjust_time (proposal->oback, proposal->time);
        obb = showtop (crawlback (proposal->oback));
        connectnodes (showtop (proposal->oback->back), proposal->osister, NULL);
        proposalack = showtop (crawlback (proposal->target));
        adjust_time (proposal->oback, proposal->time);
        proposal->oback->back = NULL;
        if (showtop (proposal->realtarget->back)->type == 'm')
            connectnodes (showtop (proposal->realtarget->back), proposal->oback,
                          NULL);
        connectnodes (proposal->oback, proposal->origin, proposal->target);
        connectnodes (proposalack, proposal->oback, proposal->tsister);
        /* either proposalack is osister or a descendent of her
           other case are disallowed by time constraints */
        connectnodes (obb, proposal->osister, NULL);
        adjust (proposalack, proposalack->tyme, 3);
        adjust (obb, obb->tyme, 3);
        set_dirty (proposal->oback);
        localevaluate (proposalack);

    }
}


void
target_oback (proposal_fmt * proposal)
{
    node *obb = showtop (crawlback (proposal->oback));
    copy_x (proposal, proposal->xf, proposal->origin->x);
#ifdef UEP

    if(proposal->world->options->uep)
        copy_uepx(proposal,proposal->uf, proposal->origin->ux);
#endif

    memcpy (proposal->mf, proposal->origin->scale,
            sizeof (MYREAL) * proposal->endsite);
    proposal->v = (proposal->time - proposal->origin->tyme);
    proposal->vs = (proposal->time - proposal->osister->tyme);
    pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
                  proposal->osister->x, proposal->osister->scale, proposal->vs);
    pseudoevaluate (proposal, proposal->xf, proposal->mf, obb, proposal->oback,
                    fabs (proposal->time - obb->tyme));
}

void
target_ocousin (proposal_fmt * proposal)
{
    node *obb;
    obb = showtop (crawlback (proposal->oback));
    copy_x (proposal, proposal->xf, proposal->origin->x);
#ifdef UEP

    if(proposal->world->options->uep)
        copy_uepx(proposal,proposal->uf, proposal->origin->ux);
#endif

    memcpy (proposal->mf, proposal->origin->scale,
            sizeof (MYREAL) * proposal->endsite);
    copy_x (proposal, proposal->xt, proposal->ocousin->x);
#ifdef UEP

    if(proposal->world->options->uep)
        copy_uepx(proposal,proposal->ut, proposal->ocousin->ux);
#endif

    memcpy (proposal->mt, proposal->ocousin->scale,
            sizeof (MYREAL) * proposal->endsite);
    proposal->v = (proposal->time - proposal->origin->tyme);
    proposal->vs = fabs (proposal->ocousin->tyme - proposal->time);
    pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
                  proposal->xt, proposal->mt, proposal->vs);
    proposal->v = fabs (proposal->time - obb->tyme);
    proposal->vs = obb->tyme - proposal->osister->tyme;
    pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
                  proposal->osister->x, proposal->osister->scale, proposal->vs);
    pseudoevaluate (proposal, proposal->xf, proposal->mf,
                    showtop (crawlback (obb)), obb,
                    (showtop (crawlback (obb))->tyme - obb->tyme));
}

void
target_rbcoa (proposal_fmt * proposal)
{
    node *obb;
    if (proposal->ocousin != NULL)
    {
        obb = showtop (crawlback (proposal->oback));
        proposal->vs = obb->tyme - proposal->osister->tyme;
        proposal->v = obb->tyme - proposal->ocousin->tyme;
        copy_x (proposal, proposal->xt, proposal->osister->x);
#ifdef UEP

        if(proposal->world->options->uep)
            copy_uepx(proposal,proposal->ut, proposal->osister->ux);
#endif

        memcpy (proposal->mt, proposal->osister->scale,
                sizeof (MYREAL) * proposal->endsite);
        pseudonuview (proposal, proposal->xt, proposal->mt, proposal->vs,
                      proposal->ocousin->x, proposal->ocousin->scale,
                      proposal->v);
        pseudoevaluate (proposal, proposal->xt, proposal->mt,
                        showtop (crawlback (obb)), obb,
                        showtop (crawlback (obb))->tyme - obb->tyme);
        proposal->v = proposal->time - proposal->origin->tyme;
        proposal->vs = proposal->time - proposal->target->tyme;
        copy_x (proposal, proposal->xf, proposal->origin->x);
#ifdef UEP

        if(proposal->world->options->uep)
            copy_uepx(proposal,proposal->uf, proposal->origin->ux);
#endif

        memcpy (proposal->mf, proposal->origin->scale,
                sizeof (MYREAL) * proposal->endsite);
        pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
                      proposal->xt, proposal->mt, proposal->vs);
    }
    else
    {
        error ("error in target_rbcoa\n");
    }
}

void
target_obbcoa (proposal_fmt * proposal)
{
    node *obb, *obbb;
    if (proposal->ocousin != NULL)
    {
        obb = showtop (crawlback (proposal->oback));
        proposal->vs = obb->tyme - proposal->osister->tyme;
        proposal->v = obb->tyme - proposal->ocousin->tyme;
        copy_x (proposal, proposal->xt, proposal->osister->x);
#ifdef UEP

        if(proposal->world->options->uep)
            copy_uepx(proposal,proposal->ut, proposal->osister->ux);
#endif

        memcpy (proposal->mt, proposal->osister->scale,
                sizeof (MYREAL) * proposal->endsite);
        pseudonuview (proposal, proposal->xt, proposal->mt, proposal->vs,
                      proposal->ocousin->x, proposal->ocousin->scale,
                      proposal->v);
        proposal->v = proposal->time - proposal->origin->tyme;
        proposal->vs = fabs (obb->tyme - proposal->time);
        copy_x (proposal, proposal->xf, proposal->origin->x);
#ifdef UEP

        if(proposal->world->options->uep)
            copy_uepx(proposal,proposal->uf, proposal->origin->ux);
#endif

        memcpy (proposal->mf, proposal->origin->scale,
                sizeof (MYREAL) * proposal->endsite);
        pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
                      proposal->xt, proposal->mt, proposal->vs);
        proposal->v = fabs (showtop (crawlback (obb))->tyme - proposal->time);
        obbb = showtop (crawlback (obb));
        pseudoevaluate (proposal, proposal->xf, proposal->mf, obbb, obb,
                        proposal->v);
    }
    else
    {
        error ("error in target_obbcoa\n");
    }
}

void
target_stancoa (proposal_fmt * proposal)
{
    node *obb, *nn, *oldnn, *d1, *d2, *proposalack, *uppern = NULL; /* o R */
    node **double_line;
    double_line = (node **) mycalloc (2, sizeof (node *));
    if (proposal->ocousin != NULL)
    {
        obb = showtop (crawlback (proposal->oback));
        gotoroot (proposal->target, proposal->line_t);
        double_line[0] = proposal->osister;
        /*findcrossing needs a last NULL element */
        if (findcrossing (double_line, proposal->line_t) != NULL)
        {
            t_s_upper (proposal, proposal->osister, obb);
        }
        else
        {
            double_line[0] = proposal->ocousin;
            /*findcrossing needs a last NULL element */
            if (findcrossing (double_line, proposal->line_t) != NULL)
            {
                t_s_upper (proposal, proposal->ocousin, obb);
            }
            else
            {
                gotoroot (obb, proposal->line_f);
                proposal->connect =
                    findcrossing (proposal->line_f, proposal->line_t);
                proposal->vs = obb->tyme - proposal->osister->tyme;
                proposal->v = obb->tyme - proposal->ocousin->tyme;
                copy_x (proposal, proposal->xt, proposal->osister->x);
#ifdef UEP

                if(proposal->world->options->uep)
                    copy_uepx(proposal,proposal->ut, proposal->osister->ux);
#endif

                memcpy (proposal->mt, proposal->osister->scale,
                        sizeof (MYREAL) * proposal->endsite);
                /*       if (obb!=proposal->connect) { */
                uppern = t_s_obbncon (proposal, obb);
                /*} */
                proposal->v = proposal->time - proposal->origin->tyme;
                copy_x (proposal, proposal->xf, proposal->origin->x);
#ifdef UEP

                if(proposal->world->options->uep)
                    copy_uepx(proposal,proposal->uf, proposal->origin->ux);
#endif

                memcpy (proposal->mf, proposal->origin->scale,
                        sizeof (MYREAL) * proposal->endsite);
                if (proposal->target != proposal->connect)
                {
                    t_s_tncon (proposal, obb);
                }
                else
                {
                    t_s_tcon (proposal, uppern);
                }
            }
        }
    }
    else
    {
        proposal->v = proposal->time - proposal->origin->tyme;
        proposal->vs = proposal->time - proposal->target->tyme;
        copy_x (proposal, proposal->xf, proposal->origin->x);
#ifdef UEP

        if(proposal->world->options->uep)
            copy_uepx(proposal,proposal->uf, proposal->origin->ux);
#endif

        memcpy (proposal->mf, proposal->origin->scale,
                sizeof (MYREAL) * proposal->endsite);
        copy_x (proposal, proposal->xt, proposal->target->x);
#ifdef UEP

        if(proposal->world->options->uep)
            copy_uepx(proposal,proposal->ut, proposal->target->ux);
#endif

        memcpy (proposal->mt, proposal->target->scale,
                sizeof (MYREAL) * proposal->endsite);
        //CHECK HEREH      memcpy(proposal->mt,proposal->osister->scale, sizeof(MYREAL)*proposal->endsite);
        pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
                      proposal->xt, proposal->mt, proposal->vs);
        proposal->v =
            fabs (((proposalack =
                        showtop (crawlback (proposal->target)))->tyme -
                   proposal->time));
        if (proposalack != proposal->oback)
        {
            children (proposalack, &d1, &d2);
            if (d1 == proposal->target)
                d1 = d2;
            pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
                          d1->x, d1->scale, d1->v);
            oldnn = proposalack;
            nn = showtop (crawlback (oldnn));
            while (nn != proposal->oback)
            {
                children (nn, &d1, &d2);
                if (d1 == oldnn)
                    d1 = d2;
                pseudonuview (proposal, proposal->xf, proposal->mf, oldnn->v,
                              d1->x, d1->scale, d1->v);
                oldnn = nn;
                nn = showtop (crawlback (nn));
            }
        }
    }
  myfree(double_line);
}

void
t_s_upper (proposal_fmt * proposal, node * connect, node * obb)
{
    node *nn, *oldnn, *d1, *d2;
    proposal->v = proposal->time - proposal->origin->tyme;
    proposal->vs = proposal->time - proposal->target->tyme;
    copy_x (proposal, proposal->xf, proposal->origin->x);
#ifdef UEP

    if(proposal->world->options->uep)
        copy_uepx(proposal,proposal->uf, proposal->origin->ux);
#endif

    memcpy (proposal->mf, proposal->origin->scale,
            sizeof (MYREAL) * proposal->endsite);
    pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
                  proposal->target->x, proposal->target->scale, proposal->vs);
    nn = showtop (crawlback (proposal->target));
    oldnn = proposal->target;
    proposal->vs = fabs (nn->tyme - proposal->time);
    while (nn != connect)
    {
        children (nn, &d1, &d2);
        if (d1 == oldnn)
            d1 = d2;
        pseudonuview (proposal, proposal->xf, proposal->mf, proposal->vs, d1->x,
                      d1->scale, d1->v);
        proposal->vs = nn->v;
        oldnn = nn;
        nn = showtop (crawlback (nn));
    }
    children (nn, &d1, &d2);
    if (d1 == oldnn)
        d1 = d2;
    pseudonuview (proposal, proposal->xf, proposal->mf, proposal->vs, d1->x,
                  d1->scale, d1->v);
    proposal->v = obb->tyme - proposal->osister->tyme;
    proposal->vs = obb->tyme - proposal->ocousin->tyme;
    if (connect == proposal->ocousin)
        pseudonuview (proposal, proposal->xf, proposal->mf, proposal->vs,
                      proposal->osister->x, proposal->osister->scale,
                      proposal->v);
    else
        pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
                      proposal->ocousin->x, proposal->ocousin->scale,
                      proposal->vs);
    pseudoevaluate (proposal, proposal->xf, proposal->mf,
                    showtop (crawlback (obb)), obb, obb->v);
}

node *
t_s_obbncon (proposal_fmt * proposal, node * obb)
{

    node *nn, *oldnn, *d1, *d2;
    /*     FPRINTF(stdout,"t_s_obbncon\n"); */
    pseudonuview (proposal, proposal->xt, proposal->mt, proposal->vs,
                  proposal->ocousin->x, proposal->ocousin->scale, proposal->v);
    nn = showtop (crawlback (obb));
    oldnn = obb;
    proposal->vs = fabs (nn->tyme - obb->tyme);
    while (nn != proposal->connect)
    {
        children (nn, &d1, &d2);
        if (d1 == oldnn)
            d1 = d2;
        pseudonuview (proposal, proposal->xt, proposal->mt, proposal->vs, d1->x,
                      d1->scale, d1->v);
        proposal->vs = nn->v;
        oldnn = nn;
        nn = showtop (crawlback (nn));
    }
    return oldnn;
}

void
t_s_tncon (proposal_fmt * proposal, node * obb)
{
    node *nn, *oldnn, *d1, *d2;
    pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
                  proposal->target->x, proposal->target->scale,
                  proposal->time - proposal->target->tyme);
    nn = showtop (crawlback (proposal->target));
    oldnn = proposal->target;
    proposal->v = fabs (nn->tyme - proposal->time);
    while (nn != proposal->connect)
    {
        children (nn, &d1, &d2);
        if (d1 == oldnn)
            d1 = d2;
        pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v, d1->x,
                      d1->scale, d1->v);
        proposal->v = nn->v;
        oldnn = nn;
        nn = showtop (crawlback (nn));
    }
    if (obb != proposal->connect)
        pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v,
                      proposal->xt, proposal->mt, proposal->vs);
    proposal->v = proposal->connect->v;
    pseudoevaluate (proposal, proposal->xf, proposal->mf,
                    showtop (crawlback (proposal->connect)), proposal->connect,
                    proposal->v);
}

void
t_s_tcon (proposal_fmt * proposal, node * uppern)
{
    node *obb, *d1, *d2;
    children (proposal->target, &d1, &d2);
    if (d1 == uppern)
        d1 = d2;
    proposal->vs = fabs (proposal->target->tyme - showtop (uppern)->tyme);
    pseudonuview (proposal, proposal->xt, proposal->mt, proposal->vs, d1->x, d1->scale, d1->v); /* eval target */
    proposal->vs = proposal->time - proposal->target->tyme;
    pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v, proposal->xt, proposal->mt, proposal->vs); /*eval oback */
    obb = showtop (crawlback (proposal->target));
    proposal->v = fabs (obb->tyme - proposal->time);
    if (obb->type == 'r')
        return;
    else
    {
        children (obb, &d1, &d2);
        if (d1 == proposal->target)
            d1 = d2;
        pseudonuview (proposal, proposal->xf, proposal->mf, proposal->v, d1->x,
                      d1->scale, d1->v);
        pseudoevaluate (proposal, proposal->xf, proposal->mf,
                        showtop (crawlback (obb)), obb, obb->v);
    }
}


/*-----------------------------------------------------------------
erase all migration nodes between proposal->origin and proposal->oback 
*/
long
erase_migr_nodes (node * up)
{
    long deleted = 0;
    node *theNode, *down, *oldNode;
    down = crawlback (up);
    theNode = up->back;

    while (theNode->type == 'm')
    {

        oldNode = theNode;
        theNode = theNode->next->back;
        if (oldNode->type == 'm')
        {
            myfree(oldNode->next);
            myfree(oldNode);
            deleted++;
        }
    }


    down->back = up;
    up->back = down;
    return deleted;
}

/*-------------------------------------------------------
Connects and adjusts three nodes, the first is the 
mother node and the next two are child nodes, if one
of the child nodes is NULL it just connects mother with
the not-NULL child.
*/
void
connectnodes (node * mother, node * brother, node * sister)
{
    node *tmp;

    if ((sister != NULL) && (brother != NULL))
    {
        if ((mother->id == brother->id) || (mother->id == sister->id)
                || (sister->id == brother->id))
        {
            error ("connectnodes() conflict");
        }
        tmp = crawl_down (brother, mother->tyme);
        mother->next->back = tmp;
        tmp->back = mother->next;
        tmp = crawl_down (sister, mother->tyme);
        mother->next->next->back = tmp;
        tmp->back = mother->next->next;
    }
    else
    {
        if (sister == NULL)
        {
            tmp = crawl_down (brother, mother->tyme);
            mother->next->back = tmp;
            tmp->back = mother->next;
        }
        else
        {

            if (brother == NULL)
            {
                if (mother->type == 'm')
                {
                    tmp = crawl_down (sister, mother->tyme);
                    mother->next->back = tmp;
                    tmp->back = mother->next;
                }
                else
                {
                    tmp = crawl_down (sister, mother->tyme);
                    mother->next->next->back = tmp;
                    tmp->back = mother->next->next;
                }
            }
            else
            {
                error ("Single, lonely rootnode detected in connectenodes()\n");
            }
        }
    }
}

void
gotoroot (node * origin, node ** ptrlist)
{
    node *theNode;
    long i = 0;

    for (theNode = origin;
            (theNode != NULL) && (crawlback (theNode)->type != 'r');
            theNode = showtop (crawlback (theNode)))
    {
        ptrlist[i++] = theNode;
    }
    ptrlist[i] = theNode;  /* adds root->back to the list */
    ptrlist[i+1] = NULL;   /* guarantees an end for the array)*/
}

void
adjust (node * theNode, MYREAL tyme, long level)
{
    if (level < 0 || theNode == NULL)
        return;

    theNode->tyme = tyme;
    if (theNode->type == 'r')
    {
        theNode->v = 1;
        theNode->length = DBL_MAX;
        if (theNode->next->back != NULL)
            adjust (crawlback (theNode->next), crawlback (theNode->next)->tyme,
                    level);
        if (theNode->next->next->back != NULL)
            adjust (crawlback (theNode->next->next),
                    crawlback (theNode->next->next)->tyme, level);
    }
    else if (theNode->type == 't')
    {
#ifndef TESTINGDATE
      theNode->tyme = 0.0;
#endif
        theNode->length = lengthof (theNode);
        ltov (theNode);
        return;
    }
    else if ((theNode->type != 't'))
    {
        if (theNode->type != 'm')
        {
            theNode->length = lengthof (theNode);
            ltov (theNode);
            if (theNode->next->back != NULL)
            {
                theNode->next->tyme = tyme;
                theNode->next->v = theNode->v;
                theNode->next->length = theNode->length;
                adjust (crawlback (theNode->next),
                        crawlback (theNode->next)->tyme, level - 1);
            }
            if (theNode->next->next->back != NULL)
            {
                theNode->next->next->tyme = tyme;
                theNode->next->next->v = theNode->v;
                theNode->next->next->length = theNode->length;
                adjust (crawlback (theNode->next->next),
                        crawlback (theNode->next->next)->tyme, level - 1);
            }
        }
        else
            adjust (crawlback (theNode->next), crawlback (theNode->next)->tyme,
                    level);
    }
}


/* calculates x-array down the tree assuming that only one line
   is affected by the change of the sub-likelihoods above that line
   BUT does NOT calculate the tree-likelihood as evaluate() does.
   THIS CHANGES THE TREE
 */
void
localevaluate (node * mother)
{
    node *nn = NULL;

    if (mother->type != 'r')
    {
        set_dirty (mother);
        for (nn = mother; crawlback (nn)->type != 'r';
                nn = showtop (crawlback (nn)))
        {
            set_dirty (nn);
        }
    }
}


///
/// copy the sequence data
void
copy_x (proposal_fmt * proposal, xarray_fmt xx1, xarray_fmt xx2)
{
    long locus        = proposal->world->locus;
    long endsite      = proposal->endsite;
    long rcategs      = proposal->world->options->rcategs;
    long *maxalleles = proposal->world->data->maxalleles;
    
#ifdef ALTIVEC            
    const size_t sitelikesize = sizeof (FloatVec) * rcategs;
#else
    const size_t sitelikesize = sizeof (sitelike) * rcategs;
#endif    
        long i;
	//        long j;

    switch (proposal->datatype)
    {
    case 'a':
    case 'b':
    case 'm':
        memcpy (xx1.a, xx2.a, sizeof (MYREAL) * maxalleles[locus]);
        break;
    case 's':
    case 'n':
    case 'u':
    case 'f':
        for (i = 0; i < endsite; i++)
        {
            //for (j = 0; j < rcategs; j++)
            //{
                memcpy (xx1.s[i], xx2.s[i], sitelikesize);
            //}
        }
        break;
    }
}

void
fix_root_pop (node * p)
{
    if (crawlback (p) != p->back)
    {
        erase_migr_nodes (p);
    }
    if (p->back->actualpop != p->actualpop)
    {
        p->back->pop = p->back->next->pop = p->back->next->next->pop = p->pop;
        p->back->actualpop = p->back->next->actualpop = p->actualpop;
        p->back->next->next->actualpop = p->actualpop;
    }
}


/* transfers an x-array down the tree assuming that only one line
   is affected by the change of the sub-likelihoods above that line
   BUT does NOT calculate the tree-likelihood as evaluate() does.
   DOES NOT CHANGE THE TREE
 */
void
pseudoevaluate (proposal_fmt * proposal, xarray_fmt x, MYREAL *lx,
                node * mother, node * newdaughter, MYREAL v)
{
    node *nn = NULL, *d1 = NULL, *d2 = NULL, *oldnn = NULL;

    if (mother->type != 'r')
    {
        children (mother, &d1, &d2);
        if (d1 == newdaughter)
            d1 = d2;
        pseudonuview (proposal, x, lx, v, d1->x, d1->scale, d1->v);
        oldnn = mother;
        nn = showtop (crawlback (mother));
        while (nn->type != 'r')
        {
            children (nn, &d1, &d2);
            if (d1 == oldnn)
                d1 = d2;
            pseudonuview (proposal, x, lx, oldnn->v, d1->x, d1->scale, d1->v);
            oldnn = nn;
            nn = showtop (crawlback (nn));
        }
    }
}

node *
findcrossing (node ** ptrl1, node ** ptrl2)
{
    long i = 0, j = 0;

    /* assumes that there is an NULL element at the end */
    for (i = 0; ptrl1[i] != NULL; j = 0, i++)
    {
        while ((ptrl2[j] != NULL) && (ptrl1[i] != ptrl2[j]))
            j++;
        if (ptrl2[j] != NULL)
        {
            break;
        }
    }
    return ptrl1[i];
}

node *
crawl_down (node * theNode, MYREAL tyme)
{
    node *otmp, *tmp = theNode->back;

    otmp = theNode;
    if (tmp == NULL)
        return otmp;
    while ((tmp->type == 'm') && showtop (tmp)->tyme < tyme)
    {
        otmp = tmp->next;
        tmp = tmp->next->back;
        if (tmp == NULL)
            return otmp;
    }
    return otmp;
}