基於OpenCV ML的Affinity Propagation原型實現

春節實現的APCluster，會在最近有空的時候和RBM一起嘗試提交到OpenCV ml中。這是一個原型代碼，缺少很多錯誤控制的流程，但已支援CvSparseMat。

mlapcluster.h

#ifndef GUARD_mlapcluster_h<br />#define GUARD_mlapcluster_h<br />#include <ml.h><br />#define CV_TYPE_NAME_ML_APCLUSTER"opencv-ml-affinity-propagation-cluster"<br />struct CV_EXPORTS CvAPCParams<br />{<br />int maxiteration;<br />int stopcriterion;<br />double lambda;<br />CvAPCParams()<br />: maxiteration(2000),<br /> stopcriterion(200),<br /> lambda(0.9)<br />{}<br />CvAPCParams( int _maxiteration, int _stopcriterion, double _lambda )<br />: maxiteration(_maxiteration),<br /> stopcriterion(_stopcriterion),<br /> lambda(_lambda)<br />{}<br />};<br />class CV_EXPORTS CvAPCluster : public CvStatModel<br />{<br />private:<br />CvAPCParams* params;<br />public:<br />CvAPCluster( CvAPCParams* _params )<br />: params(_params)<br />{}<br />virtual ~CvAPCluster()<br />{<br />clear();<br />}<br />virtual bool train( const CvMat* _train_data, const CvMat* _response );<br />virtual bool train( const CvSparseMat* _train_data, const CvMat* _response );<br />virtual void clear();<br />virtual void write( CvFileStorage* fs, const char* name );<br />virtual void read( CvFileStorage* fs, CvFileNode* root_node );<br />};<br />#endif

mlapcluster.cpp

