Skip to content

Commit

Permalink
added Jacobi method for eigenvectors / eigenvalues
Browse files Browse the repository at this point in the history
  • Loading branch information
@liuliu
liuliu committed Feb 10, 2014
1 parent 29c161b commit e84d5a7
Show file tree
Hide file tree
Showing 8 changed files with 344 additions and 113 deletions.
6 changes: 4 additions & 2 deletions bin/cwc-bench-runtime.cu
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,12 @@ extern "C" {
extern "C" void cwc_bench_runtime(ccv_convnet_t* convnet, ccv_array_t* categorizeds, ccv_convnet_train_param_t params)
{
int batch = params.mini_batch;
int i, x, y, k, c;
_cwc_convnet_alloc_reserved(convnet, batch, params.layer_params);
cwc_convnet_context_t* context = GPU(convnet)->contexts;
_cwc_convnet_batch_formation(0, categorizeds, 0, ccv_size(251, 251), convnet->rows, convnet->cols, convnet->channels, batch, 0, batch, context->host.input, context->host.c);
for (i = 0; i < convnet->rows * convnet->cols * convnet->channels; i++)
convnet->mean_activity->data.f32[i] = 128;
_cwc_convnet_batch_formation(0, categorizeds, convnet->mean_activity, 0, ccv_size(251, 251), convnet->rows, convnet->cols, convnet->channels, 0, batch, 0, batch, context->host.input, context->host.c);
cudaMemcpy(context->device.input, context->host.input, sizeof(float) * convnet->rows * convnet->cols * convnet->channels * batch, cudaMemcpyHostToDevice);

ccv_convnet_t* update_params = _ccv_convnet_update_new(convnet);
Expand Down Expand Up @@ -85,7 +88,6 @@ extern "C" void cwc_bench_runtime(ccv_convnet_t* convnet, ccv_array_t* categoriz
cudaMemcpy(first_grad, first_gpu_configuration->w, sizeof(float) * first_gpu_layer->wnum, cudaMemcpyDeviceToHost);
printf("finished backward propagate first convolutional layer on GPU\n");

int i, x, y, k, c;
for (i = 0; i < batch; i++)
{
printf("doing batch %d of %d\n", i + 1, batch);
Expand Down
62 changes: 2 additions & 60 deletions bin/image-net.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,8 @@ int main(int argc, char** argv)
.max_epoch = 100,
.mini_batch = 256,
.iterations = 5000,
.symmetric = 1,
.color_gain = 0.1,
};
int i, c;
while (getopt_long_only(argc, argv, "", image_net_options, &c) != -1)
Expand Down Expand Up @@ -399,66 +401,6 @@ int main(int argc, char** argv)
layer_params[11].dor = 0.5;
train_params.layer_params = layer_params;
train_params.size = ccv_size(257, 257);
/*
ccv_size_t size = ccv_size(257, 257);
for (i = 0; i < categorizeds->rnum; i++)
{
ccv_categorized_t* categorized = (ccv_categorized_t*)ccv_array_get(categorizeds, i);
ccv_dense_matrix_t* image = 0;
ccv_read(categorized->file.filename, &image, CCV_IO_ANY_FILE | CCV_IO_RGB_COLOR);
if (image)
{
ccv_dense_matrix_t* norm = 0;
if (image->rows > size.height && image->cols > size.width)
ccv_resample(image, &norm, 0, ccv_max(size.height, (int)(image->rows * (float)size.height / image->cols + 0.5)), ccv_max(size.width, (int)(image->cols * (float)size.width / image->rows + 0.5)), CCV_INTER_AREA);
else if (image->rows < size.height || image->cols < size.width)
ccv_resample(image, &norm, 0, ccv_max(size.height, (int)(image->rows * (float)size.height / image->cols + 0.5)), ccv_max(size.width, (int)(image->cols * (float)size.width / image->rows + 0.5)), CCV_INTER_CUBIC);
else
norm = image;
if (norm != image)
ccv_matrix_free(image);
char filename[1024];
snprintf(filename, 1024, "%s.resize.png", categorized->file.filename);
ccv_write(norm, filename, 0, CCV_IO_PNG_FILE, 0);
ccv_dense_matrix_t* patch = 0;
int x = (norm->cols - size.width) / 2;
int y = (norm->rows - size.height) / 2;
ccv_slice(norm, (ccv_matrix_t**)&patch, CCV_64F, y, x, size.width, size.height);
ccv_matrix_free(norm);
ccv_matrix_free(patch);
FLUSH("done %s, %d / %d", filename, i + 1, categorizeds->rnum);
} else {
printf("\ncannot handle %s\n", categorized->file.filename);
}
}
printf("\n");
for (i = 0; i < tests->rnum; i++)
{
ccv_categorized_t* categorized = (ccv_categorized_t*)ccv_array_get(tests, i);
ccv_dense_matrix_t* image = 0;
ccv_read(categorized->file.filename, &image, CCV_IO_ANY_FILE | CCV_IO_RGB_COLOR);
if (image)
{
ccv_dense_matrix_t* norm = 0;
if (image->rows > size.height && image->cols > size.width)
ccv_resample(image, &norm, 0, ccv_max(size.height, (int)(image->rows * (float)size.height / image->cols + 0.5)), ccv_max(size.width, (int)(image->cols * (float)size.width / image->rows + 0.5)), CCV_INTER_AREA);
else if (image->rows < size.height || image->cols < size.width)
ccv_resample(image, &norm, 0, ccv_max(size.height, (int)(image->rows * (float)size.height / image->cols + 0.5)), ccv_max(size.width, (int)(image->cols * (float)size.width / image->rows + 0.5)), CCV_INTER_CUBIC);
else
norm = image;
if (norm != image)
ccv_matrix_free(image);
char filename[1024];
snprintf(filename, 1024, "%s.resize.png", categorized->file.filename);
ccv_write(norm, filename, 0, CCV_IO_PNG_FILE, 0);
ccv_matrix_free(norm);
FLUSH("done %s, %d / %d", filename, i + 1, tests->rnum);
} else {
printf("\ncannot handle %s\n", categorized->file.filename);
}
}
printf("\n");
*/
ccv_convnet_supervised_train(convnet, categorizeds, tests, working_dir, train_params);
ccv_convnet_free(convnet);
ccv_disable_cache();
Expand Down
5 changes: 4 additions & 1 deletion lib/ccv.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -524,7 +524,7 @@ void ccv_contour_free(ccv_contour_t* contour);

void ccv_invert(ccv_matrix_t* a, ccv_matrix_t** b, int type);
void ccv_solve(ccv_matrix_t* a, ccv_matrix_t* b, ccv_matrix_t** d, int type);
void ccv_eigen(ccv_matrix_t* a, ccv_matrix_t* b, ccv_matrix_t** d, int type);
void ccv_eigen(ccv_dense_matrix_t* a, ccv_dense_matrix_t** vector, ccv_dense_matrix_t** lambda, int type, double epsilon);

typedef struct {
double interp;
Expand Down Expand Up @@ -1156,6 +1156,7 @@ typedef struct {
int channels;
// count and layer of the convnet
int count;
ccv_dense_matrix_t* mean_activity; // mean activity to subtract from
ccv_convnet_layer_t* layers; // the layer configuration
// these can be reused and we don't need to reallocate memory
ccv_dense_matrix_t** denoms; // denominators
Expand Down Expand Up @@ -1184,7 +1185,9 @@ typedef struct {
int max_epoch;
int mini_batch;
int iterations;
int symmetric;
ccv_size_t size;
float color_gain; // the gaussian value for color variations
ccv_convnet_layer_train_param_t* layer_params;
} ccv_convnet_train_param_t;

Expand Down
8 changes: 4 additions & 4 deletions lib/ccv_algebra.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -218,7 +218,7 @@ void ccv_add(ccv_matrix_t* a, ccv_matrix_t* b, ccv_matrix_t** c, int type)
{
ccv_dense_matrix_t* da = ccv_get_dense_matrix(a);
ccv_dense_matrix_t* db = ccv_get_dense_matrix(b);
assert(da->rows == db->rows && da->cols == db->cols && CCV_GET_DATA_TYPE(da->type) == CCV_GET_DATA_TYPE(db->type) && CCV_GET_CHANNEL(da->type) == CCV_GET_CHANNEL(db->type));
assert(da->rows == db->rows && da->cols == db->cols && CCV_GET_CHANNEL(da->type) == CCV_GET_CHANNEL(db->type));
ccv_declare_derived_signature(sig, da->sig != 0 && db->sig != 0, ccv_sign_with_literal("ccv_add"), da->sig, db->sig, CCV_EOF_SIGN);
int no_8u_type = (da->type & CCV_8U) ? CCV_32S : da->type;
type = (type == 0) ? CCV_GET_DATA_TYPE(no_8u_type) | CCV_GET_CHANNEL(da->type) : CCV_GET_DATA_TYPE(type) | CCV_GET_CHANNEL(da->type);
Expand All @@ -228,16 +228,16 @@ void ccv_add(ccv_matrix_t* a, ccv_matrix_t* b, ccv_matrix_t** c, int type)
unsigned char* aptr = da->data.u8;
unsigned char* bptr = db->data.u8;
unsigned char* cptr = dc->data.u8;
#define for_block(_for_get, _for_set) \
#define for_block(_for_get_a, _for_get_b, _for_set) \
for (i = 0; i < da->rows; i++) \
{ \
for (j = 0; j < da->cols * ch; j++) \
_for_set(cptr, j, _for_get(aptr, j, 0) + _for_get(bptr, j, 0), 0); \
_for_set(cptr, j, _for_get_a(aptr, j, 0) + _for_get_b(bptr, j, 0), 0); \
aptr += da->step; \
bptr += db->step; \
cptr += dc->step; \
}
ccv_matrix_getter(da->type, ccv_matrix_setter, dc->type, for_block);
ccv_matrix_getter_a(da->type, ccv_matrix_getter_b, db->type, ccv_matrix_setter, dc->type, for_block);
#undef for_block
}

Expand Down
39 changes: 38 additions & 1 deletion lib/ccv_convnet.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,8 @@ ccv_convnet_t* ccv_convnet_new(int use_cwc_accel, ccv_convnet_layer_param_t para
convnet->rows = params[0].input.matrix.rows;
convnet->cols = params[0].input.matrix.cols;
convnet->channels = params[0].input.matrix.channels;
convnet->mean_activity = ccv_dense_matrix_new(convnet->rows, convnet->cols, convnet->channels | CCV_32F, 0, 0);
ccv_zero(convnet->mean_activity);
ccv_convnet_layer_t* layers = convnet->layers;
int i, j;
for (i = 0; i < count; i++)
Expand Down Expand Up @@ -917,6 +919,7 @@ static ccv_convnet_t* _ccv_convnet_update_new(ccv_convnet_t* convnet)
update_params->cols = convnet->cols;
update_params->count = convnet->count;
update_params->channels = convnet->channels;
update_params->mean_activity = 0;
int i;
for (i = 0; i < convnet->count; i++)
{
Expand Down Expand Up @@ -1065,6 +1068,8 @@ void ccv_convnet_write(ccv_convnet_t* convnet, const char* filename, ccv_convnet
"input_matrix_rows INTEGER, input_matrix_cols INTEGER, input_matrix_channels INTEGER, input_node_count INTEGER, "
"output_rows INTEGER, output_cols INTEGER, output_channels INTEGER, output_count INTEGER, output_strides INTEGER, output_border INTEGER, "
"output_size INTEGER, output_kappa REAL, output_alpha REAL, output_beta REAL);"
"CREATE TABLE IF NOT EXISTS convnet_params "
"(convnet INTEGER PRIMARY KEY ASC, mean_activity BLOB);"
"CREATE TABLE IF NOT EXISTS layer_data "
"(layer INTEGER PRIMARY KEY ASC, weight BLOB, bias BLOB, half_precision INTEGER);";
assert(SQLITE_OK == sqlite3_exec(db, layer_create_table_qs, 0, 0, 0));
Expand Down Expand Up @@ -1147,8 +1152,24 @@ void ccv_convnet_write(ccv_convnet_t* convnet, const char* filename, ccv_convnet
sqlite3_clear_bindings(layer_data_insert_stmt);
}
}
// insert convnet related params
const char convnet_params_insert_qs[] =
"REPLACE INTO convnet_params "
"(convnet, mean_activity) VALUES (0, $mean_activity);";
sqlite3_stmt* convnet_params_insert_stmt = 0;
assert(SQLITE_OK == sqlite3_prepare_v2(db, convnet_params_insert_qs, sizeof(convnet_params_insert_qs), &convnet_params_insert_stmt, 0));
assert(convnet->mean_activity->rows == convnet->rows);
assert(convnet->mean_activity->cols == convnet->cols);
assert(CCV_GET_CHANNEL(convnet->mean_activity->type) == convnet->channels);
assert(CCV_GET_DATA_TYPE(convnet->mean_activity->type) == CCV_32F);
sqlite3_bind_blob(convnet_params_insert_stmt, 1, convnet->mean_activity->data.f32, sizeof(float) * convnet->rows * convnet->cols * convnet->channels, SQLITE_STATIC);
assert(SQLITE_DONE == sqlite3_step(convnet_params_insert_stmt));
sqlite3_reset(convnet_params_insert_stmt);
sqlite3_clear_bindings(convnet_params_insert_stmt);

sqlite3_finalize(layer_params_insert_stmt);
sqlite3_finalize(layer_data_insert_stmt);
sqlite3_finalize(convnet_params_insert_stmt);
sqlite3_close(db);
}
}
Expand Down Expand Up @@ -1214,7 +1235,7 @@ ccv_convnet_t* ccv_convnet_read(int use_cwc_accel, const char* filename)
"SELECT layer, weight, bias, half_precision FROM layer_data;";
if (SQLITE_OK == sqlite3_prepare_v2(db, layer_data_qs, sizeof(layer_data_qs), &layer_data_stmt, 0))
{
while(sqlite3_step(layer_data_stmt) == SQLITE_ROW)
while (sqlite3_step(layer_data_stmt) == SQLITE_ROW)
{
ccv_convnet_layer_t* layer = convnet->layers + sqlite3_column_int(layer_data_stmt, 0);
int half_precision = sqlite3_column_int(layer_data_stmt, 3);
Expand Down Expand Up @@ -1252,6 +1273,20 @@ ccv_convnet_t* ccv_convnet_read(int use_cwc_accel, const char* filename)
}
sqlite3_finalize(layer_data_stmt);
}
sqlite3_stmt* convnet_params_stmt = 0;
// load convnet params
const char convnet_params_qs[] =
"SELECT mean_activity FROM convnet_params WHERE convnet = 0;";
if (SQLITE_OK == sqlite3_prepare_v2(db, convnet_params_qs, sizeof(convnet_params_qs), &convnet_params_stmt, 0))
{
if (sqlite3_step(convnet_params_stmt) == SQLITE_ROW)
{
int elems = sqlite3_column_bytes(convnet_params_stmt, 0) / sizeof(float);
if (elems == convnet->rows * convnet->cols * convnet->channels)
memcpy(convnet->mean_activity->data.f32, sqlite3_column_blob(convnet_params_stmt, 0), sizeof(float) * elems);
}
sqlite3_finalize(convnet_params_stmt);
}
}
sqlite3_close(db);
return convnet;
Expand All @@ -1266,6 +1301,8 @@ void ccv_convnet_free(ccv_convnet_t* convnet)
for (i = 0; i < convnet->count; i++)
if (convnet->layers[i].w)
ccfree(convnet->layers[i].w);
if (convnet->mean_activity)
ccv_matrix_free(convnet->mean_activity);
ccfree(convnet);
}

Expand Down
89 changes: 88 additions & 1 deletion lib/ccv_numeric.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -16,8 +16,95 @@ void ccv_solve(ccv_matrix_t* a, ccv_matrix_t* b, ccv_matrix_t** d, int type)
{
}

void ccv_eigen(ccv_matrix_t* a, ccv_matrix_t* b, ccv_matrix_t** d, int type)
void ccv_eigen(ccv_dense_matrix_t* a, ccv_dense_matrix_t** vector, ccv_dense_matrix_t** lambda, int type, double epsilon)
{
ccv_declare_derived_signature(vsig, a->sig != 0, ccv_sign_with_literal("ccv_eigen(vector)"), a->sig, CCV_EOF_SIGN);
ccv_declare_derived_signature(lsig, a->sig != 0, ccv_sign_with_literal("ccv_eigen(lambda)"), a->sig, CCV_EOF_SIGN);
assert(CCV_GET_CHANNEL(a->type) == 1);
assert(a->rows == a->cols);
type = (type == 0) ? CCV_GET_DATA_TYPE(a->type) | CCV_C1 : CCV_GET_DATA_TYPE(type) | CCV_C1;
// as of now, this function only support real symmetric matrix
ccv_dense_matrix_t* dvector = *vector = ccv_dense_matrix_renew(*vector, a->rows, a->cols, CCV_32F | CCV_64F | CCV_C1, type, vsig);
ccv_dense_matrix_t* dlambda = *lambda = ccv_dense_matrix_renew(*lambda, 1, a->cols, CCV_32F | CCV_64F | CCV_C1, type, lsig);
assert(CCV_GET_DATA_TYPE(dvector->type) == CCV_GET_DATA_TYPE(dlambda->type));
ccv_object_return_if_cached(, dvector, dlambda);
double* ja = (double*)ccmalloc(sizeof(double) * a->rows * a->cols);
int i;
unsigned char* aptr = a->data.u8;
#define for_block(_, _for_get) \
for (i = 0; i < a->rows * a->cols; i++) \
ja[i] = _for_get(aptr, i, 0);
ccv_matrix_getter(a->type, for_block);
#undef for_block
ccv_zero(dvector);
ccv_zero(dlambda);
unsigned char* dvptr = dvector->data.u8;
#define for_block(_, _for_set) \
for (i = 0; i < a->rows; i++) \
_for_set(dvptr + dvector->step * i, i, 1, 0);
ccv_matrix_setter(dvector->type, for_block);
#undef for_block
double accuracy = 0;
for (i = 0; i < a->rows * a->cols; i++) \
accuracy += ja[i];
accuracy = sqrt(2 * accuracy);
int p, q;
unsigned char* dlptr = dlambda->data.u8;
int flag = 1;
#define for_block(_, _for_set, _for_get) \
do { \
if (!flag) \
accuracy = accuracy * 0.5; \
flag = 0; \
for (p = 0; p < a->rows; p++) \
{ \
for (q = p + 1; q < a->cols; q++) \
if (fabs(ja[p * a->cols + q]) > accuracy) \
{ \
double x = -ja[p * a->cols + q]; \
double y = (ja[q * a->cols + q] - ja[p * a->cols + p]) * 0.5; \
double omega = (x == 0 && y == 0) ? 1 : x / sqrt(x * x + y * y); \
if (y < 0) \
omega = -omega; \
double sn = 1.0 + sqrt(1.0 - omega * omega); \
sn = omega / sqrt(2 * sn); \
double cn = sqrt(1.0 - sn * sn); \
double fm = ja[p * a->cols + p]; \
ja[p * a->cols + p] = fm * cn * cn + ja[q * a->cols + q] * sn * sn + ja[p * a->cols + q] * omega; \
ja[q * a->cols + q] = fm * sn * sn + ja[q * a->cols + q] * cn * cn - ja[p * a->cols + q] * omega; \
ja[p * a->cols + q] = ja[q * a->cols + p] = 0; \
for (i = 0; i < a->cols; i++) \
if (i != q && i != p) \
{ \
fm = ja[p * a->cols + i]; \
ja[p * a->cols + i] = fm * cn + ja[q * a->cols + i] * sn; \
ja[q * a->cols + i] = -fm * sn + ja[q * a->cols + i] * cn; \
} \
for (i = 0; i < a->rows; i++) \
if (i != q && i != p) \
{ \
fm = ja[i * a->cols + p]; \
ja[i * a->cols + p] = fm * cn + ja[i * a->cols + q] * sn; \
ja[i * a->cols + q] = -fm * sn + ja[i * a->cols + q] * cn; \
} \
for (i = 0; i < a->cols; i++) \
{ \
fm = _for_get(dvptr + p * dvector->step, i, 0); \
_for_set(dvptr + p * dvector->step, i, fm * cn + _for_get(dvptr + q * dvector->step, i, 0) * sn, 0); \
_for_set(dvptr + q * dvector->step, i, -fm * sn + _for_get(dvptr + q * dvector->step, i, 0) * cn, 0); \
} \
for (i = 0; i < a->cols; i++) \
_for_set(dlptr, i, ja[i * a->cols + i], 0); \
flag = 1; \
break; \
} \
if (flag) \
break; \
} \
} while (accuracy > epsilon);
ccv_matrix_setter_getter(dvector->type, for_block);
#undef for_block
ccfree(ja);
}

void ccv_minimize(ccv_dense_matrix_t* x, int length, double red, ccv_minimize_f func, ccv_minimize_param_t params, void* data)
Expand Down
Loading

0 comments on commit e84d5a7

Please sign in to comment.