📄 h2m_glvq_model.c
字号:
/*
Harmonic to Minimum Generalized Learning Vector Quantization (H2M-GLVQ) classification algorithm.
Usage
------
[Wproto_est , yproto_est , E_H2MGLVQ] = h2m_glvq_model(Xtrain , ytrain , [Wproto] , [yproto] , [options]);
Inputs
-------
Xtrain Train data (d x Ntrain)
ytrain Train label (1 x Ntrain) with m different classes
Wproto Initial prototypes weights (d x Nproto) (default Wproto (d x Nproto) where Nproto=round(sqrt(Ntrain)).
Wproto ~ N(m_{ytrain=i},sigma_{ytrain=i}), where m_{ytrain=i} = E[Xtrain|ytrain=i], sigma_{ytrain=i} = E[XtrainXtrain'|ytrain=i]
yproto Intial prototypes label (1 x Nproto) where card(yproto)=card(ytrain)
options Options'structure
.epsilonk Update weight's for common class label between prototype and train data (default = 0.01).
.epsilonl Update weight's for different class label between prototype and train data (default = 0.01).
.tmax Maximum temperature for cooling schudeling (default = 1)
.tmin Minimum temperature for cooling schudeling (default = 0.6)
.xi Constant in the sigmoid function (default = 0.1)
.nb_iterations Number of iterations (default = 1000)
.shuffle Randomly permute order of train data between each iteration if shuffle=1 (default = 1)
Outputs
-------
Wproto_est Final prototypes weigths (d x Nproto)
yproto_est Final prototypes labels (1 x Nproto)
E_H2MGLVQ Energy criteria versur iteration (1 x options.nb_iterations)
To compile
----------
mex -DranSHR3 -output h2m_glvq_model.dll h2m_glvq_model.c
mex -DranSHR3 -f mexopts_intel10amd.bat -output h2m_glvq_model.dll h2m_glvq_model.c
Example 1
---------
d = 2;
Ntrain = 300;
m = 2;
M0 = [0 ; 0];
R0 = [1 0 ; 0 1];
M1 = [2 ; 3];
R1 = [0.5 0.1 ; 0.2 1];
vect_test = (-4:0.1:8);
options.epsilonk = 0.01;
options.epsilonl = 0.01;
options.tmax = 1;
options.tmin = 0.5;
options.xi = 1;
options.nb_iterations = 2000;
options.shuffle = 1;
Xtrain = [M0(: , ones(1 , Ntrain/2)) + chol(R0)'*randn(d , Ntrain/2) , M1(: , ones(1 , Ntrain/2)) + chol(R1)'*randn(d , Ntrain/2)];
ytrain = [zeros(1 , Ntrain/2) , ones(1 , Ntrain/2)];
[X , Y] = meshgrid(vect_test);
Xtest = [X(:)' ; Y(:)'];
[Wproto_est , yproto_est , E_H2MGLVQ] = h2m_glvq_model(Xtrain , ytrain , [] , [] , options);
ytest_est = glvq_predict(Xtest , Wproto_est , yproto_est);
indtrain0 = (ytrain == 0);
indtrain1 = (ytrain == 1);
indproto0 = (yproto_est == 0);
indproto1 = (yproto_est == 1);
figure(1)
imagesc(vect_test , vect_test , reshape(ytest_est , length(vect_test) , length(vect_test)) )
axis ij
hold on
plot(Xtrain(1 , indtrain0) , Xtrain(2 , indtrain0) , 'k+' , Xtrain(1 , indtrain1) , Xtrain(2 , indtrain1) , 'm+' , Wproto_est(1 , indproto0) , Wproto_est(2 , indproto0) , 'ko' , Wproto_est(1 , indproto1) , Wproto_est(2 , indproto1) , 'mo')
h = voronoi(Wproto_est(1 , :) , Wproto_est(2 , :));
set(h , 'color' , 'y' , 'linewidth' , 2)
hold off
title('E_{H2MGLVQ}(t)' , 'fontsize' , 12)
colorbar
figure(2)
plot(E_H2MGLVQ);
title('E_{H2MGLVQ}(t)' , 'fontsize' , 12)
Example 2
---------
load ionosphere
Nproto_pclass = 4*ones(1 , length(unique(y)));
options.epsilonk = 0.05;
options.epsilonl = 0.01;
options.tmax = 1;
options.tmin = 0.1;
options.xi = 1;
options.nb_iterations = 5000;
options.shuffle = 1;
[d , N] = size(X);
ON = ones(1 , N);
mindata = min(X , [] , 2);
maxdata = max(X , [] , 2);
temp = maxdata - mindata;
temp(temp==0) = 1;
X = 2*(X - mindata(: , ON))./(temp(: , ON)) - 1;
n1 = round(0.7*N);
n2 = N - n1;
ind = randperm(length(y));
ind1 = ind(1:n1);
ind2 = ind(n1+1:N);
Xtrain = X(: , ind1);
ytrain = y(ind1);
Xtest = X(: , ind2);
ytest = y(ind2);
[Wproto , yproto] = ini_glvq(Xtrain , ytrain , Nproto_pclass);
[Wproto_est , yproto_est , E_H2MGLVQ] = h2m_glvq_model(Xtrain , ytrain , Wproto , yproto , options);
ytest_est = glvq_predict(Xtest , Wproto_est , yproto_est);
Perf = sum(ytest == ytest_est)/n2;
disp(Perf)
plot(E_H2MGLVQ);
title('E_{H2MGLVQ}(t)' , 'fontsize' , 12)
Example 3
---------
load artificial
Nproto_pclass = [15 , 12 , 3];
options.epsilonk = 0.05;
options.epsilonl = 0.01;
options.tmax = 1;
options.tmin = 0.1;
options.xi = 10;
options.nb_iterations = 3000;
options.shuffle = 1;
[d , N] = size(X);
ON = ones(1 , N);
n1 = round(0.7*N);
n2 = N - n1;
ind = randperm(length(y));
ind1 = ind(1:n1);
ind2 = ind(n1+1:N);
Xtrain = X(: , ind1);
ytrain = y(ind1);
Xtest = X(: , ind2);
ytest = y(ind2);
[Wproto , yproto] = ini_glvq(Xtrain , ytrain , Nproto_pclass);
[Wproto_est , yproto_est , E_H2MGLVQ] = h2m_glvq_model(Xtrain , ytrain , Wproto , yproto , options);
ytest_est = glvq_predict(Xtest , Wproto_est , yproto_est);
Perf = sum(ytest == ytest_est)/n2;
disp(Perf)
figure(1)
plot_label(X , y);
hold on
h = voronoi(Wproto_est(1 , :) , Wproto_est(2 , :));
hold off
figure(2)
plot(E_H2MGLVQ);
title('E_{H2MGLVQ}(t)' , 'fontsize' , 12)
Example 4
---------
close all
load artificial
%load iris
options.epsilonk = 0.05;
options.epsilonl = 0.01;
options.tmax = 1;
options.tmin = 0.1;
options.xi = 10;
options.nb_iterations = 5000;
options.shuffle = 1;
options.method = 6;
options.cv.K = 10;
options.holding.rho = 0.7;
options.holding.K = 20;
options.bootstraping.rho = 0.7;
options.bootstraping.K = 20;
%Nproto_pclass = 4*ones(1 , length(unique(y)));
Nproto_pclass = [15 , 12 , 3];
[Itrain , Itest] = sampling(X , y , options);
Perf = zeros(1 , size(Itrain , 1));
for i = 1 : size(Itrain , 1)
Xtrain = X(: , Itrain(i , :));
ytrain = y(Itrain(i , :));
Xtest = X(: , Itest(i , :));
ytest = y(Itest(i , :));
[Wproto , yproto] = ini_glvq(Xtrain , ytrain , Nproto_pclass);
[Wproto_est , yproto_est] = h2m_glvq_model(Xtrain , ytrain , Wproto , yproto , options);
ytest_est = glvq_predict(Xtest , Wproto_est , yproto_est, [ ] , options);
Perf(i) = sum(ytest == ytest_est)/length(ytest);
end
disp(mean(Perf))
Author : S閎astien PARIS : sebastien.paris@lsis.org
------- Date : 04/09/2006
Reference "A new Generalized LVQ Algorithm via Harmonic to Minimumm Distance Measure Transition", A.K. Qin, P.N. Suganthan and J.J. Liang,
--------- IEEE International Conference on System, Man and Cybernetics, 2004
*/
#include <time.h>
#include <math.h>
#include <mex.h>
#define mix(a , b , c) \
{ \
a -= b; a -= c; a ^= (c>>13); \
b -= c; b -= a; b ^= (a<<8); \
c -= a; c -= b; c ^= (b>>13); \
a -= b; a -= c; a ^= (c>>12); \
b -= c; b -= a; b ^= (a<<16); \
c -= a; c -= b; c ^= (b>>5); \
a -= b; a -= c; a ^= (c>>3); \
b -= c; b -= a; b ^= (a<<10); \
c -= a; c -= b; c ^= (b>>15); \
}
#define zigstep 128 // Number of Ziggurat'Steps
#define znew (z = 36969*(z&65535) + (z>>16) )
#define wnew (w = 18000*(w&65535) + (w>>16) )
#define MWC ((znew<<16) + wnew )
#define SHR3 ( jsr ^= (jsr<<17), jsr ^= (jsr>>13), jsr ^= (jsr<<5) )
#define CONG (jcong = 69069*jcong + 1234567)
#define KISS ((MWC^CONG) + SHR3)
#ifdef ranKISS
#define randint KISS
#define rand() (randint*2.328306e-10)
#endif
#ifdef ranSHR3
#define randint SHR3
#define rand() (0.5 + (signed)randint*2.328306e-10)
#endif
typedef unsigned long UL;
static UL jsrseed = 31340134 , jsr;
#ifdef ranKISS
static UL z=362436069, w=521288629, jcong=380116160;
#endif
static UL jz , iz , kn[zigstep];
static long hz;
static double wn[zigstep] , fn[zigstep];
typedef struct OPTIONS
{
double epsilonk;
double epsilonl;
double tmax;
double tmin;
double xi;
int nb_iterations;
int shuffle;
} OPTIONS;
/* Function prototypes */
void randini(void);
void randnini(void);
double nfix(void);
double randn(void);
void qs( double * , int , int );
void qsindex( double * , int * , int , int );
void h2m_glvq_model(double * , double * , OPTIONS , int , int , int , int ,
double * , double * , double *,
double * , double * , int * , double * , double * ,
int * );
/*-------------------------------------------------------------------------------------------------------------- */
void mexFunction( int nlhs, mxArray *plhs[] , int nrhs, const mxArray *prhs[] )
{
double *Xtrain , *ytrain , *Wproto , *yproto;
OPTIONS options = {0.05 , 0.01 , 1 , 0.6 , 1 , 1000 , 1};
double *Wproto_est , *yproto_est , *E_H2MGLVQ;
int d , Ntrain , Nproto , m = 0;
int i , j , l , co , Nprotom , ld , id , indice;
double currentlabel , ind , temp;
double *tmp , *ytrainsorted , *labels , *mtemp , *stdtemp , *temp_train , *dist , *powdist;
int *Nk , *index_train;
mxArray *mxtemp;
/*Initialize Random generator */
randini();
randnini();
/* Input 1 Xtrain */
Xtrain = mxGetPr(prhs[0]);
if( mxGetNumberOfDimensions(prhs[0]) !=2 )
{
mexErrMsgTxt("Xtrain must be (d x Ntrain)");
}
d = mxGetM(prhs[0]);
Ntrain = mxGetN(prhs[0]);
Nproto = floor(sqrt(Ntrain));
/* Input 2 ytrain */
ytrain = mxGetPr(prhs[1]);
if(mxGetNumberOfDimensions(prhs[1]) !=2 || mxGetN(prhs[1]) != Ntrain)
{
mexErrMsgTxt("ytrain must be (1 x Ntrain)");
}
/* Determine unique Labels */
ytrainsorted = mxMalloc(Ntrain*sizeof(double));
for ( i = 0 ; i < Ntrain; i++ )
{
ytrainsorted[i] = ytrain[i];
}
qs( ytrainsorted , 0 , Ntrain - 1 );
labels = mxMalloc(sizeof(double));
labels[m] = ytrainsorted[0];
currentlabel = labels[0];
for (i = 0 ; i < Ntrain ; i++)
{
if (currentlabel != ytrainsorted[i])
{
labels = (double *)mxRealloc(labels , (m+2)*sizeof(double));
labels[++m] = ytrainsorted[i];
currentlabel = ytrainsorted[i];
}
}
m++;
/* Input 5 option */
if ( (nrhs > 4) && !mxIsEmpty(prhs[4]) )
{
mxtemp = mxGetField(prhs[4] , 0, "epsilonk");
if(mxtemp != NULL)
{
tmp = mxGetPr(mxtemp);
options.epsilonk = tmp[0];
}
mxtemp = mxGetField(prhs[4] , 0, "epsilonl");
if(mxtemp != NULL)
{
tmp = mxGetPr(mxtemp);
options.epsilonl = tmp[0];
}
mxtemp = mxGetField(prhs[4] , 0, "tmax");
if(mxtemp != NULL)
{
tmp = mxGetPr(mxtemp);
options.tmax = tmp[0];
}
mxtemp = mxGetField(prhs[4] , 0, "tmin");
if(mxtemp != NULL)
{
tmp = mxGetPr(mxtemp);
options.tmin = tmp[0];
}
mxtemp = mxGetField(prhs[4] , 0, "xi");
if(mxtemp != NULL)
{
tmp = mxGetPr(mxtemp);
options.xi = tmp[0];
}
mxtemp = mxGetField(prhs[4] , 0, "nb_iterations");
if(mxtemp != NULL)
{
tmp = mxGetPr(mxtemp);
options.nb_iterations = (int) tmp[0];
}
mxtemp = mxGetField(prhs[4] , 0, "shuffle");
if(mxtemp != NULL)
{
tmp = mxGetPr(mxtemp);
options.shuffle = (int) tmp[0];
}
}
/* Input 4 yproto */
if ((nrhs < 4) || mxIsEmpty(prhs[3]) )
{
plhs[1] = mxCreateDoubleMatrix(1 , Nproto, mxREAL);
yproto_est = mxGetPr(plhs[1]);
co = 0;
Nprotom = ceil((double)Nproto/(double)m);
for(i = 0 ; i < m-1 ; i++)
{
ind = labels[i];
for(j = 0 ; j < Nprotom ; j++)
{
yproto_est[co] = labels[i];
co++;
}
}
ind = labels[m-1];
for(j = (m-1)*Nprotom ; j < Nproto ; j++)
{
yproto_est[co] = ind;
co++;
}
}
else
{
yproto = mxGetPr(prhs[3]);
if(mxGetNumberOfDimensions(prhs[3]) !=2)
{
mexErrMsgTxt("yproto must be (1 x Nproto)");
}
Nproto = mxGetN(prhs[3]);
plhs[1] = mxCreateDoubleMatrix(1 , Nproto, mxREAL);
yproto_est = mxGetPr(plhs[1]);
for( i = 0 ; i < Nproto ; i++)
{
yproto_est[i] = yproto[i];
}
}
/* Input 3 Wproto */
Nk = mxMalloc(m*sizeof(int));
if ((nrhs < 3) || mxIsEmpty(prhs[2]) )
{
mtemp = mxMalloc(d*m*sizeof(double));
stdtemp = mxMalloc(d*m*sizeof(double));
for(i = 0 ; i < d*m ; i++)
{
mtemp[i] = 0.0;
stdtemp[i] = 0.0;
}
for (l = 0 ; l < m ; l++)
{
ind = labels[l];
ld = l*d;
for (i = 0 ; i < Ntrain ; i++)
{
if(ytrain[i] == ind)
{
id = i*d;
for(j = 0 ; j < d ; j++)
{
mtemp[j + ld] += Xtrain[j + id];
}
Nk[l]++;
}
}
}
for (l = 0 ; l < m ; l++)
{
ld = l*d;
temp = 1.0/Nk[l];
for(j = 0 ; j < d ; j++)
{
mtemp[j + ld] *=temp;
}
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -