?? clda.cpp
字號:
#include "stdlib.h"
#include <iostream>
#include <string.h>
#include <fstream>
#include <math.h>
using namespace std;
#include "global.h"
CLda::CLda()
{
this->numF = -1;
}
CLda::~CLda()
{
}
bool CLda::extract(char* trainFile, char* newTrainFile, char* testFile, char* newTestFile)
{
bool ret = true;
int i, j, k, index, cnum;
if (this->numF <= 0)
this->numF = 1;
ret = this->sdata.readFile(trainFile);//讀入訓練數據
if (!ret)
return false;//如果數據格式不正確，退出程序。
int numClass = this->sdata.numClass;
int numFeature = this->sdata.numFeature;
int numSample = this->sdata.numSample;
CSampleData testData;
ret = testData.readFile(testFile);
if (!ret)
return false;//如果數據格式不正確，退出程序。
DOUBLE temp1, temp2;
if (this->numF > this->sdata.numFeature)
this->numF = this->sdata.numFeature;
int newNumF = this->numF;
DOUBLE *sumX = new DOUBLE[numClass*numFeature];
DOUBLE *sumXX = new DOUBLE[numClass*numFeature*numFeature];
int* classSample = new int[numClass];
int total;
for (i=0; i<numClass; i++)
{
classSample[i]=0;
}
total = numClass * numFeature;
for (i=0; i<total; i++)
{
sumX[i]=0;
}
total = numClass * numFeature * numFeature;
for (i=0; i<total; i++)
{
sumXX[i]=0;
}
for (i=0; i<numSample; i++)
{
cnum = this->sdata.ydata[i];
classSample[cnum] ++;
for (j=0; j<numFeature; j++)
{
sumX[cnum*numFeature+j] += this->sdata.xdata[i*numFeature+j];
for (k=j; k<numFeature; k++)
{
sumXX[(cnum*numFeature+j)*numFeature+k] +=this->sdata.xdata[i*numFeature+j] * this->sdata.xdata[i*numFeature+k];
}
}
}
for (i=0; i<numClass; i++)
{
for (j=0; j<numFeature; j++)
{
sumX[i*numFeature+j] = sumX[i*numFeature+j] / ((DOUBLE) classSample[i]);
}
}
for (i=0; i<numClass; i++)
{
for (j=0; j<numFeature; j++)
{
for (k=j; k<numFeature; k++)
{
index = (i*numFeature+j)*numFeature+k;
sumXX[index] = sumXX[index]/((DOUBLE)classSample[i])-(sumX[i*numFeature+j]*sumX[i*numFeature+k]);
}
}
}
DOUBLE *U0 = new DOUBLE[numFeature];
DOUBLE *SW = new DOUBLE[numFeature*numFeature];
DOUBLE *SB = new DOUBLE[numFeature*numFeature];
for (j=0; j<numFeature; j++)
{
U0[j] = 0;
for (k=0; k<numClass; k++)
{
U0[j] += sumX[k*numFeature+j]*((DOUBLE)classSample[k]);
}
U0[j] = U0[j]/((DOUBLE)numSample);
}
for (i=0; i<numFeature; i++)
{
for (j=i; j<numFeature; j++)
{
index = i*numFeature+j;
SW[index] = 0;
for (k=0; k<numClass; k++)
{
SW[index] += sumXX[((k*numFeature)+i)*numFeature+j]*classSample[k];
}
SW[index] = SW[index]/((DOUBLE)numSample);
}
}
for (i=0; i<numFeature; i++)
{
for (j=i; j<numFeature; j++)
{
index = i*numFeature+j;
SB[index] = 0;
for (k=0; k<numClass; k++)
{
SB[index] += (sumX[k*numFeature+i]-U0[i])*(sumX[k*numFeature+j]-U0[j])*((DOUBLE)classSample[k]);
}
SB[index] = SB[index]/((DOUBLE)numSample);
}
}
for (j=1; j<numFeature; j++)
{
for (k=0; k<j; k++)
{
SW[j*numFeature+k] = SW[k*numFeature+j];
SB[j*numFeature+k] = SB[k*numFeature+j];
}
}
/*for (j=0; j<numFeature; j++)//限制協方差矩陣的對角元素的下限
{
if (SW[j*numFeature+j]<0.001)
{
SW[j*numFeature+j] = 0.001;
}
}*/
//計算特征值
CMatrix theMat;
float** Cova;
// Allocate space for the covariance matrix
Cova = theMat.allocMat(numFeature); // n: dimensionality of vector
for (i=0; i<numFeature; i++)
{
for (j=0; j<numFeature; j++)
{
Cova[i][j] = (float)SW[i*numFeature+j];
}
}
// Compute the covariance matrix......
// Diagonalization
// P: matrix with columns as the eigenvectors
float** P;
P = theMat.allocMat(numFeature);
theMat.diagonalize( Cova, numFeature, P );
//顯示計算的特征值
for (i=0; i<numFeature; i++)
{
cout<<endl;
for (j=0; j<numFeature; j++)
{
temp2 = 0;
for (k=0; k<numFeature; k++)
{
temp2 += SW[i*numFeature+k]*((DOUBLE)P[k][j]);
}
cout<<temp2/P[i][j]<<" ";
}
}
cout<<endl;
DOUBLE* featureValue = new DOUBLE[numFeature];
int* featureIndex = new int[numFeature];
for (i=0; i<numFeature; i++)
{
featureValue[i] = 0;
featureIndex[i] = i;
}
for (j=0; j<numFeature; j++)//計算特征值
{
temp1 = 0;
for (i=0; i<numFeature; i++)
{
temp2 = 0;
for (k=0; k<numFeature; k++)
{
temp2 += SW[i*numFeature+k]*((DOUBLE)P[k][j]);
}
if (P[i][j]>0)
{
temp1 += (DOUBLE)P[i][j];
featureValue[j] += temp2;
}
else
{
temp1 -= (DOUBLE)P[i][j];
featureValue[j] -= temp2;
}
}
featureValue[j] = featureValue[j]/temp1;
cout<<featureValue[j]<<endl;
}
DOUBLE* W1 = new DOUBLE[numFeature*numFeature];//計算W1
for (i=0; i<numFeature; i++)
{
if (featureValue[i]<0.001)//限制特征值最小值為0.001
{
temp1 = sqrt(0.001);
}
else
{
temp1 = sqrt(featureValue[i]);
}
for (j=0; j<numFeature; j++)
{
W1[j*numFeature+i] = ((DOUBLE)P[j][i])/temp1;
}
}
for (i=0; i<numFeature; i++)
{
for (j=0; j<numFeature; j++)
{
sumXX[i*numFeature+j] = 0;
for (k=0; k<numFeature; k++)
{
sumXX[i*numFeature+j] += W1[k*numFeature+i]*SB[k*numFeature+j];
}
}
}
for (i=0; i<numFeature; i++)
{
for (j=i; j<numFeature; j++)
{
SB[i*numFeature+j] = 0;
for (k=0; k<numFeature; k++)
{
SB[i*numFeature+j] += sumXX[i*numFeature+k]*W1[k*numFeature+j];
}
}
}
for (i=1; i<numFeature; i++)
{
for (j=0; j<i; j++)
{
SB[i*numFeature+j] = SB[j*numFeature+i];
}
}
cout<<endl;
for (i=0; i<numFeature; i++)
{
cout<<endl;
for (j=0; j<numFeature; j++)
{
Cova[i][j] = (float)SB[i*numFeature+j];
cout<<Cova[i][j]<<" ";
}
}
cout<<endl;
theMat.diagonalize( Cova, numFeature, P );//
//顯示計算的特征值
for (i=0; i<numFeature; i++)
{
cout<<endl;
for (j=0; j<numFeature; j++)
{
temp2 = 0;
for (k=0; k<numFeature; k++)
{
temp2 += SB[i*numFeature+k]*((DOUBLE)P[k][j]);
}
cout<<temp2/P[i][j]<<" ";
}
}
cout<<endl;
for (i=0; i<numFeature; i++)
{
featureValue[i] = 0;
featureIndex[i] = i;
}
for (j=0; j<numFeature; j++)//計算特征值
{
temp1 = 0;
for (i=0; i<numFeature; i++)
{
temp2 = 0;
for (k=0; k<numFeature; k++)
{
temp2 += SB[i*numFeature+k]*((DOUBLE)P[k][j]);
}
if (P[i][j]>0)
{
temp1 += (DOUBLE)P[i][j];
featureValue[j] += temp2;
}
else
{
temp1 -= (DOUBLE)P[i][j];
featureValue[j] -= temp2;
}
}
featureValue[j] = featureValue[j]/temp1;
cout<<featureValue[j]<<endl;
}
for (i=0; i<numFeature-1; i++)//特征值排序
{
for (j=numFeature-1; j>i; j--)
{
if (featureValue[j-1]<featureValue[j])
{
temp1 = featureValue[j-1];
featureValue[j-1] = featureValue[j];
featureValue[j] = temp1;
k = featureIndex[j-1];
featureIndex[j-1] = featureIndex[j];
featureIndex[j] = k;
}
}
}
for (i=0; i<numFeature; i++)
{
for (j=0; j<numFeature; j++)
{
sumXX[i*numFeature+j] = (DOUBLE)P[i][featureIndex[j]];
}
}
DOUBLE* W = new DOUBLE[numFeature*numFeature];//計算W
for (i=0; i<numFeature; i++)
{
for (j=0; j<numFeature; j++)
{
W[i*numFeature+j] = 0;
for (k=0; k<numFeature; k++)
{
W[i*numFeature+j] += W1[i*numFeature+k]*sumXX[k*numFeature+j];
}
}
}
//提取特征
DOUBLE *temp = new DOUBLE[newNumF];
for (i=0; i<numSample; i++)//中心化
{
for (j=0; j<numFeature; j++)
{
this->sdata.xdata[i*numFeature+j] -= U0[j];
}
}
for (i=0; i<numSample; i++)
{
for (j=0; j<newNumF; j++)
{
temp[j] = 0;
}
for (j=0; j<newNumF; j++)
{
for (k=0; k<numFeature; k++)
{
temp[j] += (W[k*numFeature+j])*this->sdata.xdata[i*numFeature+k];
}
}
for (j=0; j<newNumF; j++)
{
this->sdata.xdata[i*newNumF+j] = temp[j];
}
}
this->sdata.numFeature = newNumF;
this->sdata.writeFile(newTrainFile);//寫新的訓練數據文件
for (i=0; i<testData.numSample; i++)//中心化測試數據
{
for (j=0; j<numFeature; j++)
{
testData.xdata[i*numFeature+j] -= U0[j];
}
}
for (i=0; i<testData.numSample; i++)
{
for (j=0; j<newNumF; j++)
{
temp[j] = 0;
}
for (j=0; j<newNumF; j++)
{
for (k=0; k<numFeature; k++)
{
temp[j] += (W[k*numFeature+j])*testData.xdata[i*numFeature+k];
}
}
for (j=0; j<newNumF; j++)
{
testData.xdata[i*newNumF+j] = temp[j];
}
}
testData.numFeature = newNumF;
testData.writeFile(newTestFile);
delete[] Cova;
delete[] P;
delete[] temp;
delete[] featureValue;
delete[] featureIndex;
delete[] W1;
delete[] W;
delete[] sumX;
delete[] sumXX;
delete[] U0;
delete[] SW;
delete[] SB;
return ret;
}
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -