📄 srng_model.c
字号:
/*
Supervised Relevance Natural Gaz classification algorithm.
Usage
------
[Wproto_est , yproto_est , E_SRNG] = srng_model(Xtrain , ytrain , [Wproto] , [yproto] , [lambda] , [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)
lambda Weigths factor (d x 1). Default lambda = (1/d)*ones(d , 1);
options Options'structure
.epsilonk Update weight's for class label between prototype and train data (default = 0.01).
.epsilonl Update weight's for class label between prototype and train data (default = 0.01).
.epsilonlambda lambda rate (default 10e-6)
sigmastart Start sigma size neigbhourg (default 2)
.sigmaend Final sigma size
.sigmastretch Sigma step stretch;
.threshold Threshold to speed up computation (default 10e-10);
.xi Constant in the sigmoid function (default = 1)
.nb_iterations Number of iterations (default = 1000)
.metric_method Euclidean 2 or Euclidean4l1 (default = 1)
.shuffle Randomly permute order of train data between each iteration if shuffle=1 (default = 1)
.updatelambda Update Lambda if = 1 (default = 1)
Outputs
-------
Wproto_est Final prototypes weigths (d x Nproto)
yproto_est Final prototypes labels (1 x Nproto)
lambda_est Final Weigths factor estimated (d x 1)
E_SRNG Energy criteria versur iteration (1 x options.nb_iterations)
To compile
----------
mex -DranSHR3 -output srng_model.dll srng_model.c
mex -DranSHR3 -f mexopts_intel10amd.bat -output srng_model.dll srng_model.c
Example 1
---------
d = 2;
Ntrain = 100;
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.005;
options.epsilonl = 0.001;
options.epsilonlambda = 10e-8;
options.sigmastart = 2;
options.sigmaend = 10e-4;
options.sigmastretch = 10e-3;
options.threshold = 10e-10;
options.xi = 0.1;
options.nb_iterations = 3000;
options.metric_method = 1;
options.shuffle = 1;
options.updatelambda = 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(:)'];
Nproto_pclass = 4*ones(1 , length(unique(ytrain)));
[Wproto , yproto , lambda] = ini_proto(Xtrain , ytrain , Nproto_pclass);
[Wproto_est , yproto_est , lambda_est, E_SRNG] = srng_model(Xtrain , ytrain , Wproto , yproto , lambda , options);
ytest_est = NN_predict(Xtest , Wproto_est , yproto_est,lambda_est,options);
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_{SRNG}(t)' , 'fontsize' , 12)
colorbar
figure(2)
plot(E_SRNG);
title('E_{SRNG}(t)' , 'fontsize' , 12)
figure(3)
stem(lambda_est);
title('\lambda' , 'fontsize' , 12)
Example 2
---------
close all
load ionosphere
options.epsilonk = 0.005;
options.epsilonl = 0.001;
options.epsilonlambda = 10e-7;
options.sigmastart = 2;
options.sigmaend = 10e-5;
options.sigmastretch = 10e-4;
options.threshold = 10e-10;
options.xi = 10;
options.nb_iterations = 10000;
options.metric_method = 1;
options.shuffle = 1;
options.updatelambda = 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);
Nproto_pclass = 4*ones(1 , length(unique(y)));
[Wproto , yproto , lambda] = ini_proto(Xtrain , ytrain , Nproto_pclass);
[Wproto_est , yproto_est , lambda_est, E_SRNG] = srng_model(Xtrain , ytrain , Wproto , yproto , lambda, options);
ytest_est = NN_predict(Xtest , Wproto_est , yproto_est,lambda_est,options);
Perf = sum(ytest == ytest_est)/n2;
disp(Perf)
plot(E_SRNG);
title('E_{SRNG}(t)' , 'fontsize' , 12)
figure(2)
stem(lambda_est);
title('\lambda' , 'fontsize' , 12)
Example 3
---------
clear,close all hidden
load artificial
Nproto_pclass = [15 , 12 , 3];
options.epsilonk = 0.005;
options.epsilonl = 0.001;
options.epsilonlambda = 10e-8;
options.sigmastart = 4;
options.sigmaend = 10e-5;
options.sigmastretch = 10e-4;
options.threshold = 10e-10;
options.xi = 1;
options.nb_iterations = 25000;
options.metric_method = 1;
options.shuffle = 1;
options.updatelambda = 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, lambda] = ini_proto(Xtrain , ytrain , Nproto_pclass);
[Wproto_est , yproto_est , lambda_est, E_SRNG] = srng_model(Xtrain , ytrain , Wproto , yproto , lambda, options);
ytest_est = NN_predict(Xtest , Wproto_est , yproto_est,lambda_est,options);
Perf = sum(ytest == ytest_est)/n2;
disp(Perf)
figure(1)
plot_label(X , y);
hold on
h = voronoi(Wproto_est(1 , :) , Wproto_est(2 , :));
figure(2)
plot(E_SRNG);
title('E_{SRNG}(t)' , 'fontsize' , 12)
figure(3)
stem(lambda_est);
Example 4
---------
load glass
options.epsilonk = 0.005;
options.epsilonl = 0.001;
options.epsilonlambda = 10e-8;
options.sigmastart = 2;
options.sigmaend = 10e-5;
options.sigmastretch = 10e-4;
options.threshold = 10e-10;
options.xi = 10;
options.nb_iterations = 10000;
options.metric_method = 1;
options.shuffle = 1;
options.updatelambda = 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);
Nproto_pclass = 4*ones(1 , length(unique(y)));
[Wproto , yproto , lambda] = ini_proto(Xtrain , ytrain , Nproto_pclass);
[Wproto_est , yproto_est , lambda_est, E_SRNG] = srng_model(Xtrain , ytrain , Wproto , yproto , lambda , options);
ytest_est = NN_predict(Xtest , Wproto_est , yproto_est,lambda_est,options);
Perf = sum(ytest == ytest_est)/n2;
disp(Perf)
plot(E_SRNG);
title('E_{SRNG}(t)' , 'fontsize' , 12)
figure(3)
stem(lambda_est);
Author : S閎astien PARIS : sebastien.paris@lsis.org
------- Date : 04/09/2006
Reference "Margin based Active Learning for LVQ Networks", F.-M. Schleif and B. Hammer and Th. Villmann, ESANN'2006 proceedings
*/
#include <time.h>
#include <math.h>
#include <mex.h>
#define MAX(A,B) (((A) > (B)) ? (A) : (B) )
#define MIN(A,B) (((A) < (B)) ? (A) : (B) )
#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 iz , kn[zigstep];
static long hz;
static double wn[zigstep] , fn[zigstep];
typedef struct OPTIONS
{
double epsilonk;
double epsilonl;
double epsilonlambda;
double sigmastart;
double sigmaend;
double sigmastretch;
double threshold;
double xi;
int nb_iterations;
int metric_method;
int shuffle;
int updatelambda;
} 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 srng_model(double * , double * , OPTIONS , int , int , int , int ,
double * , double * , double *, double *, double * ,
double * , double * , int * , double * , int * , int *, double* , double* , int);
/*-------------------------------------------------------------------------------------------------------------- */
void mexFunction( int nlhs, mxArray *plhs[] , int nrhs, const mxArray *prhs[] )
{
double *Xtrain , *ytrain , *Wproto , *yproto;
OPTIONS options = {0.005 , 0.001 , 10e-6 , 2 , 10e-5 , 10e-4 , 10e-10 , 1 , 1000 , 1 , 1 , 1};
double *Wproto_est , *yproto_est , *E_SRNG , *lambda_est;
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 , *lambda;
double *hproto , *sproto , *distj;
int *Np ;
int *Nk , *index_train , *rank_distj , *index_distj;
int Nproto_max=0;
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 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;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -