lda3.c

LDA程式設計

	for( i=0;i<lda_sample;i++ )
	{
		gsl_matrix_add( v[0].m, class1[i].m );
		gsl_matrix_add( v[1].m, class2[i].m );
		gsl_matrix_add( v[2].m, class3[i].m );
		gsl_matrix_add( v[3].m, class4[i].m );
		gsl_matrix_add( v[4].m, class5[i].m );
		gsl_matrix_add( v[5].m, class6[i].m );
		gsl_matrix_add( v[6].m, class7[i].m );
		gsl_matrix_add( v[7].m, class8[i].m );
	}

	for( i=0;i<class_num;i++ )
	{
		gsl_matrix_memcpy( mean[i].m, v[i].m );
		gsl_matrix_scale( mean[i].m, 0.04 );
		gsl_matrix_memcpy( mean_temp[i].m, mean[i].m );
		gsl_matrix_set_zero( v[i].m );
	}

	for( j=0;j<lda_sample;j++ )
	{
		gsl_matrix_add( v[0].m, class1[j].m );
		gsl_matrix_add( v[0].m, class2[j].m );
		gsl_matrix_add( v[0].m, class3[j].m );
		gsl_matrix_add( v[0].m, class4[j].m );
		gsl_matrix_add( v[0].m, class5[j].m );
		gsl_matrix_add( v[0].m, class6[j].m );
		gsl_matrix_add( v[0].m, class7[j].m );
		gsl_matrix_add( v[0].m, class8[j].m );
	}
	gsl_matrix_memcpy( all_mean->m, v[0].m );
	gsl_matrix_scale( all_mean->m, 0.005 );
	gsl_matrix_set_zero( v[0].m );
/***********************************************************************/

/*************************向量減去平均值******************************/

	for( i=0;i<lda_sample;i++ )
	{
		gsl_matrix_sub( class1[i].m, mean[0].m );
		gsl_matrix_sub( class2[i].m, mean[1].m );
		gsl_matrix_sub( class3[i].m, mean[2].m );
		gsl_matrix_sub( class4[i].m, mean[3].m );
		gsl_matrix_sub( class5[i].m, mean[4].m );
		gsl_matrix_sub( class6[i].m, mean[5].m );
		gsl_matrix_sub( class7[i].m, mean[6].m );
		gsl_matrix_sub( class8[i].m, mean[7].m );
	}

/***********************************************************************/



/************************** 求E[(x-E[x])(x-E[x])T] *********************/

	for( i=0;i<lda_sample;i++ )
	{
		class1_T[i].m = gsl_matrix_alloc(1, sample_dim);
		gsl_matrix_transpose_memcpy ( class1_T[i].m, class1[i].m);
		class2_T[i].m = gsl_matrix_alloc(1, sample_dim);
		gsl_matrix_transpose_memcpy ( class2_T[i].m, class2[i].m);
		class3_T[i].m = gsl_matrix_alloc(1, sample_dim);
		gsl_matrix_transpose_memcpy ( class3_T[i].m, class3[i].m);
		class4_T[i].m = gsl_matrix_alloc(1, sample_dim);
		gsl_matrix_transpose_memcpy ( class4_T[i].m, class4[i].m);
		class5_T[i].m = gsl_matrix_alloc(1, sample_dim);
		gsl_matrix_transpose_memcpy ( class5_T[i].m, class5[i].m);
		class6_T[i].m = gsl_matrix_alloc(1, sample_dim);
		gsl_matrix_transpose_memcpy ( class6_T[i].m, class6[i].m);
		class7_T[i].m = gsl_matrix_alloc(1, sample_dim);
		gsl_matrix_transpose_memcpy ( class7_T[i].m, class7[i].m);
		class8_T[i].m = gsl_matrix_alloc(1, sample_dim);
		gsl_matrix_transpose_memcpy ( class8_T[i].m, class8[i].m);
	}

	for( i=0;i<lda_sample;i++ )
	{
		class1_mul[i].m = gsl_matrix_alloc(sample_dim, sample_dim);
		class2_mul[i].m = gsl_matrix_alloc(sample_dim, sample_dim);
		class3_mul[i].m = gsl_matrix_alloc(sample_dim, sample_dim);
		class4_mul[i].m = gsl_matrix_alloc(sample_dim, sample_dim);
		class5_mul[i].m = gsl_matrix_alloc(sample_dim, sample_dim);
		class6_mul[i].m = gsl_matrix_alloc(sample_dim, sample_dim);
		class7_mul[i].m = gsl_matrix_alloc(sample_dim, sample_dim);
		class8_mul[i].m = gsl_matrix_alloc(sample_dim, sample_dim);
	}

	for( i=0;i<lda_sample;i++ )
	{
		gsl_matrix_mul_matrix( class1[i].m, class1_T[i].m, class1_mul[i].m );
		gsl_matrix_mul_matrix( class2[i].m, class2_T[i].m, class2_mul[i].m );
		gsl_matrix_mul_matrix( class3[i].m, class3_T[i].m, class3_mul[i].m );
		gsl_matrix_mul_matrix( class4[i].m, class4_T[i].m, class4_mul[i].m );
		gsl_matrix_mul_matrix( class5[i].m, class5_T[i].m, class5_mul[i].m );
		gsl_matrix_mul_matrix( class6[i].m, class6_T[i].m, class6_mul[i].m );
		gsl_matrix_mul_matrix( class7[i].m, class7_T[i].m, class7_mul[i].m );
		gsl_matrix_mul_matrix( class8[i].m, class8_T[i].m, class8_mul[i].m );
	}

	class1_sw->m = gsl_matrix_calloc(sample_dim, sample_dim);
	class2_sw->m = gsl_matrix_calloc(sample_dim, sample_dim);
	class3_sw->m = gsl_matrix_calloc(sample_dim, sample_dim);
	class4_sw->m = gsl_matrix_calloc(sample_dim, sample_dim);
	class5_sw->m = gsl_matrix_calloc(sample_dim, sample_dim);
	class6_sw->m = gsl_matrix_calloc(sample_dim, sample_dim);
	class7_sw->m = gsl_matrix_calloc(sample_dim, sample_dim);
	class8_sw->m = gsl_matrix_calloc(sample_dim, sample_dim);
	m_sw->m = gsl_matrix_calloc(sample_dim, sample_dim);

	for( i=0;i<lda_sample;i++ )
	{
		gsl_matrix_add( class1_sw->m, class1_mul[i].m );
		gsl_matrix_add( class2_sw->m, class2_mul[i].m );
		gsl_matrix_add( class3_sw->m, class3_mul[i].m );
		gsl_matrix_add( class4_sw->m, class4_mul[i].m );
		gsl_matrix_add( class5_sw->m, class5_mul[i].m );
		gsl_matrix_add( class6_sw->m, class6_mul[i].m );
		gsl_matrix_add( class7_sw->m, class7_mul[i].m );
		gsl_matrix_add( class8_sw->m, class8_mul[i].m );
	}

	gsl_matrix_add( m_sw->m, class1_sw->m );
	gsl_matrix_add( m_sw->m, class2_sw->m );
	gsl_matrix_add( m_sw->m, class3_sw->m );
	gsl_matrix_add( m_sw->m, class4_sw->m );
	gsl_matrix_add( m_sw->m, class5_sw->m );
	gsl_matrix_add( m_sw->m, class6_sw->m );
	gsl_matrix_add( m_sw->m, class7_sw->m );
	gsl_matrix_add( m_sw->m, class8_sw->m );

	

/***********************************************************************/

