#include "matrix.h"
#include "utils.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <math.h>
void free_matrix(matrix m)
{
int i;
for(i = 0; i < m.rows; ++i) free(m.vals[i]);
free(m.vals);
}
float matrix_topk_accuracy(matrix truth, matrix guess, int k)
{
int* indexes = (int*)xcalloc(k, sizeof(int));
int n = truth.cols;
int i,j;
int correct = 0;
for(i = 0; i < truth.rows; ++i){
top_k(guess.vals[i], n, k, indexes);
for(j = 0; j < k; ++j){
int class_id = indexes[j];
if(truth.vals[i][class_id]){
++correct;
break;
}
}
}
free(indexes);
return (float)correct/truth.rows;
}
void scale_matrix(matrix m, float scale)
{
int i,j;
for(i = 0; i < m.rows; ++i){
for(j = 0; j < m.cols; ++j){
m.vals[i][j] *= scale;
}
}
}
matrix resize_matrix(matrix m, int size)
{
int i;
if (m.rows == size) return m;
if (m.rows < size) {
m.vals = (float**)xrealloc(m.vals, size * sizeof(float*));
for (i = m.rows; i < size; ++i) {
m.vals[i] = (float*)xcalloc(m.cols, sizeof(float));
}
} else if (m.rows > size) {
for (i = size; i < m.rows; ++i) {
free(m.vals[i]);
}
m.vals = (float**)xrealloc(m.vals, size * sizeof(float*));
}
m.rows = size;
return m;
}
void matrix_add_matrix(matrix from, matrix to)
{
assert(from.rows == to.rows && from.cols == to.cols);
int i,j;
for(i = 0; i < from.rows; ++i){
for(j = 0; j < from.cols; ++j){
to.vals[i][j] += from.vals[i][j];
}
}
}
matrix make_matrix(int rows, int cols)
{
int i;
matrix m;
m.rows = rows;
m.cols = cols;
m.vals = (float**)xcalloc(m.rows, sizeof(float*));
for(i = 0; i < m.rows; ++i){
m.vals[i] = (float*)xcalloc(m.cols, sizeof(float));
}
return m;
}
matrix hold_out_matrix(matrix *m, int n)
{
int i;
matrix h;
h.rows = n;
h.cols = m->cols;
h.vals = (float**)xcalloc(h.rows, sizeof(float*));
for(i = 0; i < n; ++i){
int index = rand()%m->rows;
h.vals[i] = m->vals[index];
m->vals[index] = m->vals[--(m->rows)];
}
return h;
}
float *pop_column(matrix *m, int c)
{
float* col = (float*)xcalloc(m->rows, sizeof(float));
int i, j;
for(i = 0; i < m->rows; ++i){
col[i] = m->vals[i][c];
for(j = c; j < m->cols-1; ++j){
m->vals[i][j] = m->vals[i][j+1];
}
}
--m->cols;
return col;
}
matrix csv_to_matrix(char *filename)
{
FILE *fp = fopen(filename, "r");
if(!fp) file_error(filename);
matrix m;
m.cols = -1;
char *line;
int n = 0;
int size = 1024;
m.vals = (float**)xcalloc(size, sizeof(float*));
while((line = fgetl(fp))){
if(m.cols == -1) m.cols = count_fields(line);
if(n == size){
size *= 2;
m.vals = (float**)xrealloc(m.vals, size * sizeof(float*));
}
m.vals[n] = parse_fields(line, m.cols);
free(line);
++n;
}
m.vals = (float**)xrealloc(m.vals, n * sizeof(float*));
m.rows = n;
return m;
}
void matrix_to_csv(matrix m)
{
int i, j;
for(i = 0; i < m.rows; ++i){
for(j = 0; j < m.cols; ++j){
if(j > 0) printf(",");
printf("%.17g", m.vals[i][j]);
}
printf("\n");
}
}
void print_matrix(matrix m)
{
int i, j;
printf("%d X %d Matrix:\n",m.rows, m.cols);
printf(" __");
for(j = 0; j < 16*m.cols-1; ++j) printf(" ");
printf("__ \n");
printf("| ");
for(j = 0; j < 16*m.cols-1; ++j) printf(" ");
printf(" |\n");
for(i = 0; i < m.rows; ++i){
printf("| ");
for(j = 0; j < m.cols; ++j){
printf("%15.7f ", m.vals[i][j]);
}
printf(" |\n");
}
printf("|__");
for(j = 0; j < 16*m.cols-1; ++j) printf(" ");
printf("__|\n");
}
matrix make_matrix(int rows, int cols);
void copy(float *x, float *y, int n);
float dist(float *x, float *y, int n);
int *sample(int n);
int closest_center(float *datum, matrix centers)
{
int j;
int best = 0;
float best_dist = dist(datum, centers.vals[best], centers.cols);
for (j = 0; j < centers.rows; ++j) {
float new_dist = dist(datum, centers.vals[j], centers.cols);
if (new_dist < best_dist) {
best_dist = new_dist;
best = j;
}
}
return best;
}
float dist_to_closest_center(float *datum, matrix centers)
{
int ci = closest_center(datum, centers);
return dist(datum, centers.vals[ci], centers.cols);
}
int kmeans_expectation(matrix data, int *assignments, matrix centers)
{
int i;
int converged = 1;
for (i = 0; i < data.rows; ++i) {
int closest = closest_center(data.vals[i], centers);
if (closest != assignments[i]) converged = 0;
assignments[i] = closest;
}
return converged;
}
void kmeans_maximization(matrix data, int *assignments, matrix centers)
{
matrix old_centers = make_matrix(centers.rows, centers.cols);
int i, j;
int *counts = (int*)xcalloc(centers.rows, sizeof(int));
for (i = 0; i < centers.rows; ++i) {
for (j = 0; j < centers.cols; ++j) {
old_centers.vals[i][j] = centers.vals[i][j];
centers.vals[i][j] = 0;
}
}
for (i = 0; i < data.rows; ++i) {
++counts[assignments[i]];
for (j = 0; j < data.cols; ++j) {
centers.vals[assignments[i]][j] += data.vals[i][j];
}
}
for (i = 0; i < centers.rows; ++i) {
if (counts[i]) {
for (j = 0; j < centers.cols; ++j) {
centers.vals[i][j] /= counts[i];
}
}
}
for (i = 0; i < centers.rows; ++i) {
for (j = 0; j < centers.cols; ++j) {
if(centers.vals[i][j] == 0) centers.vals[i][j] = old_centers.vals[i][j];
}
}
free(counts);
free_matrix(old_centers);
}
void random_centers(matrix data, matrix centers) {
int i;
int *s = sample(data.rows);
for (i = 0; i < centers.rows; ++i) {
copy(data.vals[s[i]], centers.vals[i], data.cols);
}
free(s);
}
int *sample(int n)
{
int i;
int* s = (int*)xcalloc(n, sizeof(int));
for (i = 0; i < n; ++i) s[i] = i;
for (i = n - 1; i >= 0; --i) {
int swap = s[i];
int index = rand() % (i + 1);
s[i] = s[index];
s[index] = swap;
}
return s;
}
float dist(float *x, float *y, int n)
{
float mw = (x[0] < y[0]) ? x[0] : y[0];
float mh = (x[1] < y[1]) ? x[1] : y[1];
float inter = mw*mh;
float sum = x[0] * x[1] + y[0] * y[1];
float un = sum - inter;
float iou = inter / un;
return 1 - iou;
}
void copy(float *x, float *y, int n)
{
int i;
for (i = 0; i < n; ++i) y[i] = x[i];
}
model do_kmeans(matrix data, int k)
{
matrix centers = make_matrix(k, data.cols);
int* assignments = (int*)xcalloc(data.rows, sizeof(int));
random_centers(data, centers);
if (k == 1) kmeans_maximization(data, assignments, centers);
int i;
for(i = 0; i < 1000 && !kmeans_expectation(data, assignments, centers); ++i) {
kmeans_maximization(data, assignments, centers);
}
printf("\n iterations = %d \n", i);
model m;
m.assignments = assignments;
m.centers = centers;
return m;
}