#include <vector>
#include <iostream>
#include <fstream>
#include <iomanip>
using namespace std;

#ifndef HETEROHOT_H
#define HETEROHOT_H

// GENERAL
double normrv(double mn, double var) {
  int flag;
  double u,y,x;

  flag = 1;
  while (flag) {
    u = drand48();
    y = -log(drand48());
    if ( u <= exp(-0.5*(y-1)*(y-1)) ) {
      x = y;
      flag = 0;
    }
  }
  if (drand48() < 0.5) { x = -x; }
  x = mn + sqrt(var)*x;

  return x;
}

template<class T>
T min2(T x, T y) {
  T z; 
  z = x;
  if (z > y) { z = y; }
  return z;
}

template<class T>
T max2(T x, T y) {
  T z; 
  z = x;
  if (z < y) { z = y; }
  return z;
}

int poisson(double x) {
  double temp;
  int flag, cnt;
  temp = 1;
  flag = 1;
  cnt = 0;
  while (flag) {
    cnt += 1;
    temp = temp * drand48();
    if (temp < exp(-x)) {
      flag = 0;
      cnt -= 1;
    }
  }
  return cnt;
}

// CLASS INTERVAL
class interval {
  friend class node;
  friend class vecinterval;
  friend class arg;
 public:
  interval(double = 0, double = 0);
  int ifoverlapunion(interval);
  const interval &operator=(const interval &);
  void print();
 private:
  double lend;
  double rend;
};

interval::interval(double le, double re) {
  lend = le;
  rend = re;
}

int interval::ifoverlapunion(interval i1) {
  int flag;
  double tlend, trend;

  tlend = lend;
  trend = rend;
  if (i1.lend < lend) { tlend = i1.lend; }
  if (i1.rend > rend) { trend = i1.rend; }
  flag = 0;
  if (i1.lend <= lend) {
    if (i1.rend >= lend) { flag = 1; }
  }
  else {
    if (i1.lend <= rend) { flag = 1; }
  }
  if (flag) {
    if (tlend == trend) { flag = 0; }
  }
  if (flag) {
    lend = tlend;
    rend = trend;
  }
  return flag;
}

const interval &interval::operator=(const interval &right) {
  if (&right != this) {
    lend = right.lend;
    rend = right.rend;
  }
  return *this;
}

void interval::print() {
  cout << lend << " " << rend << " ";
}

// CLASS VECINTERVAL
class vecinterval {
  friend class node;
  friend class arg;
 public:
  vecinterval();
  ~vecinterval();
  int there();
  void onion(const vecinterval &,const vecinterval &);
  int intersect(const vecinterval &,const vecinterval &);
  int complement(const vecinterval &);
  int pointin(double);
  void setfirstinterval(interval);
  void addto(interval);
  const vecinterval &operator=(const vecinterval &);
  void print();
 private:
  vector<interval> vi;
};

vecinterval::vecinterval() {
  vi.clear();
}

vecinterval::~vecinterval() {
  vi.clear();
}

int vecinterval::there() {
  int i, flag;
  flag = 0;
  for (i=0;i<vi.size();i++) {
    if (vi[i].lend < vi[i].rend) { flag = 1; break; }
  }

  return flag;
}

int vecinterval::intersect(const vecinterval &v1,const vecinterval &v2) {
  int i,j,flag;
  double tleft,tright;
  interval itemp;

  vi.clear();
  for (i=0;i<v2.vi.size();i++) {
    for (j=0;j<v1.vi.size();j++) {
      tleft  = max2(v1.vi[j].lend,v2.vi[i].lend);
      tright = min2(v1.vi[j].rend,v2.vi[i].rend);
      if (tleft < tright) {
	itemp.lend = tleft;
	itemp.rend = tright;
	vi.push_back(itemp);
      }
    }
  }

  flag = 1;
  if (vi.size() == 0) { 
    flag = 0; 
    vi.push_back(itemp);
  }
  return flag;
}

int vecinterval::complement(const vecinterval &v1) {
  vecinterval nvi;
  int i,j,flag,lo,hi;
  double lborder,rborder;
  interval itemp;

  for (i=0;i<vi.size();i++) {
    lborder = vi[i].lend;
    rborder = vi[i].rend;
    flag = 0; lo = -9;
    for (j=0;j<v1.vi.size();j++) {
      if ((lborder <= v1.vi[j].lend) && (v1.vi[j].rend <= rborder)) { 
	flag = 1; 
	if (lo == -9) { lo = j; }
	hi = j;
      }
    }
    if (flag == 1) {
      if (lborder < v1.vi[lo].lend) {
	itemp.lend = lborder;
	itemp.rend = v1.vi[lo].lend;
	nvi.vi.push_back(itemp);
      }
      for (j=lo;j<hi;j++) {
	itemp.lend = v1.vi[j].rend;
	itemp.rend = v1.vi[(j+1)].lend;
	nvi.vi.push_back(itemp);
      }
      if (v1.vi[hi].rend < rborder) {
	itemp.lend = v1.vi[hi].rend;
	itemp.rend = rborder;
	nvi.vi.push_back(itemp);
      }
    }
    else { nvi.vi.push_back(vi[i]); }
  }

  vi.clear();
  flag = 1;
  if (nvi.vi.size() == 0) { 
    flag = 0; 
    nvi.vi.push_back(itemp); 
  }
  for (i=0;i<nvi.vi.size();i++) {
    vi.push_back(nvi.vi[i]);
  }

  return flag;
}

