void ParticleWorld::updateParticles(double deltaTime)
{
    for(int i = 0; i < particleList.size(); i++) {
        double x = particleList[i].box.x;
        double y = particleList[i].box.y;
        double w = particleList[i].box.w;
        double h = particleList[i].box.h;
        
        double jitter = 20.0;
        double cx = x+w/2.0;
        double cy = y+h/2.0;
        w = clamp(w+drand(-jitter, jitter), (double)minBoxSize, (double)width );
        h = clamp(h+drand(-jitter, jitter), (double)minBoxSize, (double)height);
        
        cx = clamp(cx+drand(-jitter, jitter), w/2.0+1, ((double)width )-w/2.0-1);
        cy = clamp(cy+drand(-jitter, jitter), h/2.0+1, ((double)height)-h/2.0-1);
        particleList[i].box.x = cx-w/2.0;
        particleList[i].box.y = cy-h/2.0;
        particleList[i].box.w = w;
        particleList[i].box.h = h;
        
        if(particleList[i].box.x < 0) { std::cout << "WTF1" << std::endl; }
        if(particleList[i].box.y < 0) { std::cout << "WTF2" << std::endl; }
        if(particleList[i].box.x+particleList[i].box.w >= width) { std::cout << "WTF3" << std::endl; }
        if(particleList[i].box.y+particleList[i].box.h >= height) { std::cout << "WTF4" << std::endl; }
    }
}

int main(int argc, char **argv)
{
  int numBodies = 1000;
  bool checkErrors = true;

  // Parse custom command line args
  for (int i = 1; i < argc; ++i) {
    if (strcmp(argv[i],"-N") == 0) {
      i++;
      numBodies = atoi(argv[i]);
    } else if (strcmp(argv[i],"-nocheck") == 0) {
      checkErrors = false;
    }
  }

  // Init the FMM Kernel and options
  FMMOptions opts = get_options(argc, argv);
  typedef UnitKernel kernel_type;
  kernel_type K;

  typedef kernel_type::point_type point_type;
  typedef kernel_type::source_type source_type;
  typedef kernel_type::target_type target_type;
  typedef kernel_type::charge_type charge_type;
  typedef kernel_type::result_type result_type;

  // Init points and charges
  std::vector<source_type> points(numBodies);
  for (int k = 0; k < numBodies; ++k)
    points[k] = source_type(drand(), drand(), drand());

  std::vector<charge_type> charges(numBodies);
  for (int k = 0; k < numBodies; ++k)
    charges[k] = drand();

  // Build the FMM
  //fmm_plan plan = fmm_plan(K, bodies, opts);
  FMM_plan<kernel_type> plan = FMM_plan<kernel_type>(K, points, opts);

  // Execute the FMM
  //fmm_execute(plan, charges, target_points);
  std::vector<result_type> result = plan.execute(charges);

  // Check the result
  if (checkErrors) {
    std::vector<result_type> exact(numBodies);

    // Compute the result with a direct matrix-vector multiplication
    Direct::matvec(K, points, charges, exact);

    int wrong_results = 0;
    for (unsigned k = 0; k < result.size(); ++k) {
      printf("[%03d] exact: %lg, FMM: %lg\n", k, exact[k], result[k]);

      if ((exact[k] - result[k]) / exact[k] > 1e-13)
        ++wrong_results;
    }
    printf("Wrong counts: %d\n", wrong_results);
  }
}

// plays event group extinction and group replication
void event_group_ext_rep()
{    
  int j, g, gr;    
  double m, t;
  m = pi_g[0];
  t = -m;
  // set distribution probability    
  for(j = G; j--; ){    
    Gr_dist->p[j] = pi_g[j];                // probabilith of recolonisation of group proportional to +ve of group payoff
    Ge_dist->p[j] = -pi_g[j];             // probability of extinction of group proporttional to -ve of group payoff
    if(m < Gr_dist->p[j]) m = Gr_dist->p[j];
    if(t < Ge_dist->p[j]) t = Ge_dist->p[j];    
  }  
  exp_initdist(Gr_dist, 0);
  exp_initdist(Ge_dist, 0);
  g = drand(Ge_dist);        // select group to extinction   
  gr = drand(Gr_dist);      // select group to repopulate     
  if( g == gr) return;       // if same, then noop
  // make selected group extinct and repopulate by another
  for(j = N; j--;){    
    // copy strategies from repopulating group individuals to extincting individuals    
    x[g][j] = x[gr][j];    // replace strategy            
    dxi[g][j] = dxi[gr][j];  //replace threshold
    dsi[g][j] = dsi[gr][j];  // replace aggressiveness
    pi[g][j] = pi[gr][j];     // replace payoff     
    Api[g][j] = Api[gr][j];
  }    
  pi_g[g] = pi_g[gr];                          // replace group payoff  
}

void
TcpSrcPeriodic::doNextEvent() {
    if (_idle_time==0||_active_time==0) {
        _is_active = true;
        startflow();
        return;
    }

    if (_is_active) {
        if (eventlist().now()>=_end_active && _idle_time!=0 && _active_time!=0 ) {
            _is_active = false;

            //this clears RTOs too
            reset();
            eventlist().sourceIsPendingRel(*this,(simtime_picosec)(2*drand()*_idle_time));
        }
        else if (_rtx_timeout_pending) {
            TcpSrc::doNextEvent();
        }
        else {
            cout << "Wrong place to be in: doNextDEvent 1" << eventlist().now() << " end active " << _end_active << "timeout " << _rtx_timeout_pending << endl;
            //maybe i got a timeout here. How?
            exit(1);
        }
    }
    else {
        _is_active = true;
        ((TcpSinkPeriodic*)_sink)->reset();
        _start_active = eventlist().now();
        _end_active = _start_active + (simtime_picosec)(2*drand()*_active_time);
        eventlist().sourceIsPending(*this,_end_active);
        startflow();
    }
}

int main(int argc, char** argv)
{
  if (argc < 2) {
    std::cerr << "Usage: test_tree NUM_POINTS\n";
    exit(1);
  }

  int N = atoi(argv[1]);

  typedef Vec<3,double[3]> point_type;

  Octree<point_type> otree(BoundingBox<point_type>(point_type(0), point_type(1)));

  std::vector<point_type> points;
  for (int k = 0; k < N; ++k) {
    point_type p;
    p[0] = drand();
    p[1] = drand();
    p[2] = drand();
    points.push_back(p);
  }

  otree.construct_tree(points.begin(), points.end());

  std::cout << otree << "\n";

  return 0;
}

static int artefact_loose_to_cave (db_t *database, struct Cave *cave)
{
  db_result_t *result;
  dstring_t *query;
  int x, y;
  int minX, minY, maxX, maxY, rangeX, rangeY;

  /* number between -ALR <= n <= ALR */
  x = (int) ((ARTEFACT_LOST_RANGE * 2 + 1) * drand()) - ARTEFACT_LOST_RANGE;
  y = (int) ((ARTEFACT_LOST_RANGE * 2 + 1) * drand()) - ARTEFACT_LOST_RANGE;

  x += cave->xpos;
  y += cave->ypos;   /* these numbers may be out of range */

  if (! map_get_bounds(database, &minX, &maxX, &minY, &maxY)) {
    return 0;
  }
  rangeX = maxX - minX +1;
  rangeY = maxY - minY +1;

  x = ( (x-minX+rangeX) % (rangeX) ) + minX;
  y = ( (y-minY+rangeY) % (rangeY) ) + minY;

  query = dstring_new("SELECT caveID FROM " DB_TABLE_CAVE
          " WHERE xCoord = %d AND yCoord = %d", x, y);

  debug(DEBUG_SQL, "%s", dstring_str(query));

  result = db_query_dstring(database, query);

  return db_result_next_row(result) ? db_result_get_int(result, "caveID") : 0;
}

