/*
 * Search space for pattern closest to given pattern.
 */

#include <stdio.h>
#include <stdlib.h>
#include <memory.h>
#include <math.h>
#include <time.h>

#define RADIUS .9	/* space configuration parameter */

/* search states */
#define DISTPENDING 0	/* distance pending */
#define DISTDONE 1		/* distance computed */
#define EXPANDED 2		/* pattern expanded */
#define SRCHDONE 3		/* pattern searched */
#define EMPTY (struct SRCHWORK *)(-1)

/* pattern */
struct PATTERN
{
	struct PATTERN *childlist;	/* child pattern list */
	struct PATTERN *childlast;	/* last child */
	struct PATTERN *sibnext;	/* next sibling */
	struct PATTERN *sibback;	/* previous sibling */
	unsigned char *vals;		/* pattern values */
	double distance;			/* distance from child to parent */
};
struct PATTERN *EntPatt;	/* space entry pattern */

/* pattern search element */
struct SRCHWORK
{
	struct PATTERN *pattern;	/* pattern */
	double distance;			/* comparison distance */
	double workdist;			/* work distance */
	int state;					/* search state */
	struct SRCHWORK *childlist;	/* children */
	struct SRCHWORK *sibnext;	/* next sibling */
	struct SRCHWORK *sibback;	/* previous sibling */
	struct SRCHWORK *srchnext;	/* next on search return list */
};

/* pattern search stack */
#define STKMEM 1024			/* stack malloc increment */
static struct SRCHSTK
{
	struct SRCHWORK *currsrch;
	struct SRCHWORK *child;
	struct SRCHWORK *childnext;
} *SrchStk;
static struct SRCHSTK *SrchStkp,*SrchStkEnd;	/* stack pointers */
static int SrchStkSz;				/* current stack size */

/* SrchPatt data */
struct SRCHWORK *SrchList;	/* pattern search return list */
int SearchCount, BestSearch;

/* search work memory */
#define SRCHWORKMEM 100000		/* search work memory size */
struct SRCHWORK *SrchWork;		/* search work */
int SrchWorkUse;

/* functions */
void AddPatt(struct PATTERN *, struct PATTERN *);
void SrchPatt(struct PATTERN *, struct PATTERN *, int, int);
double Dist(struct PATTERN *, struct PATTERN *);
void FoundPatt(struct SRCHWORK *, int, int *, struct SRCHWORK **);
struct SRCHWORK *GetSrchWork();
void ClrSrchWork();

/* test parameters */
#define NUM_VALS 50
#define NUM_PATTERNS 100000
#define MAX_FIND 1
#define NUM_SEARCH NUM_PATTERNS
void PrintPattern(struct PATTERN *);

int
main(argc,argv)
int argc;
char *argv[];
{
	int i,j;
	struct PATTERN *p,*s,*t;
	double d,d2;

	srand(time(NULL));

	/* initialize search pattern */
	s = (struct PATTERN *)malloc(sizeof(struct PATTERN));
	s->sibnext = s->sibback = NULL;
	s->childlist = s->childlast = NULL;
	s->vals = (unsigned char *)malloc(NUM_VALS);
	for (j = 0; j < NUM_VALS; j++)
	{
		s->vals[j] = (rand() % 27) + 'a';
	}
	printf("Search pattern:\n");
	PrintPattern(s);

	/* add patterns to space */
	d = -1.0;
	EntPatt = NULL;
	for (i = 0; i < NUM_PATTERNS; i++)
	{
		p = (struct PATTERN *)malloc(sizeof(struct PATTERN));
		p->sibnext = p->sibback = NULL;
		p->childlist = p->childlast = NULL;
		p->vals = (unsigned char *)malloc(NUM_VALS);
		for (j = 0; j < NUM_VALS; j++)
		{
			p->vals[j] = (rand() % 26) + 'a';
		}
		if (EntPatt == NULL)
		{
			EntPatt = p;
		} else {
			AddPatt(EntPatt, p);
		}
		d2 = Dist(s,p);
		if (d < 0.0 || d2 < d) d = d2;
	}

	/* add a search target pattern */
	t = (struct PATTERN *)malloc(sizeof(struct PATTERN));
	t->sibnext = t->sibback = NULL;
	t->childlist = t->childlast = NULL;
	t->vals = (unsigned char *)malloc(NUM_VALS);
	for (j = 0; j < NUM_VALS; j++)
	{
		t->vals[j] = s->vals[j];
/*		if ((rand() % 10) == 0) t->vals[j] = (rand() % 26) + 'a'; */
		if ((rand() % 5) == 0)
		{
			if (rand() % 2) t->vals[j] += rand() % 5; else t->vals[j] -= rand() % 5;
		}
	}
	AddPatt(EntPatt, t);

	/* allocate search working memory */
	if ((SrchWork = (struct SRCHWORK *)malloc(SRCHWORKMEM*sizeof(struct SRCHWORK))) == NULL)
	{
		fprintf(stderr,"Cannot malloc search work elements");
		exit(1);
	}
	memset(SrchWork,0,SRCHWORKMEM*sizeof(struct SRCHWORK));
	SrchWorkUse = 0;
	if ((SrchStk = (struct SRCHSTK *)malloc(sizeof(struct SRCHSTK)*STKMEM)) == NULL)
	{
		fprintf(stderr,"Cannot malloc recursion memory");
		exit(1);
	}
	SrchStkSz = STKMEM;
	SrchStkEnd = (struct SRCHSTK *)((int)SrchStk+(sizeof(struct SRCHSTK)*SrchStkSz));

	/* search */
	SrchList = NULL;
	SrchPatt(s,EntPatt,MAX_FIND,NUM_SEARCH);
	
	/* print search results */
	printf("Result:\n");
	PrintPattern(SrchList->pattern);
	printf("Distance=%f\n", Dist(s,SrchList->pattern));

	return(0);
}

/* print pattern */
void PrintPattern(struct PATTERN *p)
{
	int j;
	for (j = 0; j < NUM_VALS; j++)
	{	
		printf("%c", p->vals[j]);
	}
	printf("\n");
}

/* add a pattern to the space */
void AddPatt(curr,new)
register struct PATTERN *curr;		/* current pattern */
register struct PATTERN *new;		/* new pattern */
{
	register struct PATTERN *p,*p2,*p3;
	register double dcn,dnn;

	/* clear children links */
	new->childlist = NULL;
	new->childlast = NULL;

	/* add pattern to first acceptable branch */
	dcn = Dist(curr,new);
	while (1)
	{
		for (p = curr->childlist; p != NULL; p = p->sibnext)
		{
			/* check relative distances */
			dnn = Dist(p,new);
			if (dnn <= (p->distance*RADIUS))
			{
				curr = p;	/* change current fragment */
				dcn = dnn;
				break;
			}
		}
		if (p == NULL) { break; }
	}

	/* link new as child of current pattern */
	new->distance = dcn;
	new->sibnext = NULL;
	new->sibback = curr->childlast;
	if (curr->childlast != NULL)
	{
		curr->childlast->sibnext = new;
	} else {
		curr->childlist = new;
	}
	curr->childlast = new;

	/*
	 * Check if previously added patterns should be un-linked from the
	 * current pattern and linked as children of the new pattern.
	 */
	for (p = curr->childlist; p != new && p != NULL;)
	{
		dnn = Dist(p,new);

		/* if should be linked to new pattern */
		if (dnn <= (new->distance*RADIUS))
		{
			/* re-link children */
			p2 = p->sibnext;
			if (p->sibback != NULL)
			{
				p->sibback->sibnext = p2;
			} else {
				curr->childlist = p2;
			}
			if (p2 != NULL)
			{
				p2->sibback = p->sibback;
			}

			/* convert child sub-tree to list */
			for (p2 = p3 = p, p->sibnext = NULL;; p3 = p3->sibnext)
			{
				for (; p2->sibnext != NULL; p2 = p2->sibnext) {}
				for (; p3 != NULL && p3->childlist == NULL; p3 = p3->sibnext) {}
				if (p3 == NULL) { break; }
				p2->sibnext = p3->childlist;
			}

			/* add list to current pattern */
			for (p2 = p; p2 != NULL; p2 = p3)
			{
				p3 = p2->sibnext;
				AddPatt(curr,p2);
			}

			/* restart check (since child configuration may have changed) */
			for (p = curr->childlist; p != new && p != NULL; p = p->sibnext) {}
			if (p == NULL) { break; }
			p = curr->childlist;
		} else {
			p = p->sibnext;
		}
	}
}

/* distance function */
double Dist(p1,p2)
struct PATTERN *p1,*p2;
{
	register int i;
	double dist,d;

	for (i = 0, dist = 0.0; i < NUM_VALS; i++)
	{
		d = (double)(p1->vals[i] - p2->vals[i]);
		if (d < 0.0) d = -d;
		dist += d;
	}
	return(dist);
}

/* search space for patterns closest to the given pattern */
/* put best matches on SrchList, and the best match in SrchBest */
void SrchPatt(srchpatt,entpatt,maxfind,maxsrch)
register struct PATTERN *srchpatt;	/* search for this pattern */
register struct PATTERN *entpatt;	/* search space at this entry point */
register int maxfind;				/* maximum patterns to find and return */
register int maxsrch;				/* maximum patterns to search */
{
	register int numsrch;
	register struct SRCHSTK *stkp,*stkp2;
	register struct PATTERN *p;
	register struct SRCHWORK *sw,*sw2,*bsw,*bsw2;
	struct SRCHWORK *swcut;
	int numfind;

	/* initialize */
	SearchCount = BestSearch = 0;
	ClrSrchWork();
	numsrch = numfind = 0;
	SrchList = swcut = NULL;
	if (numsrch >= maxsrch) return;
	if (entpatt == NULL) return;
	SrchStkp = SrchStk;
	SrchStkp->currsrch = GetSrchWork();
	SrchStkp->currsrch->pattern = entpatt;
	SrchStkp->currsrch->distance = Dist(SrchStkp->currsrch->pattern,srchpatt);
	SrchStkp->currsrch->state = DISTDONE;
	FoundPatt(SrchStkp->currsrch,maxfind,&numfind,&swcut);
	numsrch++;
	if (numsrch >= maxsrch) return;

	/* for each level of recursion */
	while (SrchStkp >= SrchStk)
	{
		/* find best and next best distances for current search branch */
		for (stkp = SrchStkp; stkp >= SrchStk; stkp--)
		{
			/* expand pattern? */
			if (stkp->currsrch->state == DISTDONE)
			{
				for (p = stkp->currsrch->pattern->childlist, sw2 = NULL; p != NULL; p = p->sibnext)
				{
					sw = GetSrchWork();
					sw->pattern = p;
					sw->state = DISTPENDING;
					if (sw2 == NULL)
					{
						stkp->currsrch->childlist = sw;
					} else {
						sw2->sibnext = sw;
						sw->sibback = sw2;
					}
					sw2 = sw;
				}
				stkp->child = stkp->childnext = NULL;
				stkp->currsrch->state = EXPANDED;
			}

			/* best and next best distances must be (re)computed? */
			if (stkp->child == NULL ||
			   (stkp->child != NULL &&
			    stkp->child->workdist > 0.0 &&
			    stkp->childnext == EMPTY) ||
			   (stkp->child != NULL &&
			    stkp->childnext != NULL &&
			    stkp->childnext != EMPTY &&
			    stkp->child->workdist > stkp->childnext->workdist))
			{
				bsw = bsw2 = NULL;
				for (sw = stkp->currsrch->childlist; sw != NULL; sw = sw2)
				{
					sw2 = sw->sibnext;

					/* have best possible child? */
					if (bsw != NULL && bsw->workdist == 0.0)
					{
						if (sw2 == NULL)
						{
							bsw2 = NULL;
						} else {
							bsw2 = EMPTY;
						}
						break;
					}

					/* compute distance? */
					if (sw->state == DISTPENDING)
					{
						sw->distance = Dist(sw->pattern,srchpatt);
						if ((sw->workdist = sw->distance -
						    (sw->pattern->distance*RADIUS)) < 0.0)
						{
							sw->workdist = 0.0;
						}
						sw->state = DISTDONE;

						/* save pattern on return list */
						FoundPatt(sw,maxfind,&numfind,&swcut);
						numsrch++;

						/* check for termination of search */
						if (numsrch >= maxsrch)
						{
							printf("Search %d, best at %d\n",SearchCount,BestSearch);
							return;
						}
					}

					/* cut off infeasible or finished branch */
					if (sw->state == SRCHDONE || 
					   (swcut != NULL && sw->workdist >= swcut->distance))
					{
						sw->state = SRCHDONE;	/* cut off */
						if (sw->sibback == NULL)
						{
							stkp->currsrch->childlist = sw->sibnext;
						} else {
							sw->sibback->sibnext = sw->sibnext;
						}
						if (sw->sibnext != NULL)
						{
							sw->sibnext->sibback = sw->sibback;
						}
						continue;
					}

					/* find best and next best child branches */
					if (bsw == NULL || bsw->workdist > sw->workdist)
					{
						bsw2 = bsw;
						bsw = sw;
						continue;
					}
					if (bsw2 == NULL || bsw2->workdist > sw->workdist)
					{
						bsw2 = sw;
						continue;
					}
				}

				/* change to better branch level? */
				if (stkp->child != bsw)
				{
					SrchStkp = stkp;
				}
				stkp->child = bsw;
				stkp->childnext = bsw2;
			}

			/* finished with this level? */
			if (stkp->child == NULL)
			{
				stkp->currsrch->state = SRCHDONE;
				SrchStkp = stkp-1;
				if (stkp > SrchStk)
				{
					SrchStkp->child = NULL;
				}

			} else {

				/* set new best distance for branch */
				for (stkp2 = stkp; stkp2 >= SrchStk; stkp2--)
				{
					stkp2->currsrch->workdist = stkp->child->workdist;
				}
			}
		}

		/* expand search deeper */
		if (SrchStkp >= SrchStk)
		{
			SrchStkp++;
			SrchStkp->currsrch = (SrchStkp-1)->child;
			SrchStkp->child = SrchStkp->childnext = NULL;
		}
	}
	printf("Search %d, best at %d\n",SearchCount,BestSearch);
}

/* add pattern to the found list */
void FoundPatt(swfound,maxfind,numfind,swcut)
register struct SRCHWORK *swfound;
register int maxfind;
register int *numfind;
register struct SRCHWORK **swcut;
{
	register struct SRCHWORK *sw,*sw2;

	for (sw = SrchList, sw2 = NULL; sw != NULL; sw2 = sw, sw = sw->srchnext)
	{
		if (sw->distance <= swfound->distance) { break; }
	}
	SearchCount++;
	if (sw == NULL) BestSearch = SearchCount;
	if (sw2 == NULL)
	{
		if (*numfind < maxfind)
		{
			(*numfind)++;
			swfound->srchnext = SrchList;
			SrchList = swfound;
		}
	} else {
		swfound->srchnext = sw2->srchnext;
		sw2->srchnext = swfound;
		if (*numfind < maxfind)
		{
			(*numfind)++;
		} else {
			SrchList = (SrchList)->srchnext;
		}
	}
	if (*numfind == maxfind)   /* set cut off */
	{
		*swcut = SrchList;
	}
}

/* get a pattern search element */
struct SRCHWORK *GetSrchWork()
{
	register struct SRCHWORK *sw;
	
	if (SrchWorkUse == SRCHWORKMEM)
	{
		fprintf(stderr, "Search work elements exhausted");
		exit(1);
	}
	sw = &SrchWork[SrchWorkUse];
	SrchWorkUse++;

	return(sw);
}

/* clear (free all) search work elements */
void ClrSrchWork()
{
	memset(SrchWork,0,sizeof(struct SRCHWORK)*SRCHWORKMEM);
	SrchWorkUse = 0;
}