C++ graph::get_descmatrix方法代码示例

match_result algorithm_ca::match(graph& g, graph& h,gsl_matrix* gm_P_i,gsl_matrix* _gm_ldh,double dalpha_ldh)
{
	double dhung_max=get_param_d("hungarian_max");
	bool bblast_match_end=(get_param_i("blast_match_proj")==1);	
	bool bblast_match=(get_param_i("blast_match")==1);
	double dfw_xeps=get_param_d("algo_fw_xeps");
	double bgreedy=(get_param_i("hungarian_greedy")==1);
	double dfvalue;
	if (bverbose)
		*gout<<"C-adaptation strategy"<<std::endl;	
	//some duplicate variables
	gsl_matrix* gm_Ag_d=g.get_descmatrix(cdesc_matrix);
	gsl_matrix* gm_Ah_d=h.get_descmatrix(cdesc_matrix);
	gsl_matrix * gm_C_a=gsl_matrix_alloc(N,N);
	gsl_matrix * gm_C_up=gsl_matrix_alloc(N,N);
	gsl_matrix * gm_P=gsl_matrix_alloc(N,N);
	gsl_matrix * gm_temp=gsl_matrix_alloc(N,N);
	//initialization
	gsl_matrix_memcpy(gm_C_a,gm_P_i);
	//cycle over linear combination of linear and non linear terms
	double astep=1e-30;bool bcont=true;
	if (dalpha_ldh>0)
		gsl_matrix_scale(gm_ldh,dalpha_ldh);

	while (bcont){
	bcont=(astep<1);
	//optimal permutation
	gsl_matrix_scale(gm_C_a,-dhung_max);
	gsl_matrix_transpose(gm_C_a);
	gsl_matrix_hungarian(gm_C_a,gm_P,NULL,NULL,false,(bblast_match?gm_ldh:NULL),bgreedy);
	gsl_matrix_transpose(gm_C_a);
	gsl_matrix_scale(gm_C_a,-1.0/dhung_max);
	//C-matrix update
	gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1,gm_Ag_d,gm_P,0,gm_temp);
	gsl_blas_dgemm(CblasNoTrans,CblasNoTrans,1,gm_temp,gm_Ah_d,0,gm_C_up);
	gsl_matrix_scale(gm_C_up,1-dalpha_ldh);
	if (dalpha_ldh>0) gsl_matrix_add(gm_C_up,gm_ldh);
	gsl_matrix_scale(gm_C_a,1-astep);
	gsl_matrix_scale(gm_C_up,astep);
	gsl_matrix_add(gm_C_a,gm_C_up);
	dfvalue=f_qcv(gm_Ag_d,gm_Ah_d,gm_P,gm_temp,true);
	if (bverbose) *gout<<"astep="<<astep<<", f="<<dfvalue<<std::endl;
	astep+=0.01;//0.01*astep;
	astep=(astep<1)?astep:1;
	};
	gsl_matrix_free(gm_Ag_d);
	gsl_matrix_free(gm_Ah_d);
	gsl_matrix_free(gm_temp);
	gsl_matrix_free(gm_C_up);
	gsl_matrix_free(gm_C_a);
	if (dalpha_ldh>0)
		gsl_matrix_scale(gm_ldh,1.0/dalpha_ldh);
	
	match_result mres;
	mres.gm_P=gm_P;
	
	//initial score
	mres.vd_trace.push_back(graph_dist(g,h,cscore_matrix));
	//final score
	mres.vd_trace.push_back(graph_dist(g,h,gm_P,cscore_matrix));
	//other output parameters
	mres.dres=mres.vd_trace.at(1);
	mres.inum_iteration=2;
	//transpose matrix save
	mres.gm_P=gm_P;
	mres.gm_P_exact=NULL;
	return mres;
}