int main(void) {

    for (int i = 0; i < num; ++i) {
        bool neg = rand() < mid;
        bool zrw = rand() < mid;
        //bool zrn = rand() < mid;

        if (neg) {
            if (zrw) {
                std::cout << "0 -" << rand() << "/" << drand() << std::endl;
            }
            //else if (zrn) {
            //    std::cout << "-" << rand() << " 0/" << drand() << std::endl;
            //}
            else {
                std::cout << "-" << rand() << " " << rand() << "/" << drand() << std::endl;
            }
        }
        else {
            if (zrw) {
                std::cout << "0 " << rand() << "/" << drand() << std::endl;
            }
            //else if (zrn) {
            //    std::cout << rand() << " 0/" << drand() << std::endl;
            //}
            else {
                std::cout << rand() << " " << rand() << "/" << drand() << std::endl;
            }
        }
    }

    return 0;
}

void ParticleWorld::breedParticles()
{
    std::vector<Particle> newParticleList;
    double scoreSum = 0.0;
    for(int i = 0; i < particleList.size(); i++) {
        scoreSum += particleList[i].score;
    }
    int desiredParticleCount = particleList.size();
    for(int newParticleI = 0; newParticleI < desiredParticleCount; newParticleI++) {
        Particle newParticle;
        
        //double createNewScore = 1-scoreSum;
        //std::cout << scoreSum << "; " << createNewScore << std::endl;
        //double selection = drand(0, scoreSum+createNewScore);
        //if(selection > scoreSum) {
        //    setRandParticle(newParticle);
        //} else {
        
        if(drand(0.0, 1.0) < 0.03) {
            setRandParticle(newParticle);
        } else {
            double selection = drand(0, scoreSum);
            double total = 0;
            for(int oldParticleI = 0; oldParticleI < particleList.size(); oldParticleI++){
                total += particleList[oldParticleI].score;
                if (total > selection) {
                    newParticle = particleList[oldParticleI];
                    break;
                }
            }
        }
        newParticleList.push_back(newParticle);
    }
    particleList = newParticleList;
}

void uniform_stress(int dim, SparseMatrix A, real *x, int *flag){
  UniformStressSmoother sm;
  real lambda0 = 10.1, M = 100, scaling = 1.;
  real res;
  int maxit = 300, samepoint = TRUE, i, k, n = A->m;
  SparseMatrix B = NULL;

  *flag = 0;

  /* just set random initial for now */
  for (i = 0; i < dim*n; i++) {
    x[i] = M*drand();
  }

  /* make sure x is not all at the same point */
  for (i = 1; i < n; i++){
    for (k = 0; k < dim; k++) {
      if (ABS(x[0*dim+k] - x[i*dim+k]) > MACHINEACC){
	samepoint = FALSE;
	i = n;
	break;
      }
    }
  }

  if (samepoint){
    srand(1);
#ifdef DEBUG_PRINT
    fprintf(stderr,"input coordinates to uniform_stress are the same, use random coordinates as initial input");
#endif
    for (i = 0; i < dim*n; i++) x[i] = M*drand();
  }

  B = get_distance_matrix(A, scaling);
  assert(SparseMatrix_is_symmetric(B, FALSE));

  sm = UniformStressSmoother_new(dim, B, x, 1000000*lambda0, M, flag);
  res = UniformStressSmoother_smooth(sm, dim, x, maxit);
  UniformStressSmoother_delete(sm);

  sm = UniformStressSmoother_new(dim, B, x, 10000*lambda0, M, flag);
  res = UniformStressSmoother_smooth(sm, dim, x, maxit);
  UniformStressSmoother_delete(sm);

  sm = UniformStressSmoother_new(dim, B, x, 100*lambda0, M, flag);
  res = UniformStressSmoother_smooth(sm, dim, x, maxit);
  UniformStressSmoother_delete(sm);

  sm = UniformStressSmoother_new(dim, B, x, lambda0, M, flag);
  res = UniformStressSmoother_smooth(sm, dim, x, maxit);
  UniformStressSmoother_delete(sm);

  scale_to_box(0,0,7*70,10*70,A->m,dim,x);;

  SparseMatrix_delete(B);

}

