Added displacement distribution handling, a process similar to turning angles distribution.

#include <vector>

///define to obtain some debugging auxiliary outputs

//#define DEBUG

#define DEBUG

/**

	// ================= INITIALIZATION STUFF ================= 

	RandomWalk()

		: rnd_ForceNewSeeds(-1)

		: rnd_ForceNewSeedsAng(-1), rnd_ForceNewSeedsDisp(-1)

	{

		//initiate rnd number generators

		rnd_State=gsl_rng_alloc(gsl_rng_default);

		if (!rnd_State)

			throw std::string("RandomWalk: Couldn't initialize random number generator.");

		rnd_StateAng=gsl_rng_alloc(gsl_rng_default);

		rnd_StateDisp=gsl_rng_alloc(gsl_rng_default);

		if (!rnd_StateAng || !rnd_StateDisp)

			throw std::string("RandomWalk: Couldn't initialize random number generators.");

		unsigned long seed=-1 * (int)time(NULL) * (int)getpid();

		gsl_rng_set(rnd_State,seed);

		rnd_UseCounter=0;

		gsl_rng_set(rnd_StateAng,seed);

		gsl_rng_set(rnd_StateDisp,seed+1);

		rnd_UseCounterAng=0;

		rnd_UseCounterDisp=0;

		//initiate rnd walk parameters

		wlk_displacementBias_x=wlk_displacementBias_y=wlk_displacementBias_z=0.;

	~RandomWalk()

	{

		//clear rnd number generators

		if (rnd_State) { gsl_rng_free(rnd_State); rnd_State=NULL; }

		if (rnd_StateAng)  { gsl_rng_free(rnd_StateAng);  rnd_StateAng=NULL; }

		if (rnd_StateDisp) { gsl_rng_free(rnd_StateDisp); rnd_StateDisp=NULL; }

	}

	// ================= RND WALKING STUFF ================= 

	/*

	/**

	 * The map<angle,probability> describes probability distribution function

	 * of the walk turning angles (in radians), angles must be within [-3.14; 3.14],

	 * probability values must be non-negative.

	 *

	 * This is the main input parameter of the walk. It must be, however, pre-processed

	 * with RandomWalk::PrepareRandomTurningField() for the walks to follow it.

	 *

	 * Note the (in principle) similar attribute RandomWalk::displacementStepSizeDistribution.

	 */

	std::map<FLOAT,FLOAT> turningAngleDistribution;

	/*

	/**

	 * Sets the RandomWalk::turningAngleDistribution such that generated random

	 * walks will simulate isotropic random walks.

	 */

	{

		//relays on the interpolation in the RandomWalk::PrepareRandomTurningField()

		turningAngleDistribution.clear();

		turningAngleDistribution[-3.138f]=1.f;

		turningAngleDistribution[+3.138f]=1.f;

		turningAngleDistribution[-3.138]=1.;

		turningAngleDistribution[+3.138]=1.;

	}

	/*

	/**

	 * Sets the RandomWalk::turningAngleDistribution such that the observed

	 * distribution of turning angles will be more-or-less uniform (which

	 * is unfortunatelly not the case when RandomWalk::UseUniformDistribution()

	 * is used).

	 */

	void UseCorrectedUniformDistribution(void)

	{

		turningAngleDistribution.clear();

		for (int i=0; i < 180; ++i)

		{

			turningAngleDistribution[FLOAT(i-180)/180.*3.14]=-18.0*(float(i)/180.) +31.0;

			turningAngleDistribution[    FLOAT(i)/180.*3.14]= 18.0*(float(i)/180.) +13.0;

		}

	}

	/**

	 * Sets the RandomWalk::turningAngleDistribution such that generated random

	 * walks will simulate correlated random walks. The underlying turning angle

	 * distribution will be due to the wrapped Cauchy (circular) distribution.

		}

	}

	/*

	/**

	 * Sets the RandomWalk::turningAngleDistribution such that generated random

	 * walks will simulate correlated random walks. The underlying turning angle

	 * distribution will be due to the wrapped Gauss/normal (circular) distribution.

  protected:

	/*

	/**

	 * This attribute is used for drawing random turning angles and is hence important!

	 *

	 * It is a helper array to provide turning angles: a uniform random number generator

	std::vector<FLOAT> randomTurningField;

  public:

	/*

	/**

	 * This updates the RandomWalk::randomTurningField

	 * from the RandomWalk::turningAngleDistribution.

	 *

	 * in the distribution are "guessed" using linear interpolation.

	 * It is assumed that edge values pose both zero (0) probability, unless

	 * the RandomWalk::turningAngleDistribution dictates otherwise.

	 *

	 * Note the (in principle) similar method RandomWalk::PrepareRandomStepSizeField().

	 */

	void PrepareRandomTurningField(void)

	{

		{

			if (it->first >= -3.14 && it->first <= 3.14)

			{

				distribution[(int)floor((it->first+3.14)*360./6.28)]=it->second;

				distribution[(int)round((it->first+3.14)*360./6.28)]=it->second;

#ifdef DEBUG

				file1 << it->first << "\t" << it->second << "\n";

#endif

	}

	/**

	 * The map<stepSize,probability> describes probability distribution function

	 * of the walk displacement steps (in whatever length unit, e.g. microns),

	 * the stepSize must not be negative, probability values must be non-negative too.

	 *

	 * This attribute can be typically filled from a histogram from a representant

	 * real walk, or can be filled with some theoretical distribution.

	 *

	 * The values for unspecified angles will be interpolated from their "neighbors",

	 * it is assumed zero probability for steps less than 0.0. The "right" limit

	 * (on maximum displacement size) must be explicitly given by setting it to zero

	 * probability: \e displacementStepSizeDistribution[maxStepSize]=0.f;

	 *

	 * This is the second main input parameter of the walk. It must be, however,

	 * pre-processed with RandomWalk::PrepareRandomStepSizeField() for the walks

	 * to follow it.

	 *

	 * Note the (in principle) similar attribute RandomWalk::turningAngleDistribution.

	 */

	std::map<FLOAT,FLOAT> displacementStepSizeDistribution;

	/**

	 * Attempts to create a single-peaked distribution of available displacement

	 * step sizes.

	 *

	 * \param[in] dispSize		the position of the peak

	 */

	void UseFixedDisplacementStep(const FLOAT dispSize)

	{

		displacementStepSizeDistribution.clear();

		displacementStepSizeDistribution[0.]=0.;

		displacementStepSizeDistribution[dispSize-1.0]=0.;

		displacementStepSizeDistribution[dispSize-0.2]=0.;

		displacementStepSizeDistribution[dispSize-0.1]=1.;

		displacementStepSizeDistribution[dispSize+0.1]=1.;

		displacementStepSizeDistribution[dispSize+0.2]=0.;

		displacementStepSizeDistribution[dispSize+1.0]=0.;

	}

  protected:

	/**

	 * This attribute is used for drawing random displacement step sizes and is hence important!

	 *

	 * It is a helper array to provide step sizes: a uniform random number generator

	 * chooses index to this array and the value found there is the next step size.

	 * Clearly, the histogram of step sizes values present in this array should be

	 * the same as is given in the RandomWalk::displacementStepSizeDistribution

	 * attribute. For convenience, use the RandomWalk::PrepareRandomStepSizeField()

	 * method to fill this attribute.

	 *

	 * The array should be such long that every "discrete step" (assuming user

	 * given sampling rate) whose probability is non-zero should be present

	 * in the array.

	 */

	std::vector<FLOAT> randomStepSizeField;

  public:

	/**

	 * This updates the RandomWalk::randomStepSizeField

	 * from the RandomWalk::displacementStepSizeDistribution.

	 *

	 * Since the latter may sample the interval [0.0; some_max] arbitrarily

	 * while this function considers given user stepping, the probabilities of step

	 * sizes missing in the distribution are "guessed" using linear interpolation.

	 * It is assumed that right edge value poses zero (0) probability.

	 *

	 * \param[in] stepSampling	user given granularity of step sizes

	 *

	 * The granularity should be fine such that the RandomWalk::displacementStepSizeDistribution

	 * will not tend to aggregate its sampling nodes due to overly large

	 * \e stepSampling -- this method is not accommodated for that.

	 *

	 * Note the (in principle) similar method RandomWalk::PrepareRandomTurningField().

	 */

	void PrepareRandomStepSizeField(const FLOAT stepSampling=1.)

	{

		//see the RandomWalk::PrepareRandomTurningField() for some more comments...

		if (stepSampling <= 0.)

			throw std::string("RandomWalk: stepSampling parameter must be strictly positive.");

		randomStepSizeField.clear();

		//first:

		//create and init the local resampled distribution 

		//for which we need to determine the maximum value

		typename std::map<FLOAT,FLOAT>::const_iterator it

		   =displacementStepSizeDistribution.end();

		--it;

		//if the last stored probability is not zero, inject ours...

		if (it->second != 0.)

		{

			displacementStepSizeDistribution[std::max(it->first,FLOAT(0.))+stepSampling]=0.;

			++it;

		}

		//create and init the local resampled distribution 

		const int distLength=(int)round(it->first/stepSampling) +1;

		FLOAT distribution[distLength];

		for (int i=0; i < distLength; ++i) distribution[i]=-1.;

		distribution[0]=distribution[distLength-1]=0.;

#ifdef DEBUG

		std::ofstream file1("stepSizeDist_coarse.dat");

#endif

		//fill it where we can

		it=displacementStepSizeDistribution.begin();

		while (it != displacementStepSizeDistribution.end())

		{

			if (it->first >= 0.0)

			{

				distribution[(int)round(it->first/stepSampling)]=it->second;

#ifdef DEBUG

				file1 << it->first << "\t" << it->second << "\n";

#endif

			}

			++it;

		}

#ifdef DEBUG

		file1.close();

#endif

		//search for uninitialized values and calculate them

		int lastLeftInited=0;

		int lastRightInited=1;

		while (lastRightInited < distLength && distribution[lastRightInited] == -1.) ++lastRightInited;

		//note that distribution[distLength-1] != -1. (which renders the first condition useless)

		//sum over all probabilities...

		FLOAT Sum=distribution[0]+distribution[distLength-1];

		for (int i=1; i < distLength-1; ++i)

		{

			if (distribution[i] == -1.)

			{

				//found uninitialized value, interpolate from lastLeft and lastRight ones

				distribution[i]=FLOAT(i-lastLeftInited)/FLOAT(lastRightInited-lastLeftInited)

				                *(distribution[lastRightInited]-distribution[lastLeftInited])

									 +distribution[lastLeftInited];

			}

			else

			{

				//found initialized one, move the boundaries

				lastLeftInited=i; //should be also equal to lastRightInited

				lastRightInited=i+1;

				while (lastRightInited < distLength && distribution[lastRightInited] == -1.) ++lastRightInited;

			}

			Sum+=distribution[i];

		}

#ifdef DEBUG

		//debug: print out the fine-sampled distribution

		std::ofstream file2("stepSizeDist_fine.dat");

		for (int i=0; i < distLength; ++i)

			file2 << i*stepSampling << "\t" << distribution[i] << "\n";

		file2.close();

#endif

		//second:

		const FLOAT Factor=10000. / Sum;

		//third:

#ifdef DEBUG

		std::ofstream file3("stepSizeDist_rndFieldHistogram.dat");

#endif

		for (int i=0; i < distLength; ++i)

		{

			const int length=(int)floor(distribution[i]*Factor);

			for (int j=0; j < length; ++j)

				randomStepSizeField.push_back(i*stepSampling);

#ifdef DEBUG

			file3 << i*stepSampling << "\t" << length << "\n";

#endif

		}

#ifdef DEBUG

		file3.close();

#endif

	}

	///Set bias vector to use.

	void SetBiasVector(const FLOAT x,const FLOAT y,const FLOAT z=0.)

	{ wlk_displacementBias_x=x; wlk_displacementBias_y=y; wlk_displacementBias_z=z; }

		wlk_xy_heading=newHeading;

	}

	/**

	 * Indication whether all data structures are prepared.

	 * Of course, it does not investigate the content/meaningfulness

	 * of the data.

	 */

	bool ReadyToStartSimulation(void)

	{

		return (randomTurningField.size() > 0

		     && randomStepSizeField.size() > 0);

	}

	///Do one step of the 2D random walk

	void DoOneStep2D(void)

	{

		//choose random turning angle from the underlying distribution

		//(which is represented in this helper field)

		FLOAT deltaHeading=randomTurningField[

		          (int)floor(rnd_GetValue(0.,randomTurningField.size()-0.001)) ];

		const FLOAT deltaHeading=randomTurningField[

		          (int)floor(rnd_GetValueAng(0.,randomTurningField.size()-0.001)) ];

		//const FLOAT deltaHeading=rnd_GetValueAng(-3.14,+3.14);

		//update the azimuth (absolute heading)

		wlk_xy_heading+=deltaHeading;

		//choose random step size from the underlying distribution

		const FLOAT dispSize=randomStepSizeField[

		          (int)floor(rnd_GetValueDisp(0.,randomStepSizeField.size()-0.001)) ];

		//update the (normalized) position displacement "vector"

		wlk_displacement_x=cos(wlk_xy_heading);

		wlk_displacement_y=sin(wlk_xy_heading);

		wlk_displacement_x=dispSize*cos(wlk_xy_heading);

		wlk_displacement_y=dispSize*sin(wlk_xy_heading);

		//add the bias...

		wlk_displacement_x+=wlk_displacementBias_x;

	// ================= RND GENERATORS STUFF ================= 

	///helper random number generator: Flat/uniform distribution within [\e A,\e B]

	FLOAT rnd_GetValue(const FLOAT A, const FLOAT B)

	FLOAT rnd_GetValueAng(const FLOAT A, const FLOAT B)

	{

		if (rnd_UseCounterAng == rnd_ForceNewSeedsAng)

		{

			unsigned long seed=-1 * (int)time(NULL) * (int)getpid();

			gsl_rng_set(rnd_StateAng,seed);

			rnd_UseCounterAng=0;

		}

		++rnd_UseCounterAng;

		return ( gsl_ran_flat(rnd_StateAng, A,B) );

	}

	///helper random number generator: Flat/uniform distribution within [\e A,\e B]

	FLOAT rnd_GetValueDisp(const FLOAT A, const FLOAT B)

	{

		if (rnd_UseCounter == rnd_ForceNewSeeds)

		if (rnd_UseCounterDisp == rnd_ForceNewSeedsDisp)

		{

			unsigned long seed=-1 * (int)time(NULL) * (int)getpid();

			gsl_rng_set(rnd_State,seed);

			rnd_UseCounter=0;

			gsl_rng_set(rnd_StateDisp,seed);

			rnd_UseCounterDisp=0;

		}

		++rnd_UseCounter;

		return ( gsl_ran_flat(rnd_State, A,B) );

		++rnd_UseCounterDisp;

		return ( gsl_ran_flat(rnd_StateDisp, A,B) );

	}

	///helper random number generator: generator states

	gsl_rng *rnd_State;

	gsl_rng *rnd_StateAng, *rnd_StateDisp;

	///helper random number generator: use counters (to enforce re-seeding)

	long rnd_UseCounter;

	long rnd_UseCounterAng, rnd_UseCounterDisp;

	///helper random number generator: how many rnd numbers before re-seeding

	const long rnd_ForceNewSeeds;

	const long rnd_ForceNewSeedsAng, rnd_ForceNewSeedsDisp;

};

const int TimeDeltaMinutes=29; //simulations - cell cycle = 50 frames

#define HISTOGRAM_OUTPUT

//#define INDIVIDUAL_PATHS_OUTPUT

//#define COMBINED_PATHS_OUTPUT

#define INDIVIDUAL_PATHS_OUTPUT

#define COMBINED_PATHS_OUTPUT

int main(int argc,char **argv)

	//init this walk

	RandomWalk<float> rw;

	//rw.UseWrappedCauchyDistribution(0.8); // => r=0.8

	rw.UseWrappedNormalDistribution(0.4); // => r=0.8

	//rw.UseWrappedCauchyDistribution(0.5);

	//rw.UseWrappedNormalDistribution(0.4);

	//rw.UseUniformDistribution();

	rw.UseCorrectedUniformDistribution();

	//My own testing 3-peaks distribution:

	/*

	rw.turningAngleDistribution[-1.0]=0.;

	rw.turningAngleDistribution[+1.0]=0.;

	*/

	rw.PrepareRandomTurningField();

	float speed=1.;

	//rw.UseFixedDisplacementStep(2.0f); //this is basically the speed

	rw.displacementStepSizeDistribution[0]=1979;

	rw.displacementStepSizeDistribution[0.33333]=607;

	rw.displacementStepSizeDistribution[0.66667]=533;

	rw.displacementStepSizeDistribution[1]=663;

	rw.displacementStepSizeDistribution[1.3333]=610;

	rw.displacementStepSizeDistribution[1.6667]=660;

	rw.displacementStepSizeDistribution[2]=573;

	rw.displacementStepSizeDistribution[2.3333]=510;

	rw.displacementStepSizeDistribution[2.6667]=361;

	rw.displacementStepSizeDistribution[3]=361;

	rw.displacementStepSizeDistribution[3.3333]=231;

	rw.displacementStepSizeDistribution[3.6667]=229;

	rw.displacementStepSizeDistribution[4]=173;

	rw.displacementStepSizeDistribution[4.3333]=157;

	rw.displacementStepSizeDistribution[4.6667]=87;

	rw.displacementStepSizeDistribution[5]=138;

	rw.displacementStepSizeDistribution[5.3333]=68;

	rw.displacementStepSizeDistribution[5.6667]=73;

	rw.displacementStepSizeDistribution[6]=39;

	rw.displacementStepSizeDistribution[6.3333]=64;

	rw.displacementStepSizeDistribution[6.6667]=44;

	rw.displacementStepSizeDistribution[7]=34;

	rw.displacementStepSizeDistribution[7.3333]=26;

	rw.displacementStepSizeDistribution[7.6667]=27;

	rw.displacementStepSizeDistribution[8]=20;

	rw.displacementStepSizeDistribution[8.3333]=15;

	rw.displacementStepSizeDistribution[8.6667]=5;

	rw.displacementStepSizeDistribution[9]=10;

	rw.displacementStepSizeDistribution[9.3333]=4;

	rw.displacementStepSizeDistribution[9.6667]=16;

	rw.PrepareRandomStepSizeField(0.1f);

	/*

	std::cout << "Ready for simulation? "

	          << ((rw.ReadyToStartSimulation())? "yes":"no")

				 << "\n";

	*/

	//several random walks...

	for (int attemptNo=0; attemptNo < 10000; ++attemptNo)

	for (int attemptNo=0; attemptNo < 200; ++attemptNo)

	{

		//std::cout << "========== " << attemptNo << "==========\n";

		pathsI.SetVoxel(round(posX),round(posY),0, 100);

#endif

#ifdef COMBINED_PATHS_OUTPUT

		if (attemptNo < 50)

		if (attemptNo < 10)

			pathsC.SetVoxel(round(posX),round(posY),0, attemptNo+1);

#endif

		rw.ResetXYHeading(0.);

		//do the walk...

		for (int step=0; step < 60; ++step)

		for (int step=0; step < 40; ++step)

		{

			//update heading

			float dx,dy;

			rw.DoOneStep2D(dx,dy);

			//update speed...

			//speed += GetRandomGauss(0.f,1.f, &randState);

			//update current position

			posX += speed * dx;

			posY += speed * dy;

			posX += dx;

			posY += dy;

			//draw current position

#ifdef INDIVIDUAL_PATHS_OUTPUT

			pathsI.SetVoxel(round(posX),round(posY),0, 100+5*step);

#endif

#ifdef COMBINED_PATHS_OUTPUT

			if (attemptNo < 50)

			if (attemptNo < 10)

				pathsC.SetVoxel(round(posX),round(posY),0, attemptNo+1);

#endif