Aujuuuuuhhhh hey guys

someshit.txt

+#include "RIVtools.h"
+#include <float.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_eigen.h>
+
+#define ROUNDOFF 0.000000001
+typedef struct sparseRIVKernel{
+	int count;
+	int* locations;
+	double* values;
+}sparseDRIV;
+
+typedef struct cluster{
+	char name[100];
+	double* centroid;
+	sparseRIV* primaryMap;
+	sparseRIV* secondaryMap;
+	double** rotation;
+	double* eigenValues;
+}RIVcluster;
+void printVector(double* vec, int size){
+	for(int i=0; i<size; i++){
+		printf("%8.3lf,", vec[i]);
+	}puts("");
+}
+double** allocateDoubleBlock(int x, int y, int spaghetti){
+	double** output = malloc(x*sizeof(double*));
+	if(spaghetti){
+		for(int i=0; i<x; i++){
+			output[i] = calloc(y, sizeof(double));
+		}
+	}else{
+		output[0] = calloc(x*y, sizeof(double));
+		for(int i=1; i<x; i++){
+			output[i] = output[0] + y*i;
+		}
+	}
+	return output;
+}
+
+void addDoubleVec(double* a, double* b, int size){
+	for(int i=0; i<size; i++){
+		a[i] += b[i];
+	}
+}
+
+void addSetToVec(double* out, double** vecSet, int setSize, int dimensionality){
+	for(int i=0; i<setSize; i++){
+		addDoubleVec(out, vecSet[i], dimensionality);
+	}
+}
+
+void copyVectorSet(double** dest, double** src, int x, int y){
+	for(int i=0; i<x; i++){
+		memcpy(dest[i], src[i], y*sizeof(double));
+	}
+}
+
+double** convertSetToDouble(sparseRIV** set, int setSize){
+	double** output = allocateDoubleBlock(setSize, RIVSIZE, 0);
+	for(int i=0; i<setSize; i++){
+		for(int j=0; j<set[i]->count; j++){
+			output[i][set[i]->locations[j]] += set[i]->values[j];
+		}
+	}
+	return output;
+}
+double magnitude(double* vec1, int size){
+	double mag = 0;
+	for(int i=0; i<size; i++){
+		mag += vec1[i]*vec1[i];
+	}
+	return sqrt(mag);
+}
+double dotProduct(double* vec1, double* vec2, int size){
+	double dot = 0;
+	for(int i=0; i<size; i++){
+		dot+= vec1[i]*vec2[i];
+	}
+	return dot;
+}
+void normalizeVec(double* vec, int size, double divisor){
+	if(!divisor){
+		divisor = magnitude(vec, size);
+	}
+	for(int i=0; i<size; i++){
+		vec[i] /= divisor;
+	}
+}
+double doubleCosine(double* a, double* b, int size){
+	double maga = 0;
+	double magb = 0;
+	double dot = 0;
+	for(int i=0; i<size; i++){
+		maga += a[i]*a[i];
+		magb += b[i]*b[i];
+		dot += a[i]*b[i];
+	}
+	return dot/(sqrt(maga)*sqrt(magb));
+}
+
+int tangentMap(double** set, int setSize, double* centroid, double** rotatingSet, int dimensionality){
+	int netted = 0;
+	memset(*rotatingSet, 0, setSize*dimensionality*sizeof(double));
+	for(int i=0; i<setSize; i++){
+		double dot = dotProduct(centroid, set[i], dimensionality);
+		if(dot <= 0){
+			continue;
+		}
+		addDoubleVec(rotatingSet[netted], set[i], dimensionality);
+		normalizeVec(rotatingSet[netted], dimensionality, dot);
+		
+		netted ++;
+	}
+	return netted;
+}
+
+double* iterativelyCenter(double** cluster, int setSize, int dimensionality){
+	double centroid[dimensionality];
+	memset(centroid, 0, dimensionality*sizeof(double));
+	addSetToVec(centroid, cluster, setSize, dimensionality);
+	normalizeVec(centroid, dimensionality, 0);
+	
+	double** mappedSet = allocateDoubleBlock(setSize, dimensionality, 0);
+
+	int netted = tangentMap(cluster, setSize, centroid, mappedSet, dimensionality);
+	double* newCentroid = calloc(dimensionality, sizeof(double));
+	addSetToVec(newCentroid, mappedSet, netted, dimensionality);
+	normalizeVec(newCentroid, dimensionality, 0);
+	while( doubleCosine(centroid, newCentroid, dimensionality) < 0.99){
+		addDoubleVec(centroid, newCentroid, dimensionality);
+		memset(newCentroid, 0, dimensionality*sizeof(double));
+		netted = tangentMap(cluster, setSize, centroid, mappedSet, dimensionality);
+		addSetToVec(newCentroid, mappedSet, netted, dimensionality);
+		normalizeVec(newCentroid, dimensionality, 0);
+	}
+	free(*mappedSet);
+	free(mappedSet);
+	return newCentroid;
+}
+
+void transposeSquare(void* vsquare, int size){
+	long long int** square = (long long int**)vsquare;
+	for(int i=0; i<size; i++){
+		for(int j=0; j<i; j++){
+			square[j][i] = square[j][i]^square[i][j];
+			square[i][j] = square[j][i]^square[i][j];
+			square[j][i] = square[j][i]^square[i][j];
+		}
+	}
+}
+double** generateRotation(double* centroid){
+	double** rotation = allocateDoubleBlock(RIVSIZE, RIVSIZE, 0);
+	double temp[RIVSIZE] = {0};
+	addDoubleVec(temp, centroid, RIVSIZE);
+	addDoubleVec(rotation[0], centroid, RIVSIZE);
+	for(int i=1; i<RIVSIZE; i++){
+		int back = i-1;
+		if(!temp[back]){
+			rotation[i][back] = 1;
+			continue;
+		}
+
+		double SS = dotProduct(temp + i, temp+i, RIVSIZE- i);
+		if(!SS){
+			rotation[i][i] = 1;
+			continue;
+		}
+		double factor = SS*-1;
+		temp[back] = factor/temp[back];
+		addDoubleVec(rotation[i]+back, temp + back, RIVSIZE - back);
+		temp[back] = 0;
+	}
+	for(int i=0; i<RIVSIZE; i++){
+		normalizeVec(rotation[i], RIVSIZE, 0);
+	}
+	return rotation;
+}
+void rhoMapVec(double* dest, double* src, double** rotation, double rho){
+	double xAxis = dotProduct(src, rotation[0], RIVSIZE);
+	for(int i=1; i<RIVSIZE; i++){
+		dest[i-1] = dotProduct(src, rotation[i], RIVSIZE);
+		dest[i-1] *= rho/xAxis;
+	}
+}
+double tanDist(double cos){
+	if(cos <= 0.0001)return __FLT_MAX__;
+	return sqrt(1-cos*cos)/cos;
+}
+
+double** rhoMapSet(double** rotation, double* centroid, double** doubleSet, int setSize, double (*f)(double cos), int* nettedCount){
+	int netted = 0;
+	double** rhoMapped = allocateDoubleBlock(setSize, RIVSIZE - 1, 0);
+	for(int i=0; i<setSize; i++){
+		double cos = doubleCosine(centroid, doubleSet[i], RIVSIZE);
+		if(cos < 0) continue;
+		double rho = f(cos)/tanDist(cos);
+		rhoMapVec(rhoMapped[netted], doubleSet[i], rotation, rho);
+		netted++;
+	}
+	*nettedCount = netted;
+	return rhoMapped;
+}
+
+double secDist(double cos){
+	if(cos <= 0.00001) return __FLT_MAX__;
+	return (1.0/cos) -1; 
+}
+void addVectorToCovariance(double** covariance, double* rhoMapped, int size){
+	for(int i=0; i<size; i++){
+		for(int j=0; j<i+1; j++){
+			covariance[i][j] += rhoMapped[i]*rhoMapped[j];
+		}
+	}
+}
+double** generateCovarianceWithCentroid(double** tanMapped, int setSize, double* centroid, int dim){
+	double** covariance = allocateDoubleBlock(dim, dim, 0);
+	for(int i=0; i<setSize; i++){
+		double tempVec[dim];
+		for(int j=0; j<dim; j++){
+			tempVec[j] = tanMapped[i][j] - centroid[j];
+		}
+		addVectorToCovariance(covariance, tempVec, dim);
+	}
+	for(int i=0; i<dim; i++){
+		normalizeVec(covariance[i], i+1, setSize);
+	}
+	for(int i=0; i<dim; i++){
+		for(int j=0; j<i; j++){
+			covariance[j][i] = covariance[i][j];
+		}
+	}
+	return covariance;
+}
+double** generateCovariance(double** rhoMapped, int setSize){
+	double** covariance = allocateDoubleBlock(RIVSIZE-1, RIVSIZE-1, 0);
+	for(int i=0; i<setSize; i++){
+		addVectorToCovariance(covariance, rhoMapped[i], RIVSIZE-1);
+	}
+	for(int i=0; i<RIVSIZE-1; i++){
+		normalizeVec(covariance[i], i+1, setSize);
+	}
+	for(int i=0; i<RIVSIZE-1; i++){
+		for(int j=0; j<i; j++){
+			covariance[j][i] = covariance[i][j];
+		}
+	}
+
+	return covariance;
+}
+sparseRIV* getMatrixMap(double** covariance){
+	sparseRIV* out = sparseAllocate(RIVSIZE-1);
+	out->count = 0;
+	for(int i=0; i<RIVSIZE-1; i++){
+		if(covariance[i][i] > 0.000000001){
+			out->locations[out->count] = i;
+			out->count++;
+		}
+	}
+	return out;
+}
+double** compressCovariance(double** covariance, sparseRIV* map){
+	double** out = allocateDoubleBlock(map->count, map->count, 0);
+	for(int i=0; i<map->count; i++){
+		for(int j=0; j < map->count; j++){
+			out[i][j] = covariance[map->locations[i]][map->locations[j]];
+		}
+	}
+	
+
+	free(*covariance);
+	free(covariance);
+	return out;
+}
+void rotateVector(double* vec, double** rotation, double* out, int dim){
+	for(int i=0; i<dim; i++){
+		out[i] = dotProduct(vec, rotation[i], dim);
+	}
+}
+double** rotateVectors(double** set, double** rotation, int count, int dim){
+	double** rotated = allocateDoubleBlock(count, dim, 0);
+	for(int i=0; i<count; i++){
+		rotateVector(set[i], rotation, rotated[i], dim);
+	}
+	return rotated;
+}
+#define CULLBYMAP(vec, map) for(int i=0,j=0;i<map->count;i++)vec[j++]=vec[map->locations[i]];
+
+sparseRIV* primaryCull(double** set, int setSize, int threshold){
+	int counts[RIVSIZE] = {0};
+	for(int i=0; i<setSize; i++){
+		for(int j=0; j<RIVSIZE; j++){
+			if(set[i][j]) counts[j]++;
+		}
+	}
+	for(int i=0; i<RIVSIZE; i++){
+		if(counts[i] < threshold) counts[i] = 0;
+	}
+	sparseRIV* map = consolidateD2S(counts);
+	for(int i=0; i<setSize; i++){
+		CULLBYMAP(set[i], map);
+	}
+	return map;
+}
+sparseRIV* secondaryCull(double* eval, int length, double threshold){
+	int counts[RIVSIZE] = {0};
+	for(int i=0; i<length; i++){
+		if(eval[i] > threshold)counts[i] = 1;
+	}
+	sparseRIV* map = consolidateD2S(counts);
+	CULLBYMAP(eval, map);
+	
+	return map;
+}
+RIVcluster gaussianCluster(sparseRIV** cluster, int setSize){
+	puts("converting to double");
+	double** doubleSet = convertSetToDouble(cluster, setSize);
+	sparseRIV* primaryMap = primaryCull(doubleSet, setSize, 15);
+	printf("culled: %d\n", primaryMap->count);
+	puts("finding Center");
+	int dimension = primaryMap->count;
+	double* centroid = iterativelyCenter(doubleSet, setSize, dimension);
+	printf("mag: %lf\n", dotProduct(centroid, centroid, dimension));
+	
+	double** tanMapped = allocateDoubleBlock(setSize, RIVSIZE, 0);
+	int nettedCount = tangentMap(doubleSet, setSize, centroid, tanMapped, dimension);
+	puts("making covariance");
+	double** covariance = generateCovarianceWithCentroid(tanMapped, nettedCount, centroid, dimension);
+	
+	
+	gsl_matrix_view matrix = gsl_matrix_view_array(*covariance, dimension, dimension);
+	gsl_vector *eval = gsl_vector_alloc(dimension);
+	gsl_matrix *evec = gsl_matrix_alloc(dimension, dimension);
+
+	gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(dimension);
+	gsl_eigen_symmv(&matrix.matrix, eval, evec, w);
+	// gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_DESC); 
+	printf("mag: %lf\n", dotProduct(centroid, centroid, dimension));
+	double** rotation = malloc(dimension*sizeof(double*));
+	for(int i=0; i<dimension; i++){
+		rotation[i] = &(evec->data[i*dimension]);
+	}
+	double* eigenValues = eval->data;
+	printVector(eigenValues, dimension);
+	sparseRIV* secondaryMap = secondaryCull(eigenValues, dimension, 0.0001);
+	transposeSquare(rotation, dimension);
+	double rotated[dimension];
+	rotateVector(centroid, rotation, rotated, dimension);
+	double** a = rotateVectors(tanMapped, rotation, nettedCount, dimension);
+	covariance = generateCovarianceWithCentroid(a, nettedCount, rotated, dimension);
+	for(int i=0; i<dimension; i++){
+		printVector(covariance[i], dimension);
+	}
+	printf("culled: %d \n", secondaryMap->count);
+	RIVcluster out;
+	centroid = calloc(dimension, sizeof(double));
+	for(int i=0; i<dimension; i++){
+		addDoubleVec(centroid, doubleSet[i] , dimension);
+	}
+	out.centroid = centroid;
+	out.rotation = rotation;
+	out.primaryMap = primaryMap;
+	out.secondaryMap = secondaryMap;
+	out.eigenValues = eigenValues;
+	printVector(eigenValues, secondaryMap->count);
+	printVector(rotated, dimension);
+	free(*covariance);
+	free(covariance);
+	free(*tanMapped);
+	free(tanMapped);
+	
+	return out;
+}
+int tanMapVector(double* vec, double* centroid, int size){
+	double dot = dotProduct(vec, centroid, size);
+	printf("dot:%lf\n", dot);
+	if(dot <= 0){
+		return 0;
+	}
+	normalizeVec(vec, size, dot);
+	// for(int i=0; i<size; i++){
+	// 	vec[i] -= centroid[i];
+	// }
+	return 1;
+}
+double distFromCluster(RIVcluster cluster, sparseRIV* vector){
+	//puts("");
+	double doublevec[RIVSIZE] = {0};
+	for(int i=0; i<vector->count; i++){
+		doublevec[vector->locations[i]] = vector->values[i];
+	}
+	
+	CULLBYMAP(doublevec, cluster.primaryMap);
+	return 1-doubleCosine(cluster.centroid, doublevec, cluster.primaryMap->count);
+	double rotatedCenter[RIVSIZE] = {0};
+	double cos1 = doubleCosine(cluster.centroid, doublevec, cluster.primaryMap->count);
+	rotateVector(cluster.centroid, cluster.rotation, rotatedCenter, cluster.primaryMap->count);
+	
+	if(!tanMapVector(doublevec, cluster.centroid, cluster.primaryMap->count)) return DBL_MAX;
+
+
+	double rho = (1-cos1)/tanDist(cos1);
+	double rotated[RIVSIZE];
+	rotateVector(doublevec, cluster.rotation, rotated, cluster.primaryMap->count);
+	double cos2 = doubleCosine(rotatedCenter, rotated, cluster.primaryMap->count);
+	printf("%lf, %lf\n", cos1, cos2);
+	printVector(rotated, cluster.primaryMap->count);
+	CULLBYMAP(rotated, cluster.secondaryMap);
+	double out = 0;
+	for(int i=0; i<cluster.secondaryMap->count; i++){
+		out += rho*rho*rotated[i]*rotated[i];
+	}
+	printf("out: %lf\n", sqrt(out));
+	return sqrt(out);
+}
\ No newline at end of file