static inline void
random_move (gint *x,
             gint *y,
             gint  max)
{
  gdouble angle  = drand () * G_PI;
  gdouble radius = drand () * (gdouble) max;
  *x = (gint) (radius * cos (angle));
  *y = (gint) (radius * sin (angle));
}

GLvoid InitGL(GLvoid) {
int i;
glClearColor(0.0f, 0.0f, 0.0f, 0.0f);
// This Will Clear The Background Color To Black
glClearDepth(1.0);
// Enables Clearing Of The Depth Buffer
glDepthFunc(GL_LESS);
// The Type Of Depth Test To Do
glEnable(GL_DEPTH_TEST);
// Enables Depth Testing
glShadeModel(GL_SMOOTH);
// Enables Smooth Color Shading
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
// Reset The Projection Matrix
gluPerspective(45.0f,
 (GLfloat)kWindowWidth/(GLfloat)kWindowHeight,0.1f,100.0f);
// Calculate The Aspect Ratio Of The Window
glMatrixMode(GL_MODELVIEW);
for (i=0; i<N; i++) {
mr[i].r = rect_open();
rect_init(mr[i].r, RECT5whxyz
,drand()/5.0+0.1,drand()/5.0+0.1,
drand(),drand(),drand());
rect_setcolor(mr[i].r, drand(), drand(), drand());
mr[i].dx=drand()/10;
mr[i].dy=mr[i].dx;
} /*for*/
}

int main()
{
  typedef LaplaceSpherical kernel_type;
  kernel_type K(5);
  typedef kernel_type::point_type point_type;
  typedef kernel_type::charge_type charge_type;
  typedef kernel_type::result_type result_type;

  FMMOptions opts;
  opts.set_mac_theta(.5);
  opts.set_max_per_box(125);    // optimal ncrit

  std::vector<std::pair<unsigned,double>> times;
  double tic, toc;

  
  int numBodies = 1000000;
  // initialize points
  std::vector<point_type> points(numBodies);
  for (int k=0; k<numBodies; ++k){
          points[k] = point_type(drand(), drand(), drand());
  }

  // initialize charges
  std::vector<charge_type> charges(numBodies);
  for (int k=0; k<numBodies; ++k){
          charges[k] = drand();
  }
  

  // loop in ncrit
  for (int ncrit=50; ncrit<=400; ncrit+=50){
  
  opts.set_max_per_box(ncrit);
  // create FMM plan
  FMM_plan<kernel_type> plan = FMM_plan<kernel_type>(K, points, opts);
  // execute FMM
  // run 3 times and make an average
  int nt = 3;    // number of identical runs for timing
  std::vector<double> timings(nt);
  std::vector<result_type> result(numBodies);
  for (int i=0; i<nt; i++){
	tic = get_time();
  	result = plan.execute(charges);
  	toc = get_time();
	timings[i] = toc-tic;
  }
  
  double FMM_time = std::accumulate(timings.begin(), timings.end(), 0.0) / timings.size();
  std::cout << ncrit << " " << FMM_time << std::endl;
 }
  return 0;

}

int fill8bitData(void)
{
	int len, i;

	len = 1 + (int)drand((MAX_LENGTH - 2));
	for(i=0; i<len; i++) {
		data[i] = (unsigned char)drand(256);
	}
	data[len] = '\0';

	return len;
}

int fillANData(void)
{
	int len, i;

	len = 1 + (int)drand((MAX_LENGTH - 2));
	for(i=0; i<len; i++) {
		data[i] = AN[(int)drand(45)];
	}
	data[len] = '\0';

	return len;
}

Compute::Compute(void) {
  current_time = GetTickCount();
  double speed = 0.25*GetSystemMetrics(SM_CXSCREEN);
  double angle = 8.0*std::atan(1.0)*drand();
  velocity_x = speed*std::cos(angle);
  velocity_y = speed*std::sin(angle);
  int screen_width = GetSystemMetrics(SM_CXSCREEN);
  int screen_height =  GetSystemMetrics(SM_CYSCREEN);
  position_x = screen_width*drand();
  position_y =screen_height*drand();
  offset = int(std::sqrt(0.01*screen_width*screen_height));
}

// Box-Muller optimized
static double gaussRandom() 
{
	double u,v,r,c;
	u = 2.0*drand()-1;
	v = 2.0*drand()-1;
	r = u*u + v*v;
	if (r == 0.0 || r > 1.0) 
		return gaussRandom();

	c = sqrtl(-2*logl(r)/logl(M_E)/r);
	return u*c;
}

/* Fullscreen antialiasing. Ultra-slow. */
Ray camera_ray_aa(Camera *cam, int i, int j, int sample, double near)
{
	float offx, offy;
	int p, q;
	const float n = (float) config->aa_samples;

	p = sample % config->aa_samples;
	q = sample / config->aa_samples;
	offx = (p + drand()) / n;
	offy = (q + drand()) / n;

	return cam_ray_internal(cam, i, j, offx, offy, near);
}

int main(int argc, char* argv){
  uint64_t seed;
  Compl z, c, tmp;
  int i, j;
  int N1, N0;
  double etime0, etime1, cptime;
  double A; 
  double r;
  int n;
  // boundary
  lowerLeft.real = -2.0;
  lowerLeft.imag = 0;
  upperRight.real = 0.5;
  upperRight.imag = 1.125;
  omp_set_dynamic(0);
  #pragma omp parallel
  dsrand(12345);

  N1 = N0 = 0;
  timing(&etime0, &cptime);

  #pragma omp parallel firstprivate(j, n, c, z, tmp, stride) \
                      shared(i, N0, N1, lowerLeft, upperRight, threshold, MAX_ITERATE)
  #pragma omp for schedule(static, 10) collapse(2)
  for(i = 0; i < (int)((upperRight.real - lowerLeft.real) / stride); i++){
    for(j = 0; j < (int)((upperRight.imag - lowerLeft.imag) / stride); j++){
      if(i == 0 && j == 0 && omp_get_thread_num() == 0) printf("Threads: %d\n", omp_get_num_threads());
      c.real = lowerLeft.real + (drand() + i) * stride;
      c.imag = lowerLeft.imag + (drand() + j) * stride;
      z = c;
      for(n = 0; n < MAX_ITERATE; n++){
        multiply(&z, &c, &tmp);
        z = tmp;
        if(lengthsq(&z) > threshold * threshold){
            break;
        }        
      }
      if(n == MAX_ITERATE){
        #pragma omp critical(N1)
        N1++;
      }else{
        #pragma omp critical(N0)
        N0++;
      }
    }
  }
  timing(&etime1, &cptime);
  A = 2.0 * N1 / (N0 + N1) * (upperRight.real - lowerLeft.real) * (upperRight.imag - lowerLeft.imag);
  printf("area is %f, time elapsed is %f, total cell: %d\n", A, etime1 - etime0, N1 + N0);
}

