ctmq_batch2.c

统计模式识别算法包

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "memoryutilf.h"
#define RAND_MAX 2147483647




/* CTM code */



void presskey()

{
printf("press any key:");
getchar();
putchar('\n');

}


/*	given a coordinate index, this function places the value of the next coordinate
	in the standard ordering.	*/
	
void next_coord(coords,dim,units_p_dim,index)

int **coords;
int dim,units_p_dim,index;

{

int i,carry;

carry=1;
for (i = 1;i<=dim; i++)
	{
	
	coords[index+1][i]=coords[index][i]+carry;
	carry=coords[index+1][i] / (units_p_dim+1);
	if(carry>0)
		coords[index+1][i]=1;
		
	}
	
}


/* this function fills a coord matrix using the specified parameters */

void create_coords(coords, dim, units_p_dim,total_units)

int **coords;
int dim,units_p_dim,total_units;

{

int i,j;
	
		for (j = 1; j <= dim; j++)
			coords[1][j]=1;	/* initialize first unit coord */
		
	for (i = 1; i<=total_units-1; i++)
		{
		next_coord(coords,dim,units_p_dim,i);
		}
}


void create_map(map, dim, units_p_dim, ivars,total_units, coords)

int dim,units_p_dim,ivars,total_units;
int **coords;
float **map;


{

	int i,j;
for (i = 1; i<= total_units; i++)
	{
	
	for ( j = 1; j<= dim; j++)
		map[i][j]=2*(coords[i][j]-(.5*units_p_dim)-.5)/(units_p_dim-1);
	
	for (j=dim+1; j<=ivars; j++)
		map[i][j]=0;

	/*
	for(j=1; j<=ivars; j++)
		printf("%f,",map[i][j]);
	putchar('\n');
	*/
	
	}
}

void create_a(a, ivars,total_units)

int ivars,total_units;
float **a;


{

	int i,j;
for (i = 1; i<= total_units; i++)
	for (j=1; j<=ivars+1; j++)
		a[i][j]=0;
	
	
}

void create_win_list(win_list,map,data,ivars,total_units,samples)
int *win_list,total_units,samples,ivars;
float **map,**data;

{
int i,find_winner();

for (i=1; i<=samples; i++)
	win_list[i]=find_winner(map,data,i,ivars,total_units);
}





/* this function finds the winning unit its a little faster */

int find_winner(map,data,point,ivars,total_units)

float **map,**data;
int point,ivars, total_units;

{
	int i,j,best;
	float sum,diff,bestd;
	
	best=1;
	bestd=0.0;
	
	for (j = 1; j <= ivars ; j++)
			{
			diff=map[1][j]-data[point][j];
			bestd=bestd+diff*diff;
			}
	for (i = 2; i<=total_units; i++)
		{
		sum=0;
		j=1;
		while(j <= ivars && sum <= bestd)
			{
			diff=map[i][j]-data[point][j];
			sum=sum+diff*diff;
			j++;
			}
		if(sum<bestd)
			{
			bestd=sum;
			best=i;
			}
		}
		
	return(best);
}



/* this function finds the winning unit with variable weighing its a little faster */

int find_winnerw(map,data,weights,point,ivars,total_units)

float **map,**data,*weights;
int point,ivars, total_units;

{
	int i,j,best;
	float sum,diff,bestd;
	
	best=1;
	bestd=0.0;
	
	for (j = 1; j <= ivars ; j++)
			{
			diff=map[1][j]-data[point][j];
			bestd=bestd+diff*diff*weights[j];
			}
	for (i = 2; i<=total_units; i++)
		{
		sum=0;
		j=1;
		while(j <= ivars && sum <= bestd)
			{
			diff=map[i][j]-data[point][j];
			sum=sum+diff*diff*weights[j];
			j++;
			}
		if(sum<bestd)
			{
			bestd=sum;
			best=i;
			}
		}
		
	return(best);
}





/* this routine finds the winner searching in neighborhood of map_winner from
last time */

int find_winner2(map,coords,data,map_winner,point,ivars,total_units,map_dim,nb)

float **map, **data;
int ivars,total_units,point,map_winner,**coords,nb,map_dim;

{
float bestd,diff,sum;

int i,j,k,ma=ivars+1,best,disti;


	
best=map_winner;
bestd=0.0;
	
for (j = 1; j <= ivars ; j++)
		{
		diff=map[best][j]-data[point][j];
		bestd=bestd+diff*diff;
		}

for(i=1; i<=total_units; i++)
	{
	
	/* check if neighbor to map_winner */
	disti=0;
	k=1;
	while(k <=map_dim && disti<=nb)
		{
		disti += abs(coords[map_winner][k]-coords[i][k]);
		k++;
		}
		
		
	if(disti<=nb)
		{
		sum=0;
		j=1;
		while(j <= ivars && sum <= bestd)
			{
			diff=map[i][j]-data[point][j];
			sum=sum+diff*diff;
			j++;
			}
		if(sum<bestd)
			{
			bestd=sum;
			best=i;
			}
		}
	}

return(best);

}







/* this routine finds the winner based on pred error in neighborhood of map_winner */

int find_winner_pe2(map,a,coords,data,point,ivars,map_winner,total_units,map_dim)

float **map,**a, **data;
int ivars,total_units,point,map_winner,**coords;

{
float d,bestf,diff,goodf,dist_a,dist_b,dist_c;

int i,j,k,ma=ivars+1,besti,diffi,disti;


diff=a[map_winner][1]-data[point][ma];
dist_b=0.0;
for(j=1; j<=ivars; j++)
	{
	diff += data[point][j]*a[map_winner][j+1];
	d=data[point][j]-map[map_winner][j];
	dist_b += d*d;
	}
	

goodf=bestf=fabs(diff);
besti=map_winner;


for(i=1; i<=total_units; i++)
	{
	
	/* check if neighbor to map_winner */
	disti=0;
	k=1;
	while(k <=map_dim && disti<2)
		{
		disti += abs(coords[map_winner][k]-coords[i][k]);
		k++;
		}
		
		
	if(disti<2)
		{
		dist_a=0.0;
		dist_c=0.0;
		diff=a[i][1]-data[point][ma];
		for(j=1; j<=ivars; j++)
			{
			diff += data[point][j]*a[i][j+1];
			d=map[map_winner][j]-map[i][j];
			dist_a += d*d;
			d=data[point][j]-map[i][j];
			dist_c += d*d;
			}
		diff=fabs(diff);
		
		if(diff<bestf && dist_b<dist_a && dist_c<dist_a)
			
			{
			bestf=diff;
			besti=i;
			/*printf("better\n");*/
			}
		}
	}

return(besti);

}




/* this function gives a weighing for a data point based on measuring the
second derivative */

float curve_weight(map,a,coords,data,point,ivars,map_winner,total_units,map_dim)

float **map,**a,**data;
int **coords,point,ivars,map_winner,total_units,map_dim;

{
float bestf,diffx,diffy,goodf,xsum,ysum,slopem;

int i,j,k,count,ma=ivars+1,besti,diffi,disti;



slopem=0;
count=0;
for(i=1; i<=total_units; i++)
	{
	
	/* check if neighbor to map_winner */
	disti=0;
	k=1;
	while(k <=map_dim && disti<=1)
		{
		disti += abs(coords[map_winner][k]-coords[i][k]);
		k++;
		}
		
		
	if(disti==1)
		{
		xsum=ysum=0;
		count++;
		for(j=2;j<=ma;j++)
			{
			diffy= a[i][j]-a[map_winner][j];
			/*diffx= map[i][j-1]-map[map_winner][j-1];
			xsum += diffx*diffx; */
			ysum += diffy*diffy;
			}
		/*slopem += ysum*ysum;*/
		slopem += sqrt(ysum);

		}
	}
/*slopem=atan(slopem/((float)count));*/
slopem=slopem/((float)count);
slopem=slopem*slopem;
return(slopem);

}



float error_weight(data,point,ivars,a,unit,dd,diag,atb,workspace)
float **data,**a,**dd,*diag,*atb,*workspace;
int point,ivars,unit;

{
float fact,fit1,e;

int i,j,k,ma=ivars+1;

atb[1]=1;
fit1=a[unit][1]-data[point][ma];
for(j=1; j<=ivars; j++)
	{
	fit1 += data[point][j]*a[unit][j+1];
	atb[j+1]=data[point][j];
	}
/*print_mat(dd,1,ma,1,ma);*/
cholsl(dd,ma,diag,atb,workspace);
fact=1;
for(j=1;j<=ma;j++)
	{
	fact -= workspace[j]*atb[j];
	/*printf("%f,",workspace[j]);*/
	}
/*putchar('\n');*/
if(fact<=0.0)
	e=-1.0;
else
	e=fit1/fact;
	
return(e);	

}



/* this function makes predictions, piecewise lin, original y in data is replaced
   and a weighing is used to find winning units */

float predict_plw(map,a,weights,data,ivars,total_units,samples)

float **map,**a, **data, *weights;
int ivars,total_units,samples;

{

int i,j,winner,ma=ivars+1;
float mean,std,rms_error,tmp;

mean=0;
for(i=1; i<=samples; i++)
	mean += data[i][ma];
mean=mean/((float)samples);

rms_error=0;
std=0;
for(i=1; i<=samples; i++)
	{
	tmp=data[i][ma];
	std += (tmp-mean)*(tmp-mean);
	winner=find_winnerw(map,data,weights,i,ivars,total_units);
	data[i][ma]=a[winner][1];
	for(j=1; j<=ivars; j++)
		data[i][ma] += data[i][j]*a[winner][j+1];
	rms_error += (data[i][ma]-tmp)*(data[i][ma]-tmp);
	}
std=sqrt(std/((float)samples));
rms_error=sqrt(rms_error/((float)samples));
printf("mean(y)=%f std(y)=%f rms error=%f ",mean,std,rms_error);
rms_error=rms_error/std;
printf("normalized rms err=%f\n",rms_error);
return(rms_error);

}

float rms_plw(map,a,weights,data,ivars,total_units,samples)

float **map,**a, **data, *weights;
int ivars,total_units,samples;

{

int i,j,winner,ma=ivars+1;
float rms_error,tmp;

rms_error=0;
for(i=1; i<=samples; i++)
        {
        winner=find_winnerw(map,data,weights,i,ivars,total_units);
        tmp=a[winner][1]-data[i][ma];
        for(j=1; j<=ivars; j++)
                tmp += data[i][j]*a[winner][j+1];
        rms_error += tmp*tmp;
        }
rms_error=sqrt(rms_error/((float)samples));
return(rms_error);

}

/* this function makes predictions, piecewise lin, original y in data is replaced */

float predict_pl(map,a,data,ivars,total_units,samples)

float **map,**a, **data;
int ivars,total_units,samples;

{

int i,j,winner,ma=ivars+1;
float mean,std,rms_error,tmp;

mean=0;
for(i=1; i<=samples; i++)
	mean += data[i][ma];
mean=mean/((float)samples);

rms_error=0;
std=0;
for(i=1; i<=samples; i++)
	{
	tmp=data[i][ma];
	std += (tmp-mean)*(tmp-mean);
	winner=find_winner(map,data,i,ivars,total_units);
	data[i][ma]=a[winner][1];
	for(j=1; j<=ivars; j++)
		data[i][ma] += data[i][j]*a[winner][j+1];
	rms_error += (data[i][ma]-tmp)*(data[i][ma]-tmp);
	}
std=sqrt(std/((float)samples));
rms_error=sqrt(rms_error/((float)samples));
printf("mean=%f std=%f err=%f ",mean,std,rms_error);
rms_error=rms_error/std;
printf("serr=%f\n",rms_error);
return(rms_error);

}







/* this function makes predictions, piecewise const, original y in data is replaced */

void predict_pc(map,a,data,ivars,total_units,samples)

float **map,**a, **data;
int ivars,total_units,samples;

{

int i,j,winner,ma=ivars+1;





for(i=1; i<=samples; i++)
	{
	winner=find_winner(map,data,i,ivars,total_units);
	data[i][ma]=a[winner][1];
	for(j=1; j<=ivars; j++)
		data[i][ma] += map[winner][j]*a[winner][j+1];
	}



}






void som(map,coords,data,win_list,sigmax,dim,ivars,u_p_dim,total_units,samples)

float **map,**data,sigmax;
int **coords,*win_list;
int dim,ivars,total_units,samples,u_p_dim;


{
	int k,j,i;
	int winner,dist,diff,ndist;
	float weightx,weight_sumx,factorx;
		
	
	factorx=-.5/(sigmax*sigmax);	
	ndist=2*sigmax*sigmax;
	if (ndist<2)
		ndist=2;
	if(ndist<u_p_dim)
		for(k=1; k<=samples; k++)
			win_list[k]=find_winner2(map,coords,data,win_list[k],
				k,ivars,total_units,dim,ndist);
	else
		for(k=1; k<=samples; k++)
			win_list[k]=find_winner(map,data,k,ivars,total_units);
	
	for(j=1; j<=total_units; j++)
		{
		
		/* zero unit */
		for(k=1; k<=ivars; k++)
			map[j][k]=0;
		
		weight_sumx=0;
		
		for(i=1; i<=samples; i++)
			{
			
			dist=0;
			winner=win_list[i];
			for (k=1; k<=dim; k++)
				{
				diff=coords[winner][k]-coords[j][k];
				dist=dist+diff*diff;
				}
			weightx=exp(dist*factorx);
			if(weightx>.01)
				{
				weight_sumx += weightx;
				for(k=1; k<=ivars; k++)
					map[j][k] += data[i][k]*weightx;
				}
				
			}
		/* normalize map unit */
		weight_sumx =1.0/weight_sumx;
		for(k=1; k<=ivars; k++)
			map[j][k] = map[j][k]*weight_sumx;
			
		
		}	
}




/*      this function performs 2 step CTM in batch mode with local lin. fits
        The map has units x1,x2,..,xn  and
        the 0th and 1st order coef are stored in a1,a2,..,an+1
*/

        /* For fixed map positions, find optimal linear fits with neighborhood
           cutoff  uses faster Cholesky decomp */


float update_fit3(map,a,coords,data,win_list,sigmay,
        dim,ivars,total_units,samples,u,b,dd,diag,workspace,atb)

float **map,**data,**a,*b,**u,**dd,*diag,*workspace,*atb,sigmay;
int **coords,*win_list;
int dim,ivars,total_units,samples;


{

        int i,j,k,p,winner,dist,diff,s_count;
        float factory,c,distf,sum,d;
        int d_count;
        float weight,rms_error,error_weight();
        float wmax,thresh,press;
        int ma=ivars+1;

        factory=-.5/(sigmay*sigmay);

        /* For fixed map positions, find optimal linear fits */
        /* get winner list */

        for(i=1; i<=samples; i++)
                win_list[i]=winner=find_winner(map,data,
                        i,ivars,total_units);

        /*determine weights for regression for each unit */

        press=0;
        s_count=0;
        for(j=1; j<=total_units; j++)