vectorsimplex.txt

vectorsimplex in c, please enjoy it.

/* 驴誊弄润俐妨呵努步缄恕 Vector Simplex 恕 */
/* Proposed and Coded by T.Takahama, 1998 */

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

/* 誊弄簇眶 y=f(x) の年盗: xは N 肌傅、yは M 肌傅 */
#define N		2	/* 掐蜗の肌傅 */
#define M		2	/* 叫蜗の肌傅 */
#define square(x)       ((x)*(x))

/* 誊弄簇眶 */
void f(double x[], double y[])
{
	/* f1=x_1^2+x_2^2 */
	y[0]=square(x[0])+square(x[1]);

	/* f2=(x_1-1)^2+(x_2-1)^2 */
	y[1]=square(x[0]-1.0)+square(x[1]-1.0);
}

/* Vector Simplex 恕のパラメ〖タ */
#define N_VERTEX	26	/* 暮爬の眶 */
#define S		5	/* 介袋帽挛の络きさ */
#define Alpha		1.0
#define Beta		0.5
#define Gamma		2.0

#ifndef RAND_MAX
#define RAND_MAX 0x7fffffff
#endif
#define rnd(low, high)  ((low)+(double)rand()/RAND_MAX*((high)-(low)))

/* 暮爬の年盗 */
typedef struct {
        double x[N];
        double y[M];
        int set;        /* 掳する礁圭: Uh, Us, Ul */
} Vertex;

enum { All, Uh='h', Us='s', Ul='l' };

Vertex Vs[N_VERTEX];

main()
{
	Initialize(Vs, N_VERTEX);
	vsimplex(Vs, N_VERTEX);
	Print(Vs, N_VERTEX);
}

/* 介袋暮爬の菇喇: (0,0)を奶る染仿S/2の亩靛惧の爬 */
Initialize(Vertex *V, int n)
{
	register int i, j;

	double r, s;

	s=sqrt((double)N);
	for(i=0; i<n; i++) {
		r=0;
		for(j=0; j<N; j++) {
			V[i].x[j]=rnd(-1.0, 1.0);
			r+=square(V[i].x[j]);
		}
		r=sqrt(r);
		for(j=0; j<N; j++) {
			V[i].x[j]*=0.5*S/r;
			V[i].x[j]+=0.5*S/s;
		}
	}
	for(i=0; i<n; i++)
		Evaluate(&V[i]);
}

/* ベクトル粗の络井簇犯 */
enum relation { LT, EQ, GT, Unknown, LE, GE };

/* Vector Simplex恕の悸乖 */
vsimplex(Vertex *V, int n)
{
	int h, l;
	double x0[N];
	Vertex Vr, Ve, Vc;

	for(;;) {
		PartitionHSL(V, n);
		h=SelectH(V, n);
		if(h<0) break;
		CenterU0(V, n, h, x0);
	/* reflection */
		Operation(Vr.x, N, 1+Alpha, x0, -Alpha, V[h].x);
		Evaluate(&Vr);
		if(Exists(V, n, Ul, GT, Vr.y)) {
	/* expansion */
			Operation(Ve.x, N, Gamma, Vr.x, 1-Gamma, x0);
			Evaluate(&Ve);
			if(Exists(V, n, Ul, GT, Ve.y))
				Copy(&V[h], &Ve);
			else
				Copy(&V[h], &Vr);
		}
		else if(!Exists(V, n, Ul, LT, Vr.y) ||
			Exists(V, n, Us, GE, Vr.y))
			Copy(&V[h], &Vr);
		else {
			if(Exists(V, n, Uh, GT, Vr.y))
				Copy(&V[h], &Vr);
	/* contraction */
			Operation(Vc.x, N, Beta, V[h].x, 1-Beta, x0);
			Evaluate(&Vc);
			if(!Exists(V, n, Ul, LT, Vc.y) ||
			   Exists(V, n, Uh, GT, Vc.y))
				Copy(&V[h], &Vc);
			else {
	/* reduction */
				l=SelectL(V, n, V[h].y);
				if(l<0)
					fprintf(stderr, "SelectL failed\n");
				Operation(V[h].x, N, 0.5, V[h].x, 0.5, V[l].x);
				Evaluate(&V[h]);
			}
		}
	}
}

/* 礁圭Ul, Us, Uh を疯年する */
PartitionHSL(Vertex *V, int n)
{
	register int i;

	for(i=0; i<n; i++)
		V[i].set=Us;
	for(i=0; i<n; i++)
		if(!Exists(V, n, All, LT, V[i].y))
			V[i].set=Ul;
	for(i=0; i<n; i++)
		if(V[i].set==Us && !Exists(V, n, Us, GT, V[i].y))
			V[i].set=Uh;
}

/* 回年した簇犯を塔颅する妥燎があるかどうかを冉年する */
int Exists(Vertex *V, int n, int set, int relation, double y[])
{
	register int i;
	int cmp;

	for(i=0; i<n; i++)
		if(set==All || V[i].set==set) {
			cmp=Compare(V[i].y, y, M);
			switch(relation) {
			case LT: if(cmp==LT) return(1); break;
			case GT: if(cmp==GT) return(1); break;
			case LE: if(cmp==LT || cmp==EQ) return(1); break;
			case GE: if(cmp==GT || cmp==EQ) return(1); break;
			case EQ: if(cmp==EQ) return(1); break;
			}
		}
	return(0);
}

/* ベクトルを孺秤し、>(GT), <(LT), =(EQ) or 稍汤(Unknown) を手す */
int Compare(double y1[], double y2[], int n)
{
	register int i;
	int status;

	status=EQ;
	for(i=0; i<n; i++)
		if(status==EQ) {
			if(y1[i]>y2[i]) status=GT;
			else if(y1[i]<y2[i]) status=LT;
		}
		else if(status==GT) {
			if(y1[i]<y2[i]) return(Unknown);      /* unknown */
		}
		else if(y1[i]>y2[i]) return(Unknown); /* unknown */
	return(status);
}

/* Uh から Xh を联买する。Uh が鄂ならば -1 を手す */
int SelectH(Vertex *V, int n)
{
	static int h=0;
	int h0;

	if(h>=n) h=0;
	h0=h;
	do {
		if(V[h].set==Uh)
			return(h++);
		h++;
		if(h>=n) h=0;
	} while(h!=h0);
	return(-1);
}

/* U0 を疯年し、哭看を滇める */
CenterU0(Vertex *V, int n, int h, double x0[])
{
	static int i=0;
	register int j, k;

	for(j=0; j<N; j++)
		x0[j]=0.0;
	for(k=0; k<N; k++) {
		if(i>=n) i=0;
		if(i!=h)
			for(j=0; j<N; j++) x0[j]+=V[i].x[j];
		i++;
	}
	for(j=0; j<N; j++)
		x0[j]/=N;
}

/* 称硷拎侯を悸乖する */
Operation(double x[], int n, double w1, double x1[], double w2, double x2[])
{
	register int j;

	for(j=0; j<n; j++)
		x[j]=w1*x1[j]+w2*x2[j];
}

/* 暮爬を删擦する */
Evaluate(Vertex *V)
{
	f(V->x, V->y);
}

/* 暮爬を弥き垂える */
Copy(Vertex *dest, Vertex *src)
{
	register int j;

	for(j=0; j<N; j++)
		dest->x[j]=src->x[j];
	for(j=0; j<M; j++)
		dest->y[j]=src->y[j];
}

/* Xh よりの惧疤の Xl を手す */
int SelectL(Vertex *V, int n, double y[])
{
	static int l=0;
	int l0;

	if(l>=n) l=0;
	l0=l;
	do {
		if(V[l].set==Ul && Compare(V[l].y, y, M)==LT)
			return(l++);
		l++;
		if(l>=n) l=0;
	} while(l!=l0);
	return(-1);
}

Print(Vertex *V, int n)
{
	register int i, j;

	for(i=0; i<n; i++) {
		 printf("%2d %c: ", i, V[i].set);
		for(j=0; j<N; j++)
			printf("%lf ", V[i].x[j]);
		printf(": ");
		for(j=0; j<M; j++)
			printf("%lf ", V[i].y[j]);
		printf("\n");
	}
}