#include "mlapcluster.h"<br />#include <limits><br />bool CvAPCluster::train( const CvMat* similarity,<br /> const CvMat* labels )<br />{<br />int stop = 0;<br />CV_FUNCNAME( "CvAPCluster::train" );<br />__BEGIN__;<br />CvMat* rsp;<br />CvMat* avl;<br />int* rspmaxidx;<br />double* rspmaxval;<br />double* rspmaxval2;<br />double* diagrsp;<br />double& lam = params->lambda;<br />double _lam = 1.-lam;<br />int num;<br />if ( !CV_IS_MAT(similarity) || similarity->rows != similarity->cols || CV_MAT_TYPE( similarity->type ) != CV_64FC1 )<br />CV_ERROR( CV_StsBadArg,<br />"similarity matrix is a double-point matrix with equal rows and cols" );</p><p>num = similarity->rows;<br />if ( !CV_IS_MAT(labels) || labels->cols != num || CV_MAT_TYPE( labels->type ) != CV_32SC1 )<br />CV_ERROR( CV_StsBadArg,<br />"labels array (when passed) must be a valid 1d integer vector of <sample_count> elements" );<br />rsp = cvCreateMat( num, num, CV_64FC1 );<br />avl = cvCreateMat( num, num, CV_64FC1 );<br />cvZero( rsp );<br />cvZero( avl );</p><p>rspmaxidx = (int*)cvAlloc( num*sizeof(rspmaxidx[0]) );<br />rspmaxval = (double*)cvAlloc( num*sizeof(rspmaxval[0]) );<br />rspmaxval2 = (double*)cvAlloc( num*sizeof(rspmaxval2[0]) );<br />diagrsp = (double*)cvAlloc( num*sizeof(diagrsp[0]) );<br />for ( int k = 0; k < params->maxiteration; k++ )<br />{<br />double* rsp_vec = rsp->data.db;<br />double* avl_vec = avl->data.db;<br />double* sim_vec = similarity->data.db;<br />double* sim_vec2 = sim_vec;<br />for ( int i = 0; i < num; i++ )<br />{<br />int yk = -1;<br />double y = -std::numeric_limits<double>::max(), y2 = -std::numeric_limits<double>::max();<br />for ( int j = 0; j < num; j++ )<br />{<br />double t = *sim_vec2 + *avl_vec;<br />if ( t > y )<br />{<br />y = t;<br />yk = j;<br />} else if ( t > y2 )<br />y2 = t;<br />sim_vec2++;<br />avl_vec++;<br />}<br />for ( int j = 0; j < num; j++ )<br />{<br />if ( j != yk )<br />*rsp_vec = *rsp_vec*lam+(*sim_vec-y)*_lam;<br />else<br />*rsp_vec = *rsp_vec*lam+(*sim_vec-y2)*_lam;<br />rsp_vec++;<br />sim_vec++;<br />}<br />}<br />int* rspmaxidx_vec = rspmaxidx;<br />double* rspmaxval_vec = rspmaxval;<br />double* rspmaxval2_vec = rspmaxval2;<br />for ( int i = 0; i < num; i++, rspmaxidx_vec++, rspmaxval_vec++, rspmaxval2_vec++ )<br />{<br />*rspmaxidx_vec = -1;<br />*rspmaxval_vec = 0;<br />*rspmaxval2_vec = 0;<br />}<br />rsp_vec = rsp->data.db;<br />double* diagrsp_vec = diagrsp;<br />for ( int i = 0; i < num; i++, diagrsp_vec++ )<br />{<br />rspmaxidx_vec = rspmaxidx;<br />rspmaxval_vec = rspmaxval;<br />rspmaxval2_vec = rspmaxval2;<br />for ( int j = 0; j < num; j++, rspmaxidx_vec++, rspmaxval_vec++, rspmaxval2_vec++ )<br />{<br />if ( i != j )<br />{<br />if ( *rsp_vec > *rspmaxval_vec )<br />{<br />*rspmaxidx_vec = i;<br />*rspmaxval_vec = *rsp_vec;<br />} else if ( *rsp_vec > *rspmaxval2_vec )<br />*rspmaxval2_vec = *rsp_vec;<br />} else<br />*diagrsp_vec = *rsp_vec;<br />rsp_vec++;<br />}<br />}<br />avl_vec = avl->data.db;<br />diagrsp_vec = diagrsp;<br />for ( int i = 0; i < num; i++, diagrsp_vec++ )<br />{<br />rspmaxidx_vec = rspmaxidx;<br />rspmaxval_vec = rspmaxval;<br />rspmaxval2_vec = rspmaxval2;<br />for ( int j = 0; j < num; j++, rspmaxidx_vec++, rspmaxval_vec++, rspmaxval2_vec++ )<br />{<br />double tmp;<br />if ( i != j )<br />{<br />if ( i != *rspmaxidx_vec )<br />tmp = *diagrsp_vec + *rspmaxval_vec;<br />else<br />tmp = *diagrsp_vec + *rspmaxval2_vec;<br />if ( tmp > 0 )<br />tmp = 0;<br />} else<br />tmp = *rspmaxval_vec;<br />*avl_vec = *avl_vec*lam+tmp*_lam;<br />avl_vec++;<br />}<br />}<br />stop++;<br />int* cls_vec = labels->data.i;<br />rsp_vec = rsp->data.db;<br />avl_vec = avl->data.db;<br />for ( int i = 0; i < num; i++, cls_vec++ )<br />{<br />int maxidx = i;<br />double maxval = -std::numeric_limits<double>::max();<br />for ( int j = 0; j < num; j++ )<br />{<br />double t = *rsp_vec + *avl_vec;<br />if ( t > maxval )<br />{<br />maxval = t;<br />maxidx = j;<br />}<br />rsp_vec++;<br />avl_vec++;<br />}<br />if ( *cls_vec != maxidx )<br />stop = 0;<br />*cls_vec = maxidx;<br />}<br />if ( stop > params->stopcriterion )<br />break;<br />}<br />cvFree( &diagrsp );<br />cvFree( &rspmaxval2 );<br />cvFree( &rspmaxval );<br />cvFree( &rspmaxidx );<br />cvReleaseMat( &avl );<br />cvReleaseMat( &rsp );<br />__END__;<br />return stop > params->stopcriterion;<br />}<br />struct CvSparseNode2D<br />{<br />int i;<br />int k;<br />double val;<br />};<br />bool<br />CvAPCluster::train( const CvSparseMat* similarity,<br /> const CvMat* labels )<br />{<br />int stop = 0;<br />CV_FUNCNAME( "CvAPCluster::train" );<br />__BEGIN__;<br />int* psize;<br />int* psize_vec;<br />int num;<br />double* rsp;<br />double* avl;<br />double* rsp_vec;<br />double* avl_vec;<br />CvSparseNode2D* nodes;<br />CvSparseNode2D** node_entries;<br />int* rspmaxidx;<br />double* rspmaxval;<br />double* rspmaxval2;<br />double* diagrsp;<br />int total = 0;<br />double& lam = params->lambda;<br />double _lam = 1.-lam;<br />if ( !CV_IS_SPARSE_MAT(similarity) || similarity->dims != 2 || similarity->size[0] != similarity->size[1] || CV_MAT_TYPE( similarity->type ) != CV_64FC1 )<br />CV_ERROR( CV_StsBadArg,<br />"similarity matrix is a double-point sparse matrix with equal rows and cols" );<br />num = similarity->size[0];<br />if ( !CV_IS_MAT(labels) || labels->cols != num || CV_MAT_TYPE( labels->type ) != CV_32SC1 )<br />CV_ERROR( CV_StsBadArg,<br />"labels array (when passed) must be a valid 1d integer vector of <sample_count> elements" );<br />psize = (int*)cvAlloc( num*sizeof(psize[0]) );<br />psize_vec = psize;<br />for ( int i = 0; i < num; i++, psize_vec++ )<br />*psize_vec = 0;<br />CvSparseMatIterator mat_iterator;<br />CvSparseNode* node;<br />node = cvInitSparseMatIterator( similarity, &mat_iterator );<br />for ( ; node != 0; node = cvGetNextSparseNode( &mat_iterator ) )<br />{<br />const int* idx = CV_NODE_IDX( similarity, node );<br />psize[idx[0]]++;<br />total++;<br />}<br />rsp = (double*)cvAlloc( total*sizeof(rsp[0]) );<br />avl = (double*)cvAlloc( total*sizeof(avl[0]) );<br />rsp_vec = rsp;<br />avl_vec = avl;<br />nodes = (CvSparseNode2D*)cvAlloc( total*sizeof(nodes[0]) );<br />node_entries = (CvSparseNode2D**)cvAlloc( num*sizeof(node_entries[0]) );<br />psize_vec = psize;<br />node_entries[0] = nodes;<br />for ( int i = 1; i < num; i++, psize_vec++ )<br />node_entries[i] = node_entries[i-1] + *psize_vec;<br />node = cvInitSparseMatIterator( similarity, &mat_iterator );<br />for ( ; node != 0; node = cvGetNextSparseNode( &mat_iterator ), rsp_vec++, avl_vec++ )<br />{<br />*rsp_vec = *avl_vec = 0;<br />const int* idx = CV_NODE_IDX( similarity, node );<br />CvSparseNode2D*& node_entry = node_entries[idx[0]];<br />node_entry->i = idx[0];<br />node_entry->k = idx[1];<br />node_entry->val = *(double*)CV_NODE_VAL( similarity, node );<br />node_entry++;<br />}<br />cvFree( &node_entries );<br />rspmaxidx = (int*)cvAlloc( num*sizeof(rspmaxidx[0]) );<br />rspmaxval = (double*)cvAlloc( num*sizeof(rspmaxval[0]) );<br />rspmaxval2 = (double*)cvAlloc( num*sizeof(rspmaxval2[0]) );<br />diagrsp = (double*)cvAlloc( num*sizeof(diagrsp[0]) );<br />for ( int k = 0; k < params->maxiteration; k++ )<br />{<br />CvSparseNode2D* nodes_vec = nodes;<br />psize_vec = psize;<br />CvSparseNode2D* nodes_vec2 = nodes;<br />avl_vec = avl;<br />rsp_vec = rsp;<br />for ( int i = 0; i < num; i++, psize_vec++ )<br />{<br />int yk = -1;<br />double y = -std::numeric_limits<double>::max(), y2 = -std::numeric_limits<double>::max();<br />for ( int j = 0; j < *psize_vec; j++ )<br />{<br />double t = nodes_vec2->val + *avl_vec;<br />if ( t > y )<br />{<br />y = t;<br />yk = nodes_vec2->k;<br />} else if ( t > y2 )<br />y2 = t;<br />nodes_vec2++;<br />avl_vec++;<br />}<br />for ( int j = 0; j < *psize_vec; j++ )<br />{<br />if ( nodes_vec->k != yk )<br />*rsp_vec = *rsp_vec*lam+(nodes_vec->val-y)*_lam;<br />else<br />*rsp_vec = *rsp_vec*lam+(nodes_vec->val-y2)*_lam;<br />nodes_vec++;<br />rsp_vec++;<br />}<br />}<br />int* rspmaxidx_vec = rspmaxidx;<br />double* rspmaxval_vec = rspmaxval;<br />double* rspmaxval2_vec = rspmaxval2;<br />for ( int i = 0; i < num; i++, rspmaxidx_vec++, rspmaxval_vec++, rspmaxval2_vec++ )<br />{<br />*rspmaxidx_vec = -1;<br />*rspmaxval_vec = 0;<br />*rspmaxval2_vec = 0;<br />}<br />rsp_vec = rsp;<br />psize_vec = psize;<br />nodes_vec = nodes;<br />double* diagrsp_vec = diagrsp;<br />for ( int i = 0; i < num; i++, psize_vec++, diagrsp_vec++ )<br />for ( int j = 0; j < *psize_vec; j++ )<br />{<br />if ( i != nodes_vec->k )<br />{<br />if ( *rsp_vec > rspmaxval[nodes_vec->k] )<br />{<br />rspmaxidx[nodes_vec->k] = i;<br />rspmaxval[nodes_vec->k] = *rsp_vec;<br />} else if ( *rsp_vec > rspmaxval2[nodes_vec->k] )<br />rspmaxval2[nodes_vec->k] = *rsp_vec;<br />} else<br />*diagrsp_vec = *rsp_vec;<br />nodes_vec++;<br />rsp_vec++;<br />}<br />avl_vec = avl;<br />diagrsp_vec = diagrsp;<br />psize_vec = psize;<br />nodes_vec = nodes;<br />for ( int i = 0; i < num; i++, psize_vec++, diagrsp_vec++ )<br />for ( int j = 0; j < *psize_vec; j++ )<br />{<br />double tmp;<br />if ( i != nodes_vec->k )<br />{<br />if ( i != rspmaxidx[nodes_vec->k] )<br />tmp = *diagrsp_vec + rspmaxval[nodes_vec->k];<br />else<br />tmp = *diagrsp_vec + rspmaxval2[nodes_vec->k];<br />if ( tmp > 0 )<br />tmp = 0;<br />} else<br />tmp = rspmaxval[nodes_vec->k];<br />*avl_vec = *avl_vec*lam+tmp*_lam;<br />nodes_vec++;<br />avl_vec++;<br />}<br />stop++;<br />int* cls_vec = labels->data.i;<br />rsp_vec = rsp;<br />avl_vec = avl;<br />psize_vec = psize;<br />nodes_vec = nodes;<br />for ( int i = 0; i < num; i++, psize_vec++, cls_vec++ )<br />{<br />int maxidx = i;<br />double maxval = -std::numeric_limits<double>::max();<br />for ( int j = 0; j < *psize_vec; j++ )<br />{<br />double t = *rsp_vec + *avl_vec;<br />if ( t > maxval )<br />{<br />maxval = t;<br />maxidx = nodes_vec->k;<br />}<br />nodes_vec++;<br />rsp_vec++;<br />avl_vec++;<br />}<br />if ( *cls_vec != maxidx )<br />stop = 0;<br />*cls_vec = maxidx;<br />}<br />if ( stop > params->stopcriterion )<br />break;<br />}<br />cvFree( &diagrsp );<br />cvFree( &rspmaxval2 );<br />cvFree( &rspmaxval );<br />cvFree( &rspmaxidx );</p><p>cvFree( &rsp );<br />cvFree( &avl );<br />cvFree( &psize );<br />cvFree( &nodes );<br />__END__;<br />return stop > params->stopcriterion;<br />}<br />void<br />CvAPCluster::clear()<br />{<br />}<br />void<br />CvAPCluster::write( CvFileStorage* fs,<br /> const char* name )<br />{<br />}<br />void<br />CvAPCluster::read( CvFileStorage* fs,<br /> CvFileNode* root_node )<br />{<br />}

測試例子

#include "mlapcluster.h"<br />#include <iostream><br />int main()<br />{<br />CvAPCParams* params = new CvAPCParams(2000, 200, 0.5);<br />CvAPCluster* apcluster = new CvAPCluster(params);<br />//CvMat* dist = cvCreateMat( 25, 25, CV_64FC1 );<br />int sizes[] = {25, 25};<br />CvSparseMat* dist = cvCreateSparseMat( 2, sizes, CV_64FC1 );<br />cvZero( dist );<br />freopen( "ToyProblemSimilarities.txt", "r", stdin );<br />double total = 0;<br />for ( int k = 0; k < 25*24; k++ )<br />{<br />int i, j;<br />double d;<br />std::cin>>i>>j>>d;<br />cvSetReal2D( dist, i-1, j-1, d );<br />total += d;<br />}<br />total = -15.561256;<br />//total = total / (25*24);<br />printf("preference: %f/n", total);<br />for ( int k = 0; k < 25; k++ )<br />cvSetReal2D( dist, k, k, total );<br />CvMat* response = cvCreateMat( 1, 25, CV_32SC1 );<br /> long long t;<br />t = cvGetTickCount();<br />apcluster->train( dist, response );<br /> t = cvGetTickCount() - t;<br />printf( "apcluster trained in %lld ms./n", t/1000000 );<br />for ( int i = 0; i < 25; i++ )<br />printf("%d ", response->data.i[i]);<br />printf("/n");<br />return 0;<br />}