/*  Ant-based Clustering    Copyright (C) 2004 Julia Handl    Email: Julia.Handl@gmx.de    This program is free software; you can redistribute it and/or modify    it under the terms of the GNU General Public License as published by    the Free Software Foundation; either version 2 of the License, or    (at your option) any later version.    This program is distributed in the hope that it will be useful,    but WITHOUT ANY WARRANTY; without even the implied warranty of    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the    GNU General Public License for more details.    You should have received a copy of the GNU General Public License    along with this program; if not, write to the Free Software    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA*//***************************************************date: 7.4.2003author: Julia Handl (julia.handl@gmx.de)description: implementation of a the two-dimensional gridused for ant-based clusteringThe grid is a simple two-dimensional array thatcan store integers (these index into the data collection). It contains an additionalindex structure that keeps track of the dataitems on the grid and their respective positions(such that ants can directly access them and do not have to search). The grid is toroidaland the update of grid positions thereforerequires to "wrap around" corrdinates at the borders. ***************************************************/#include "grid.h"#include <iostream>#include "random.h"extern long idum;/* Constructor:    memory allocation and initialisation*/grid::grid(conf * c, databin<USED_DATA_TYPE> * b) {       /* storage of the pointer to data and current parameter settings */    par = c;    bin = b;    /* preparation of the grid-cells */    cells = new int* [par->imax];    for (int i=0; i<par->imax; i++) {	cells[i] = new int [par->jmax];	for (int j=0; j<par->jmax; j++) {	    cells[i][j] = FREE;	}    }    /* preparation of the index of free data items */    index = new position<int>[par->binsize];    last = 0;    items = 0;}/* Access to index */position<int> & grid::getposition(int i) {    return index[i];}/* Destructor: delete all allocated memory */grid::~grid() {    /* delete grid cells */    for (int i=0; i<par->imax; i++) {	delete [] cells[i];    }    delete [] cells;    /* delete index */    delete [] index;}/* Print the current spatial distribution to a file;   required by the visualisation routine */void grid::print() {    ofstream to("grid.dat");    for (int i=0; i<par->imax; i++) {	for (int j=0; j<par->jmax; j++) {	    if (cells[i][j] != FREE) {		to << cells[i][j]+1 << " ";	    }	    else {		to << 0 << " ";		    }	}	to << endl;    }}void grid::init() {    // generation of the initial random data distribution    // update of both grid cells and index      last = 0;    for (int i=0; i<par->binsize; ) {	int ipos = int(ran0(&idum)*par->imax);	int jpos = int(ran0(&idum)*par->jmax);       	if (cells[ipos][jpos] == FREE) {	    cells[ipos][jpos] = i;	    index[i].set(ipos, jpos);	    i++;	}     }    items = par->binsize;    }/* write-read access to grid cells  */int &grid::operator()(int i, int j) {    return cells[i][j];}/* return the position of a random element    from the index of free data items */const position<int> & grid::next() {    last = int(ran0(&idum) * items);        return index[last];}/* advise data item <data> to the grid at position <pos> */void grid::add(position<int> & pos, const int data) {    /* copy data identifier to grid */    cells[int(pos[0])][int(pos[1])] = data;    /* update index */    index[items++] = pos;}/* same as add, without index update */const void grid::replace(position<int> & pos, int data, int replaced) {    cells[int(pos[0])][int(pos[1])] = data;}const void grid::remove(position<int> & pos) {        // update index (last points to the last element that    // was returned by next(). This is the one to be replaced.    // This isn't very elegant, I know :-(        index[last] = index[--items];    // mark cell as free    cells[int(pos[0])][int(pos[1])] = FREE;    return;}const int grid::remove_init(position<int> * pos) {     // notice that the index must be treated different    // from a simple remove    (*pos) = index[--items];        int temp = cells[int((*pos)[0])][int((*pos)[1])];    // mark cell as free    cells[int((*pos)[0])][int((*pos)[1])] = FREE;    return temp;}// neighborhood function f providing environment informationconst double grid::fold(const int data, position<int> & pos, const int radius, const int nsize, double scalefactor) {    double sum = 0.0;    for (int i=int(pos[0]-radius); i<=pos[0]+radius; i++) {	for (int j=int(pos[1]-radius); j<=pos[1]+radius; j++) {	    int ii = checkborder(i, par->imax);	    int jj = checkborder(j, par->jmax);	 	    if ((cells[ii][jj] != FREE) && (cells[ii][jj] != data)) {		double d = dissimilarity(data, cells[ii][jj]) / scalefactor;		sum += 1.0 - d;	    }	}    }        if (sum > 0) return sum/nsize;    else return 0.0;} /* neighborhood function f providing environment information */const double grid::f(const int data, position<int> & pos, const int radius, const int nsize, double scalefactor, int clusterphase) {    double sum = 0.0;     double ctr = 0.0;    double weight = 0.0;    int dx, dy;    int incx, incy;    incy = -1;    dy = radius;    for (int i=int(pos[0]-radius); i<=pos[0]+radius; i++) {	dx = radius;	incx = -1;	for (int j=int(pos[1]-radius); j<=pos[1]+radius; j++) {	    int ii = checkborder(i, par->imax);	    int jj = checkborder(j, par->jmax);	 	    if ((cells[ii][jj] != FREE) && (cells[ii][jj] != data)) {		double d;		if (par->weighting == FALSE) {		    d = (1.0 - dissimilarity(data, cells[ii][jj]) / scalefactor);		    if (d <= 0 ) return 0.0;		    ctr += 1.0; 			}		else {		    weight = exp(-double(max(dx, dy))/double(radius));		    d = (1.0-dissimilarity(data, cells[ii][jj]) / scalefactor) / weight; 		    ctr += weight;		}		sum += d;	    }	    if (dx == 0) incx = 1;	    dx += incx;	}	if (dy == 0) incy = 1;	dy += incy;    }   if (clusterphase == TRUE) {       if (sum > 0) return sum/double(nsize);	 else return 0.0;    }    else {	if (sum > 0) return sum/double(ctr);	else return 0.0;	}}/* access to precomputed distance matrix */const double grid::dissimilarity(const int data1, const int data2) {    return bin->precomputed_d(data1, data2);}/* generate random position within a specified radius of the current position */void grid::generate(position<int> & pos, const int radius, position<int> & temp) {    int dx, dy;    double xt = (ran0(&idum)*2.0-1.0)*radius;    double yt = (ran0(&idum)*2.0-1.0)*radius;    dx = int((xt < 0) ? floor(xt) : ceil(xt));    dy = int((yt < 0) ? floor(yt) : ceil(yt));    temp = pos;    temp.add(dx,dy,par->imax,par->jmax);}void grid::generate(position<double> & pos, const int radius, position<double> & temp) {    int dx, dy;    double xt = (ran0(&idum)*2.0-1.0)*radius;    double yt = (ran0(&idum)*2.0-1.0)*radius;    dx = int((xt < 0) ? ceil(xt) : floor(xt));    dy = int((yt < 0) ? ceil(yt) : floor(yt));    temp = pos;    temp.add((int)dx,(int)dy,par->imax,par->jmax);}/* search for dropping position near <pos> */const position<int> grid::searchdroppingsite(position<int> & pos, const int data) {    double radius = 1.0;    position<int> temp = pos;    int ctr = 0;        /* try each radius five times */    while (cells[int(temp[0])][int(temp[1])] != FREE) {	generate(pos, (int)radius, temp);	radius += 0.2;	ctr++;    }    return temp;}/* Search for a free dropping site near the grid position <pos> */const position<double> grid::searchdroppingsite(position<double> & pos, const int data) {      double radius = 1.0;    position<double> temp = pos;    int ctr = 0;        /* try each radius five times */    while (cells[int(temp[0])][int(temp[1])] != FREE) {	generate(pos, (int)radius, temp);	radius += 0.2;	ctr++;    }        return temp;}/* border checks to stay on toroidal grid */const int grid::checkborder(const int p, const int lim) {    if (p >= lim) {	return p % lim;    }    if (p < 0) {	int temp = p%lim;	if (temp != 0) return lim + (p%lim);	else return 0;    }    return p;}double grid::square(double x) {    return x*x;}/* Euclidean distance for integer positions */double grid::spatialdistance(position<int> * p1, position<int> * p2, int torus) {    int dx = abs((*p1)[0]-(*p2)[0]);    if (torus == TRUE) dx = min(dx, par->imax - dx);    int dy = abs((*p1)[1]-(*p2)[1]);    if (torus == TRUE) dy = min(dy, par->jmax - dy);    return sqrt(double(square(dx)+square(dy)));}/* Euclidean distance for double positions */double grid::spatialdistance(position<double> * p1, position<double> * p2, int torus) {    double dx = abs((*p1)[0]-(*p2)[0]);    if (torus == TRUE) dx = min(dx, par->imax - dx);    double dy = abs((*p1)[1]-(*p2)[1]);    if (torus == TRUE) dy = min(dy, par->jmax - dy);    return sqrt(double(square(dx)+square(dy)));}double max(double i, double j) {    if (i > j) return i;    else return j;}double min(double i, double j) {    if (i < j) return i;    else return j;}int max(int i, int j) {    if (i > j) return i;    else return j;}int min(int i, int j) {    if (i < j) return i;    else return j;}double abs(double x) {    if (x < 0) return -x;    else return x;}