const vecinterval &vecinterval::operator=(const vecinterval &right) {
  int i;
  if (&right != this) {
    vi.clear();
    for (i=0;i<right.vi.size();i++) {
      vi.push_back(right.vi[i]);
    }
  }
  return *this;
}

int vecinterval::pointin(double x) {
	int i,flag;
	flag = 0;
	for (i=0;i<vi.size();i++) {
		if ((vi[i].lend<x) && (x<=vi[i].rend)) { flag = 1; break; }
	}
	return flag;
}

void vecinterval::setfirstinterval(interval x) {
  vi.clear();
  vi.push_back(x);
}

void vecinterval::addto(interval x) {
  vi.push_back(x);
}

void vecinterval::onion(const vecinterval &v1, const vecinterval &v2) {
  vector<int> o1,o2;
  int tot,cnt,flag,temp,i,iflag;
  double tflag;
  interval tinterval;

  vi.clear();
  tot = v1.vi.size() + v2.vi.size();
  for (i=0;i<v1.vi.size();i++) {
    if (v1.vi[i].lend == v1.vi[i].rend) { o1.push_back(1); tot -= 1; }
    else { o1.push_back(0); }
  }
  for (i=0;i<v2.vi.size();i++) {
    if (v2.vi[i].lend == v2.vi[i].rend) { o2.push_back(1); tot -= 1; }
    else { o2.push_back(0); }
  }
  cnt = -1;

  while (tot > 0) {
    flag = 0;
    for (i=0;i<v1.vi.size();i++) {
      if (o1[i] == 0) {
	if (flag == 0) {
	  flag = 1;
	  iflag = i;
	  tflag = v1.vi[i].lend;
	}
	else {
	  if (v1.vi[i].lend < tflag) {
	    flag = 1;
	    iflag = i;
	    tflag = v1.vi[i].lend;
	  }
	}
      }
    }
    for (i=0;i<v2.vi.size();i++) {
      if (o2[i] == 0) {
	if (flag == 0) {
	  flag = 2;
	  iflag = i;
	  tflag = v2.vi[i].lend;
	}
	else {
	  if (v2.vi[i].lend < tflag) {
	    flag = 2;
	    iflag = i;
	    tflag = v2.vi[i].lend;
	  }
	}
      }
    }
    if (flag == 1) {
      o1[iflag] = 1;
      tot -= 1;
      cnt += 1;
      tinterval.lend = v1.vi[iflag].lend;
      tinterval.rend = v1.vi[iflag].rend;
      vi.push_back(tinterval);
    }
    if (flag == 2) {
      o2[iflag] = 1;
      tot -= 1;
      cnt += 1;
      tinterval.lend = v2.vi[iflag].lend;
      tinterval.rend = v2.vi[iflag].rend;
      vi.push_back(tinterval);
    }
    while (flag > 0) {
      flag = 0;
      for (i=0;i<v1.vi.size();i++) {
	if (o1[i] == 0) {
	  temp = vi[cnt].ifoverlapunion(v1.vi[i]);
	  if (temp == 1) { o1[i] = 1; tot -= 1; }
	  flag += temp;
	}
      }
      for (i=0;i<v2.vi.size();i++) {
	if (o2[i] == 0) {
	  temp = vi[cnt].ifoverlapunion(v2.vi[i]);
	  if (temp == 1) { o2[i] = 1; tot -= 1; }
	  flag += temp;
	}
      }
    }
  }

  if (vi.size() == 0) {
    tinterval.lend = 0;
    tinterval.rend = 0;
    vi.push_back(tinterval);
  }
}

void vecinterval::print() {
  int i;
  cout << "SIZE: " << vi.size() << endl;
  for (i=0;i<vi.size();i++) {
    vi[i].print();
  }
  cout << endl;
}

// CLASS NODE
class node {
  friend class arg;
 public:
  node();
  ~node();
  void setss(int);
  const node &operator=(const node &);
  void print();
  void errorcatch(int);
 private:
  int thru;
  node *replace;
  vector<vecinterval> origsampleregions;
  double time;
  int hot;
  vector<double> snps;
  vecinterval regions;
  node *child1;
  node *child2;
  node *parent1;
  node *parent2;
  vecinterval tochild1;
  vecinterval tochild2;
};

