SHOGUN
v3.2.0
|
00001 /* 00002 * This program is free software; you can redistribute it and/or modify 00003 * it under the terms of the GNU General Public License as published by 00004 * the Free Software Foundation; either version 3 of the License, or 00005 * (at your option) any later version. 00006 * 00007 * Copyright (W) 2012 Sergey Lisitsyn 00008 */ 00009 00010 #include <shogun/statistics/KernelMeanMatching.h> 00011 #include <shogun/lib/external/libqp.h> 00012 00013 00014 static float64_t* kmm_K = NULL; 00015 static int32_t kmm_K_ld = 0; 00016 00017 static const float64_t* kmm_get_col(uint32_t i) 00018 { 00019 return kmm_K + kmm_K_ld*i; 00020 } 00021 00022 namespace shogun 00023 { 00024 CKernelMeanMatching::CKernelMeanMatching() : 00025 CSGObject(), m_kernel(NULL) 00026 { 00027 } 00028 00029 CKernelMeanMatching::CKernelMeanMatching(CKernel* kernel, SGVector<index_t> training_indices, 00030 SGVector<index_t> test_indices) : 00031 CSGObject(), m_kernel(NULL) 00032 { 00033 set_kernel(kernel); 00034 set_training_indices(training_indices); 00035 set_test_indices(test_indices); 00036 } 00037 00038 SGVector<float64_t> CKernelMeanMatching::compute_weights() 00039 { 00040 int32_t i,j; 00041 ASSERT(m_kernel) 00042 ASSERT(m_training_indices.vlen) 00043 ASSERT(m_test_indices.vlen) 00044 00045 int32_t n_tr = m_training_indices.vlen; 00046 int32_t n_te = m_test_indices.vlen; 00047 00048 SGVector<float64_t> weights(n_tr); 00049 weights.zero(); 00050 00051 kmm_K = SG_MALLOC(float64_t, n_tr*n_tr); 00052 kmm_K_ld = n_tr; 00053 float64_t* diag_K = SG_MALLOC(float64_t, n_tr); 00054 for (i=0; i<n_tr; i++) 00055 { 00056 float64_t d = m_kernel->kernel(m_training_indices[i], m_training_indices[i]); 00057 diag_K[i] = d; 00058 kmm_K[i*n_tr+i] = d; 00059 for (j=i+1; j<n_tr; j++) 00060 { 00061 d = m_kernel->kernel(m_training_indices[i],m_training_indices[j]); 00062 kmm_K[i*n_tr+j] = d; 00063 kmm_K[j*n_tr+i] = d; 00064 } 00065 } 00066 float64_t* kappa = SG_MALLOC(float64_t, n_tr); 00067 for (i=0; i<n_tr; i++) 00068 { 00069 float64_t avg = 0.0; 00070 for (j=0; j<n_te; j++) 00071 avg+= m_kernel->kernel(m_training_indices[i],m_test_indices[j]); 00072 00073 avg *= float64_t(n_tr)/n_te; 00074 kappa[i] = -avg; 00075 } 00076 float64_t* a = SG_MALLOC(float64_t, n_tr); 00077 for (i=0; i<n_tr; i++) a[i] = 1.0; 00078 float64_t* LB = SG_MALLOC(float64_t, n_tr); 00079 float64_t* UB = SG_MALLOC(float64_t, n_tr); 00080 float64_t B = 2.0; 00081 for (i=0; i<n_tr; i++) 00082 { 00083 LB[i] = 0.0; 00084 UB[i] = B; 00085 } 00086 for (i=0; i<n_tr; i++) 00087 weights[i] = 1.0/float64_t(n_tr); 00088 00089 libqp_state_T result = 00090 libqp_gsmo_solver(&kmm_get_col,diag_K,kappa,a,1.0,LB,UB,weights,n_tr,1000,1e-9,NULL); 00091 00092 SG_DEBUG("libqp exitflag=%d, %d iterations passed, primal objective=%f\n", 00093 result.exitflag,result.nIter,result.QP); 00094 00095 SG_FREE(kappa); 00096 SG_FREE(a); 00097 SG_FREE(LB); 00098 SG_FREE(UB); 00099 SG_FREE(diag_K); 00100 SG_FREE(kmm_K); 00101 00102 return weights; 00103 } 00104 00105 }