match_result algorithm_umeyama::match(graph& g, graph& h,gsl_matrix* gm_P_i,gsl_matrix* gm_ldh,double dalpha_ldh)
{
	if (bverbose)
		*gout<<"Umeyama algorithm"<<std::endl;	
	bool bblast_match_end=(get_param_i("blast_match_proj")==1);
	//some duplicate variables
	gsl_matrix* gm_Ag_d=g.get_descmatrix(cdesc_matrix);
	gsl_matrix* gm_Ah_d=h.get_descmatrix(cdesc_matrix);
	if (pdebug.ivalue) gsl_matrix_printout(gm_Ag_d,"Ag",pdebug.strvalue); 
	if (pdebug.ivalue) gsl_matrix_printout(gm_Ah_d,"Ah",pdebug.strvalue);
	//memory allocation
	gsl_eigen_symmv_workspace * gesw= gsl_eigen_symmv_alloc (N);
	gsl_vector* eval_g=gsl_vector_alloc(N);
	gsl_vector* eval_h=gsl_vector_alloc(N);
	gsl_matrix* evec_g=gsl_matrix_alloc(N,N);
	gsl_matrix* evec_h=gsl_matrix_alloc(N,N);
	if (bverbose) *gout<<"Memory allocation finished"<<std::endl;
	//eigenvalues and eigenvectors for both matrices
	gsl_eigen_symmv (gm_Ag_d, eval_g,evec_g,gesw);
	if (bverbose) *gout<<"Ag eigen vectors"<<std::endl;
	gsl_eigen_symmv (gm_Ah_d, eval_h,evec_h,gesw);
	
	gsl_matrix_free(gm_Ag_d);
	gsl_matrix_free(gm_Ah_d);
	
	if (bverbose) *gout<<"Ah eigen vectors"<<std::endl;
	gsl_eigen_symmv_sort (eval_g, evec_g, GSL_EIGEN_SORT_VAL_DESC);
	gsl_eigen_symmv_sort (eval_h, evec_h, GSL_EIGEN_SORT_VAL_DESC);
	
	if (pdebug.ivalue){ gsl_matrix_printout(eval_g,"eval_g",pdebug.strvalue); 
			    gsl_matrix_printout(eval_h,"eval_h",pdebug.strvalue); 
	 		    gsl_matrix_printout(evec_g,"evec_g",pdebug.strvalue);
			    gsl_matrix_printout(evec_h,"evec_h",pdebug.strvalue);};
	gsl_matrix_abs(evec_g);
	gsl_matrix_abs(evec_h);
	if (pdebug.ivalue) gsl_matrix_printout(evec_g,"abs(evec_g)",pdebug.strvalue);
	if (pdebug.ivalue) gsl_matrix_printout(evec_h,"abs(evec_h)",pdebug.strvalue);  
	//loss matrix construction
	gsl_matrix* C=gsl_matrix_alloc(N,N);
	if (bverbose) *gout<<"Loss function matrix allocation"<<std::endl;
	gsl_blas_dgemm(CblasNoTrans,CblasTrans,1,evec_g,evec_h,0,C);
	if (pdebug.ivalue) gsl_matrix_printout(C,"C=abs(evec_g)*abs(evec_h')",pdebug.strvalue); 
	//label cost matrix
	update_C_hungarian(C,-1);
	//scaling for hungarian
	double dscale_factor =gsl_matrix_max_abs(C);
	dscale_factor=(dscale_factor>EPSILON)?dscale_factor:EPSILON;
	dscale_factor=10000/dscale_factor;
	gsl_matrix_scale(C,-dscale_factor);
	gsl_matrix_transpose(C);
	if (pdebug.ivalue) gsl_matrix_printout(C,"scale(C)",pdebug.strvalue); 
	gsl_matrix* gm_P=gsl_matrix_alloc(N,N);
        gsl_matrix_hungarian(C,gm_P,NULL,NULL,false,(bblast_match_end?gm_ldh:NULL),false);
	if (pdebug.ivalue) gsl_matrix_printout(gm_P,"gm_P",pdebug.strvalue); 
	if (bverbose) *gout<<"Hungarian solved"<<std::endl;
	match_result mres;
	mres.gm_P=gm_P;
	

	//initial score
	mres.vd_trace.push_back(graph_dist(g,h,cscore_matrix));
	//final score
	mres.vd_trace.push_back(graph_dist(g,h,gm_P,cscore_matrix));
	//other output parameters
	mres.dres=mres.vd_trace.at(1);
	mres.inum_iteration=2;
	//transpose matrix save
	mres.gm_P=gm_P;
	mres.gm_P_exact=NULL;
	return mres;
}