void node::errorcatch(int type) {
  int flag, flag2, k, ecatch;
  vecinterval v1;

    ecatch = 0;
    if (regions.there() ==0) { cerr << "PROBLEM X " << type << endl; ecatch = 1; }
    flag = 0;
    flag2 = 0;
    for (k=0;k<origsampleregions.size();k++) {
      if (origsampleregions[k].there() == 1) { flag = 1; }
      v1 = origsampleregions[k];
    }
    if (flag == 0) { cerr << "PROBLEM Y " << type << endl; ecatch = 1; }
    if (ecatch == 1) { exit(1); }
}

void node::setss(int ss) {
  int i;
  interval itemp(0,0);
  vecinterval vitemp;
  vitemp.setfirstinterval(itemp);
  origsampleregions.clear();
  for (i=0;i<ss;i++) {
    origsampleregions.push_back(vitemp);
  }
}

node::node() {
  int i;
  interval itemp(0,0);
  vecinterval vitemp;
  vitemp.setfirstinterval(itemp);

  origsampleregions.clear();
  thru = 0;
  replace = 0;
  time = 0;
  hot = 0;
  snps.clear();
  regions = vitemp;
  child1 = 0;
  child2 = 0;
  parent1 = 0;
  parent2 = 0;
  tochild1 = vitemp;
  tochild2 = vitemp;
}

node::~node() {
  origsampleregions.clear();
  snps.clear();
}

const node &node::operator=(const node &right) {
  int k;
  if (&right != this) {
    for (k=0;k<origsampleregions.size();k++) {
      origsampleregions[k] = right.origsampleregions[k];
    }
    thru = right.thru;
    replace = right.replace;
    time = right.time;
    hot = right.hot;
    snps = right.snps;
    regions = right.regions;
    child1 = right.child1;
    child2 = right.child2;
    parent1 = right.parent1;
    parent2 = right.parent2;
    tochild1 = right.tochild1;
    tochild2 = right.tochild2;
  }
  return *this;
}

void node::print() {
  int i;
  cout << "ORIGSAMPLE:" << endl;
  for (i=0;i<origsampleregions.size();i++) {
    if (origsampleregions[i].there()) {
	cout << i << " ";
	origsampleregions[i].print();
    }
  }
  cout << "TIME: " << time << endl;
  cout << "HOT:  " << hot  << endl;
  cout << "SNPs: ";
  for (i=0;i<snps.size();i++) {
    cout << snps[i] << " ";
  }
  cout << endl;
  cout << "REGIONS:" << endl;
  regions.print();
  cout << "CHILD1:  " << child1 << endl;
  cout << "CHILD2:  " << child2 << endl;
  cout << "PARENT1: " << parent1 << endl;
  cout << "PARENT2: " << parent2 << endl;
  cout << "TO CHILD1:" << endl;
  tochild1.print();
  cout << "TO CHILD2:" << endl;
  tochild2.print();
  cout << "THRU: " << thru << endl;
  cout << "REPLACE: " << replace << endl;
}

// CLASS ARG
class arg {
 public:
  arg(int argc,char *argv[]);
  void print();
  void print2();
  void print3(int);
  void print4();
 private:
  void getpars(int argc, char *argv[]);
  int rseed;
  int catch2();
  double coarse; 
  int samplesize; // HAPLOIDS
  double popsize; // DIPLOIDS
  double inithotspotfreq; // FREQ
  int initsamplefreq;     // COUNT
  double mutrate;
  double recrate; // PER COPY
  double hotrate;
  double pcrossover;
  double lochew; // ASSUME SYMMETRIC
  double hichew;
  double seqpos;
  double alpha;
  double getalpha();
  vector<node *> pops;
  vector<node *> sample;
  double time; // TEMPORARY: next 4
  double hsf;
  vector<node *> currhot;
  vector<node *> currcold;
  vector<node *> bigone;
  void eventprob(vector<double> &);
  void event(vector<double>);
  void coalesce(int);
  void recombine(int);
  void mutate();
  void mutatehelper(node *);
  void mutatework(node *,node *,vecinterval);
  double hotspotfreq(double);
  double consthotspotfreq(double);
  vector<vecinterval> origsampleregions;
};

double arg::getalpha() {
  double x,f;

  f = (seqpos - 0.5 - lochew) / (hichew - lochew);
  if (f > 1) { f = 0; }
  else {
    if (f <= 0) { f = 1; }
    else { f = 1 - f; }
  }
  x = 2*((double)(popsize))*hotrate*f;

  return x;
}

int arg::catch2() {
  int i,k,flag;
  vecinterval v1,v2;
  interval itemp;

  flag = 0;
  for (k=0;k<samplesize;k++) {
    v1.setfirstinterval(itemp);
    for (i=0;i<pops.size();i++) {
      v2 = v1;
      v1.onion(v2,pops[i]->origsampleregions[k]);
    }
    for (i=0;i<currcold.size();i++) {
      v2 = v1;
      v1.onion(v2,currcold[i]->origsampleregions[k]);
    }
    for (i=0;i<currhot.size();i++) {
      v2 = v1;
      v1.onion(v2,currhot[i]->origsampleregions[k]);
    }
  
    if (v1.vi.size()==1) {
      if ((v1.vi[0].lend == 0) && (v1.vi[0].rend == 1)) {
	flag += 1;
      }
      else { cout << "NOO1 " << k << endl; v1.print(); }
    }
    else { cout << "NOO2 " << k << endl; v1.print(); }
  }

  if (flag != samplesize) { flag = 1; }
  else {flag = 0; }

  return flag;
}