int mutation(individual *ind)
{
     int i;

     if (ind == NULL)
     {
	 return (1);
     }
     
     for (i = 0; i < ind->n; i++)
     {
	 if (drand(1) <= variable_mutation_probability)
	 {
	     double eta = eta_mutation;
	     double u = drand(1.0);
	     double delta = 0;
	     double x = ind->x[i];
	     double lb = 0;    /* lower bound of variable i */
	     double ub = 1;    /* upper bound of variable i */
	     double diff = ub - lb;  /* range of variable i */
	     double maxmut0 = x - lb;
	     double maxmut = ub - x;
	     double delta_max = maxmut0 / diff;
	     if (maxmut0 > maxmut)
	     {
		 delta_max = maxmut / diff;
	     }
	     
	     if (u < 0.5)
	     {
		 double b =  2*u + (1-2*u)*(pow(1-delta_max,(eta+1)));
		 delta = pow(b,(1.0/(eta+1))) - 1.0;
	     }
	     else
	     {
		 double b = 2*(1-u) + 2*(u-0.5)*(pow(1-delta_max,(eta+1)));
		 delta = 1.0 - pow(b,(1.0/(eta+1)));
	     }
	     if (delta > delta_max)  /* machine accuracy problem */
		 delta = delta_max;
	     else if (delta < -delta_max)
		 delta = -delta_max;
	     
	     ind->x[i] = x + delta * diff;
	 }
     }
     
     return (0);
}

void CStaticStructure::ProcessFrame(unsigned int dwCurrentTime,double dTimeFraction)
{
	CEntityBase::ProcessFrame(dwCurrentTime,dTimeFraction);
	m_dwNextProcessFrame=dwCurrentTime+10;
	
	if(GetState()==eStaticStructureState_Destroyed){return;}
	if(GetState()==eStaticStructureState_Normal)
	{
		size_t nAnimationToSet=0;
		double dMaxHealth=GetMaxHealth();
		unsigned int nAnimations=m_pTypeBase->GetStateAnimations(ENTITY_STATE_BASE);
		nAnimationToSet=(size_t)(((dMaxHealth-m_dHealth)/dMaxHealth)*((double)nAnimations));
		if(nAnimationToSet>nAnimations-1){nAnimationToSet=nAnimations-1;}
		SetState(ENTITY_STATE_BASE,(int)nAnimationToSet);
	}
	
	bool bAllChildDead=true;
	for(unsigned int x=0;x<m_vChildren.size();x++){if(m_vChildren[x].piEntity->GetHealth()>0){bAllChildDead=false;}}
	m_dwDamageType=(bAllChildDead?m_nConfiguredDamageType:DAMAGE_TYPE_NONE);
	
	if(m_dwAlignment==ENTITY_ALIGNMENT_ENEMIES)
	{
		GetTarget();
	}
	
	if(m_piTarget && dwCurrentTime>m_dwNextShotTime && m_vWeapons.size())
	{
		bool bVisible=g_PlayAreaManagerWrapper.m_piInterface && g_PlayAreaManagerWrapper.m_piInterface->IsVisible(m_PhysicInfo.vPosition,0);
		if(bVisible)
		{
			double dDifficulty=g_PlayerManagerWrapper.m_piInterface->GetEffectiveDifficulty();
			double dTimeFirstShotMin=m_pType->m_dTimeFirstShotMin/dDifficulty;
			double dTimeFirstShotMax=m_pType->m_dTimeFirstShotMax/dDifficulty;
			double dTimeBetweenShotsMin=m_pType->m_dTimeBetweenShotsMin/dDifficulty;
			double dTimeBetweenShotsMax=m_pType->m_dTimeBetweenShotsMax/dDifficulty;
			
			if(m_bFirstTimeVisible)
			{
				m_bFirstTimeVisible=false;
				m_dwNextShotTime=dwCurrentTime+drand()*(dTimeFirstShotMax-dTimeFirstShotMin)+dTimeFirstShotMin;
			}
			else
			{
				for(unsigned int x=0;x<m_vWeapons.size();x++){FireWeapon(x,dwCurrentTime);}
				m_dwNextShotTime=dwCurrentTime+drand()*(dTimeBetweenShotsMax-dTimeBetweenShotsMin)+dTimeBetweenShotsMin;
			}
		}
	}
}