/***********************    SB    ******************************/
	matrixPtr mean_T;
	mean_T = malloc( class_num*sizeof( matrix ) );

	class1_sb->m = gsl_matrix_calloc(sample_dim, sample_dim);
	class2_sb->m = gsl_matrix_calloc(sample_dim, sample_dim);
	class3_sb->m = gsl_matrix_calloc(sample_dim, sample_dim);
	class4_sb->m = gsl_matrix_calloc(sample_dim, sample_dim);
	class5_sb->m = gsl_matrix_calloc(sample_dim, sample_dim);
	class6_sb->m = gsl_matrix_calloc(sample_dim, sample_dim);
	class7_sb->m = gsl_matrix_calloc(sample_dim, sample_dim);
	class8_sb->m = gsl_matrix_calloc(sample_dim, sample_dim);
	m_sb->m = gsl_matrix_calloc(sample_dim, sample_dim);

	for( i=0;i<class_num;i++ )
	{
		mean_T[i].m = gsl_matrix_alloc(1, sample_dim);
		gsl_matrix_sub( mean[i].m, all_mean->m );
		gsl_matrix_transpose_memcpy ( mean_T[i].m, mean[i].m);
	}

	gsl_matrix_mul_matrix( mean[0].m, mean_T[0].m, class1_sb->m );
	gsl_matrix_mul_matrix( mean[1].m, mean_T[1].m, class2_sb->m );
	gsl_matrix_mul_matrix( mean[2].m, mean_T[2].m, class3_sb->m );
	gsl_matrix_mul_matrix( mean[3].m, mean_T[3].m, class4_sb->m );
	gsl_matrix_mul_matrix( mean[4].m, mean_T[4].m, class5_sb->m );
	gsl_matrix_mul_matrix( mean[5].m, mean_T[5].m, class6_sb->m );
	gsl_matrix_mul_matrix( mean[6].m, mean_T[6].m, class7_sb->m );
	gsl_matrix_mul_matrix( mean[7].m, mean_T[7].m, class8_sb->m );

	gsl_matrix_add( m_sb->m, class1_sb->m );
	gsl_matrix_add( m_sb->m, class2_sb->m );
	gsl_matrix_add( m_sb->m, class3_sb->m );
	gsl_matrix_add( m_sb->m, class4_sb->m );
	gsl_matrix_add( m_sb->m, class5_sb->m );
	gsl_matrix_add( m_sb->m, class6_sb->m );
	gsl_matrix_add( m_sb->m, class7_sb->m );
	gsl_matrix_add( m_sb->m, class8_sb->m );

/*************************************************************/

/***********************    Compute SW Inverse    ******************************/

	gsl_matrix * m_sw_c = gsl_matrix_alloc(18, 18);
	gsl_matrix_memcpy( m_sw_c, m_sw->m );

	gsl_permutation * p = gsl_permutation_alloc( 18 );
	gsl_linalg_LU_decomp( m_sw_c, p, &sign );
	
	gsl_matrix * m_sw_inverse = gsl_matrix_alloc( 18, 18 );
	gsl_linalg_LU_invert(m_sw_c, p, m_sw_inverse);
/*******************************************************************************/

/***********************    Compute SW Inverse * SB    ******************************/
	gsl_matrix * m_sw_sb = gsl_matrix_alloc(18, 18);
	gsl_matrix_mul_matrix( m_sw_inverse, m_sb->m, m_sw_sb );
/************************************************************************************/