void arg::getpars(int argc, char *argv[]) {
  
  if (argc < 12) {
    cerr << "PROBLEM IN GETPARS" << endl;
    exit(1);
  }

  rseed = atoi(argv[1]);
  srand48(rseed);
  samplesize = atoi(argv[2]);
  popsize = atoi(argv[3]);
  inithotspotfreq = atof(argv[4]);
  initsamplefreq = atoi(argv[5]);
  mutrate = atof(argv[6]);
  recrate = atof(argv[7]);
  hotrate = atof(argv[8]);
  pcrossover = atof(argv[9]);
  lochew = atof(argv[10])/coarse;
  hichew = atof(argv[11])/coarse;
  seqpos = atof(argv[12])/coarse;

}

arg::arg(int argc,char *argv[]) {
  int i,j,k,zz;
  interval itemp(0,1);
  interval ztemp(0,0);
  vecinterval vitemp,zvitemp;
  vector<double> ep(5);

  origsampleregions.clear();
  pops.clear();

  coarse = 100000;
  getpars(argc,argv);

  alpha = getalpha();
  hsf = inithotspotfreq;
  time = 0;
  node *ntemp = new node[samplesize];
  sample.clear();

  // SET-UP
  vitemp.setfirstinterval(itemp);
  zvitemp.setfirstinterval(ztemp);
  for (i=0;i<samplesize;i++) {
    ntemp[i].setss(samplesize);
    ntemp[i].origsampleregions[i] = vitemp;
    ntemp[i].regions = vitemp;   
    if (i < initsamplefreq) { ntemp[i].hot = 1; }
    else { ntemp[i].hot = 0; }
    sample.push_back(&ntemp[i]);
    if (sample[i]->hot == 1) { currhot.push_back(sample[i]); }
    else { currcold.push_back(sample[i]); }
    bigone.push_back(sample[i]);
  }

  // LIFE
  while ((currhot.size()+currcold.size()) > 1) {
    hsf = hotspotfreq(hsf); 
    time += 1;
    eventprob(ep);
    if (drand48() < ep[1]) {
      recombine(0);
    }
    if ((drand48() < ep[2]) && (currhot.size() > 0)) {
      recombine(1);
    }
    if ((drand48() < ep[3]) && (currcold.size() > 1)) {
      coalesce(0);
    }
    if ((drand48() < ep[4]) && (currhot.size() > 1)) {
      coalesce(1);
    }

  }
  mutate();
}

void arg::eventprob(vector<double> &ep) {
  double coldrec,hotrec,hotcoal,coldcoal;

  coldrec = ((double)(currcold.size())) * (2*recrate + hsf*hotrate);
  hotrec = ((double)(currhot.size())) * (2*recrate + (1+hsf)*hotrate);
  if (hsf == 0) {
    if (currhot.size() > 1) { hotcoal = 1; }
    else { hotcoal = 0; }
  }
  else {
    hotcoal = 0.25 * ((double)(currhot.size())) * (((double)(currhot.size()))-1) / (hsf*popsize);
  }
  if (hsf == 1) {
    if (currcold.size() > 1) { coldcoal = 1; }
    else { coldcoal = 0; }
  }
  else {
    coldcoal = 0.25 * ((double)(currcold.size())) * (((double)(currcold.size()))-1) / ((1-hsf)*popsize);
  }
  if ((hotcoal == 1)&&(coldcoal == 1)) { cerr << "ONES AWAY" << endl; exit(1);}
  if (hotcoal == 1) { coldrec = 0; hotrec = 0; coldcoal = 0; }
  if (coldcoal == 1) { coldrec = 0; hotrec = 0; hotcoal = 0; }
  ep[1] = coldrec;
  ep[2] = hotrec;
  ep[3] = coldcoal;
  ep[4] = hotcoal;
  ep[0] = ep[1]+ep[2]+ep[3]+ep[4];
}

void arg::event(vector<double> ep) {
  double r;
  r = drand48();
  if (r < (ep[1]/ep[0])) { // COLDREC
    recombine(0);
  }
  else {
    if (r < ((ep[1]+ep[2])/ep[0])) { // HOTREC
      recombine(1);
    }
    else {
      if (r < ((ep[1]+ep[2]+ep[3])/ep[0])) { // COLDCOAL
	coalesce(0);
      }
      else { // HOTCOAL
	coalesce(1);
      }
    }
  }
}

void arg::coalesce(int hotpop) {
  int i,j,k,mi,mj,flag,flag2,othsize;
  node *ca = new node; ca->setss(samplesize);
  node *ca2 = new node; ca2->setss(samplesize);
  vecinterval v1; vecinterval v2;

  vector<node *> curr;
  if (hotpop == 1) { curr = currhot; othsize = currcold.size(); }
  else { curr = currcold; othsize = currhot.size(); }
  
  i = ((int)(floor(drand48()*((double)(curr.size())))));
  flag = 1;
  while (flag) {
    j = ((int)(floor(drand48()*((double)(curr.size())))));
    if (i != j) { flag = 0; }
  }

  ca->time = time;
  ca->hot = curr[0]->hot;
  ca->regions.onion(curr[i]->regions,curr[j]->regions);
  ca->child1 = curr[i];
  ca->child2 = curr[j];
  curr[i]->parent1 = ca;
  curr[j]->parent1 = ca;
  ca->tochild1 = curr[i]->regions;   
  ca->tochild2 = curr[j]->regions;   

  //ORIGSAMPLEREGIONSHERE
  for (k=0;k<samplesize;k++) {
    ca->origsampleregions[k].onion(curr[i]->origsampleregions[k],curr[j]->origsampleregions[k]);
  }

  k = curr.size() - 1;
  mi = max2(i,j);
  mj = min2(i,j);
  if (mi == k) { curr[mj] = curr[(k-1)]; }
  else {
	  if (mi == (k-1)) { curr[mj] = curr[k]; }
	  else {
		curr[mi] = curr[k];
		curr[mj] = curr[(k-1)];
	  }
  }
  curr.pop_back();
  curr.pop_back(); 
 
  v1 = ca->origsampleregions[0];
  for (k=1;k<samplesize;k++) {
    flag = 1;
    flag = v2.intersect(v1,ca->origsampleregions[k]);
    v1 = v2;
    if (flag == 0) { break; }
  }
  if (flag == 1) {
    *ca2 = *ca;
    flag2 = 0;
    ca->regions.complement(v1);
    ca2->regions = v1;
    for (k=0;k<samplesize;k++) {
      flag2 += ca->origsampleregions[k].complement(v1);
      ca2->origsampleregions[k] = v1;
    }
    ca2->tochild1 = v1;
    ca2->tochild2 = v1;
    ca->tochild1.complement(v1);
    ca->tochild2.complement(v1);

    ca->replace = ca2;

    bigone.push_back(ca2);
    pops.push_back(ca2);

    if (flag2 > 0) { 
      bigone.push_back(ca);
      curr.push_back(ca); 
    }
    else { 
      delete ca;
    }
  }
  else {
    bigone.push_back(ca);
    curr.push_back(ca);

    delete ca2;
  }
  
  if (hotpop == 1) { currhot = curr; }
  else { currcold = curr; }

}

