two new utilities for encoding and computing codelengths

parent babc6339
......@@ -631,12 +631,6 @@ binary_matrix& mul(const binary_matrix& A, const bool At, const binary_matrix& B
//====================================================================
static idx_t grid_width = 10;
void set_grid_width(idx_t g) {
grid_width = g;
}
// slow: I don't think anyone cares about fast dumping, since most of the time
// will be I/O anyway.
std::ostream& operator<<(std::ostream& out, const binary_matrix& A) {
......
......@@ -331,6 +331,4 @@ private:
block_t trail_mask; // mask trailing bits
};
void set_grid_width(idx_t g);
#endif
......@@ -6,7 +6,42 @@
],
"settings": {
"files.associations": {
"iostream": "cpp"
"iostream": "cpp",
"array": "cpp",
"*.tcc": "cpp",
"bitset": "cpp",
"cctype": "cpp",
"clocale": "cpp",
"cmath": "cpp",
"cstdint": "cpp",
"cstdio": "cpp",
"cstdlib": "cpp",
"cstring": "cpp",
"ctime": "cpp",
"cwchar": "cpp",
"cwctype": "cpp",
"deque": "cpp",
"unordered_map": "cpp",
"vector": "cpp",
"exception": "cpp",
"fstream": "cpp",
"initializer_list": "cpp",
"iomanip": "cpp",
"iosfwd": "cpp",
"istream": "cpp",
"limits": "cpp",
"new": "cpp",
"optional": "cpp",
"ostream": "cpp",
"sstream": "cpp",
"stdexcept": "cpp",
"streambuf": "cpp",
"string_view": "cpp",
"system_error": "cpp",
"type_traits": "cpp",
"tuple": "cpp",
"typeinfo": "cpp",
"utility": "cpp"
}
}
}
\ No newline at end of file
......@@ -95,7 +95,6 @@ int main(int argc, char **argv) {
parse_args(argc,argv);
learn_model_setup(mi_algo,es_algo,du_algo,lm_algo,lmi_algo,0);
fX = fopen(Xname,"r");
set_grid_width(W);
if (!fX) return -1;
res = read_pbm_header(fX,rows,cols);
std::cout << "rows=" << rows << " cols=" << cols << std::endl;
......
/**
* \file bmf_codelength_tool.cpp
* \brief Computes the theoretical codelength of the rows of a binary matrix using a two-pass coding*
*/
#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include "binmat.h"
#include "pbm.h"
#include "bmf_util.h"
#include "bmf_config.h"
#include "bmf_coding.h"
const char* iname = NULL;
void parse_args(int argc, char **argv) {
for (int i = 1; i < argc; ++i) {
if (argv[i][0] == '-') {
if (argv[i][1] == 'h') {
show_help(argv[0]);
exit(0);
} else if (argv[i][1] == 'v') {
inc_verbosity();
continue;
} else if (i == (argc-1)) {
std::cerr << "Missing argument for " << argv[i] << std::endl;
exit(-1);
}
const char* val = argv[i+1];
switch (argv[i][1]) {
default:
std::cerr << "Invalid option " << argv[i] << std::endl;
exit(-1);
}
i++;
} else {
if (!strlen(iname)) {
iname = argv[i];
std::cout << "input file name " << iname << std::endl;
} else {
std::cerr << "Too many arguments. " << std::endl;
exit(-1);
}
}
}
}
int main(int argc, char **argv) {
int res;
FILE* fhandle = 0;
idx_t N,K;
binary_matrix A;
parse_args(argc,argv);
//
// indirect paramters
//
if (!iname) {
std::cerr << "Missing input data file." << std::endl;
std::exit(1);
}
//
// input data
//
fhandle = fopen(iname,"r");
if (!fhandle) {
std::cerr << "Error opening data file " << iname << std::endl;
std::exit(1);
}
res = read_pbm_header(fhandle,N,K);
std::cout << "rows=" << N << " cols=" << K << std::endl;
A.allocate(N,K);
res = read_pbm_data(fhandle,A);
if (res !=PBM_OK) {
std::cerr << "Error " << res << " reading data." << std::endl;
A.destroy();
std::exit(1);
}
fclose(fhandle);
//
// compute codelength
//
double* codelengths = new double[N];
coeffs_codelength_colwise_bernoulli(A,codelengths);
//
// bye
//
for (idx_t i = 0; i < N; i++) {
std::cout << codelengths[i] << std::endl;
}
delete[] codelengths;
A.destroy();
std::exit(0);
}
......@@ -66,3 +66,40 @@ codelength model_codelength(const binary_matrix& E,
}
return L;
}
/*
* L(A_{i:}) = \sum_j -log P_j(A_{ij}|\theta_j), \theta_j = <A_{:j}>
*/
void coeffs_codelength_colwise_bernoulli(const binary_matrix& A, double*& codelengths) {
const idx_t N = A.get_rows();
const idx_t K = A.get_cols();
double* theta = new double[K];
double* L1 = new double[K];
double* L0 = new double[K];
memset(theta,0,sizeof(double));
if (!codelengths) {
codelengths = new double[N];
}
//
// compute K-T estimator for each Bernoulli column model
//
for (idx_t j = 0; j < K; j++) {
theta[j] = (A.col_weight(j) + 0.5)/ ((double) N + 0.5);
L1[j] = -log2(theta[j]);
L0[j] = -log2(1.0-theta[j]);
}
//
// compute each codelength
//
for (idx_t i = 0; i < N; i++) {
codelengths[i] = 0;
for (idx_t j = 0; j < K; j++) {
codelengths[i] += A.get(i,j) ? L1[j] : L0[j];
}
}
delete[] theta;
delete[] L1;
delete[] L0;
}
......@@ -42,4 +42,11 @@ codelength model_codelength(const binary_matrix& E,
const binary_matrix& D,
const binary_matrix& A);
/**
* Computes the codelength of a coefficients matrix
* by assuming a different Bernoulli model for each *column*
* and then encoding each row as L(A_i) = \sum_j -log P_j(A_ij)
*/
void coeffs_codelength_colwise_bernoulli(const binary_matrix& A, double*& codelengths);
#endif
......@@ -108,19 +108,34 @@ static int count_options(const char* names[]) {
return n;
}
void coding_setup(int es_algo, int pv_algo) {
int n;
n = count_options(coeffs_update_algorithm_names);
if (es_algo > n) {
std::cerr << "Invalid coefficients update algorithm (0-" << n << ')' << std::endl;
exit(-1);
}
n = count_options(preview_iter_names);
if (pv_algo > n) {
std::cerr << "Invalid preview algorithm (0-" << n << ')' << std::endl;
exit(-1);
}
coefficients_update = coeffs_update_algorithm_catalog[es_algo];
std::cout << "Using " << coeffs_update_algorithm_names[es_algo] << " for coefficients update." << std::endl;
preview_iter = preview_iter_catalog[pv_algo];
std::cout << "Using " << preview_iter_names[pv_algo] << " for previewing iterations." << std::endl;
}
void learn_model_setup(int mi_algo, int es_algo, int du_algo, int lmo_algo, int lmi_algo, int pv_algo) {
int n;
coding_setup(es_algo,pv_algo);
n = count_options(model_init_algorithm_names);
if (mi_algo > n) {
std::cerr << "Invalid model initialization algorithm (0-" << n << ')' << std::endl;
exit(-1);
}
n = count_options(coeffs_update_algorithm_names);
if (es_algo > n) {
std::cerr << "Invalid coefficients update algorithm (0-" << n << ')' << std::endl;
exit(-1);
}
n = count_options(dict_update_algorithm_names);
if (du_algo > n) {
std::cerr << "Invalid dictionary update algorithm (0-" << n << ')' << std::endl;
......@@ -145,22 +160,12 @@ void learn_model_setup(int mi_algo, int es_algo, int du_algo, int lmo_algo, int
initialize_dictionary = initialize_dictionary_dummy;
}
n = count_options(preview_iter_names);
if (pv_algo > n) {
std::cerr << "Invalid preview algorithm (0-" << n << ')' << std::endl;
exit(-1);
}
coefficients_update = coeffs_update_algorithm_catalog[es_algo];
std::cout << "Using " << coeffs_update_algorithm_names[es_algo] << " for coefficients update." << std::endl;
update_dictionary = dict_update_algorithm_catalog[du_algo];
std::cout << "Using " << dict_update_algorithm_names[du_algo] << " for dictionary update." << std::endl;
learn_model_outer = learn_model_algorithm_catalog[lmo_algo];
std::cout << "Using " << model_learn_algorithm_names[lmo_algo] << " for outer learning loop." << std::endl;
learn_model_inner = learn_model_algorithm_catalog[lmi_algo];
std::cout << "Using " << model_learn_algorithm_names[lmi_algo] << " for inner learning." << std::endl;
preview_iter = preview_iter_catalog[pv_algo];
std::cout << "Using " << preview_iter_names[pv_algo] << " for previewing iterations." << std::endl;
}
void set_max_err_weight(size_t _me) {
......
......@@ -45,6 +45,8 @@ extern model_learn_algorithm_t learn_model_inner;
extern preview_iter_t preview_iter;
void coding_setup(int es_algo, int pv_algo);
void learn_model_setup(int mi_algo, int es_algo, int du_algo, int lm_algo, int lmi_algo, int pv_algo);
void list_choices(const char* prefix, const char* opts[]);
......
/**
* \file bmf_learn_tool.cpp
* \brief Learns or updates a Binary Matrix Factorization model of the form $X = DA + E$ where $X$ is the input matrix, $D$ is a dictionary, $A$ are linear coefficients, and $E$ is a residual.
*/
#include <cstdio>
#include <cstdlib>
#include <iomanip>
#include "binmat.h"
#include "pbm.h"
#include "bmf.h"
#include "bmf_rng.h"
#include "bmf_util.h"
#include "bmf_config.h"
#include "binimage.h"
#include "preview.h"
int es_algo = 2; // default now is BMP
int pv_algo = 0;
idx_t W = 16;
idx_t stride = 0;
idx_t K = 512;
idx_t maxK = 512;
idx_t batch_size = 0;
bool image_mode = false;
bool force_mosaic = true;
bool force_residual_mosaic = true;
const char* dname = 0;
const char* iname = 0;
const char* oname = 0;
const char* mname = 0; // optional
void parse_args(int argc, char **argv) {
for (int i = 1; i < argc; ++i) {
if (argv[i][0] == '-') {
if (argv[i][1] == 'h') {
show_help(argv[0]);
exit(0);
} else if (argv[i][1] == 'v') {
inc_verbosity();
continue;
} else if (i == (argc-1)) {
std::cerr << "Missing argument for " << argv[i] << std::endl;
exit(-1);
}
const char* val = argv[i+1];
switch (argv[i][1]) {
case 'c':
es_algo = atoi(val);
break;
case 'w':
W = (idx_t) atoi(val);
break;
case 's':
stride = (idx_t) atoi(val);
break;
case 'r':
random_seed = atol(val);
break;
case 'I':
image_mode = (atoi(val) > 0);
break;
case 'B':
batch_size = (atoi(val) > 0);
break;
case 'M':
force_residual_mosaic = (atoi(val) > 0);
break;
case 'V':
pv_algo = atoi(val);
break;
default:
std::cerr << "Invalid option " << argv[i] << std::endl;
exit(-1);
}
i++;
} else {
if (!dname) {
dname = argv[i];
std::cout << "dictionary file name " << dname << std::endl;
} else if (!iname) {
iname = argv[i];
std::cout << "input file name " << iname << std::endl;
} else if (!oname) {
mname = argv[i];
std::cout << "output file name " << oname << std::endl;
} else if (!mname) {
mname = argv[i];
std::cout << "mask file name " << mname << std::endl;
} else {
std::cerr << "Too many arguments. " << std::endl;
exit(-1);
}
}
}
}
/**
* KSVD-like binary dictionary learning algorithm applied to
* image patches.
*/
int main(int argc, char **argv) {
idx_t rows,cols;
int res;
FILE* fhandle = 0;
idx_t M,N,K;
binary_matrix X,E,A,D,H;
parse_args(argc,argv);
//
// indirect paramters
//
if (stride == 0) {
stride = W;
}
if (!dname) {
std::cerr << "Error missing dictionary file " << dname << std::endl;
std::exit(1);
}
if (!iname) {
std::cerr << "Missing input data file." << std::endl;
std::exit(1);
}
//
// load dictionary
//
fhandle = fopen(dname,"r");
if (!dname) {
std::cerr << "Error opening dictionary file " << dname << std::endl;
std::exit(2);
}
res = read_pbm_header(fhandle,K,M);
D.allocate(K,M);
res = read_pbm_data(fhandle,D);
if (res !=PBM_OK) {
std::cerr << "Error " << res << " reading dictionary." << std::endl;
D.destroy();
std::exit(2);
}
fclose(fhandle);
//
// input data
//
fhandle = fopen(iname,"r");
if (!fhandle) {
std::cerr << "Error opening data file " << iname << std::endl;
D.destroy();
A.destroy();
std::exit(1);
}
res = read_pbm_header(fhandle,rows,cols);
std::cout << "rows=" << rows << " cols=" << cols << std::endl;
binary_matrix I(rows,cols);
res = read_pbm_data(fhandle,I);
if (res !=PBM_OK) {
std::cerr << "Error " << res << " reading data." << std::endl;
D.destroy();
A.destroy();
I.destroy();
std::exit(1);
}
fclose(fhandle);
//
// optional mask
//
bool have_mask = false;
if (strlen(mname) > 0) {
have_mask = true;
fhandle = fopen(mname,"r");
{
std::cerr << "Error opening mask file " << mname << std::endl;
D.destroy();
I.destroy();
std::exit(2);
}
idx_t rows2, cols2;
res = read_pbm_header(fhandle,rows2,cols2);
if (rows2 != rows) {
std::cerr << "# rows of mask and data file do not coincide. " << rows << " != " << rows2 << std::endl;
D.destroy();
I.destroy();
std::exit(2);
}
if (rows2 != rows) {
std::cerr << "# columns of mask and data file do not coincide. " << cols << " != " << cols2 << std::endl;
D.destroy();
I.destroy();
std::exit(2);
}
}
//
// mask
//
binary_matrix J;
if (have_mask) {
J.allocate(rows,cols);
res = read_pbm_data(fhandle,J);
if (res !=PBM_OK) {
std::cerr << "Error " << res << " reading mask." << std::endl;
D.destroy();
A.destroy();
I.destroy();
J.destroy();
std::exit(3);
}
;
}
if (image_mode) {
std::cout << "==== INPUT DATA TREATED AS IMAGE, VECTORS ARE PATCHES =====\n" << std::endl;
idx_t Ny = compute_grid_size(rows,W,stride,EXTRACT_FULL);
idx_t Nx = compute_grid_size(cols,W,stride,EXTRACT_FULL);
M = W*W;
std::cout << "Nx=" << Nx << " Ny=" << Ny << std::endl;
N = Nx*Ny;
X.allocate(N,M);
H.allocate(N,M);
extract_patches(I,W,stride,EXTRACT_FULL,X);
extract_patches(J,W,stride,EXTRACT_FULL,H);
J.destroy();
I.destroy();
} else {
std::cout << "==== INPUT DATA TREATED AS MATRIX, VECTORS ARE ROWS =====\n" << std::endl;
X = I;
if (have_mask) {
H = J;
}
M = I.get_cols();
N = I.get_rows();
}
//
// set up different aspects of learning algorithm
// the overall algorithm delegates parts to other algorithms.
//
coding_setup(es_algo,pv_algo);
E.allocate(N,M);
A.allocate(N,K);
X.copy_to(E);
A.clear();
//
// learn model
//
set_prefix("./");
idx_t L = coefficients_update(E,H,D,A,K,0);
//mul(A,false,D,false,X);
std::cout << "FINAL: h(E)=" << E.weight()
<< " h(A)=" << A.weight()
<< " h(D)=" << D.weight()
<< " L*(X)=" << L
<< " K*=" << D.get_rows() << std::endl;
//
// write output
//
write_pbm(A,oname);
if (image_mode) {
fhandle = fopen("__residual__.pbm","w");
if (!fhandle) std::exit(4);
integer_matrix accu,count;
I.allocate(rows,cols);
accu.allocate(rows,cols);
count.allocate(rows,cols);
stitch_patches(E,I.get_rows(),I.get_cols(),stride,EXTRACT_FULL,I,accu,count);
count.destroy();
accu.destroy();
write_pbm(I,fhandle);
I.destroy();
fclose(fhandle);
} else {
write_pbm(E,"__residual__.pbm");
}
if (force_residual_mosaic) {
render_mosaic(E,"__residual_mosaic__.pbm");
}
A.destroy();
E.destroy();
D.destroy();
X.destroy();
H.destroy();
std::exit(0);
}
......@@ -260,10 +260,6 @@ int main(int argc, char **argv) {
write_pbm(A,"__initial_coeffs__.pbm");
binary_matrix E(N,M);
//
// setup learning algorithm
//
set_grid_width(W);
//
// learn model
//
set_prefix("./");
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment