26 using namespace shogun;
30 #define MAX_LOG_TABLE_SIZE 10*1024*1024
31 #define LOG_TABLE_PRECISION 1e-6
33 #define MAX_LOG_TABLE_SIZE 123*1024*1024
34 #define LOG_TABLE_PRECISION 1e-15
36 int32_t CMath::LOGACCURACY = 0;
37 #endif // USE_LOGCACHE
62 LOGRANGE=CMath::determine_logrange();
63 LOGACCURACY=CMath::determine_logaccuracy(
LOGRANGE);
92 for (int32_t i=0; i<n; i++)
93 SG_SPRINT(
"%d%s", vector[i], i==n-1?
"" :
",");
102 for (int32_t i=0; i<n; i++)
103 SG_SPRINT(
"%d%s", vector[i], i==n-1?
"" :
",");
112 for (int32_t i=0; i<n; i++)
113 SG_SPRINT(
"%lld%s", vector[i], i==n-1?
"" :
",");
122 for (int32_t i=0; i<n; i++)
123 SG_SPRINT(
"%llu%s", vector[i], i==n-1?
"" :
",");
132 for (int32_t i=0; i<n; i++)
133 SG_SPRINT(
"%g%s", vector[i], i==n-1?
"" :
",");
142 for (int32_t i=0; i<n; i++)
143 SG_SPRINT(
"%.18g%s", vector[i], i==n-1?
"" :
",");
152 for (int32_t i=0; i<n; i++)
153 SG_SPRINT(
"%.36Lg%s", (
long double) vector[i], i==n-1?
"" :
",");
159 const int32_t* matrix, int32_t rows, int32_t cols,
const char* name)
161 ASSERT(rows>=0 && cols>=0);
163 for (int32_t i=0; i<rows; i++)
166 for (int32_t j=0; j<cols; j++)
168 j==cols-1?
"" :
",");
176 const float64_t* matrix, int32_t rows, int32_t cols,
const char* name)
178 ASSERT(rows>=0 && cols>=0);
180 for (int32_t i=0; i<rows; i++)
183 for (int32_t j=0; j<cols; j++)
184 SG_SPRINT(
"\t%.18g%s", (
double) matrix[j*rows+i],
185 j==cols-1?
"" :
",");
193 const float32_t* matrix, int32_t rows, int32_t cols,
const char* name)
195 ASSERT(rows>=0 && cols>=0);
197 for (int32_t i=0; i<rows; i++)
200 for (int32_t j=0; j<cols; j++)
201 SG_SPRINT(
"\t%.18g%s", (
float) matrix[j*rows+i],
202 j==cols-1?
"" :
",");
216 for (int32_t i=0; i<len; i++)
243 for (int32_t i=0; i<3+2; i++)
250 for (int32_t i=0; i<3*2; i++)
253 log_denomf+=gamma(table.
matrix[i]+1);
256 floatmax_t prob_table_log=log_nom - log_denom;
262 for (int32_t k=0; k<=dim1; k++)
266 x[0 + 0*2 + counter*2*3] = k;
267 x[0 + 1*2 + counter*2*3] = l;
268 x[0 + 2*2 + counter*2*3] = m[0] - x[0 + 0*2 + counter*2*3] - x[0 + 1*2 + counter*2*3];
269 x[1 + 0*2 + counter*2*3] = m[2] - x[0 + 0*2 + counter*2*3];
270 x[1 + 1*2 + counter*2*3] = m[3] - x[0 + 1*2 + counter*2*3];
271 x[1 + 2*2 + counter*2*3] = m[4] - x[0 + 2*2 + counter*2*3];
278 #ifdef DEBUG_FISHER_TABLE
283 SG_SPRINT(
"prob_table_log=%.18Lg\n", prob_table_log);
284 SG_SPRINT(
"log_denomf=%.18g\n", log_denomf);
285 SG_SPRINT(
"log_denom=%.18Lg\n", log_denom);
289 #endif // DEBUG_FISHER_TABLE
295 for (int32_t k=0; k<counter; k++)
297 for (int32_t j=0; j<3; j++)
299 for (int32_t i=0; i<2; i++)
304 for (int32_t i=0; i<counter; i++)
305 log_denom_vec[i]=log_nom-log_denom_vec[i];
307 #ifdef DEBUG_FISHER_TABLE
309 #endif // DEBUG_FISHER_TABLE
313 for (int32_t i=0; i<counter; i++)
315 if (log_denom_vec[i]<=prob_table_log)
319 #ifdef DEBUG_FISHER_TABLE
320 SG_SPRINT(
"nonrand_p=%.18g\n", nonrand_p);
322 #endif // DEBUG_FISHER_TABLE
335 int32_t CMath::determine_logrange()
346 SG_SINFO(
"determined range for x in table log(1+exp(-x)) is:%d (error:%G)\n",i,acc);
350 int32_t CMath::determine_logaccuracy(int32_t range)
352 range=MAX_LOG_TABLE_SIZE/range/((int)
sizeof(
float64_t));
353 SG_SINFO(
"determined accuracy for x in table log(1+exp(-x)) is:%d (error:%G)\n",range,1.0/(
double) range);
358 void CMath::init_log_table()
360 for (int32_t i=0; i< LOGACCURACY*
LOGRANGE; i++)
368 if (a[0]==-1) return ;
371 changed=0; int32_t i=0 ;
372 while ((a[(i+1)*cols]!=-1) && (a[(i+1)*cols+1]!=-1))
374 if (a[i*cols+sort_col]>a[(i+1)*cols+sort_col])
376 for (int32_t j=0; j<cols; j++)
391 for (int32_t i=0; i<N-1; i++)
396 swap(idx[i],idx[i+1]) ;
405 char* seq1,
char* seq2, int32_t l1, int32_t l2,
float64_t gapCost)
415 for( i1 = 0; i1 < l1; ++i1 ) {
416 gapCosts1[ i1 ] = gapCost * i1;
419 for( i2 = 0; i2 < l2; ++i2 ) {
420 gapCosts2[ i2 ] = gapCost * i2;
421 costs2_1[ i2+1 ] = costs2_1[ i2 ] + gapCosts2[ i2 ];
424 for( i1 = 0; i1 < l1; ++i1 ) {
425 swap( costs2_0, costs2_1 );
426 actCost = costs2_0[ 0 ] + gapCosts1[ i1 ];
427 costs2_1[ 0 ] = actCost;
428 for( i2 = 0; i2 < l2; ++i2 ) {
429 const float64_t actMatch = costs2_0[ i2 ] + ( seq1[i1] == seq2[i2] );
430 const float64_t actGap1 = costs2_0[ i2+1 ] + gapCosts1[ i1 ];
431 const float64_t actGap2 = actCost + gapCosts2[ i2 ];
433 actCost =
min( actMatch, actGap );
434 costs2_1[ i2+1 ] = actCost;
451 for (int32_t i=0; i<len; i++)
452 for (int32_t j=0; j<len; j++)
453 e+=
exp(p2[j*len+i])*(p2[j*len+i]-p1[i]-p1[j]);
462 for (int32_t i=0; i<len; i++)
463 e+=
exp(p[i])*(p[i]-q[i]);
472 for (int32_t i=0; i<len; i++)
503 int32_t lsize=
CMath::min((int32_t) m, (int32_t) n);
508 wrap_dgesvd(jobu, jobvt, m, n, matrix, lda, s, u, ldu, vt, ldvt, &info);
511 for (int32_t i=0; i<n; i++)
513 for (int32_t j=0; j<lsize; j++)
514 vt[i*n+j]=vt[i*n+j]/s[j];
517 cblas_dgemm(CblasColMajor, CblasTrans, CblasTrans, m, n, m, 1.0, vt, ldvt, u, ldu, 0, target, m);