/**********************  Compute Eigen Problem  **********************************/
	//gsl_eigen_nonsymmv_workspace * w = gsl_eigen_nonsymmv_alloc ( 18 );
	gsl_vector_complex *eval = gsl_vector_complex_alloc (18);
	gsl_matrix_complex *evec = gsl_matrix_complex_alloc (18, 18);
     
	gsl_eigen_nonsymmv_workspace * w = gsl_eigen_nonsymmv_alloc (18);
       
	gsl_eigen_nonsymmv (m_sw_sb, eval, evec, w);
     
	gsl_eigen_nonsymmv_free (w);
     
	gsl_eigen_nonsymmv_sort (eval, evec, GSL_EIGEN_SORT_ABS_DESC);

	gsl_matrix *W = gsl_matrix_alloc (7, 18);
	for (i = 0; i < 18; i++)
	{
		gsl_complex eval_i = gsl_vector_complex_get (eval, i);
		gsl_vector_complex_view evec_i = gsl_matrix_complex_column (evec, i);
     
		//printf ("eigenvalue = %g + %gi\n", GSL_REAL(eval_i), GSL_IMAG(eval_i));
		printf ("eigenvector %d = \n", i);
		for (j = 0; j < 18; ++j)
		{
			gsl_complex z = gsl_vector_complex_get(&evec_i.vector, j);
			printf("%g + %gi\n", GSL_REAL(z), GSL_IMAG(z));

			if(i < 7)
			{
				gsl_matrix_set (W, i, j, GSL_REAL(z)/*gsl_complex_abs (z)*/);
			}
		}
	}
/**********************   求m prime  ***********************************************/
	
	gsl_matrix_mul_matrix(W, mean_temp[0].m, class1_mean_prime);
	gsl_matrix_mul_matrix(W, mean_temp[1].m, class2_mean_prime);
	gsl_matrix_mul_matrix(W, mean_temp[2].m, class3_mean_prime);
	gsl_matrix_mul_matrix(W, mean_temp[3].m, class4_mean_prime);
	gsl_matrix_mul_matrix(W, mean_temp[4].m, class5_mean_prime);
	gsl_matrix_mul_matrix(W, mean_temp[5].m, class6_mean_prime);
	gsl_matrix_mul_matrix(W, mean_temp[6].m, class7_mean_prime);
	gsl_matrix_mul_matrix(W, mean_temp[7].m, class8_mean_prime);



	gsl_matrix_fprintf (mean_1, class1_mean_prime, "%f");
	gsl_matrix_fprintf (mean_2, class2_mean_prime, "%f");
	gsl_matrix_fprintf (mean_3, class3_mean_prime, "%f");
	gsl_matrix_fprintf (mean_4, class4_mean_prime, "%f");
	gsl_matrix_fprintf (mean_5, class5_mean_prime, "%f");
	gsl_matrix_fprintf (mean_6, class6_mean_prime, "%f");
	gsl_matrix_fprintf (mean_7, class7_mean_prime, "%f");
	gsl_matrix_fprintf (mean_8, class8_mean_prime, "%f");


	
/****************************************************************************************************/