void arg::recombine(int hotpop) {
  int i,j,othhot,ours,crossover;
  double r,tot,pos,poslo,poshi;
  node *ca1 = new node; ca1->setss(samplesize);
  node *ca2 = new node; ca2->setss(samplesize);
  interval temp;
  vecinterval vi1,vi2;
  vector<node *>::iterator p;
  vector<node *> curr1; vector<node *> curr2;

  if (hotpop == 1) { curr1 = currhot; curr2 = currcold; }
  else { curr1 = currcold; curr2 = currhot; }

  i = ((int)(floor(drand48()*((double)(curr1.size())))));
  if (drand48() < hsf) { othhot = 1; }
  else { othhot = 0; }
  r = drand48();
  tot = 2*recrate + (curr1[0]->hot+othhot)*hotrate;
  if (r < (recrate/tot)) { //UNIFORM, OURS
    pos = drand48();
    ours = 1;
  }
  else {
    if (r < (2*recrate/tot)) { //UNIFORM, OTHER
      pos = drand48();
      ours = 0;
    }
	else {
	    if (r < ((2*recrate+othhot)/tot)) { //HOT, OTHER
		  pos = 0.5;
		  ours = 0;
		}
		else { //HOT, OURS
			pos = 0.5;
			ours = 1;
		}
	}
  }
  r = drand48();
  poslo = pos - (lochew + r*(hichew-lochew));
  poshi = pos + lochew + r*(hichew-lochew);
  if (drand48() < pcrossover) { crossover = 1; }
  else { crossover = 0; }

  if ((ours>0) || (crossover>0)) { 
	ca1->time = time;
	ca2->time = time;
	ca1->child1 = curr1[i];
	ca2->child1 = curr1[i];
	if (crossover == 0) {
		temp.lend = 0;
		temp.rend = poslo;
		vi1.setfirstinterval(temp);
		temp.lend = poshi;
		temp.rend = 1;
		vi1.addto(temp);
		temp.lend = poslo;
		temp.rend = poshi;
		vi2.setfirstinterval(temp);
	}
	else {
		if (ours == 1) {
			if (drand48() < 0.5) {
				temp.lend = 0;
				temp.rend = poslo;
				vi1.setfirstinterval(temp);
				temp.lend = poslo;
				temp.rend = 1;
				vi2.setfirstinterval(temp);
			}
			else {
				temp.lend = poshi;
				temp.rend = 1;
				vi1.setfirstinterval(temp);
				temp.lend = 0;
				temp.rend = poshi;
				vi2.setfirstinterval(temp);
			}
		}
		else {
			if (drand48() < 0.5) {
				temp.lend = 0;
				temp.rend = poshi;
				vi1.setfirstinterval(temp);
				temp.lend = poshi;
				temp.rend = 1;
				vi2.setfirstinterval(temp);
			}
			else {
				temp.lend = poslo;
				temp.rend = 1;
				vi1.setfirstinterval(temp);
				temp.lend = 0;
				temp.rend = poslo;
				vi2.setfirstinterval(temp);
			}
		}
	}
	ca1->regions.intersect(vi1,curr1[i]->regions);
	ca2->regions.intersect(vi2,curr1[i]->regions);
	ca1->tochild1 = ca1->regions;
	ca2->tochild1 = ca2->regions;
	
	//ORIGSAMPLEREGIONS
	for (j=0;j<samplesize;j++) {
	  ca1->origsampleregions[j].intersect(ca1->regions,curr1[i]->origsampleregions[j]);
	  ca2->origsampleregions[j].intersect(ca2->regions,curr1[i]->origsampleregions[j]);
	}

	// DETERMINE HOT
	if (vi1.pointin(seqpos)) { ca1->hot = curr1[i]->hot; }
	else { ca1->hot = othhot; }
	if (vi2.pointin(seqpos)) { ca2->hot = othhot; }
	else { ca2-> hot = curr1[i]->hot; }

	if ((ca1->regions.there() > 0) && (ca2->regions.there() > 0)) {
	  curr1[i]->parent1 = ca1;
	  curr1[i]->parent2 = ca2;
	  if (ca1->hot == curr1[i]->hot) { curr1.push_back(ca1); }
	  else { curr2.push_back(ca1); }
	  if (ca2->hot == curr1[i]->hot) { curr1.push_back(ca2); }
	  else { curr2.push_back(ca2); }
	  p = curr1.begin();
	  for (j=0;j<i;j++) { p++; }
	  curr1.erase(p); 
	  bigone.push_back(ca1);
	  bigone.push_back(ca2);
	}
	else {
	  if (ca1->regions.there()) { 
	    if (curr1[i]->hot != ca1->hot) {
	      curr1[i]->hot = ca1->hot;
	      curr2.push_back(curr1[i]);
	      p = curr1.begin();
	      for (j=0;j<i;j++) { p++; }
	      curr1.erase(p);
	    }     
	  }
	  if (ca2->regions.there()) { 
	    if (curr1[i]->hot != ca2->hot) {
	      curr1[i]->hot = ca2->hot;
	      curr2.push_back(curr1[i]);
	      p = curr1.begin();
	      for (j=0;j<i;j++) { p++; }
	      curr1.erase(p);
	    } 
	  }
	  delete ca1;
	  delete ca2;
	}
  }
  else {
    delete ca1;
    delete ca2;
  }

  if (hotpop == 1) { currhot = curr1; currcold = curr2; }
  else { currcold = curr1; currhot = curr2; }

}

void arg::mutatework(node *x,node *child,vecinterval tochild) {
  vector<double>::iterator p;
  int i,j,k,flag,no;
  double pos;

  // COPY EXISTING SNPs
  for (i=0;i<x->snps.size();i++) {
    pos = x->snps[i];
    flag = 0;
    if ( (tochild.pointin(pos)) ) { flag = 1; }
    if (flag == 1) {
      if (child->snps.empty()) {
	child->snps.push_back(pos);
      }
      else {
	p = child->snps.begin();
	for (j=0;j<child->snps.size();j++) {
	  if (pos < *p) {
	    child->snps.insert(p,pos);
	    break;
	  }
	  if (j == (child->snps.size() - 1)) {
	    child->snps.push_back(pos);
	    break;
	  }
	  p++;
	}
      }
    }
  }

  // POSSIBLY ADD NEW SNPs
  no = poisson((x->time - child->time)*mutrate);
  for (i=1;i<=no;i++) {
    pos = drand48();
    flag = 0;
    if ( (tochild.pointin(pos)) ) { flag = 1; }
    if (flag == 1) {
      if (child->snps.empty()) {
	child->snps.push_back(pos);
      }
      else {
	p = child->snps.begin();
	for (j=0;j<child->snps.size();j++) {
	  if (pos < *p) {
	    child->snps.insert(p,pos);
	    break;
	  }
	  if (j == (child->snps.size() - 1)) {
	    child->snps.push_back(pos);
	    break;
	  }
	  p++;
	}
      }
    }
  }
}

void arg::mutatehelper(node *x) {
  int flag, k;

  if (x != 0) {
    if (x->replace != 0) {
      mutatework(x,x->child1,x->tochild1);
      mutatework(x,x->child2,x->tochild2);
      x->replace->thru += 1;
      mutatehelper(x->replace);
    }
    else {
      if (x->child1 != 0) {
	mutatework(x,x->child1,x->tochild1);
	x->child1->thru += 1;
	if (x->child1->parent2 != 0) {
	  if (x->child1->thru > 2) { cerr << "WHOOA 1" << endl; exit(1); }
	  if (x->child1->thru == 2) { mutatehelper(x->child1); }
	}
	else { 
	  if (x->child1->thru > 1) { cerr << "WHOOA 2" << endl; exit(1); }
	  mutatehelper(x->child1); 
	}
      }

      if (x->child2 != 0) {
	mutatework(x,x->child2,x->tochild2);
	x->child2->thru += 1;
	if (x->child2->parent2 != 0) {
	  if (x->child2->thru > 2) { cerr << "WHOOA 3" << endl; exit(1); }
	  if (x->child2->thru == 2) { mutatehelper(x->child2); }
	}
	else { 
	  if (x->child2->thru > 1) { cerr << "WHOOA 4" << endl; exit(1); }
	  mutatehelper(x->child2); 
	}
      }
    }
  }
}

void arg::mutate() {
  int i;
  for (i=(pops.size()-1);i>=0;i--) {
    if (pops[i]->thru == 0) {
      //cerr << i << endl;
      mutatehelper(pops[i]);
    }
  }
}

double arg::consthotspotfreq(double x) {
  return x;
}

double arg::hotspotfreq(double x) {
  double s2,mu,var,xp,timeint;

  if (x == 1) { xp = 1; }
  else {
    if (x == 0) { xp = 0; }
    else {
      s2 = x*(1-x);
      timeint = 0.5/((double)(popsize));
  
      mu = (-alpha*s2 + s2*2*alpha/(1-exp(2*alpha*(1-x)))) * timeint;
      var = s2*normrv(0,timeint);

      xp = x + mu + var;
      if (xp > 1) { xp = 1; }
      if (xp < 0) { xp = 0; }
    }
  }

  return xp;
}


void arg::print2() {
vector<double>::iterator p;
  vector<double> totsnps;
  int i,j,k,flag;
  double pos;

  for (i=0;i<sample.size();i++) {
    for (k=0;k<sample[i]->snps.size();k++) {
      pos = sample[i]->snps[k];
      if (totsnps.empty()) {
	totsnps.push_back(pos);
      }
      else {
	flag = 0;
	p = totsnps.begin();
	while (p != totsnps.end()) {
	  if (pos == *p) { 
	    flag = 1;
	    break;
	  }
	  if (pos < *p) {
		flag = 1;
		totsnps.insert(p,pos);
	    break;
	  }
	  p++;
	}
	if (flag == 0) {
	    totsnps.push_back(pos);
	}
	  }
	}
  }

  cout << ((int)(((double)(samplesize))/2)) << endl;
  cout << totsnps.size() << endl;
  cout << "P ";
  for (i=0;i<totsnps.size();i++) {
    cout << ((int)(floor(100000*totsnps[i]))) << " ";
  }
  cout << endl;
  for (i=0;i<totsnps.size();i++) {
    cout << "S";
  }
  cout << endl;
  for (i=0;i<sample.size();i+=2) {
    cout << "#" << i << endl;
    for (j=0;j<totsnps.size();j++) {
      flag = 0;
      for (k=0;k<sample[i]->snps.size();k++) {
	pos = sample[i]->snps[k];
	if (pos == totsnps[j]) {
	  flag = 1;
	  break;
	}
      }
      if (flag == 1) { cout << "1"; }
      else { cout << "0"; }
    }
    cout << endl;
    for (j=0;j<totsnps.size();j++) {
      flag = 0;
      for (k=0;k<sample[i+1]->snps.size();k++) {
	pos = sample[i+1]->snps[k];
	if (pos == totsnps[j]) {
	  flag = 1;
	  break;
	}
      }
      if (flag == 1) { cout << "1"; }
      else { cout << "0"; }
    }
    cout << endl;
    cout << endl;
  }
}