/**********************   求Sw prime  ***********************************************/
	
	gsl_matrix * temp_matrix1 = gsl_matrix_alloc (7, 18);
	gsl_matrix * temp_matrix2 = gsl_matrix_alloc (7, 18);
	gsl_matrix * temp_matrix3 = gsl_matrix_alloc (7, 18);
	gsl_matrix * temp_matrix4 = gsl_matrix_alloc (7, 18);
	gsl_matrix * temp_matrix5 = gsl_matrix_alloc (7, 18);
	gsl_matrix * temp_matrix6 = gsl_matrix_alloc (7, 18);
	gsl_matrix * temp_matrix7 = gsl_matrix_alloc (7, 18);
	gsl_matrix * temp_matrix8 = gsl_matrix_alloc (7, 18);

	gsl_matrix_mul_matrix(W, class1_sw->m, temp_matrix1);
	gsl_matrix_mul_matrix(W, class2_sw->m, temp_matrix2);
	gsl_matrix_mul_matrix(W, class3_sw->m, temp_matrix3);
	gsl_matrix_mul_matrix(W, class4_sw->m, temp_matrix4);
	gsl_matrix_mul_matrix(W, class5_sw->m, temp_matrix5);
	gsl_matrix_mul_matrix(W, class6_sw->m, temp_matrix6);
	gsl_matrix_mul_matrix(W, class7_sw->m, temp_matrix7);
	gsl_matrix_mul_matrix(W, class8_sw->m, temp_matrix8);

	gsl_matrix * W_tran = gsl_matrix_alloc (18, 7);
	gsl_matrix_transpose_memcpy (W_tran, W);

	gsl_matrix_mul_matrix(temp_matrix1, W_tran, sw_class1_prime);
	gsl_matrix_mul_matrix(temp_matrix2, W_tran, sw_class2_prime);
	gsl_matrix_mul_matrix(temp_matrix3, W_tran, sw_class3_prime);
	gsl_matrix_mul_matrix(temp_matrix4, W_tran, sw_class4_prime);
	gsl_matrix_mul_matrix(temp_matrix5, W_tran, sw_class5_prime);
	gsl_matrix_mul_matrix(temp_matrix6, W_tran, sw_class6_prime);
	gsl_matrix_mul_matrix(temp_matrix7, W_tran, sw_class7_prime);
	gsl_matrix_mul_matrix(temp_matrix8, W_tran, sw_class8_prime);

	gsl_matrix_scale (sw_class1_prime, 0.04);
	gsl_matrix_scale (sw_class2_prime, 0.04);
	gsl_matrix_scale (sw_class3_prime, 0.04);
	gsl_matrix_scale (sw_class4_prime, 0.04);
	gsl_matrix_scale (sw_class5_prime, 0.04);
	gsl_matrix_scale (sw_class6_prime, 0.04);
	gsl_matrix_scale (sw_class7_prime, 0.04);
	gsl_matrix_scale (sw_class8_prime, 0.04);



	for(i=0;i<7;i++)
	{
		for(j=0;j<7;j++)
		{
			fprintf(sw_1, "%f   ", gsl_matrix_get(sw_class1_prime, i, j));
			fprintf(sw_2, "%f   ", gsl_matrix_get(sw_class2_prime, i, j));
			fprintf(sw_3, "%f   ", gsl_matrix_get(sw_class3_prime, i, j));
			fprintf(sw_4, "%f   ", gsl_matrix_get(sw_class4_prime, i, j));
			fprintf(sw_5, "%f   ", gsl_matrix_get(sw_class5_prime, i, j));
			fprintf(sw_6, "%f   ", gsl_matrix_get(sw_class6_prime, i, j));
			fprintf(sw_7, "%f   ", gsl_matrix_get(sw_class7_prime, i, j));
			fprintf(sw_8, "%f   ", gsl_matrix_get(sw_class8_prime, i, j));
		}
		fprintf(sw_1, "\n");
		fprintf(sw_1, "\n");
		fprintf(sw_2, "\n");
		fprintf(sw_2, "\n");
		fprintf(sw_3, "\n");
		fprintf(sw_3, "\n");
		fprintf(sw_4, "\n");
		fprintf(sw_4, "\n");
		fprintf(sw_5, "\n");
		fprintf(sw_5, "\n");
		fprintf(sw_6, "\n");
		fprintf(sw_6, "\n");
		fprintf(sw_7, "\n");
		fprintf(sw_7, "\n");
		fprintf(sw_8, "\n");
		fprintf(sw_8, "\n");
	}


	
/****************************************************************************************************/


	
	return evec;
	//gsl_vector_complex_free(eval);
	//gsl_matrix_complex_free(evec);
/************************************************************************************/

}

/************************** 矩陣乘法 *********************/
void gsl_matrix_mul_matrix(const gsl_matrix * a, const gsl_matrix * b, gsl_matrix * c)
{
	assert(a->size2 == b->size1 && a->size1 == c->size1 && b->size2 == c->size2);

	/* set all elements in c to zero */
	gsl_matrix_set_zero(c);

	/* calculate mul */
	int i, j, k;
	for(i = 0; i < a->size1; i++)
	{
   		for(j = 0; j < b->size2; j++)
   		{
       			for(k = 0; k < a->size2; k++)
       			{
           			gsl_matrix_set(c, i, j, gsl_matrix_get(c, i, j) + gsl_matrix_get(a, i, k) * gsl_matrix_get(b, k, j));
       			}
   		}
	}
}