void arg::print3(int numshow) {
vector<double>::iterator p;
  vector<double> totsnps;
  int i,j,k,flag,xdiv,xlo,xhi;
  double pos;

  for (i=0;i<sample.size();i++) {
    for (k=0;k<sample[i]->snps.size();k++) {
      pos = sample[i]->snps[k];
      if (totsnps.empty()) {
	totsnps.push_back(pos);
      }
      else {
	flag = 0;
	p = totsnps.begin();
	while (p != totsnps.end()) {
	  if (pos == *p) { 
	    flag = 1;
	    break;
	  }
	  if (pos < *p) {
		flag = 1;
		totsnps.insert(p,pos);
	    break;
	  }
	  p++;
	}
	if (flag == 0) {
	    totsnps.push_back(pos);
	}
	  }
	}
  }

  xdiv = -9;
  for (i=0;i<totsnps.size();i++) {
    if (totsnps[i] > 0.5) { 
      xdiv = i; 
      break;
    }
  }
  if (xdiv == -9) {
    cerr << "PROBLEM IN PRINT3" << endl;
    exit(1);
  }
  xlo = xdiv - numshow; if (xlo < 0) { xlo = 0; }
  xhi = xdiv + numshow; if (xlo > totsnps.size()) { xhi = totsnps.size(); }
  xdiv = xhi - xlo;


  cout << ((int)(((double)(samplesize))/2)) << endl;
  cout << xdiv << endl;
  cout << "P ";
  for (i=xlo;i<xhi;i++) {
    cout << ((int)(floor(100000*totsnps[i]))) << " ";
  }
  cout << endl;
  for (i=xlo;i<xhi;i++) {
    cout << "S";
  }
  cout << endl;
  for (i=0;i<sample.size();i+=2) {
    cout << "#" << i << endl;
    for (j=xlo;j<xhi;j++) {
      flag = 0;
      for (k=0;k<sample[i]->snps.size();k++) {
	pos = sample[i]->snps[k];
	if (pos == totsnps[j]) {
	  flag = 1;
	  break;
	}
      }
      if (flag == 1) { cout << "1"; }
      else { cout << "0"; }
    }
    cout << endl;
    for (j=xlo;j<xhi;j++) {
      flag = 0;
      for (k=0;k<sample[i+1]->snps.size();k++) {
	pos = sample[i+1]->snps[k];
	if (pos == totsnps[j]) {
	  flag = 1;
	  break;
	}
      }
      if (flag == 1) { cout << "1"; }
      else { cout << "0"; }
    }
    cout << endl;
    cout << endl;
  }
}






void arg::print() {
  vector<double>::iterator p;
  vector<double> totsnps;
  int i,j,k,flag;
  double pos;

  for (i=0;i<sample.size();i++) {
    for (k=0;k<sample[i]->snps.size();k++) {
      pos = sample[i]->snps[k];
      if (totsnps.empty()) {
	totsnps.push_back(pos);
      }
      else {
	flag = 0;
	p = totsnps.begin();
	while (p != totsnps.end()) {
	  if (pos == *p) { 
	    flag = 1;
	    break;
	  }
	  if (pos < *p) {
		flag = 1;
		totsnps.insert(p,pos);
	    break;
	  }
	  p++;
	}
	if (flag == 0) {
	    totsnps.push_back(pos);
	}
	  }
	}
  }

  cout << "./ms " << samplesize << " 1 -t 1 -r 1 " << coarse << endl;
  cout << rseed << " " << rseed <<" " << rseed << endl << endl << "//" << endl;
  cout << "segsites: " << totsnps.size() << endl;
  cout << "positions: ";
  for (i=0;i<totsnps.size();i++) {
    cout << totsnps[i] << " ";
  }
  cout << endl;
  for (i=0;i<sample.size();i++) {
    for (j=0;j<totsnps.size();j++) {
      flag = 0;
      for (k=0;k<sample[i]->snps.size();k++) {
	pos = sample[i]->snps[k];
	if (pos == totsnps[j]) {
	  flag = 1;
	  break;
	}
      }
      if (flag == 1) { cout << "1"; }
      else { cout << "0"; }
    }
    cout << endl;
  }
}

void arg::print4() {
  vector<double>::iterator p;
  vector<double> totsnps;
  int i,j,k,flag;
  double pos;

  for (i=0;i<sample.size();i++) {
    for (k=0;k<sample[i]->snps.size();k++) {
      pos = sample[i]->snps[k];
      if (totsnps.empty()) {
	totsnps.push_back(pos);
      }
      else {
	flag = 0;
	p = totsnps.begin();
	while (p != totsnps.end()) {
	  if (pos == *p) { 
	    flag = 1;
	    break;
	  }
	  if (pos < *p) {
		flag = 1;
		totsnps.insert(p,pos);
	    break;
	  }
	  p++;
	}
	if (flag == 0) {
	    totsnps.push_back(pos);
	}
	  }
	}
  }

  //cout << "./ms " << samplesize << " 1 -t 1 -r 1 " << coarse << endl;
  //cout << rseed <<" "<< rseed <<" " << rseed << endl << endl << "//" << endl;
  cout << "segsites: " << totsnps.size() << endl;
  cout << "positions: ";
  for (i=0;i<totsnps.size();i++) {
    cout << ((int)(round(coarse*totsnps[i]))) << " ";
  }
  cout << endl;
  for (i=0;i<sample.size();i++) {
    for (j=0;j<totsnps.size();j++) {
      flag = 0;
      for (k=0;k<sample[i]->snps.size();k++) {
	pos = sample[i]->snps[k];
	if (pos == totsnps[j]) {
	  flag = 1;
	  break;
	}
      }
      if (flag == 1) { cout << "1"; }
      else { cout << "0"; }
    }
    cout << endl;
  }
}


#endif