return strncmp(s1, s2, strlen(s2)) == 0;
}
-static void init(int argc, char *argv[], real *angle, real *accuracy, int *check_edges_with_same_endpoint, int *seed, char **color_scheme, char **lightness){
+static void init(int argc, char *argv[], double *angle, double *accuracy, int *check_edges_with_same_endpoint, int *seed, char **color_scheme, char **lightness){
char* cmd = argv[0];
outfile = NULL;
}
-static int clarify(Agraph_t* g, real angle, real accuracy, int check_edges_with_same_endpoint, int seed, char *color_scheme, char *lightness){
+static int clarify(Agraph_t* g, double angle, double accuracy, int check_edges_with_same_endpoint, int seed, char *color_scheme, char *lightness){
if (checkG(g)) {
agerr (AGERR, "Graph %s contains loops or multiedges\n", agnameof(g));
int main(int argc, char *argv[])
{
- real accuracy;
- real angle;
+ double accuracy;
+ double angle;
int check_edges_with_same_endpoint, seed;
char *color_scheme = NULL;
char *lightness = NULL;
#include <time.h>
static void get_local_12_norm(int n, int i, const int *ia, const int *ja,
- const int *p, real *norm){
+ const int *p, double *norm){
int j, nz = 0;
norm[0] = n; norm[1] = 0;
for (j = ia[i]; j < ia[i+1]; j++){
}
if (nz > 0) norm[1] /= nz;
}
-static void get_12_norm(int n, int *ia, int *ja, int *p, real *norm){
+static void get_12_norm(int n, int *ia, int *ja, int *p, double *norm){
/* norm[0] := antibandwidth
norm[1] := (\sum_{{i,j}\in E} |p[i] - p[j]|)/|E|
norm[2] := (\sum_{i\in V} (Min_{{j,i}\in E} |p[i] - p[j]|)/|V|
*/
int i, j, nz = 0;
- real tmp;
+ double tmp;
norm[0] = n; norm[1] = 0; norm[2] = 0;
for (i = 0; i < n; i++){
tmp = n;
void improve_antibandwidth_by_swapping(SparseMatrix A, int *p){
bool improved = true;
int cnt = 1, n = A->m, i, j, *ia = A->ia, *ja = A->ja;
- real norm1[3], norm2[3], norm11[3], norm22[3];
- real pi, pj;
+ double norm1[3], norm2[3], norm11[3], norm22[3];
+ double pi, pj;
clock_t start = clock();
FILE *fp = NULL;
SparseMatrix L, A2;
int *ia = A->ia, *ja = A->ja;
int a = -1;
- real nrow;
- real *v = NULL;
- real norm1[3];
+ double nrow;
+ double *v = NULL;
+ double norm1[3];
clock_t start, start2;
start = clock();
/* largest eigen vector */
{
- real eig;
+ double eig;
power_method(L, L->n, 1, seed, &v, &eig);
}
char **infiles;
FILE* outfile;
int dim;
- real shore_depth_tol;
+ double shore_depth_tol;
int nrandom;
int show_points;
- real bbox_margin[2];
+ double bbox_margin[2];
int useClusters;
int clusterMethod;
bool plotedges;
int color_scheme;
- real line_width;
+ double line_width;
char *color_scheme_str;
const char *opacity;
char *plot_label;
- real *bg_color;
+ double *bg_color;
int improve_contiguity_n;
int nart;
bool color_optimize;
{
char* cmd = argv[0];
int c;
- real s;
+ double s;
int v, r;
char stmp[3]; /* two character string plus '\0' */
case 'g': {
gvcolor_t color;
if (colorxlate(optarg, &color, RGBA_DOUBLE) == COLOR_OK) {
- if (!pm->bg_color) pm->bg_color = N_NEW(3,real);
+ if (!pm->bg_color) pm->bg_color = N_NEW(3,double);
pm->bg_color[0] = color.u.RGBA[0];
pm->bg_color[1] = color.u.RGBA[1];
pm->bg_color[2] = color.u.RGBA[2];
}
static void
-makeMap (SparseMatrix graph, int n, real* x, real* width, int* grouping,
+makeMap (SparseMatrix graph, int n, double* x, double* width, int* grouping,
char** labels, float* fsz, float* rgb_r, float* rgb_g, float* rgb_b, params_t* pm, Agraph_t* g )
{
int dim = pm->dim;
int i, flag = 0;
SparseMatrix poly_lines, polys, poly_point_map;
- real edge_bridge_tol = 0.;
+ double edge_bridge_tol = 0.;
int npolys, nverts, *polys_groups, exclude_random;
- real *x_poly, *xcombined;
+ double *x_poly, *xcombined;
SparseMatrix country_graph;
int improve_contiguity_n = pm->improve_contiguity_n;
#ifdef TIME
}
#ifdef TIME
- fprintf(stderr, "map making time = %f\n",((real) (clock() - cpu)) / CLOCKS_PER_SEC);
+ fprintf(stderr, "map making time = %f\n",((double) (clock() - cpu)) / CLOCKS_PER_SEC);
#endif
{
SparseMatrix graph;
int n;
- real* width = NULL;
- real* x;
+ double* width = NULL;
+ double* x;
char** labels = NULL;
int* grouping;
float* rgb_r;
#include <edgepaint/lab.h>
#include <edgepaint/node_distinct_coloring.h>
-void map_palette_optimal_coloring(char *color_scheme, char *lightness, SparseMatrix A0, real accuracy, int seed,
+void map_palette_optimal_coloring(char *color_scheme, char *lightness, SparseMatrix A0, double accuracy, int seed,
float **rgb_r, float **rgb_g, float **rgb_b){
/*
for a graph A, get a distinctive color of its nodes so that the color distanmce among all nodes are maximized. Here
/*color: On input an array of size n*cdim, if NULL, will be allocated. On exit the final color assignment for node i is [cdim*i,cdim*(i+1)), in RGB (between 0 to 1)
*/
- real *colors = NULL;
+ double *colors = NULL;
int n = A0->m, i, cdim;
SparseMatrix A;
int weightedQ = TRUE;
int iter_max = 100;
- {real *dist = NULL;
+ {double *dist = NULL;
A = SparseMatrix_symmetrize(A0, FALSE);
SparseMatrix_distance_matrix(A, 0, &dist);
SparseMatrix_delete(A);
return point_poly_map->ja[point_poly_map->ia[ip]];
}
-void improve_contiguity(int n, int dim, int *grouping, SparseMatrix poly_point_map, real *x, SparseMatrix graph){
+void improve_contiguity(int n, int dim, int *grouping, SparseMatrix poly_point_map, double *x, SparseMatrix graph){
/*
grouping: which group each of the vertex belongs to
poly_point_map: a matrix of dimension npolys x (n + nrandom), poly_point_map[i,j] != 0 if polygon i contains the point j.
. If j < n, it is the original point, otherwise it is artificial point (forming the rectangle around a label) or random points.
*/
int i, j, *ia, *ja, u, v;
- real *a;
+ double *a;
SparseMatrix point_poly_map, D;
- real dist;
+ double dist;
int nbad = 0, flag;
int maxit = 10;
- real tol = 0.001;
+ double tol = 0.001;
D = SparseMatrix_get_real_adjacency_matrix_symmetrized(graph);
assert(graph->m == n);
ia = D->ia; ja = D->ja;
- a = (real*) D->a;
+ a = (double*) D->a;
/* point_poly_map: each row i has only 1 entry at column j, which says that point i is in polygon j */
point_poly_map = SparseMatrix_transpose(poly_point_map);
}
}
- if (Verbose) fprintf(stderr,"ratio (edges among discontiguous regions vs total edges)=%f\n",((real) nbad)/ia[n]);
+ if (Verbose) fprintf(stderr,"ratio (edges among discontiguous regions vs total edges)=%f\n",((double) nbad)/ia[n]);
stress_model(dim, D, D, &x, FALSE, maxit, tol, &flag);
assert(!flag);
struct Triangle {
int vertices[3];/* 3 points */
- real center[2]; /* center of the triangle */
+ double center[2]; /* center of the triangle */
};
-static void normal(real v[], real normal[]){
+static void normal(double v[], double normal[]){
if (v[0] == 0){
normal[0] = 1; normal[1] = 0;
} else {
}
}
-static void triangle_center(real x[], real y[], real z[], real c[]){
+static void triangle_center(double x[], double y[], double z[], double c[]){
/* find the "center" c, which is the intersection of the 3 vectors that are normal to each
of the edges respectively, and which passes through the center of the edges respectively
center[{x_, y_, z_}] := Module[
]
*/
- real xy[2], yz[2], nxy[2], nyz[2], ymx[2], ymz[2], beta, bot;
+ double xy[2], yz[2], nxy[2], nyz[2], ymx[2], ymz[2], beta, bot;
int i;
for (i = 0; i < 2; i++) ymx[i] = y[i] - x[i];
}
}
-static void plot_dot_labels(FILE *f, int n, int dim, real *x, char **labels, float *fsz){
+static void plot_dot_labels(FILE *f, int n, int dim, double *x, char **labels, float *fsz){
int i;
for (i = 0; i < n; i++){
}
-static void dot_polygon(agxbuf *sbuff, int np, float *xp, float *yp, real line_width,
+static void dot_polygon(agxbuf *sbuff, int np, float *xp, float *yp, double line_width,
int fill, const char *cstring){
if (np > 0){
}
}
-static void dot_one_poly(agxbuf *sbuff, real line_width, int fill, int np,
+static void dot_one_poly(agxbuf *sbuff, double line_width, int fill, int np,
float *xp, float *yp, const char *cstring) {
dot_polygon(sbuff, np, xp, yp, line_width, fill, cstring);
}
-static void plot_dot_polygons(agxbuf *sbuff, real line_width,
+static void plot_dot_polygons(agxbuf *sbuff, double line_width,
const char *line_color, SparseMatrix polys,
- real *x_poly, int *polys_groups, float *r,
+ double *x_poly, int *polys_groups, float *r,
float *g, float *b, const char *opacity) {
int i, j, *ia = polys->ia, *ja = polys->ja, *a = (int*) polys->a, npolys = polys->m, nverts = polys->n, ipoly,first;
int np = 0, maxlen = 0;
free(yp);
}
-void plot_dot_map(Agraph_t* gr, int n, int dim, real *x, SparseMatrix polys,
- SparseMatrix poly_lines, real line_width,
- const char *line_color, real *x_poly, int *polys_groups,
+void plot_dot_map(Agraph_t* gr, int n, int dim, double *x, SparseMatrix polys,
+ SparseMatrix poly_lines, double line_width,
+ const char *line_color, double *x_poly, int *polys_groups,
char **labels, float *fsz, float *r, float *g, float *b,
const char* opacity, SparseMatrix A, FILE* f) {
/* if graph object exist, we just modify some attributes, otherwise we dump the whole graph */
agxbfree(&sbuff);
}
-static void get_tri(int n, int dim, real *x, int *nt, struct Triangle **T, SparseMatrix *E) {
+static void get_tri(int n, int dim, double *x, int *nt, struct Triangle **T, SparseMatrix *E) {
/* always contains a self edge and is symmetric.
input:
n: number of points
free(elist);
}
-static void plot_cycle(int head, int *cycle, int *edge_table, real *x){
+static void plot_cycle(int head, int *cycle, int *edge_table, double *x){
int cur, next;
- real x1, x2, y1, y2;
+ double x1, x2, y1, y2;
printf("Graphics[{");
cur = head;
}
static void get_polygon_solids(int nt, SparseMatrix E, int ncomps, int *comps_ptr, int *comps,
- int *mask, real *x_poly, SparseMatrix *polys){
+ int *mask, double *x_poly, SparseMatrix *polys){
/*============================================================
polygon slids that will be colored
}
static void get_polygons(int exclude_random, int n, int nrandom, int dim, SparseMatrix graph, int *grouping,
- int nt, struct Triangle *Tp, SparseMatrix E, int *nverts, real **x_poly,
+ int nt, struct Triangle *Tp, SparseMatrix E, int *nverts, double **x_poly,
int *npolys, SparseMatrix *poly_lines, SparseMatrix *polys, int **polys_groups, SparseMatrix *poly_point_map, SparseMatrix *country_graph){
int i, j;
int *mask;
}
*npolys = ncomps;
- *x_poly = MALLOC(sizeof(real)*dim*nt);
+ *x_poly = MALLOC(sizeof(double)*dim*nt);
for (i = 0; i < nt; i++){
for (j = 0; j < dim; j++){
(*x_poly)[i*dim+j] = Tp[i].center[j];
}
static int make_map_internal(int exclude_random, int include_OK_points,
- int n, int dim, real *x0, int *grouping0, SparseMatrix graph, real bounding_box_margin[], int *nrandom, int nedgep,
- real shore_depth_tol, real **xcombined, int *nverts, real **x_poly,
+ int n, int dim, double *x0, int *grouping0, SparseMatrix graph, double bounding_box_margin[], int *nrandom, int nedgep,
+ double shore_depth_tol, double **xcombined, int *nverts, double **x_poly,
int *npolys, SparseMatrix *poly_lines, SparseMatrix *polys, int **polys_groups, SparseMatrix *poly_point_map,
SparseMatrix *country_graph, int highlight_cluster, int *flag){
- real xmax[2], xmin[2], area, *x = x0;
+ double xmax[2], xmin[2], area, *x = x0;
int i, j;
QuadTree qt;
int dim2 = 2, nn = 0;
int max_qtree_level = 10;
- real ymin[2], min;
+ double ymin[2], min;
int imin, nz, nzok = 0, nzok0 = 0, nt;
- real *xran, point[2];
+ double *xran, point[2];
struct Triangle *Tp;
SparseMatrix E;
- real boxsize[2];
+ double boxsize[2];
int INCLUDE_OK_POINTS = include_OK_points;/* OK points are random points inserted and found to be within shore_depth_tol of real/artificial points,
including them instead of throwing away increase realism of boundary */
int *grouping = grouping0;
*nrandom -= 4;
}
- if (shore_depth_tol < 0) shore_depth_tol = sqrt(area/(real) n); /* set to average distance for random distribution */
+ if (shore_depth_tol < 0) shore_depth_tol = sqrt(area/(double) n); /* set to average distance for random distribution */
if (Verbose) fprintf(stderr,"nrandom=%d shore_depth_tol=%.08f\n",*nrandom, shore_depth_tol);
two connected components be separated due to small shore depth */
{
int nz;
- real *y;
+ double *y;
int k, t, np=nedgep;
if (graph && np){
fprintf(stderr,"add art np = %d\n",np);
nz = graph->nz;
- y = MALLOC(sizeof(real)*(dim*n + dim*nz*np));
+ y = MALLOC(sizeof(double)*(dim*n + dim*nz*np));
for (i = 0; i < n*dim; i++) y[i] = x[i];
grouping = MALLOC(sizeof(int)*(n + nz*np));
for (i = 0; i < n; i++) grouping[i] = grouping0[i];
if (!HIGHLIGHT_SET || (grouping[i] == grouping[graph->ja[j]] && grouping[i] == HIGHLIGHT_SET)){
for (t = 0; t < np; t++){
for (k = 0; k < dim; k++){
- y[nz*dim+k] = t/((real) np)*x[i*dim+k] + (1-t/((real) np))*x[(graph->ja[j])*dim + k];
+ y[nz*dim+k] = t/((double) np)*x[i*dim+k] + (1-t/((double) np))*x[(graph->ja[j])*dim + k];
}
assert(n + (nz-n)*np + t < n + nz*np && n + (nz-n)*np + t >= 0);
- if (t/((real) np) > 0.5){
+ if (t/((double) np) > 0.5){
grouping[nz] = grouping[i];
} else {
grouping[nz] = grouping[graph->ja[j]];
}
}
if (Verbose) {
- real bbm0 = bounding_box_margin[0];
- real bbm1 = bounding_box_margin[1];
+ double bbm0 = bounding_box_margin[0];
+ double bbm1 = bounding_box_margin[1];
if (bbm0 > 0)
fprintf (stderr, "bounding box margin: %.06f", bbm0);
else if (bbm0 < 0)
fprintf (stderr, " %.06f\n", MAX(boxsize[1]*0.2, 2*shore_depth_tol));
}
if (*nrandom < 0) {
- real n1, n2, area2;
+ double n1, n2, area2;
area2 = (xmax[1] - xmin[1])*(xmax[0] - xmin[0]);
n1 = (int) area2/(shore_depth_tol*shore_depth_tol);
n2 = n*((int) area2/area);
*nrandom = MAX(n1, n2);
}
srand(123);
- xran = MALLOC(sizeof(real)*(*nrandom + 4)*dim2);
+ xran = MALLOC(sizeof(double)*(*nrandom + 4)*dim2);
nz = 0;
if (INCLUDE_OK_POINTS){
nzok0 = nzok = *nrandom - 1;/* points that are within tolerance of real or artificial points */
*nrandom = nz;
if (Verbose) fprintf( stderr, "nn nrandom=%d\n",*nrandom);
} else {
- xran = MALLOC(sizeof(real)*4*dim2);
+ xran = MALLOC(sizeof(double)*4*dim2);
}
if (INCLUDE_OK_POINTS){
- *xcombined = MALLOC(sizeof(real)*(nn+*nrandom)*dim2);
+ *xcombined = MALLOC(sizeof(double)*(nn+*nrandom)*dim2);
} else {
- *xcombined = MALLOC(sizeof(real)*(n+*nrandom)*dim2);
+ *xcombined = MALLOC(sizeof(double)*(n+*nrandom)*dim2);
}
for (i = 0; i < n; i++) {
for (j = 0; j < dim2; j++) (*xcombined)[i*dim2+j] = x[i*dim+j];
{
int nz, nh = 0;/* the set to highlight */
- real *xtemp;
+ double *xtemp;
if (HIGHLIGHT_SET){
if (Verbose) fprintf(stderr," highlight cluster %d, n = %d\n",HIGHLIGHT_SET, n);
- xtemp = MALLOC(sizeof(real)*n*dim);
+ xtemp = MALLOC(sizeof(double)*n*dim);
/* shift set to the beginning */
nz = 0;
for (i = 0; i < n; i++){
for (i = nh; i < n; i++){
grouping[i] = 2;
}
- memcpy(*xcombined, xtemp, n*dim*sizeof(real));
+ memcpy(*xcombined, xtemp, n*dim*sizeof(double));
*nrandom += n - nh;/* count everything except cluster HIGHLIGHT_SET as random */
n = nh;
if (Verbose) fprintf(stderr,"nh = %d\n",nh);
return 0;
}
-static void add_point(int *n, int igrp, real **x, int *nmax, real point[], int **groups){
+static void add_point(int *n, int igrp, double **x, int *nmax, double point[], int **groups){
if (*n >= *nmax){
*nmax = MAX((int) 0.2*(*n), 20) + *n;
- *x = REALLOC(*x, sizeof(real)*2*(*nmax));
+ *x = REALLOC(*x, sizeof(double)*2*(*nmax));
*groups = REALLOC(*groups, sizeof(int)*(*nmax));
}
(*n)++;
}
-static void get_boundingbox(int n, int dim, real *x, real *width, real *bbox){
+static void get_boundingbox(int n, int dim, double *x, double *width, double *bbox){
int i;
bbox[0] = bbox[1] = x[0];
bbox[2] = bbox[3] = x[1];
}
int make_map_from_rectangle_groups(int exclude_random, int include_OK_points,
- int n, int dim, real *x, real *sizes,
- int *grouping, SparseMatrix graph0, real bounding_box_margin[], int *nrandom, int *nart, int nedgep,
- real shore_depth_tol, real edge_bridge_tol,
- real **xcombined, int *nverts, real **x_poly,
+ int n, int dim, double *x, double *sizes,
+ int *grouping, SparseMatrix graph0, double bounding_box_margin[], int *nrandom, int *nart, int nedgep,
+ double shore_depth_tol, double edge_bridge_tol,
+ double **xcombined, int *nverts, double **x_poly,
int *npolys, SparseMatrix *poly_lines, SparseMatrix *polys, int **polys_groups, SparseMatrix *poly_point_map,
SparseMatrix *country_graph, int highlight_cluster, int *flag){
*/
- real dist, avgdist;
- real *X;
+ double dist, avgdist;
+ double *X;
int N, nmax, i, j, k, igrp;
int *groups, K = *nart;/* average number of points added per side of rectangle */
- real avgsize[2], avgsz, h[2], p1, p0;
- real point[2];
+ double avgsize[2], avgsz, h[2], p1, p0;
+ double point[2];
int nadded[2];
int res;
- real delta[2];
+ double delta[2];
SparseMatrix graph = graph0;
- real bbox[4];
+ double bbox[4];
if (K < 0){
K = (int) 10/(1+n/400.);/* 0 if n > 3600*/
if (Verbose) fprintf(stderr, "avgsize = {%f, %f}\n",avgsize[0], avgsize[1]);
nmax = 2*n;
- X = MALLOC(sizeof(real)*dim*(n+nmax));
+ X = MALLOC(sizeof(double)*dim*(n+nmax));
groups = MALLOC(sizeof(int)*(n+nmax));
for (i = 0; i < n; i++) {
groups[i] = grouping[i];
if (shore_depth_tol < 0) {
shore_depth_tol = -(shore_depth_tol)*avgsz;
} else if (shore_depth_tol == 0){
- real area;
+ double area;
get_boundingbox(n, dim, x, sizes, bbox);
area = (bbox[1] - bbox[0])*(bbox[3] - bbox[2]);
- shore_depth_tol = sqrt(area/(real) n);
+ shore_depth_tol = sqrt(area/(double) n);
if (Verbose) fprintf(stderr,"setting shore length ======%f\n",shore_depth_tol);
} else {
}
#include <sparse/SparseMatrix.h>
#include <cgraph/cgraph.h>
-int make_map_from_rectangle_groups(int exclude_random, int include_OK_points, int n, int dim, real *x, real *sizes,
- int *grouping, SparseMatrix graph, real bounding_box_margin[], int *nrandom,int *nart, int nedgep,
- real shore_depth_tol, real edge_bridge_tol, real **xcombined, int *nverts, real **x_poly,
+int make_map_from_rectangle_groups(int exclude_random, int include_OK_points, int n, int dim, double *x, double *sizes,
+ int *grouping, SparseMatrix graph, double bounding_box_margin[], int *nrandom,int *nart, int nedgep,
+ double shore_depth_tol, double edge_bridge_tol, double **xcombined, int *nverts, double **x_poly,
int *npolys, SparseMatrix *poly_lines, SparseMatrix *polys, int **polys_groups, SparseMatrix *poly_point_map,
SparseMatrix *country_graph, int highlight_cluster, int *flag);
-void improve_contiguity(int n, int dim, int *grouping, SparseMatrix poly_point_map, real *x, SparseMatrix graph);
+void improve_contiguity(int n, int dim, int *grouping, SparseMatrix poly_point_map, double *x, SparseMatrix graph);
-void plot_dot_map(Agraph_t* gr, int n, int dim, real *x, SparseMatrix polys,
- SparseMatrix poly_lines, real line_width,
- const char *line_color, real *x_poly, int *polys_groups,
+void plot_dot_map(Agraph_t* gr, int n, int dim, double *x, SparseMatrix polys,
+ SparseMatrix poly_lines, double line_width,
+ const char *line_color, double *x_poly, int *polys_groups,
char **labels, float *fsz, float *r, float *g, float *b,
const char* opacity, SparseMatrix A, FILE*);
void map_optimal_coloring(int seed, SparseMatrix A, float *rgb_r, float *rgb_g, float *rgb_b);
-void map_palette_optimal_coloring(char *color_scheme, char *lightness, SparseMatrix A, real accuracy, int seed, float **rgb_r, float **rgb_g, float **rgb_b);
+void map_palette_optimal_coloring(char *color_scheme, char *lightness, SparseMatrix A, double accuracy, int seed, float **rgb_r, float **rgb_g, float **rgb_b);
enum {POLY_LINE_REAL_EDGE, POLY_LINE_NOT_REAL_EDGE};
#define neighbor(t, i, edim, elist) elist[(edim)*(t)+i]
static const int maxit = 100;
// Accuracy control (convergence criterion) for power_method
-static const real tolerance = 0.00001;
+static const double tolerance = 0.00001;
void power_method(void *A, int n, int K, int random_seed,
- real **eigv, real *eigs){
+ double **eigv, double *eigs){
/* find k-largest eigenvectors of a matrix A. Result in eigv. if eigv == NULL; memory will be allocated.
This converges only if the largest eigenvectors/values are real (e.g., if A is symmetric) and the
end for
return v1,v2,...
*/
- real **v, *u, *vv;
+ double **v, *u, *vv;
int iter = 0;
- real res, unorm;
+ double res, unorm;
int i, j, k;
- real uij;
+ double uij;
K = MAX(0, MIN(n, K));
assert(K <= n && K > 0);
- if (!(*eigv)) *eigv = MALLOC(sizeof(real)*n*K);
- v = MALLOC(sizeof(real*)*K);
+ if (!(*eigv)) *eigv = MALLOC(sizeof(double)*n*K);
+ v = MALLOC(sizeof(double*)*K);
- vv = MALLOC(sizeof(real)*n);
- u = MALLOC(sizeof(real)*n);
+ vv = MALLOC(sizeof(double)*n);
+ u = MALLOC(sizeof(double)*n);
srand(random_seed);
#include <sparse/general.h>
void power_method(void *A, int n, int K, int random_seed,
- real **eigv, real *eigs);
+ double **eigv, double *eigs);
int outer_iter;
int method;
int compatibility_method;
- real K;
+ double K;
fmt_t fmt;
int nneighbors;
int max_recursion;
- real angle_param;
- real angle;
+ double angle_param;
+ double angle;
} opts_t;
static char *fname;
{
int c;
char* cmd = argv[0];
- real s;
+ double s;
int i;
opterr = 0;
{
int k, j, mm, kk;
int dim = edge->dim;
- real* x = edge->x;
- real tt1[3]={0.15,0.5,0.85};
- real tt2[4]={0.15,0.4,0.6,0.85};
- real t, *tt;
+ double* x = edge->x;
+ double tt1[3]={0.15,0.5,0.85};
+ double tt2[4]={0.15,0.4,0.6,0.85};
+ double t, *tt;
for (j = 0; j < edge->npoints; j++){
if (j != 0) {
{
int k, j;
int dim = edge->dim;
- real* x = edge->x;
+ double* x = edge->x;
for (j = 0; j < edge->npoints; j++){
if (j != 0) agxbputc(xb, ':');
}
static void
-genBundleColors (pedge edge, agxbuf* xb, real maxwgt)
+genBundleColors (pedge edge, agxbuf* xb, double maxwgt)
{
int k, j, r, g, b;
- real len, t, len_total0 = 0;
+ double len, t, len_total0 = 0;
int dim = edge->dim;
- real* x = edge->x;
- real* lens = MALLOC(sizeof(real)*edge->npoints);
+ double* x = edge->x;
+ double* lens = MALLOC(sizeof(double)*edge->npoints);
for (j = 0; j < edge->npoints - 1; j++){
len = 0;
Agnode_t* n;
Agedge_t* e;
int i, j;
- real maxwgt = 0;
+ double maxwgt = 0;
pedge edge;
agxbuf xbuf;
unsigned char buf[BUFSIZ];
static int
bundle (Agraph_t* g, opts_t* opts)
{
- real *x = NULL;
- real *label_sizes = NULL;
+ double *x = NULL;
+ double *label_sizes = NULL;
int n_edge_label_nodes;
int dim = 2;
SparseMatrix A;
SparseMatrix B;
pedge* edges;
- real *xx, eps = 0.;
+ double *xx, eps = 0.;
int nz = 0;
int *ia, *ja, i, j, k;
int rv = 0;
ia = A->ia; ja = A->ja;
nz = A->nz;
- xx = MALLOC(sizeof(real)*nz*4);
+ xx = MALLOC(sizeof(double)*nz*4);
nz = 0;
dim = 4;
for (i = 0; i < A->m; i++){
{
int ret_code, type;
MM_typecode matcode;
- real *val = NULL, *v;
+ double *val = NULL, *v;
int *vali = NULL, i, m, n, *I = NULL, *J = NULL, nz;
void *vp = NULL;
SparseMatrix A = NULL;
type = mm_get_type(matcode);
switch (type) {
case MATRIX_TYPE_REAL:
- val = malloc(nz * sizeof(real));
+ val = malloc(nz * sizeof(double));
for (i = 0; i < nz; i++) {
fscanf(f, "%d %d %lg\n", &I[i], &J[i], &val[i]);
I[i]--; /* adjust from 1-based to 0-based */
if (mm_is_symmetric(matcode)) {
I = REALLOC(I, 2 * sizeof(int) * nz);
J = REALLOC(J, 2 * sizeof(int) * nz);
- val = REALLOC(val, 2 * sizeof(real) * nz);
+ val = REALLOC(val, 2 * sizeof(double) * nz);
nzold = nz;
for (i = 0; i < nzold; i++) {
if (I[i] != J[i]) {
} else if (mm_is_skew(matcode)) {
I = REALLOC(I, 2 * sizeof(int) * nz);
J = REALLOC(J, 2 * sizeof(int) * nz);
- val = REALLOC(val, 2 * sizeof(real) * nz);
+ val = REALLOC(val, 2 * sizeof(double) * nz);
nzold = nz;
for (i = 0; i < nzold; i++) {
assert(I[i] != J[i]); /* skew symm has no diag */
}
break;
case MATRIX_TYPE_COMPLEX:
- val = malloc(2 * nz * sizeof(real));
+ val = malloc(2 * nz * sizeof(double));
v = val;
for (i = 0; i < nz; i++) {
fscanf(f, "%d %d %lg %lg\n", &I[i], &J[i], &v[0], &v[1]);
if (mm_is_symmetric(matcode)) {
I = REALLOC(I, 2 * sizeof(int) * nz);
J = REALLOC(J, 2 * sizeof(int) * nz);
- val = REALLOC(val, 4 * sizeof(real) * nz);
+ val = REALLOC(val, 4 * sizeof(double) * nz);
nzold = nz;
for (i = 0; i < nzold; i++) {
if (I[i] != J[i]) {
} else if (mm_is_skew(matcode)) {
I = REALLOC(I, 2 * sizeof(int) * nz);
J = REALLOC(J, 2 * sizeof(int) * nz);
- val = REALLOC(val, 4 * sizeof(real) * nz);
+ val = REALLOC(val, 4 * sizeof(double) * nz);
nzold = nz;
for (i = 0; i < nzold; i++) {
assert(I[i] != J[i]); /* skew symm has no diag */
} else if (mm_is_hermitian(matcode)) {
I = REALLOC(I, 2 * sizeof(int) * nz);
J = REALLOC(J, 2 * sizeof(int) * nz);
- val = REALLOC(val, 4 * sizeof(real) * nz);
+ val = REALLOC(val, 4 * sizeof(double) * nz);
nzold = nz;
for (i = 0; i < nzold; i++) {
if (I[i] != J[i]) {
if (format == FORMAT_CSR) {
A = SparseMatrix_from_coordinate_arrays(nz, m, n, I, J, vp,
- type, sizeof(real));
+ type, sizeof(double));
} else {
A = SparseMatrix_new(m, n, 1, type, FORMAT_COORD);
A = SparseMatrix_coordinate_form_add_entries(A, nz, I, J, vp);
#define MALLOC malloc
#define test_flag(a, flag) ((a)&(flag))
-#define real double
#define BUFS 1024
typedef struct {
static char *cmd;
-static real Hue2RGB(real v1, real v2, real H)
+static double Hue2RGB(double v1, double v2, double H)
{
if (H < 0.0)
H += 1.0;
return v1;
}
-static char *hue2rgb(real hue, char *color)
+static char *hue2rgb(double hue, char *color)
{
- real v1, v2, lightness = .5, saturation = 1;
+ double v1, v2, lightness = .5, saturation = 1;
int red, blue, green;
if (lightness < 0.5)
}
static Agraph_t *makeDotGraph(SparseMatrix A, char *name, int dim,
- real * x, int with_color, int with_label, int with_val)
+ double * x, int with_color, int with_label, int with_val)
{
Agraph_t *g;
Agnode_t *n;
Agsym_t *sym, *sym2 = NULL, *sym3 = NULL;
int *ia = A->ia;
int *ja = A->ja;
- real *val = (real *) (A->a);
+ double *val = (double *) (A->a);
Agnode_t **arr = N_NEW(A->m, Agnode_t *);
- real *color = NULL;
+ double *color = NULL;
char cstring[8];
name = strip_dir(name);
}
if (with_color) {
- real maxdist = 0.;
- real mindist = 0.;
+ double maxdist = 0.;
+ double mindist = 0.;
int first = TRUE;
sym2 = agattr(g, AGEDGE, "color", "");
sym3 = agattr(g, AGEDGE, "wt", "");
agattr(g, AGRAPH, "bgcolor", "black");
- color = N_NEW(A->nz, real);
+ color = N_NEW(A->nz, double);
for (n = agfstnode(g); n; n = agnxtnode(g, n)) {
i = ND_id(n);
if (A->type != MATRIX_TYPE_REAL) {
#include <sparse/QuadTree.h>
static int splines_intersect(int dim, int u1, int v1, int u2, int v2,
- real cos_critical, int check_edges_with_same_endpoint,
+ double cos_critical, int check_edges_with_same_endpoint,
char *xsplines1, char *xsplines2){
/* u1, v2 an u2, v2: the node index of the two edn points of two edges.
cos_critical: cos of critical angle
*/
int itmp;
int len1 = 100, len2 = 100;
- real *x1, *x2;
+ double *x1, *x2;
int ns1 = 0, ns2 = 0;
int i, j, iter1 = 0, iter2 = 0;
- real cos_a, tmp[2];
+ double cos_a, tmp[2];
int endp1 = 0, endp2 = 0;
tmp[0] = tmp[1] = 0;
- x1 = MALLOC(sizeof(real)*len1);
- x2 = MALLOC(sizeof(real)*len2);
+ x1 = MALLOC(sizeof(double)*len1);
+ x2 = MALLOC(sizeof(double)*len2);
assert(dim <= 3);
/* if two end points are the same, make sure they are the first in each edge */
xsplines1++;
if (ns1*dim >= len1){
len1 = ns1*dim + (int)MAX(10, 0.2*ns1*dim);
- x1 = REALLOC(x1, sizeof(real)*len1);
+ x1 = REALLOC(x1, sizeof(double)*len1);
}
}
if (endp1){/* pad the end point at the last position */
ns1++;
if (ns1*dim >= len1){
len1 = ns1*dim + (int)MAX(10, 0.2*ns1*dim);
- x1 = REALLOC(x1, sizeof(real)*len1);
+ x1 = REALLOC(x1, sizeof(double)*len1);
}
x1[(ns1-1)*dim] = tmp[0]; x1[(ns1-1)*dim + 1] = tmp[1];
}
xsplines2++;
if (ns2*dim >= len2){
len2 = ns2*dim + (int)MAX(10, 0.2*ns2*dim);
- x2 = REALLOC(x2, sizeof(real)*len2);
+ x2 = REALLOC(x2, sizeof(double)*len2);
}
}
if (endp2){/* pad the end point at the last position */
ns2++;
if (ns2*dim >= len2){
len2 = ns2*dim + (int)MAX(10, 0.2*ns2*dim);
- x2 = REALLOC(x2, sizeof(real)*len2);
+ x2 = REALLOC(x2, sizeof(double)*len2);
}
x2[(ns2-1)*dim] = tmp[0]; x2[(ns2-1)*dim + 1] = tmp[1];
}
}
-Agraph_t* edge_distinct_coloring(char *color_scheme, char *lightness, Agraph_t* g, real angle, real accuracy, int check_edges_with_same_endpoint, int seed){
+Agraph_t* edge_distinct_coloring(char *color_scheme, char *lightness, Agraph_t* g, double angle, double accuracy, int check_edges_with_same_endpoint, int seed){
/* color the edges of a graph so that conflicting edges are as dinstrinct in color as possibl.
color_scheme: rgb, lab, gray, or a list of comma separaterd RGB colors in hex, like #ff0000,#00ff00
lightness: of the form 0,70, specifying the range of lightness of LAB color. Ignored if scheme is not COLOR_LAB.
. are not consider conflict.
seed: random_seed. If negative, consider -seed as the number of random start iterations
*/
- real *x = NULL;
+ double *x = NULL;
int dim = 2;
SparseMatrix A, B, C;
int *irn, *jcn, nz, nz2 = 0;
- real cos_critical = cos(angle/180*3.14159), cos_a;
+ double cos_critical = cos(angle/180*3.14159), cos_a;
int u1, v1, u2, v2, i, j;
- real *colors = NULL;
+ double *colors = NULL;
int flag, ne;
char **xsplines = NULL;
int cdim;
}
}
#ifdef TIME
- fprintf(stderr, "cpu for dual graph =%10.3f", ((real) (clock() - start))/CLOCKS_PER_SEC);
+ fprintf(stderr, "cpu for dual graph =%10.3f", ((double) (clock() - start))/CLOCKS_PER_SEC);
#endif
} else {
}
}
#ifdef TIME
- fprintf(stderr, "cpu for dual graph (splines) =%10.3f\n", ((real) (clock() - start))/CLOCKS_PER_SEC);
+ fprintf(stderr, "cpu for dual graph (splines) =%10.3f\n", ((double) (clock() - start))/CLOCKS_PER_SEC);
#endif
}
C = SparseMatrix_from_coordinate_format(B);
flag = node_distinct_coloring(color_scheme, lightness, weightedQ, C, accuracy, iter_max, seed, &cdim, &colors);
if (flag) goto RETURN;
#ifdef TIME
- fprintf(stderr, "cpu for color assignmment =%10.3f\n", ((real) (clock() - start))/CLOCKS_PER_SEC);
+ fprintf(stderr, "cpu for color assignmment =%10.3f\n", ((double) (clock() - start))/CLOCKS_PER_SEC);
#endif
}
#pragma once
-Agraph_t* edge_distinct_coloring(char *color_scheme, char *lightness, Agraph_t* g, real angle, real accuracy, int check_edges_with_same_endpoint, int seed);
+Agraph_t* edge_distinct_coloring(char *color_scheme, char *lightness, Agraph_t* g, double angle, double accuracy, int check_edges_with_same_endpoint, int seed);
#include <edgepaint/furtherest_point.h>
#include <string.h>
-static real dist(int dim, real *x, real *y){
+static double dist(int dim, double *x, double *y){
int k;
- real d = 0;
+ double d = 0;
for (k = 0; k < dim; k++) d += (x[k] - y[k])*(x[k]-y[k]);
return sqrt(d);
}
-static real distance_to_group(int k, int dim, real *wgt, real *pts, real *center, real (*usr_dist)(int, real*, real*)){
+static double distance_to_group(int k, int dim, double *wgt, double *pts, double *center, double (*usr_dist)(int, double*, double*)){
int i;
- real distance = -1, dist_min = 0;
+ double distance = -1, dist_min = 0;
if (!wgt){
for (i = 0; i < k; i++){
distance = (*usr_dist)(dim, &(pts[i*dim]), center);
return dist_min;
}
-void furtherest_point(int k, int dim, real *wgt, real *pts, real *center, real width, int max_level, real (*usr_dist)(int, real*, real*), real *dist_max, real **argmax){
+void furtherest_point(int k, int dim, double *wgt, double *pts, double *center, double width, int max_level, double (*usr_dist)(int, double*, double*), double *dist_max, double **argmax){
/* Assume that in the box defined by {center, width} are feasible;
find, in the, a point "furtherest_point" that is furtherest away from a group of k-points pts, using quadtree,
with up to max_level. Here the distance of a point to a group of point is defined as the minimum
QuadTree *candidates, *ctmp;/* a cadidate array of quadtrees */
int ncandidates2 = 10, ncandidates2_max = 10;
QuadTree *candidates2;/* a cadidate array of quadtrees */
- real distance;
+ double distance;
int level = 0;
int i, ii, j, pruned;
- real wmax = 0;
+ double wmax = 0;
if (!usr_dist) usr_dist = dist;
if (wgt){
qt0 = qt = QuadTree_new(dim, center, width, max_level);
qt->total_weight = *dist_max = distance_to_group(k, dim, wgt, pts, center, usr_dist);/* store distance in total_weight */
- if (!(*argmax)) *argmax = MALLOC(sizeof(real)*dim);
- memcpy(*argmax, center, sizeof(real)*dim);
+ if (!(*argmax)) *argmax = MALLOC(sizeof(double)*dim);
+ memcpy(*argmax, center, sizeof(double)*dim);
candidates = MALLOC(sizeof(qt)*ncandidates_max);
candidates2 = MALLOC(sizeof(qt)*ncandidates2_max);
}
distance = qt->total_weight;/* total_weight is used to store the distance from the center to the group */
- if (distance + wmax*sqrt(((real) dim))*qt->width < *dist_max) continue;/* this could happen if this candidate was entered into the list earlier than a better one later in the list */
+ if (distance + wmax*sqrt(((double) dim))*qt->width < *dist_max) continue;/* this could happen if this candidate was entered into the list earlier than a better one later in the list */
qt->qts = MALLOC(sizeof(QuadTree)*(1<<dim));
for (ii = 0; ii < 1<<dim; ii++) {
qt->qts[ii] = QuadTree_new_in_quadrant(qt->dim, qt->center, (qt->width)/2, max_level, ii);
for (j = 0; j < dim; j++) fprintf(stderr,"%f, ", qt->qts[ii]->center[j]);
fprintf(stderr,"}\n");
}
- memcpy(*argmax, qt->qts[ii]->center, sizeof(real)*dim);
- } else if (distance + wmax*sqrt(((real) dim))*(qt->width)/2 < *dist_max){
+ memcpy(*argmax, qt->qts[ii]->center, sizeof(double)*dim);
+ } else if (distance + wmax*sqrt(((double) dim))*(qt->width)/2 < *dist_max){
pruned = TRUE;
}
if (!pruned){
}
-void furtherest_point_in_list(int k, int dim, real *wgt, real *pts, QuadTree qt, int max_level,
- real (*usr_dist)(int, real*, real*), real *dist_max, real **argmax){
+void furtherest_point_in_list(int k, int dim, double *wgt, double *pts, QuadTree qt, int max_level,
+ double (*usr_dist)(int, double*, double*), double *dist_max, double **argmax){
/* Given a list of points in the Euclidean space contained in the quadtree qt (called the feasible points), find among them one that
is closest to another list of point {dim, k, pts}.
QuadTree *candidates, *ctmp;/* a candidate array of quadtrees */
int ncandidates2 = 10, ncandidates2_max = 10;
QuadTree *candidates2;/* a candidate array of quadtrees */
- real distance;
+ double distance;
int level = 0;
int i, ii, j, pruned;
- real *average;
- real wmax = 0.;
+ double *average;
+ double wmax = 0.;
if (!usr_dist) usr_dist = dist;
average = qt->average;
qt->total_weight = *dist_max = distance_to_group(k, dim, wgt, pts, average, usr_dist);/* store distance in total_weight */
- if (!(*argmax)) *argmax = MALLOC(sizeof(real)*dim);
- memcpy(*argmax, average, sizeof(real)*dim);
+ if (!(*argmax)) *argmax = MALLOC(sizeof(double)*dim);
+ memcpy(*argmax, average, sizeof(double)*dim);
candidates = MALLOC(sizeof(qt)*ncandidates_max);
candidates2 = MALLOC(sizeof(qt)*ncandidates2_max);
}
distance = qt->total_weight;/* total_weight is used to store the distance from average feasible points to the group */
- if (qt->n == 1 || distance + wmax*2*sqrt(((real) dim))*qt->width < *dist_max) continue;/* this could happen if this candidate was entered into the list earlier than a better one later in the list. Since the distance
+ if (qt->n == 1 || distance + wmax*2*sqrt(((double) dim))*qt->width < *dist_max) continue;/* this could happen if this candidate was entered into the list earlier than a better one later in the list. Since the distance
is from the average of the feasible points in the square which may not be at the center */
if (!(qt->qts)) continue;
for (j = 0; j < dim; j++) fprintf(stderr,"%f, ", qt->qts[ii]->average[j]);
fprintf(stderr,"}\n");
}
- memcpy(*argmax, qt->qts[ii]->average, sizeof(real)*dim);
- } else if (distance + wmax*sqrt(((real) dim))*(qt->width) < *dist_max){/* average feasible point in this square is too close to the point set */
+ memcpy(*argmax, qt->qts[ii]->average, sizeof(double)*dim);
+ } else if (distance + wmax*sqrt(((double) dim))*(qt->width) < *dist_max){/* average feasible point in this square is too close to the point set */
pruned = TRUE;
}
if (!pruned){
#pragma once
-void furtherest_point(int k, int dim, real *wgt, real *pts, real *center, real width, int max_level, real (*usr_dist)(int, real*, real*), real *dist_max, real **argmax);
-void furtherest_point_in_list(int k, int dim, real *wgt, real *pts, QuadTree qt, int max_level,
- real (*usr_dist)(int, real*, real*), real *dist_max, real **argmax);
+void furtherest_point(int k, int dim, double *wgt, double *pts, double *center, double width, int max_level, double (*usr_dist)(int, double*, double*), double *dist_max, double **argmax);
+void furtherest_point_in_list(int k, int dim, double *wgt, double *pts, QuadTree qt, int max_level,
+ double (*usr_dist)(int, double*, double*), double *dist_max, double **argmax);
#include <sparse/general.h>
#include <math.h>
-static real cross(real *u, real *v){
+static double cross(double *u, double *v){
return u[0]*v[1] - u[1]*v[0];
}
*/
-static real dist(int dim, real *x, real *y){
+static double dist(int dim, double *x, double *y){
int k;
- real d = 0;
+ double d = 0;
for (k = 0; k < dim; k++) d += (x[k] - y[k])*(x[k]-y[k]);
return sqrt(d);
}
-static real point_line_distance(real *p, real *q, real *r){
+static double point_line_distance(double *p, double *q, double *r){
/* distance between point p and line q--r */
enum {dim = 2};
- real t = 0, b = 0;
+ double t = 0, b = 0;
int i;
- real tmp;
+ double tmp;
/* t = ((p - q).(r - q))/((r - q).(r - q)) gives the position of the project of p on line r--q */
for (i = 0; i < dim; i++){
}
-static real line_segments_distance(real *p1, real *p2, real *q1, real *q2){
+static double line_segments_distance(double *p1, double *p2, double *q1, double *q2){
/* distance between line segments p1--p2 and q1--q2 */
- real t1, t2, t3, t4;
+ double t1, t2, t3, t4;
t1 = point_line_distance(p1, q1, q2);
t2 = point_line_distance(p2, q1, q2);
t3 = point_line_distance(q1, p1, p2);
}
-real intersection_angle(real *p1, real *p2, real *q1, real *q2){
+double intersection_angle(double *p1, double *p2, double *q1, double *q2){
/* give two lines p1--p2 and q1--q2, find their intersection agle
and return Abs[Cos[theta]] of that angle.
. from -1 to 1.
*/
enum {dim = 2};
- real r[dim], s[dim], qp[dim];
- real rnorm = 0, snorm = 0, b, t, u;
- // real epsilon = sqrt(MACHINEACC), close = 0.01;
- //this may be better. Apply to ngk10_4 and look at double edge between 28 and 43. real epsilon = sin(10/180.), close = 0.1;
- real epsilon = sin(1/180.), close = 0.01;
+ double r[dim], s[dim], qp[dim];
+ double rnorm = 0, snorm = 0, b, t, u;
+ // double epsilon = sqrt(MACHINEACC), close = 0.01;
+ //this may be better. Apply to ngk10_4 and look at double edge between 28 and 43. double epsilon = sin(10/180.), close = 0.1;
+ double epsilon = sin(1/180.), close = 0.01;
int line_dist_close;
int i;
- real res;
+ double res;
for (i = 0; i < dim; i++) {
r[i] = p2[i] - p1[i];
u = cross(qp, r)/b;
if ((t >= 0 && t <= 1 && u >= 0 && u <= 1) /* they intersect */
|| line_dist_close){/* or lines are close */
- real rs = 0;
+ double rs = 0;
if (rnorm*snorm < MACHINEACC) return 0;
for (i = 0; i < dim; i++){
rs += r[i]*s[i];
#pragma once
-real intersection_angle(real *p1, real *p2, real *q1, real *q2);
+double intersection_angle(double *p1, double *p2, double *q1, double *q2);
return color_lab_init(L, A, B);
}
-void LAB2RGB_real_01(real *color){
+void LAB2RGB_real_01(double *color){
/* convert an array[3] of LAB colors to RGB between 0 to 1, in place */
color_rgb rgb;
color_lab lab;
color_xyz RGB2XYZ(color_rgb color);
color_rgb XYZ2RGB(color_xyz color);
color_lab RGB2LAB(color_rgb color);
-void LAB2RGB_real_01(real *color); /* convert an array[3] of LAB colors to RGB between 0 to 1, in place */
+void LAB2RGB_real_01(double *color); /* convert an array[3] of LAB colors to RGB between 0 to 1, in place */
color_rgb LAB2RGB(color_lab color);
color_rgb color_rgb_init(double r, double g, double b);
color_xyz color_xyz_init(double x, double y, double z);
#include <sparse/color_palette.h>
#include <string.h>
-static real mydist(int dim, real *x, real *y){
+static double mydist(int dim, double *x, double *y){
int k;
- real d = 0;
+ double d = 0;
for (k = 0; k < dim; k++) d += (x[k] - y[k])*(x[k]-y[k]);
return sqrt(d);
}
-static void node_distinct_coloring_internal2(int scheme, QuadTree qt, int weightedQ, SparseMatrix A, int cdim, real accuracy, int iter_max, int seed, real *colors,
- real *color_diff0, real *color_diff_sum0){
+static void node_distinct_coloring_internal2(int scheme, QuadTree qt, int weightedQ, SparseMatrix A, int cdim, double accuracy, int iter_max, int seed, double *colors,
+ double *color_diff0, double *color_diff_sum0){
/* here we assume the graph is connected. And that the matrix is symmetric */
int i, j, *ia, *ja, n, k = 0;
int max_level;
- real center[3];
- real width;
- real *a = NULL;
- real *x;
- real dist_max;
- real color_diff = 0, color_diff_old;
- real color_diff_sum = 0, color_diff_sum_old, *cc;
+ double center[3];
+ double width;
+ double *a = NULL;
+ double *x;
+ double dist_max;
+ double color_diff = 0, color_diff_old;
+ double color_diff_sum = 0, color_diff_sum_old, *cc;
int iter = 0;
- real cspace_size = 0.7;
- real red[3], black[3], min;
+ double cspace_size = 0.7;
+ double red[3], black[3], min;
int flag = 0, imin;
- real *wgt = NULL;
+ double *wgt = NULL;
assert(accuracy > 0);
max_level = MAX(1, -log(accuracy)/log(2.));
ia = A->ia;
ja = A->ja;
if (A->type == MATRIX_TYPE_REAL && A->a){
- a = (real*) A->a;
+ a = (double*) A->a;
}
/* cube [0, cspace_size]^3: only uised if not LAB */
srand(seed);
for (i = 0; i < n*cdim; i++) colors[i] = cspace_size*drand();
- x = MALLOC(sizeof(real)*cdim*n);
- if (weightedQ) wgt = MALLOC(sizeof(real)*n);
+ x = MALLOC(sizeof(double)*cdim*n);
+ if (weightedQ) wgt = MALLOC(sizeof(double)*n);
color_diff = 0; color_diff_old = -1;
color_diff_sum = 0; color_diff_sum_old = -1;
k = 0;
for (j = ia[i]; j < ia[i+1]; j++){
if (ja[j] == i) continue;
- memcpy(&(x[k*cdim]), &(colors[ja[j]*cdim]), sizeof(real)*cdim);
+ memcpy(&(x[k*cdim]), &(colors[ja[j]*cdim]), sizeof(double)*cdim);
if (wgt && a) wgt[k] = a[j];
k++;
}
free(x);
}
-static void node_distinct_coloring_internal(int scheme, QuadTree qt, int weightedQ, SparseMatrix A, int cdim, real accuracy, int iter_max, int seed, real *colors){
+static void node_distinct_coloring_internal(int scheme, QuadTree qt, int weightedQ, SparseMatrix A, int cdim, double accuracy, int iter_max, int seed, double *colors){
int i;
- real color_diff;
- real color_diff_sum;
+ double color_diff;
+ double color_diff_sum;
if (seed < 0) {
/* do multiple iterations and pick the best */
int iter, seed_max = -1;
- real color_diff_max = -1;
+ double color_diff_max = -1;
srand(123);
iter = -seed;
for (i = 0; i < iter; i++){
}
-int node_distinct_coloring(char *color_scheme, char *lightness, int weightedQ, SparseMatrix A0, real accuracy, int iter_max, int seed, int *cdim0, real **colors){
+int node_distinct_coloring(char *color_scheme, char *lightness, int weightedQ, SparseMatrix A0, double accuracy, int iter_max, int seed, int *cdim0, double **colors){
/*
for a graph A, get a distinctive color of its nodes so that the color distance among all neighboring nodes are maximized. Here
color distance on a node is defined as the minimum of color differences between a node and its neighbors (or the minimum of weighted color differences if weightedQ = true,
SparseMatrix B, A = A0;
int ncomps, *comps = NULL, *comps_ptr = NULL;
int nn, n;
- real *ctmp;
+ double *ctmp;
int i, j, jj, nnodes = 0;
QuadTree qt = NULL;
int cdim;
}
if (!(*colors)) {
- *colors = MALLOC(sizeof(real)*cdim*n);
+ *colors = MALLOC(sizeof(double)*cdim*n);
}
- ctmp = MALLOC(sizeof(real)*cdim*n);
+ ctmp = MALLOC(sizeof(double)*cdim*n);
B = SparseMatrix_symmetrize(A, FALSE);
A = B;
for (j = comps_ptr[i]; j < comps_ptr[i+1]; j++){
jj = j - comps_ptr[i];
- memcpy(&((*colors)[comps[j]*cdim]), &(ctmp[jj*cdim]), cdim*sizeof(real));
+ memcpy(&((*colors)[comps[j]*cdim]), &(ctmp[jj*cdim]), cdim*sizeof(double));
}
SparseMatrix_delete(B);
}
enum {COLOR_RGB, COLOR_GRAY, COLOR_LAB};
enum {ERROR_BAD_COLOR_SCHEME = -9};
-int node_distinct_coloring(char *color_scheme, char *lightness, int weightedQ, SparseMatrix A, real accuracy, int iter_max, int seed, int *cdim, real **colors);
+int node_distinct_coloring(char *color_scheme, char *lightness, int weightedQ, SparseMatrix A, double accuracy, int iter_max, int seed, int *cdim, double **colors);
grid->R = NULL;
grid->next = NULL;
grid->prev = NULL;
- grid->inks = MALLOC(sizeof(real)*(A->m));
+ grid->inks = MALLOC(sizeof(double)*(A->m));
grid->edges = edges;
grid->delete_top_level_A = 0;
grid->total_ink = -1;
if (level == 0){
- real total_ink = 0;
+ double total_ink = 0;
for (i = 0; i < n; i++) {
(grid->inks)[i] = ink1(edges[i]);
total_ink += (grid->inks)[i];
free(grid);
}
-static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_establish(Agglomerative_Ink_Bundling grid, int *pick, real angle_param, real angle){
+static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_establish(Agglomerative_Ink_Bundling grid, int *pick, double angle_param, double angle){
/* pick is a work array of dimension n, with n the total number of original edges */
int *matching;
SparseMatrix A = grid->A;
int i, j, k, jj, jc, jmax, ni, nj, npicks;
int *mask;
pedge *edges = grid->edges;
- real *inks = grid->inks, *cinks, inki, inkj;
- real gain, maxgain, minink, total_gain = 0;
+ double *inks = grid->inks, *cinks, inki, inkj;
+ double gain, maxgain, minink, total_gain = 0;
int *ip = NULL, *jp = NULL, ie;
Vector *cedges;/* a table listing the content of bundled edges in the coarsen grid.
cedges[i] contain the list of origonal edges that make up the bundle i in the next level */
- real ink0, ink1, grand_total_ink = 0, grand_total_gain = 0;
+ double ink0, ink1, grand_total_ink = 0, grand_total_gain = 0;
point_t meet1, meet2;
if (Verbose > 1) fprintf(stderr,"level ===================== %d, n = %d\n",grid->level, n);
cedges = MALLOC(sizeof(Vector)*n);
- cinks = MALLOC(sizeof(real)*n);
+ cinks = MALLOC(sizeof(double)*n);
for (i = 0; i < n; i++) cedges[i] = Vector_new(1, sizeof(int), NULL);
if (grid->level > 0){
}
matching = MALLOC(sizeof(int)*n);
- mask = MALLOC(sizeof(real)*n);
+ mask = MALLOC(sizeof(double)*n);
for (i = 0; i < n; i++) mask[i] = -1;
assert(n == A->n);
if (nc >= 1 && total_gain > 0){
/* now set up restriction and prolongation operator */
SparseMatrix P, R, R1, R0, B, cA;
- real one = 1.;
+ double one = 1.;
Agglomerative_Ink_Bundling cgrid;
R1 = SparseMatrix_new(nc, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD);
return grid;
}
-static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_new(SparseMatrix A0, pedge *edges, real angle_param, real angle){
+static Agglomerative_Ink_Bundling Agglomerative_Ink_Bundling_new(SparseMatrix A0, pedge *edges, double angle_param, double angle){
/* give a link of edges and their nearest neighbor graph, return a multilevel of edge bundling based on ink saving */
Agglomerative_Ink_Bundling grid;
int *pick;
return grid;
}
-static pedge* agglomerative_ink_bundling_internal(int dim, SparseMatrix A, pedge* edges, int nneighbors, int *recurse_level, int MAX_RECURSE_LEVEL, real angle_param, real angle, int open_gl, real *current_ink, real *ink00, int *flag){
+static pedge* agglomerative_ink_bundling_internal(int dim, SparseMatrix A, pedge* edges, int nneighbors, int *recurse_level, int MAX_RECURSE_LEVEL, double angle_param, double angle, int open_gl, double *current_ink, double *ink00, int *flag){
int i, j, jj, k;
int *ia, *ja;
int *pick;
Agglomerative_Ink_Bundling grid, cgrid;
SparseMatrix R;
- real ink0, ink1;
+ double ink0, ink1;
point_t meet1, meet2;
pedge e;
- real TOL = 0.0001, wgt_all;
+ double TOL = 0.0001, wgt_all;
clock_t start;
(*recurse_level)++;
start = clock();
grid = Agglomerative_Ink_Bundling_new(A, edges, angle_param, angle);
if (Verbose > 1)
- fprintf(stderr, "CPU for agglomerative bundling %f\n", ((real) (clock() - start))/CLOCKS_PER_SEC);
+ fprintf(stderr, "CPU for agglomerative bundling %f\n", ((double) (clock() - start))/CLOCKS_PER_SEC);
ink0 = grid->total_ink;
/* find coarsest */
improved later */
e = edges[jj] = pedge_double(edges[jj]);
- e->wgts = REALLOC(e->wgts, sizeof(real)*4);
+ e->wgts = REALLOC(e->wgts, sizeof(double)*4);
e->x[1*dim] = meet1.x;
e->x[1*dim+1] = meet1.y;
e->x[2*dim] = meet2.x;
}
} else {
pedge *mid_edges, midedge;/* middle section of edges that will be bundled again */
- real *xx;
+ double *xx;
int ne, npp, l;
SparseMatrix A_mid;
- real eps = 0., wgt, total_wgt = 0;
+ double eps = 0., wgt, total_wgt = 0;
/* make new edges using meet1 and meet2.
ja = R->ja;
ne = R->m;
mid_edges = MALLOC(sizeof(pedge)*ne);
- xx = MALLOC(sizeof(real)*4*ne);
+ xx = MALLOC(sizeof(double)*4*ne);
for (i = 0; i < R->m; i++){
pick = &(ja[ia[i]]);
wgt = 0.;
}
-pedge* agglomerative_ink_bundling(int dim, SparseMatrix A, pedge* edges, int nneighbor, int MAX_RECURSE_LEVEL, real angle_param, real angle, int open_gl, int *flag){
+pedge* agglomerative_ink_bundling(int dim, SparseMatrix A, pedge* edges, int nneighbor, int MAX_RECURSE_LEVEL, double angle_param, double angle, int open_gl, int *flag){
int recurse_level = 0;
- real current_ink = -1, ink0;
+ double current_ink = -1, ink0;
pedge *edges2;
ink_count = 0;
if (Verbose > 1)
- fprintf(stderr,"initial total ink = %f, final total ink = %f, inksaving = %f percent, total ink_calc = %f, avg ink_calc per edge = %f\n", ink0, current_ink, (ink0-current_ink)/ink0, ink_count, ink_count/(real) A->m);
+ fprintf(stderr,"initial total ink = %f, final total ink = %f, inksaving = %f percent, total ink_calc = %f, avg ink_calc per edge = %f\n", ink0, current_ink, (ink0-current_ink)/ink0, ink_count, ink_count/(double) A->m);
return edges2;
}
SparseMatrix R;/* striction mtrix from level to level + 1*/
Agglomerative_Ink_Bundling next;
Agglomerative_Ink_Bundling prev;
- real *inks; /* amount of ink needed to draw this edge/bundle. Dimension n. */
- real total_ink; /* amount of ink needed to draw this edge/bundle. Dimension n. */
+ double *inks; /* amount of ink needed to draw this edge/bundle. Dimension n. */
+ double total_ink; /* amount of ink needed to draw this edge/bundle. Dimension n. */
pedge* edges; /* the original edge info. This does not vary level to level and is of dimenion n0, where n0 is the number of original edges */
int delete_top_level_A;/*whether the top level matrix should be deleted on garbage collecting the grid */
};
-pedge* agglomerative_ink_bundling(int dim, SparseMatrix A, pedge* edges, int nneighbor, int max_recursion, real angle_param, real angle, int open_gl, int *flag);
+pedge* agglomerative_ink_bundling(int dim, SparseMatrix A, pedge* edges, int nneighbor, int max_recursion, double angle_param, double angle, int open_gl, int *flag);
extern int *clusters_global;
#endif
-static real norm(int n, real *x){
- real res = 0;
+static double norm(int n, double *x){
+ double res = 0;
int i;
for (i = 0; i < n; i++) res += x[i]*x[i];
return sqrt(res);
}
-static real sqr_dist(int dim, real *x, real *y){
+static double sqr_dist(int dim, double *x, double *y){
int i;
- real res = 0;
+ double res = 0;
for (i = 0; i < dim; i++) res += (x[i] - y[i])*(x[i] - y[i]);
return res;
}
-static real dist(int dim, real *x, real *y){
+static double dist(int dim, double *x, double *y){
return sqrt(sqr_dist(dim,x,y));
}
-pedge pedge_new(int np, int dim, real *x){
+pedge pedge_new(int np, int dim, double *x){
pedge e;
e = MALLOC(sizeof(struct pedge_struct));
e->npoints = np;
e->dim = dim;
e->len = np;
- e->x = MALLOC(dim*e->len*sizeof(real));
- memcpy(e->x, x, dim*e->len*sizeof(real));
+ e->x = MALLOC(dim*e->len*sizeof(double));
+ memcpy(e->x, x, dim*e->len*sizeof(double));
e->edge_length = dist(dim, &x[0*dim], &x[(np-1)*dim]);
e->wgt = 1.;
e->wgts = NULL;
return e;
}
-pedge pedge_wgt_new(int np, int dim, real *x, real wgt){
+pedge pedge_wgt_new(int np, int dim, double *x, double wgt){
pedge e;
int i;
e->npoints = np;
e->dim = dim;
e->len = np;
- e->x = MALLOC(dim*e->len*sizeof(real));
- memcpy(e->x, x, dim*e->len*sizeof(real));
+ e->x = MALLOC(dim*e->len*sizeof(double));
+ memcpy(e->x, x, dim*e->len*sizeof(double));
e->edge_length = dist(dim, &x[0*dim], &x[(np-1)*dim]);
e->wgt = wgt;
- e->wgts = MALLOC(sizeof(real)*(np - 1));
+ e->wgts = MALLOC(sizeof(double)*(np - 1));
for (i = 0; i < np - 1; i++) e->wgts[i] = wgt;
return e;
pedge pedge_flip(pedge e){
/* flip the polyline so that last point becomes the first, second last the second, etc*/
- real *y;
- real *x = e->x;
+ double *y;
+ double *x = e->x;
int i, dim = e->dim;
int n = e->npoints;
- y = MALLOC(sizeof(real)*e->dim);
+ y = MALLOC(sizeof(double)*e->dim);
for (i = 0; i < e->npoints/2; i++){
- memcpy(y, &x[i*dim], sizeof(real)*dim);
- memcpy(&x[(n-1-i)*dim], &x[i*dim], sizeof(real)*dim);
- memcpy(&x[i*dim], y, sizeof(real)*dim);
+ memcpy(y, &x[i*dim], sizeof(double)*dim);
+ memcpy(&x[(n-1-i)*dim], &x[i*dim], sizeof(double)*dim);
+ memcpy(&x[i*dim], y, sizeof(double)*dim);
}
free(y);
return e;
}
-static real edge_compatibility(pedge e1, pedge e2){
+static double edge_compatibility(pedge e1, pedge e2){
/* two edges are u1->v1, u2->v2.
return 1 if two edges are exactly the same, 0 if they are very different.
*/
- real *u1, *v1, *u2, *v2, *u, dist1, dist2, len1, len2;
+ double *u1, *v1, *u2, *v2, *u, dist1, dist2, len1, len2;
int dim = e1->dim, flipped = FALSE;
u1 = e1->x;
}
}
-static real edge_compatibility_full(pedge e1, pedge e2){
+static double edge_compatibility_full(pedge e1, pedge e2){
/* two edges are u1->v1, u2->v2.
return 1 if two edges are exactly the same, 0 if they are very different.
This is based on Holten and van Wijk's paper
*/
- real *u1, *v1, *u2, *v2, *u, dist1, dist2, len1, len2, len;
- real tmp, ca, cp, cs;
+ double *u1, *v1, *u2, *v2, *u, dist1, dist2, len1, len2, len;
+ double tmp, ca, cp, cs;
int dim = e1->dim, flipped = FALSE, i;
u1 = e1->x;
void pedge_export_gv(FILE *fp, int ne, pedge *edges){
pedge edge;
- real *x, t;
+ double *x, t;
int i, j, k, kk, dim, sta, sto;
- real maxwgt = 0, len, len_total, len_total0;
+ double maxwgt = 0, len, len_total, len_total0;
int r, g, b;
- real tt1[3]={0.15,0.5,0.85};
- real tt2[4]={0.15,0.4,0.6,0.85};
- real *tt;
+ double tt1[3]={0.15,0.5,0.85};
+ double tt2[4]={0.15,0.4,0.6,0.85};
+ double *tt;
fprintf(fp,"strict graph{\n");
/* points */
void pedge_export_mma(FILE *fp, int ne, pedge *edges){
pedge edge;
- real *x;
+ double *x;
int i, j, k, dim;
fprintf(fp,"Graphics[{");
pedge pedge_realloc(pedge e, int n){
if (n <= e->npoints) return e;
- e->x = REALLOC(e->x, e->dim*n*sizeof(real));
- if (e->wgts) e->wgts = REALLOC(e->wgts, (n-1)*sizeof(real));
+ e->x = REALLOC(e->x, e->dim*n*sizeof(double));
+ if (e->wgts) e->wgts = REALLOC(e->wgts, (n-1)*sizeof(double));
e->len = n;
return e;
}
/* diff from pedge_alloc: allocate wgts if do not exist and initialize to wgt */
int i;
if (n <= e->npoints) return e;
- e->x = REALLOC(e->x, e->dim*n*sizeof(real));
+ e->x = REALLOC(e->x, e->dim*n*sizeof(double));
if (!(e->wgts)){
- e->wgts = REALLOC(e->wgts, (n-1)*sizeof(real));
+ e->wgts = REALLOC(e->wgts, (n-1)*sizeof(double));
for (i = 0; i < e->npoints; i++) e->wgts[i] = e->wgt;
} else {
- e->wgts = REALLOC(e->wgts, (n-1)*sizeof(real));
+ e->wgts = REALLOC(e->wgts, (n-1)*sizeof(double));
}
e->len = n;
return e;
pedge pedge_double(pedge e){
/* double the number of points (more precisely, add a point between two points in the polyline */
int npoints = e->npoints, len = e->len, i, dim = e->dim;
- real *x;
+ double *x;
int j, ii, ii2, np;
assert(npoints >= 2);
if (npoints*2-1 > len){
len = 3*npoints;
- e->x = REALLOC(e->x, dim*len*sizeof(real));
+ e->x = REALLOC(e->x, dim*len*sizeof(double));
}
x = e->x;
return e;
}
-static void edge_tension_force(real *force, pedge e){
- real *x = e->x;
+static void edge_tension_force(double *force, pedge e){
+ double *x = e->x;
int dim = e->dim;
int np = e->npoints;
int i, left, right, j;
- real s;
+ double s;
/* tention force = ((np-1)*||2x-xleft-xright||)/||e||, so the force is norminal and unitless
*/
}
}
-static void edge_attraction_force(real similarity, pedge e1, pedge e2, real *force){
+static void edge_attraction_force(double similarity, pedge e1, pedge e2, double *force){
/* attrractive force from x2 applied to x1 */
- real *x1 = e1->x, *x2 = e2->x;
+ double *x1 = e1->x, *x2 = e2->x;
int dim = e1->dim;
int np = e1->npoints;
int i, j;
- real dist, s, ss;
- real edge_length = e1->edge_length;
+ double dist, s, ss;
+ double edge_length = e1->edge_length;
assert(e1->npoints == e2->npoints);
}
-static pedge* force_directed_edge_bundling(SparseMatrix A, pedge* edges, int maxit, real step0, real K, int open_gl){
+static pedge* force_directed_edge_bundling(SparseMatrix A, pedge* edges, int maxit, double step0, double K, int open_gl){
int i, j, ne = A->n, k;
int *ia = A->ia, *ja = A->ja, iter = 0;
- real *a = (real*) A->a;
+ double *a = (double*) A->a;
pedge e1, e2;
- real *force_t, *force_a;
+ double *force_t, *force_a;
int np = edges[0]->npoints, dim = edges[0]->dim;
- real *x;
- real step = step0;
- real fnorm_a, fnorm_t, edge_length, start;
+ double *x;
+ double step = step0;
+ double fnorm_a, fnorm_t, edge_length, start;
if (Verbose > 1)
- fprintf(stderr, "total interaction pairs = %d out of %d, avg neighbors per edge = %f\n",A->nz, A->m*A->m, A->nz/(real) A->m);
+ fprintf(stderr, "total interaction pairs = %d out of %d, avg neighbors per edge = %f\n",A->nz, A->m*A->m, A->nz/(double) A->m);
- force_t = MALLOC(sizeof(real)*dim*np);
- force_a = MALLOC(sizeof(real)*dim*np);
+ force_t = MALLOC(sizeof(double)*dim*np);
+ force_a = MALLOC(sizeof(double)*dim*np);
while (step > 0.001 && iter < maxit){
start = clock();
iter++;
}
step = step*0.9;
if (Verbose > 1)
- fprintf(stderr, "iter ==== %d cpu = %f npoints = %d\n",iter, ((real) (clock() - start))/CLOCKS_PER_SEC, np - 2);
+ fprintf(stderr, "iter ==== %d cpu = %f npoints = %d\n",iter, ((double) (clock() - start))/CLOCKS_PER_SEC, np - 2);
#ifdef OPENGL
if (open_gl){
return edges;
}
-static pedge* modularity_ink_bundling(int dim, int ne, SparseMatrix B, pedge* edges, real angle_param, real angle){
+static pedge* modularity_ink_bundling(int dim, int ne, SparseMatrix B, pedge* edges, double angle_param, double angle){
int *assignment = NULL, flag, nclusters;
- real modularity;
+ double modularity;
int *clusterp, *clusters;
SparseMatrix D, C;
point_t meet1, meet2;
- real ink0, ink1;
+ double ink0, ink1;
pedge e;
int i, j, jj;
int use_value_for_clustering = TRUE;
return edges;
}
-static SparseMatrix check_compatibility(SparseMatrix A, int ne, pedge *edges, int compatibility_method, real tol){
+static SparseMatrix check_compatibility(SparseMatrix A, int ne, pedge *edges, int compatibility_method, double tol){
/* go through the links and make sure edges are compatible */
SparseMatrix B, C;
int *ia, *ja, i, j, jj;
- real start;
- real dist;
+ double start;
+ double dist;
B = SparseMatrix_new(1, 1, 1, MATRIX_TYPE_REAL, FORMAT_COORD);
ia = A->ia; ja = A->ja;
SparseMatrix_delete(B);
B = C;
if (Verbose > 1)
- fprintf(stderr, "edge compatibilitu time = %f\n",((real) (clock() - start))/CLOCKS_PER_SEC);
+ fprintf(stderr, "edge compatibilitu time = %f\n",((double) (clock() - start))/CLOCKS_PER_SEC);
return B;
}
-pedge* edge_bundling(SparseMatrix A0, int dim, real *x, int maxit_outer, real K, int method, int nneighbor, int compatibility_method,
- int max_recursion, real angle_param, real angle, int open_gl){
+pedge* edge_bundling(SparseMatrix A0, int dim, double *x, int maxit_outer, double K, int method, int nneighbor, int compatibility_method,
+ int max_recursion, double angle_param, double angle, int open_gl){
/* bundle edges.
A: edge graph
x: edge i is at {p,q},
pedge *edges;
SparseMatrix A = A0, B = NULL;
int i;
- real tol = 0.001;
+ double tol = 0.001;
int k;
- real step0 = 0.1, start = 0.0;
+ double step0 = 0.1, start = 0.0;
int maxit = 10;
int flag;
assert(0);
}
if (Verbose)
- fprintf(stderr, "total edge bundling cpu = %f\n",((real) (clock() - start))/CLOCKS_PER_SEC);
+ fprintf(stderr, "total edge bundling cpu = %f\n",((double) (clock() - start))/CLOCKS_PER_SEC);
if (B != A) SparseMatrix_delete(B);
if (A != A0) SparseMatrix_delete(A);
#include <sparse/SparseMatrix.h>
struct pedge_struct {
- real wgt; /* weight, telling how many original edges this edge represent. If this edge consists of multiple sections of different weights then this is a lower bound. This only applied for agglomerative bundling */
+ double wgt; /* weight, telling how many original edges this edge represent. If this edge consists of multiple sections of different weights then this is a lower bound. This only applied for agglomerative bundling */
int npoints;/* number of poly points */
int len;/* length of arra x. len >= npoints */
int dim;/* dim >= 2. Point i is stored from x[i*dim]*/
- real edge_length;
- real *x;/* coordinates of the npoints poly points. Dimension dim*npoints */
- real *wgts;/* number of original edges each section represnet. Dimension npoint - 1. This only applied for agglomerative bundling Null for other methods */
+ double edge_length;
+ double *x;/* coordinates of the npoints poly points. Dimension dim*npoints */
+ double *wgts;/* number of original edges each section represnet. Dimension npoint - 1. This only applied for agglomerative bundling Null for other methods */
};
typedef struct pedge_struct* pedge;
-pedge* edge_bundling(SparseMatrix A, int dim, real *x, int maxit_outer, real K, int method, int nneighbor, int compatibility_method, int max_recursion, real angle_param, real angle, int open_gl);
+pedge* edge_bundling(SparseMatrix A, int dim, double *x, int maxit_outer, double K, int method, int nneighbor, int compatibility_method, int max_recursion, double angle_param, double angle, int open_gl);
void pedge_delete(pedge e);
pedge pedge_realloc(pedge e, int np);
pedge pedge_wgts_realloc(pedge e, int n);
void pedge_export_gv(FILE *fp, int ne, pedge *edges);
enum {METHOD_NONE = -1, METHOD_FD, METHOD_INK_AGGLOMERATE, METHOD_INK};
enum {COMPATIBILITY_DIST = 0, COMPATIBILITY_FULL};
-pedge pedge_new(int np, int dim, real *x);
-pedge pedge_wgt_new(int np, int dim, real *x, real wgt);
+pedge pedge_new(int np, int dim, double *x);
+pedge pedge_wgt_new(int np, int dim, double *x, double wgt);
pedge pedge_double(pedge e);
/* flip the polyline so that last point becomes the first, second last the second, etc*/
/* sumLengths:
*/
-static double sumLengths_avoid_bad_angle(point_t* points, int npoints, point_t end, point_t meeting, real angle_param)
+static double sumLengths_avoid_bad_angle(point_t* points, int npoints, point_t end, point_t meeting, double angle_param)
{
/* avoid sharp turns, we want cos_theta to be as close to -1 as possible */
int i;
/* bestInk:
*/
-static double bestInk(point_t* points, int npoints, point_t begin, point_t end, double prec, point_t *meet, real angle_param)
+static double bestInk(point_t* points, int npoints, point_t begin, point_t end, double prec, point_t *meet, double angle_param)
{
point_t first, second, third, fourth, diff, meeting;
double value1, value2, value3, value4;
}
-static double project_to_line(point_t pt, point_t left, point_t right, real angle){
+static double project_to_line(point_t pt, point_t left, point_t right, double angle){
/* pt
^ ^
. \ \
point_t b, a;
- real bnorm, dnorm;
- real alpha, ccord;
+ double bnorm, dnorm;
+ double alpha, ccord;
if (angle <=0 || angle >= M_PI) return 2;/* return outside of the interval which should be handled as a sign of infeasible turning angle */
alpha = angle;
* Compute minimal ink used the input edges are bundled.
* Assumes tails all occur on one side and heads on the other.
*/
-double ink(pedge* edges, int numEdges, int *pick, double *ink0, point_t *meet1, point_t *meet2, real angle_param, real angle)
+double ink(pedge* edges, int numEdges, int *pick, double *ink0, point_t *meet1, point_t *meet2, double angle_param, double angle)
{
int i;
point_t begin, end, mid, diff;
pedge e;
- real *x;
+ double *x;
point_t* sources = N_NEW(numEdges, point_t);
point_t* targets = N_NEW(numEdges, point_t);
double inkUsed;
}
double ink1(pedge e){
- real *x, xx, yy;
+ double *x, xx, yy;
- real ink0 = 0;
+ double ink0 = 0;
x = e->x;
xx = x[0] - x[e->dim*e->npoints - e->dim];
meet1, meet2: meeting point
return: best ink needed if bundled.
*/
-double ink(pedge* edges, int numEdges, int *pick, double *ink0, point_t *meet1, point_t *meet2, real angle_param, real angle);
+double ink(pedge* edges, int numEdges, int *pick, double *ink0, point_t *meet1, point_t *meet2, double angle_param, double angle);
double ink1(pedge e);
extern double ink_count;
*/
int *irn = NULL, *jcn = NULL, nz;
- real *val = NULL;
+ double *val = NULL;
SparseMatrix A;
int k = num_neigbors;
#ifdef HAVE_ANN
nearest_neighbor_graph_ann(nPts, num_neigbors, eps, x, &nz, &irn, &jcn, &val);
- A = SparseMatrix_from_coordinate_arrays(nz, nPts, nPts, irn, jcn, val, MATRIX_TYPE_REAL, sizeof(real));
+ A = SparseMatrix_from_coordinate_arrays(nz, nPts, nPts, irn, jcn, val, MATRIX_TYPE_REAL, sizeof(double));
#else
A = NULL;
#endif
double *getSizes(Agraph_t * g, pointf pad, int* n_elabels, int** elabels)
{
Agnode_t *n;
- real *sizes = N_GNEW(Ndim * agnnodes(g), real);
+ double *sizes = N_GNEW(Ndim * agnnodes(g), double);
int i, nedge_nodes = 0;
int* elabs;
int i, row;
int *I;
int *J;
- real *val;
- real v;
+ double *val;
+ double v;
int type = MATRIX_TYPE_REAL;
Agsym_t* symD = NULL;
- real* valD = NULL;
+ double* valD = NULL;
if (!g)
return NULL;
I = N_GNEW(nedges, int);
J = N_GNEW(nedges, int);
- val = N_GNEW(nedges, real);
+ val = N_GNEW(nedges, double);
sym = agfindedgeattr(g, "weight");
if (D) {
symD = agfindedgeattr(g, "len");
- valD = N_NEW(nedges, real);
+ valD = N_NEW(nedges, double);
}
i = 0;
}
A = SparseMatrix_from_coordinate_arrays(nedges, nnodes, nnodes, I, J,
- val, type, sizeof(real));
+ val, type, sizeof(double));
- if (D) *D = SparseMatrix_from_coordinate_arrays(nedges, nnodes, nnodes, I, J, valD, type, sizeof(real));
+ if (D) *D = SparseMatrix_from_coordinate_arrays(nedges, nnodes, nnodes, I, J, valD, type, sizeof(double));
free(I);
free(J);
{
SparseMatrix A0 = makeMatrix(g, NULL);
SparseMatrix A = A0;
- real *sizes;
- real *pos = N_NEW(Ndim * agnnodes(g), real);
+ double *sizes;
+ double *pos = N_NEW(Ndim * agnnodes(g), double);
Agnode_t *n;
int flag = 0, i;
expand_t sep = sepFactor(g);
sizes = getSizes(g, pad, NULL, NULL);
for (n = agfstnode(g); n; n = agnxtnode(g, n)) {
- real* npos = pos + Ndim * ND_id(n);
+ double* npos = pos + Ndim * ND_id(n);
for (i = 0; i < Ndim; i++) {
npos[i] = ND_pos(n)[i];
}
ELSCHEME_NONE, 0, NULL, NULL, mapBool (agget(g, "overlap_shrink"), TRUE));
for (n = agfstnode(g); n; n = agnxtnode(g, n)) {
- real *npos = pos + Ndim * ND_id(n);
+ double *npos = pos + Ndim * ND_id(n);
for (i = 0; i < Ndim; i++) {
ND_pos(n)[i] = npos[i];
}
#include <neatogen/delaunay.h>
#include <stddef.h>
-SparseMatrix call_tri(int n, int dim, real * x)
+SparseMatrix call_tri(int n, int dim, double * x)
{
- real one = 1;
+ double one = 1;
int i, ii, jj;
SparseMatrix A;
SparseMatrix B;
int* edgelist = NULL;
- real* xv = N_GNEW(n, real);
- real* yv = N_GNEW(n, real);
+ double* xv = N_GNEW(n, double);
+ double* yv = N_GNEW(n, double);
int numberofedges = 0;
for (i = 0; i < n; i++) {
return B;
}
-SparseMatrix call_tri2(int n, int dim, real * xx)
+SparseMatrix call_tri2(int n, int dim, double * xx)
{
- real *x, *y;
+ double *x, *y;
v_data *delaunay;
int i, j;
SparseMatrix A;
SparseMatrix B;
- real one = 1;
- x = N_GNEW(n, real);
- y = N_GNEW(n, real);
+ double one = 1;
+ x = N_GNEW(n, double);
+ y = N_GNEW(n, double);
for (i = 0; i < n; i++) {
x[i] = xx[dim * i];
#pragma once
-SparseMatrix call_tri(int n, int dim, real * x);
-SparseMatrix call_tri2(int n, int dim, real * x);
+SparseMatrix call_tri(int n, int dim, double * x);
+SparseMatrix call_tri2(int n, int dim, double * x);
#include <common/globals.h>
#include <time.h>
-static void ideal_distance_avoid_overlap(int dim, SparseMatrix A, real *x, real *width, real *ideal_distance, real *tmax, real *tmin){
+static void ideal_distance_avoid_overlap(int dim, SparseMatrix A, double *x, double *width, double *ideal_distance, double *tmax, double *tmin){
/* if (x1>x2 && y1 > y2) we want either x1 + t (x1-x2) - x2 > (width1+width2), or y1 + t (y1-y2) - y2 > (height1+height2),
hence t = MAX(expandmin, MIN(expandmax, (width1+width2)/(x1-x2) - 1, (height1+height2)/(y1-y2) - 1)), and
new ideal distance = (1+t) old_distance. t can be negative sometimes.
*/
int i, j, jj;
int *ia = A->ia, *ja = A->ja;
- real dist, dx, dy, wx, wy, t;
- real expandmax = 1.5, expandmin = 1;
+ double dist, dx, dy, wx, wy, t;
+ double expandmax = 1.5, expandmin = 1;
*tmax = 0;
*tmin = 1.e10;
struct scan_point_struct{
int node;
- real x;
+ double x;
int status;
};
static void InfoDest(void *a){
}
-static SparseMatrix get_overlap_graph(int dim, int n, real *x, real *width, int check_overlap_only){
+static SparseMatrix get_overlap_graph(int dim, int n, double *x, double *width, int check_overlap_only){
/* if check_overlap_only = TRUE, we only check whether there is one overlap */
scan_point *scanpointsx, *scanpointsy;
int i, k, neighbor;
SparseMatrix A = NULL, B = NULL;
rb_red_blk_node *newNode, *newNode0, *newNode2 = NULL;
rb_red_blk_tree* treey;
- real one = 1;
+ double one = 1;
A = SparseMatrix_new(n, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD);
RBTreeInsert(treey, &(scanpointsy[k+n]), NULL);
} else {
- real bsta, bbsta, bsto, bbsto; int ii;
+ double bsta, bbsta, bsto, bbsto; int ii;
assert(scanpointsx[i].node >= n);
return data;
}
-static void scale_coord(int dim, int m, real *x, real scale){
+static void scale_coord(int dim, int m, double *x, double scale){
int i;
for (i = 0; i < dim*m; i++) {
x[i] *= scale;
}
}
-real overlap_scaling(int dim, int m, real *x, real *width, real scale_sta, real scale_sto, real epsilon, int maxiter){
+double overlap_scaling(int dim, int m, double *x, double *width, double scale_sta, double scale_sto, double epsilon, int maxiter){
/* do a bisection between scale_sta and scale_sto, up to maxiter iterations or till interval <= epsilon, to find the best scaling to avoid overlap
m: number of points
x: the coordinates
- for shrinking down a layout to reduce white space, we will assume scale_sta and scale_sto are both given and positive, and scale_sta is the current guess.
- for scaling up, we assume scale_sta, scale_sto <= 0
*/
- real scale = -1, scale_best = -1;
+ double scale = -1, scale_best = -1;
SparseMatrix C = NULL;
int check_overlap_only = 1;
int overlap = 0;
- real two = 2;
+ double two = 2;
int iter = 0;
assert(epsilon > 0);
}
OverlapSmoother OverlapSmoother_new(SparseMatrix A, int m,
- int dim, real lambda0, real *x, real *width, int include_original_graph, int neighborhood_only,
- real *max_overlap, real *min_overlap,
+ int dim, double lambda0, double *x, double *width, int include_original_graph, int neighborhood_only,
+ double *max_overlap, double *min_overlap,
int edge_labeling_scheme, int n_constr_nodes, int *constr_nodes, SparseMatrix A_constr, int shrink
){
OverlapSmoother sm;
int i, j, k, *iw, *jw, jdiag;
SparseMatrix B;
- real *lambda, *d, *w, diag_d, diag_w, dist;
+ double *lambda, *d, *w, diag_d, diag_w, dist;
assert((!A) || SparseMatrix_is_symmetric(A, FALSE));
sm->tol_cg = 0.01;
sm->maxit_cg = sqrt((double) A->m);
- lambda = sm->lambda = N_GNEW(m,real);
+ lambda = sm->lambda = N_GNEW(m,double);
for (i = 0; i < m; i++) sm->lambda[i] = lambda0;
B= call_tri(m, dim, x);
assert((sm->Lwd)->type == MATRIX_TYPE_REAL);
- ideal_distance_avoid_overlap(dim, sm->Lwd, x, width, (real*) (sm->Lwd->a), max_overlap, min_overlap);
+ ideal_distance_avoid_overlap(dim, sm->Lwd, x, width, (double*) (sm->Lwd->a), max_overlap, min_overlap);
/* no overlap at all! */
if (*max_overlap < 1 && shrink){
- real scale_sta = MIN(1, *max_overlap*1.0001), scale_sto = 1;
+ double scale_sta = MIN(1, *max_overlap*1.0001), scale_sto = 1;
if (Verbose) fprintf(stderr," no overlap (overlap = %f), rescale to shrink\n", *max_overlap - 1);
}
iw = sm->Lw->ia; jw = sm->Lw->ja;
- w = (real*) sm->Lw->a; d = (real*) sm->Lwd->a;
+ w = (double*) sm->Lw->a; d = (double*) sm->Lwd->a;
for (i = 0; i < m; i++){
diag_d = diag_w = 0;
}
-real OverlapSmoother_smooth(OverlapSmoother sm, int dim, real *x){
+double OverlapSmoother_smooth(OverlapSmoother sm, int dim, double *x){
int maxit_sm = 1;/* only using 1 iteration of stress majorization
is found to give better results and save time! */
- real res = StressMajorizationSmoother_smooth(sm, dim, x, maxit_sm, 0.001);
+ double res = StressMajorizationSmoother_smooth(sm, dim, x, maxit_sm, 0.001);
#ifdef DEBUG
{FILE *fp;
fp = fopen("/tmp/222","w");
/*================================= end OverlapSmoother =============*/
-static void scale_to_edge_length(int dim, SparseMatrix A, real *x, real avg_label_size){
- real dist;
+static void scale_to_edge_length(int dim, SparseMatrix A, double *x, double avg_label_size){
+ double dist;
int i;
if (!A) return;
for (i = 0; i < dim*A->m; i++) x[i] *= dist;
}
-static void print_bounding_box(int n, int dim, real *x){
- real *xmin, *xmax;
+static void print_bounding_box(int n, int dim, double *x){
+ double *xmin, *xmax;
int i, k;
- xmin = N_GNEW(dim,real);
- xmax = N_GNEW(dim,real);
+ xmin = N_GNEW(dim,double);
+ xmax = N_GNEW(dim,double);
for (i = 0; i < dim; i++) xmin[i]=xmax[i] = x[i];
free(xmax);
}
-static int check_convergence(real max_overlap, real res, int has_penalty_terms, real epsilon){
+static int check_convergence(double max_overlap, double res, int has_penalty_terms, double epsilon){
if (!has_penalty_terms) return (max_overlap <= 1);
return res < epsilon;
}
-void remove_overlap(int dim, SparseMatrix A, real *x, real *label_sizes, int ntry, real initial_scaling,
+void remove_overlap(int dim, SparseMatrix A, double *x, double *label_sizes, int ntry, double initial_scaling,
int edge_labeling_scheme, int n_constr_nodes, int *constr_nodes, SparseMatrix A_constr, int do_shrinking){
/*
edge_labeling_scheme: if ELSCHEME_NONE, n_constr_nodes/constr_nodes/A_constr are not used
*/
- real lambda = 0.00;
+ double lambda = 0.00;
OverlapSmoother sm;
int include_original_graph = 0, i;
- real LARGE = 100000;
- real avg_label_size, res = LARGE;
- real max_overlap = 0, min_overlap = 999;
+ double LARGE = 100000;
+ double avg_label_size, res = LARGE;
+ double max_overlap = 0, min_overlap = 999;
int neighborhood_only = TRUE;
int has_penalty_terms = FALSE;
- real epsilon = 0.005;
+ double epsilon = 0.005;
int shrink = 0;
#ifdef TIME
}
#endif
#ifdef TIME
- fprintf(stderr, "post processing %f\n",((real) (clock() - cpu)) / CLOCKS_PER_SEC);
+ fprintf(stderr, "post processing %f\n",((double) (clock() - cpu)) / CLOCKS_PER_SEC);
#endif
}
#else
#include <common/types.h>
#include <sparse/SparseMatrix.h>
-void remove_overlap(int dim, SparseMatrix A, real *x, real *label_sizes, int ntry, real initial_scaling,
+void remove_overlap(int dim, SparseMatrix A, double *x, double *label_sizes, int ntry, double initial_scaling,
int edge_labeling_scheme, int n_constr_nodes, int *constr_nodes, SparseMatrix A_constr, int do_shrinking)
{
static int once;
void OverlapSmoother_delete(OverlapSmoother sm);
OverlapSmoother OverlapSmoother_new(SparseMatrix A, int m,
- int dim, real lambda0, real *x, real *width, int include_original_graph, int neighborhood_only,
- real *max_overlap, real *min_overlap,
+ int dim, double lambda0, double *x, double *width, int include_original_graph, int neighborhood_only,
+ double *max_overlap, double *min_overlap,
int edge_labeling_scheme, int n_constr_nodes, int *constr_nodes, SparseMatrix A_constr, int shrink
);
enum {ELSCHEME_NONE = 0, ELSCHEME_PENALTY, ELSCHEME_PENALTY2, ELSCHEME_STRAIGHTLINE_PENALTY, ELSCHEME_STRAIGHTLINE_PENALTY2};
struct relative_position_constraints_struct{
- real constr_penalty; /* penalty parameter used in making edge labels as much on the line as possible */
+ double constr_penalty; /* penalty parameter used in making edge labels as much on the line as possible */
int edge_labeling_scheme;/* specifying whether to treat node of the form |edgelabel|* as a special node representing an edge label.
0 (no action, default), 1 (penalty based method to make that kind of node close to the center of its neighbor),
2 (penalty based method to make that kind of node close to the "old" center of its neighbor),
int *constr_nodes;/*constr_nodes: a list of nodes that need to be constrained. If NULL, unused.*/
int *irn;/* working arrays to hold the Laplacian of the constrain graph */
int *jcn;
- real *val;
+ double *val;
SparseMatrix A_constr; /*A_constr: neighbors of node i are in the row i of this matrix. i needs to sit
in between these neighbors as much as possible. this must not be NULL
if constr_nodes != NULL.*/
typedef struct relative_position_constraints_struct* relative_position_constraints;
-real OverlapSmoother_smooth(OverlapSmoother sm, int dim, real *x);
+double OverlapSmoother_smooth(OverlapSmoother sm, int dim, double *x);
-void remove_overlap(int dim, SparseMatrix A, real *x, real *label_sizes, int ntry, real initial_scaling,
+void remove_overlap(int dim, SparseMatrix A, double *x, double *label_sizes, int ntry, double initial_scaling,
int edge_labeling_scheme, int n_constr_nodes, int *constr_nodes, SparseMatrix A_constr, int doShrink);
-real overlap_scaling(int dim, int m, real *x, real *width, real scale_sta, real scale_sto, real epsilon, int maxiter);
+double overlap_scaling(int dim, int m, double *x, double *width, double scale_sta, double scale_sto, double epsilon, int maxiter);
#include <common/render.h>
#include <patchwork/tree_map.h>
-static void squarify(int n, real *area, rectangle *recs, int nadded, real maxarea, real minarea, real totalarea,
- real asp, rectangle fillrec){
+static void squarify(int n, double *area, rectangle *recs, int nadded, double maxarea, double minarea, double totalarea,
+ double asp, rectangle fillrec){
/* add a list of area in fillrec using squarified treemap alg.
n: number of items to add
area: area of these items, Sum to 1 (?).
asp: current worst aspect ratio of the already added items so far
fillrec: the rectangle to be filled in.
*/
- real w = MIN(fillrec.size[0], fillrec.size[1]);
+ double w = MIN(fillrec.size[0], fillrec.size[1]);
int i;
if (n <= 0) return;
totalarea = area[0];
squarify(n, area, recs, nadded, maxarea, minarea, totalarea, asp, fillrec);
} else {
- real newmaxarea, newminarea, s, h, maxw, minw, newasp, hh, ww, xx, yy;
+ double newmaxarea, newminarea, s, h, maxw, minw, newasp, hh, ww, xx, yy;
if (nadded < n){
newmaxarea = MAX(maxarea, area[nadded]);
newminarea = MIN(minarea, area[nadded]);
* fillred - rectangle to be filled
* return array of rectangles
*/
-rectangle* tree_map(int n, real *area, rectangle fillrec){
+rectangle* tree_map(int n, double *area, rectangle fillrec){
/* fill a rectangle rec with n items, each item i has area[i] area. */
rectangle *recs;
int i;
- real total = 0, minarea = 1., maxarea = 0., asp = 1, totalarea = 0;
+ double total = 0, minarea = 1., maxarea = 0., asp = 1, totalarea = 0;
int nadded = 0;
for (i = 0; i < n; i++) total += area[i];
/* rectangle_new:
* Create and initialize a new rectangle structure
*/
-rectangle rectangle_new(real x, real y, real width, real height){
+rectangle rectangle_new(double x, double y, double width, double height){
rectangle r;
r.x[0] = x;
r.x[1] = y;
#include <sparse/SparseMatrix.h>
typedef struct rectangle_struct {
- real x[2];/* center */
- real size[2]; /* total width/height*/
+ double x[2];/* center */
+ double size[2]; /* total width/height*/
} rectangle;
-extern rectangle* tree_map(int n, real *area, rectangle fillrec);
+extern rectangle* tree_map(int n, double *area, rectangle fillrec);
-extern rectangle rectangle_new(real x, real y, real width, real height);
+extern rectangle rectangle_new(double x, double y, double width, double height);
free(ctrl);
}
-static Multilevel Multilevel_init(SparseMatrix A, SparseMatrix D, real *node_weights){
+static Multilevel Multilevel_init(SparseMatrix A, SparseMatrix D, double *node_weights){
Multilevel grid;
if (!A) return NULL;
assert(A->m == A->n);
static void maximal_independent_edge_set_heavest_edge_pernode(SparseMatrix A, int randomize, int **matching, int *nmatch){
int i, ii, j, *ia, *ja, m, n, *p = NULL;
- real *a, amax = 0;
+ double *a, amax = 0;
int first = TRUE, jamax = 0;
assert(A);
assert(SparseMatrix_is_symmetric(A, FALSE));
assert(A->type == MATRIX_TYPE_REAL);
- a = (real*) A->a;
+ a = (double*) A->a;
if (!randomize){
for (i = 0; i < m; i++){
first = TRUE;
static void maximal_independent_edge_set_heavest_edge_pernode_leaves_first(SparseMatrix A, int randomize, int **cluster, int **clusterp, int *ncluster){
int i, ii, j, *ia, *ja, m, n, *p = NULL, q;
NOTUSED(n);
- real *a, amax = 0;
+ double *a, amax = 0;
int first = TRUE, jamax = 0;
int *matched, nz, ncmax = 0, nz0, nzz,k ;
enum {UNMATCHED = -2, MATCHED = -1};
*ncluster = 0;
(*clusterp)[0] = 0;
nz = 0;
- a = (real*) A->a;
+ a = (double*) A->a;
if (!randomize){
for (i = 0; i < m; i++){
if (matched[i] == MATCHED || node_degree(i) != 1) continue;
static void maximal_independent_edge_set_heavest_edge_pernode_supernodes_first(SparseMatrix A, int randomize, int **cluster, int **clusterp, int *ncluster){
int i, ii, j, *ia, *ja, m, n, *p = NULL;
NOTUSED(n);
- real *a, amax = 0;
+ double *a, amax = 0;
int first = TRUE, jamax = 0;
int *matched, nz, nz0;
enum {UNMATCHED = -2, MATCHED = -1};
*ncluster = 0;
(*clusterp)[0] = 0;
nz = 0;
- a = (real*) A->a;
+ a = (double*) A->a;
for (i = 0; i < nsuper; i++){
if (superp[i+1] - superp[i] <= 1) continue;
}
static int scomp(const void *s1, const void *s2){
- const real *ss1, *ss2;
- ss1 = (const real*) s1;
- ss2 = (const real*) s2;
+ const double *ss1, *ss2;
+ ss1 = (const double*) s1;
+ ss2 = (const double*) s2;
if ((ss1)[1] > (ss2)[1]){
return -1;
int randomize, int **cluster, int **clusterp, int *ncluster){
int i, ii, j, *ia, *ja, m, n, *p = NULL, q, iv;
NOTUSED(n);
- real *a;
+ double *a;
int *matched, nz, nz0, nzz,k, nv;
enum {UNMATCHED = -2, MATCHED = -1};
- real *vlist;
+ double *vlist;
assert(A);
assert(SparseMatrix_known_strucural_symmetric(A));
*cluster = N_GNEW(m,int);
*clusterp = N_GNEW((m+1),int);
matched = N_GNEW(m,int);
- vlist = N_GNEW(2*m,real);
+ vlist = N_GNEW(2*m,double);
for (i = 0; i < m; i++) matched[i] = i;
*ncluster = 0;
(*clusterp)[0] = 0;
nz = 0;
- a = (real*) A->a;
+ a = (double*) A->a;
p = random_permutation(m);
for (ii = 0; ii < m; ii++){
}
}
if (nv > 0){
- qsort(vlist, nv, sizeof(real)*2, scomp);
+ qsort(vlist, nv, sizeof(double)*2, scomp);
for (j = 0; j < MIN(csize - 1, nv); j++){
iv = (int) vlist[2*j];
matched[iv] = MATCHED;
}
static void maximal_independent_edge_set_heavest_edge_pernode_scaled(SparseMatrix A, int randomize, int **matching, int *nmatch){
int i, ii, j, *ia, *ja, m, n, *p = NULL;
- real *a, amax = 0;
+ double *a, amax = 0;
int first = TRUE, jamax = 0;
assert(A);
assert(SparseMatrix_is_symmetric(A, FALSE));
assert(A->type == MATRIX_TYPE_REAL);
- a = (real*) A->a;
+ a = (double*) A->a;
if (!randomize){
for (i = 0; i < m; i++){
first = TRUE;
}
static void Multilevel_coarsen_internal(SparseMatrix A, SparseMatrix *cA, SparseMatrix D, SparseMatrix *cD,
- real *node_wgt, real **cnode_wgt,
+ double *node_wgt, double **cnode_wgt,
SparseMatrix *P, SparseMatrix *R, Multilevel_control ctrl, int *coarsen_scheme_used){
int *matching = NULL, nmatch = 0, nc, nzc, n, i;
int *irn = NULL, *jcn = NULL, *ia = NULL, *ja = NULL;
- real *val = NULL;
+ double *val = NULL;
SparseMatrix B = NULL;
int *vset = NULL, nvset, ncov, j;
int *cluster=NULL, *clusterp=NULL, ncluster;
}
irn = N_GNEW(n,int);
jcn = N_GNEW(n,int);
- val = N_GNEW(n,real);
+ val = N_GNEW(n,double);
nzc = 0;
for (i = 0; i < ncluster; i++){
for (j = clusterp[i]; j < clusterp[i+1]; j++){
}
}
assert(nzc == n);
- *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn, (void *) val, MATRIX_TYPE_REAL, sizeof(real));
+ *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn, (void *) val, MATRIX_TYPE_REAL, sizeof(double));
*R = SparseMatrix_transpose(*P);
*cD = DistanceMatrix_restrict_cluster(ncluster, clusterp, cluster, *P, *R, D);
}
irn = N_GNEW(n,int);
jcn = N_GNEW(n,int);
- val = N_GNEW(n,real);
+ val = N_GNEW(n,double);
nzc = 0; nc = 0;
for (i = 0; i < n; i++){
if (matching[i] >= 0){
}
assert(nc == nmatch);
assert(nzc == n);
- *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn, (void *) val, MATRIX_TYPE_REAL, sizeof(real));
+ *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn, (void *) val, MATRIX_TYPE_REAL, sizeof(double));
*R = SparseMatrix_transpose(*P);
*cA = SparseMatrix_multiply3(*R, A, *P);
/*
}
irn = N_GNEW(nzc,int);
jcn = N_GNEW(nzc,int);
- val = N_GNEW(nzc,real);
+ val = N_GNEW(nzc,double);
nzc = 0;
for (i = 0; i < n; i++){
if (vset[i] == MAX_IND_VTX_SET_F){
}
}
- *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn, (void *) val, MATRIX_TYPE_REAL, sizeof(real));
+ *P = SparseMatrix_from_coordinate_arrays(nzc, n, nc, irn, jcn, (void *) val, MATRIX_TYPE_REAL, sizeof(double));
*R = SparseMatrix_transpose(*P);
*cA = SparseMatrix_multiply3(*R, A, *P);
if (!*cA) goto RETURN;
free(clusterp);
}
-void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, SparseMatrix D, SparseMatrix *cD, real *node_wgt, real **cnode_wgt,
+void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, SparseMatrix D, SparseMatrix *cD, double *node_wgt, double **cnode_wgt,
SparseMatrix *P, SparseMatrix *R, Multilevel_control ctrl, int *coarsen_scheme_used){
SparseMatrix cA0 = A, cD0 = NULL, P0 = NULL, R0 = NULL, M;
- real *cnode_wgt0 = NULL;
+ double *cnode_wgt0 = NULL;
int nc = 0, n;
*P = NULL; *R = NULL; *cA = NULL; *cnode_wgt = NULL, *cD = NULL;
static Multilevel Multilevel_establish(Multilevel grid, Multilevel_control ctrl){
Multilevel cgrid;
int coarsen_scheme_used;
- real *cnode_weights = NULL;
+ double *cnode_weights = NULL;
SparseMatrix P, R, A, cA, D, cD;
#ifdef DEBUG_PRINT
}
-Multilevel Multilevel_new(SparseMatrix A0, SparseMatrix D0, real *node_weights, Multilevel_control ctrl){
+Multilevel Multilevel_new(SparseMatrix A0, SparseMatrix D0, double *node_weights, Multilevel_control ctrl){
/* A: the weighting matrix. D: the distance matrix, could be NULL. If not null, the two matrices must have the same sparsity pattern */
Multilevel grid;
SparseMatrix A = A0, D = D0;
but different entry values. For spring-electrical method, D = NULL. */
SparseMatrix P;
SparseMatrix R;
- real *node_weights;
+ double *node_weights;
Multilevel next;
Multilevel prev;
int delete_top_level_A;
struct Multilevel_control_struct {
int minsize;
- real min_coarsen_factor;
+ double min_coarsen_factor;
int maxlevel;
int randomize;
int coarsen_scheme;
void Multilevel_delete(Multilevel grid);
-Multilevel Multilevel_new(SparseMatrix A, SparseMatrix D, real *node_weights, Multilevel_control ctrl);
+Multilevel Multilevel_new(SparseMatrix A, SparseMatrix D, double *node_weights, Multilevel_control ctrl);
Multilevel Multilevel_get_coarsest(Multilevel grid);
#define Multilevel_is_finest(grid) (!((grid)->prev))
#define Multilevel_is_coarsest(grid) (!((grid)->next))
-void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, SparseMatrix D, SparseMatrix *cD, real *node_wgt, real **cnode_wgt,
+void Multilevel_coarsen(SparseMatrix A, SparseMatrix *cA, SparseMatrix D, SparseMatrix *cD, double *node_wgt, double **cnode_wgt,
SparseMatrix *P, SparseMatrix *R, Multilevel_control ctrl, int *coarsen_scheme_used);
#define node_degree(i) (ia[(i)+1] - ia[(i)])
-static SparseMatrix ideal_distance_matrix(SparseMatrix A, int dim, real *x){
+static SparseMatrix ideal_distance_matrix(SparseMatrix A, int dim, double *x){
/* find the ideal distance between edges, either 1, or |N[i] \Union N[j]| - |N[i] \Intersection N[j]|
*/
SparseMatrix D;
int *ia, *ja, i, j, k, l, nz;
- real *d;
+ double *d;
int *mask = NULL;
- real len, di, sum, sumd;
+ double len, di, sum, sumd;
assert(SparseMatrix_is_symmetric(A, FALSE));
if (D->type != MATRIX_TYPE_REAL){
free(D->a);
D->type = MATRIX_TYPE_REAL;
- D->a = N_GNEW(D->nz,real);
+ D->a = N_GNEW(D->nz,double);
}
- d = (real*) D->a;
+ d = (double*) D->a;
mask = N_GNEW(D->m,int);
for (i = 0; i < D->m; i++) mask[i] = -1;
}
-StressMajorizationSmoother StressMajorizationSmoother2_new(SparseMatrix A, int dim, real lambda0, real *x,
+StressMajorizationSmoother StressMajorizationSmoother2_new(SparseMatrix A, int dim, double lambda0, double *x,
int ideal_dist_scheme){
/* use up to dist 2 neighbor */
/* use up to dist 2 neighbor. This is used in overcoming pherical effect with ideal distance of
StressMajorizationSmoother sm;
int i, j, k, l, m = A->m, *ia = A->ia, *ja = A->ja, *iw, *jw, *id, *jd;
int *mask, nz;
- real *d, *w, *lambda;
- real *avg_dist, diag_d, diag_w, dist, s = 0, stop = 0, sbot = 0;
+ double *d, *w, *lambda;
+ double *avg_dist, diag_d, diag_w, dist, s = 0, stop = 0, sbot = 0;
SparseMatrix ID;
assert(SparseMatrix_is_symmetric(A, FALSE));
sm->tol_cg = 0.01;
sm->maxit_cg = (int)sqrt((double) A->m);
- lambda = sm->lambda = N_GNEW(m,real);
+ lambda = sm->lambda = N_GNEW(m,double);
for (i = 0; i < m; i++) sm->lambda[i] = lambda0;
mask = N_GNEW(m,int);
- avg_dist = N_GNEW(m,real);
+ avg_dist = N_GNEW(m,double);
for (i = 0; i < m ;i++){
avg_dist[i] = 0;
iw = sm->Lw->ia; jw = sm->Lw->ja;
- w = (real*) sm->Lw->a; d = (real*) sm->Lwd->a;
+ w = (double*) sm->Lw->a; d = (double*) sm->Lwd->a;
id = sm->Lwd->ia; jd = sm->Lwd->ja;
iw[0] = id[0] = 0;
return sm;
}
-StressMajorizationSmoother SparseStressMajorizationSmoother_new(SparseMatrix A, int dim, real lambda0, real *x,
+StressMajorizationSmoother SparseStressMajorizationSmoother_new(SparseMatrix A, int dim, double lambda0, double *x,
int weighting_scheme, int scale_initial_coord){
/* solve a stress model to achieve the ideal distance among a sparse set of edges recorded in A.
A must be a real matrix.
StressMajorizationSmoother sm;
int i, j, k, m = A->m, *ia, *ja, *iw, *jw, *id, *jd;
int nz;
- real *d, *w, *lambda;
- real diag_d, diag_w, *a, dist, s = 0, stop = 0, sbot = 0;
- real xdot = 0;
+ double *d, *w, *lambda;
+ double diag_d, diag_w, *a, dist, s = 0, stop = 0, sbot = 0;
+ double xdot = 0;
assert(SparseMatrix_is_symmetric(A, FALSE) && A->type == MATRIX_TYPE_REAL);
ia = A->ia;
ja = A->ja;
- a = (real*) A->a;
+ a = (double*) A->a;
sm = MALLOC(sizeof(struct StressMajorizationSmoother_struct));
sm->tol_cg = 0.01;
sm->maxit_cg = (int)sqrt((double) A->m);
- lambda = sm->lambda = MALLOC(sizeof(real)*m);
+ lambda = sm->lambda = MALLOC(sizeof(double)*m);
for (i = 0; i < m; i++) sm->lambda[i] = lambda0;
nz = A->nz;
iw = sm->Lw->ia; jw = sm->Lw->ja;
id = sm->Lwd->ia; jd = sm->Lwd->ja;
- w = (real*) sm->Lw->a; d = (real*) sm->Lwd->a;
+ w = (double*) sm->Lw->a; d = (double*) sm->Lwd->a;
iw[0] = id[0] = 0;
nz = 0;
return sm;
}
-static real total_distance(int m, int dim, real* x, real* y){
- real total = 0, dist = 0;
+static double total_distance(int m, int dim, double* x, double* y){
+ double total = 0, dist = 0;
int i, j;
for (i = 0; i < m; i++){
}
-real SparseStressMajorizationSmoother_smooth(SparseStressMajorizationSmoother sm, int dim, real *x, int maxit_sm, real tol){
+double SparseStressMajorizationSmoother_smooth(SparseStressMajorizationSmoother sm, int dim, double *x, int maxit_sm, double tol){
return StressMajorizationSmoother_smooth(sm, dim, x, maxit_sm, tol);
}
-static void get_edge_label_matrix(relative_position_constraints data, int m, int dim, real *x, SparseMatrix *LL, real **rhs){
+static void get_edge_label_matrix(relative_position_constraints data, int m, int dim, double *x, SparseMatrix *LL, double **rhs){
int edge_labeling_scheme = data->edge_labeling_scheme;
int n_constr_nodes = data->n_constr_nodes;
int *constr_nodes = data->constr_nodes;
SparseMatrix A_constr = data->A_constr;
int *ia = A_constr->ia, *ja = A_constr->ja, ii, jj, nz, l, ll, i, j;
int *irn = data->irn, *jcn = data->jcn;
- real *val = data->val, dist, kk, k;
- real *x00 = NULL;
+ double *val = data->val, dist, kk, k;
+ double *x00 = NULL;
SparseMatrix Lc = NULL;
- real constr_penalty = data->constr_penalty;
+ double constr_penalty = data->constr_penalty;
if (edge_labeling_scheme == ELSCHEME_PENALTY || edge_labeling_scheme == ELSCHEME_STRAIGHTLINE_PENALTY){
/* for an node with two neighbors j--i--k, and assume i needs to be between j and k, then the contribution to P is
}
}
}
- Lc = SparseMatrix_from_coordinate_arrays(nz, m, m, irn, jcn, val, MATRIX_TYPE_REAL, sizeof(real));
+ Lc = SparseMatrix_from_coordinate_arrays(nz, m, m, irn, jcn, val, MATRIX_TYPE_REAL, sizeof(double));
} else if (edge_labeling_scheme == ELSCHEME_PENALTY2 || edge_labeling_scheme == ELSCHEME_STRAIGHTLINE_PENALTY2){
/* for an node with two neighbors j--i--k, and assume i needs to be between the old position of j and k, then the contribution to P is
1/d_jk, and to the right hand side: {0,...,average_position_of_i's neighbor if i is an edge node,...}
jcn = data->jcn = MALLOC(sizeof(int)*nz);
val = data->val = MALLOC(sizeof(double)*nz);
}
- x00 = MALLOC(sizeof(real)*m*dim);
+ x00 = MALLOC(sizeof(double)*m*dim);
for (i = 0; i < m*dim; i++) x00[i] = 0;
nz = 0;
for (i = 0; i < n_constr_nodes; i++){
x00[ii*dim+l] *= constr_penalty/(dist)/(ia[ii+1] - ia[ii]);
}
}
- Lc = SparseMatrix_from_coordinate_arrays(nz, m, m, irn, jcn, val, MATRIX_TYPE_REAL, sizeof(real));
+ Lc = SparseMatrix_from_coordinate_arrays(nz, m, m, irn, jcn, val, MATRIX_TYPE_REAL, sizeof(double));
}
*LL = Lc;
*rhs = x00;
}
-real get_stress(int m, int dim, int *iw, int *jw, real *w, real *d, real *x, real scaling, void *data, int weighted){
+double get_stress(int m, int dim, int *iw, int *jw, double *w, double *d, double *x, double scaling, void *data, int weighted){
int i, j;
- real res = 0., dist;
+ double res = 0., dist;
/* we use the fact that d_ij = w_ij*graph_dist(i,j). Also, d_ij and x are scalinged by *scaling, so divide by it to get actual unscaled streee. */
for (i = 0; i < m; i++){
for (j = iw[i]; j < iw[i+1]; j++){
}
-static void uniform_stress_augment_rhs(int m, int dim, real *x, real *y, real alpha, real M){
+static void uniform_stress_augment_rhs(int m, int dim, double *x, double *y, double alpha, double M){
int i, j, k;
- real dist, distij;
+ double dist, distij;
for (i = 0; i < m; i++){
for (j = i+1; j < m; j++){
dist = distance_cropped(x, dim, i, j);
}
}
-static real uniform_stress_solve(SparseMatrix Lw, real alpha, int dim, real *x0, real *rhs, real tol, int maxit){
+static double uniform_stress_solve(SparseMatrix Lw, double alpha, int dim, double *x0, double *rhs, double tol, int maxit){
Operator Ax;
Operator Precon;
}
-real StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, real *x, int maxit_sm, real tol) {
+double StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, double *x, int maxit_sm, double tol) {
SparseMatrix Lw = sm->Lw, Lwd = sm->Lwd, Lwdd = NULL;
int i, j, k, m, *id, *jd, *iw, *jw, idiag, flag = 0, iter = 0;
- real *w, *dd, *d, *y = NULL, *x0 = NULL, *x00 = NULL, diag, diff = 1, *lambda = sm->lambda, alpha = 0., M = 0.;
+ double *w, *dd, *d, *y = NULL, *x0 = NULL, *x00 = NULL, diag, diff = 1, *lambda = sm->lambda, alpha = 0., M = 0.;
SparseMatrix Lc = NULL;
- real dij, dist;
+ double dij, dist;
Lwdd = SparseMatrix_copy(Lwd);
m = Lw->m;
- x0 = calloc(dim*m, sizeof(real));
+ x0 = calloc(dim*m, sizeof(double));
if (!x0) goto RETURN;
- memcpy(x0, x, sizeof(real)*dim*m);
- y = calloc(dim*m, sizeof(real));
+ memcpy(x0, x, sizeof(double)*dim*m);
+ y = calloc(dim*m, sizeof(double));
if (!y) goto RETURN;
id = Lwd->ia; jd = Lwd->ja;
- d = (real*) Lwd->a;
- dd = (real*) Lwdd->a;
- w = (real*) Lw->a;
+ d = (double*) Lwd->a;
+ dd = (double*) Lwdd->a;
+ w = (double*) Lw->a;
iw = Lw->ia; jw = Lw->ja;
#ifdef DEBUG_PRINT
get_edge_label_matrix(sm->data, m, dim, x, &Lc, &x00);
if (Lc) Lw = SparseMatrix_add(Lw, Lc);
} else if (sm->scheme == SM_SCHEME_UNIFORM_STRESS){
- alpha = ((real*) (sm->data))[0];
- M = ((real*) (sm->data))[1];
+ alpha = ((double*) (sm->data))[0];
+ M = ((double*) (sm->data))[1];
}
while (iter++ < maxit_sm && diff > tol){
#endif
- memcpy(x, y, sizeof(real)*m*dim);
+ memcpy(x, y, sizeof(double)*m*dim);
}
#ifdef DEBUG
}
-TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, real lambda0, real *x, int use_triangularization){
+TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, double lambda0, double *x, int use_triangularization){
TriangleSmoother sm;
int i, j, k, m = A->m, *ia = A->ia, *ja = A->ja, *iw, *jw, jdiag, nz;
SparseMatrix B;
- real *avg_dist, *lambda, *d, *w, diag_d, diag_w, dist;
- real s = 0, stop = 0, sbot = 0;
+ double *avg_dist, *lambda, *d, *w, diag_d, diag_w, dist;
+ double s = 0, stop = 0, sbot = 0;
assert(SparseMatrix_is_symmetric(A, FALSE));
- avg_dist = N_GNEW(m,real);
+ avg_dist = N_GNEW(m,double);
for (i = 0; i < m ;i++){
avg_dist[i] = 0;
sm->tol_cg = 0.01;
sm->maxit_cg = (int)sqrt((double) A->m);
- lambda = sm->lambda = N_GNEW(m,real);
+ lambda = sm->lambda = N_GNEW(m,double);
for (i = 0; i < m; i++) sm->lambda[i] = lambda0;
if (m > 2){
iw = sm->Lw->ia; jw = sm->Lw->ja;
- w = (real*) sm->Lw->a; d = (real*) sm->Lwd->a;
+ w = (double*) sm->Lw->a; d = (double*) sm->Lwd->a;
for (i = 0; i < m; i++){
diag_d = diag_w = 0;
}
-void TriangleSmoother_smooth(TriangleSmoother sm, int dim, real *x){
+void TriangleSmoother_smooth(TriangleSmoother sm, int dim, double *x){
StressMajorizationSmoother_smooth(sm, dim, x, 50, 0.001);
}
/* ================================ spring and spring-electrical based smoother ================ */
-SpringSmoother SpringSmoother_new(SparseMatrix A, int dim, spring_electrical_control ctrl, real *x){
+SpringSmoother SpringSmoother_new(SparseMatrix A, int dim, spring_electrical_control ctrl, double *x){
SpringSmoother sm;
int i, j, k, l, m = A->m, *ia = A->ia, *ja = A->ja, *id, *jd;
int *mask, nz;
- real *d, *dd;
- real *avg_dist;
+ double *d, *dd;
+ double *avg_dist;
SparseMatrix ID = NULL;
assert(SparseMatrix_is_symmetric(A, FALSE));
ID = ideal_distance_matrix(A, dim, x);
- dd = (real*) ID->a;
+ dd = (double*) ID->a;
sm = N_GNEW(1,struct SpringSmoother_struct);
mask = N_GNEW(m,int);
- avg_dist = N_GNEW(m,real);
+ avg_dist = N_GNEW(m,double);
for (i = 0; i < m ;i++){
avg_dist[i] = 0;
}
id = sm->D->ia; jd = sm->D->ja;
- d = (real*) sm->D->a;
+ d = (double*) sm->D->a;
id[0] = 0;
nz = 0;
-void SpringSmoother_smooth(SpringSmoother sm, SparseMatrix A, real *node_weights, int dim, real *x){
+void SpringSmoother_smooth(SpringSmoother sm, SparseMatrix A, double *node_weights, int dim, double *x){
int flag = 0;
spring_electrical_spring_embedding(dim, A, sm->D, sm->ctrl, node_weights, x, &flag);
/*=============================== end of spring and spring-electrical based smoother =========== */
-void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control ctrl, real *node_weights, real *x, int *flag){
+void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control ctrl, double *node_weights, double *x, int *flag){
#ifdef TIME
clock_t cpu;
#endif
}/* end switch between smoothing methods */
#ifdef TIME
- if (Verbose) fprintf(stderr, "post processing %f\n",((real) (clock() - cpu)) / CLOCKS_PER_SEC);
+ if (Verbose) fprintf(stderr, "post processing %f\n",((double) (clock() - cpu)) / CLOCKS_PER_SEC);
#endif
}
SparseMatrix D;/* distance matrix. The diagonal is removed hence the ia, ja structure is different from Lw and Lwd!! */
SparseMatrix Lw;/* the weighted laplacian. with offdiag = -1/w_ij */
SparseMatrix Lwd;/* the laplacian like matrix with offdiag = -scaling*d_ij/w_ij. RHS in stress majorization = Lwd.x */
- real* lambda;
+ double* lambda;
void (*data_deallocator)(void*);
void *data;
int scheme;
- real scaling;/* scaling. It is multiplied to Lwd. need to divide coordinate x at the end of the stress majorization process */
- real tol_cg;/* tolerance and maxit for conjugate gradient that solves the Laplacian system.
+ double scaling;/* scaling. It is multiplied to Lwd. need to divide coordinate x at the end of the stress majorization process */
+ double tol_cg;/* tolerance and maxit for conjugate gradient that solves the Laplacian system.
typically the Laplacian only needs to be solved very crudely as it is part of an
outer iteration.*/
int maxit_cg;
void StressMajorizationSmoother_delete(StressMajorizationSmoother sm);
enum {IDEAL_GRAPH_DIST, IDEAL_AVG_DIST, IDEAL_POWER_DIST};
-StressMajorizationSmoother StressMajorizationSmoother2_new(SparseMatrix A, int dim, real lambda, real *x, int ideal_dist_scheme);
+StressMajorizationSmoother StressMajorizationSmoother2_new(SparseMatrix A, int dim, double lambda, double *x, int ideal_dist_scheme);
-real StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, real *x, int maxit, real tol);
+double StressMajorizationSmoother_smooth(StressMajorizationSmoother sm, int dim, double *x, int maxit, double tol);
/*-------------------- triangle/neirhborhood graph based smoother ------------------- */
typedef StressMajorizationSmoother TriangleSmoother;
void TriangleSmoother_delete(TriangleSmoother sm);
-TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, real lambda, real *x, int use_triangularization);
+TriangleSmoother TriangleSmoother_new(SparseMatrix A, int dim, double lambda, double *x, int use_triangularization);
-void TriangleSmoother_smooth(TriangleSmoother sm, int dim, real *x);
+void TriangleSmoother_smooth(TriangleSmoother sm, int dim, double *x);
typedef struct SpringSmoother_struct *SpringSmoother;
-SpringSmoother SpringSmoother_new(SparseMatrix A, int dim, spring_electrical_control ctrl, real *x);
+SpringSmoother SpringSmoother_new(SparseMatrix A, int dim, spring_electrical_control ctrl, double *x);
void SpringSmoother_delete(SpringSmoother sm);
-void SpringSmoother_smooth(SpringSmoother sm, SparseMatrix A, real *node_weights, int dim, real *x);
+void SpringSmoother_smooth(SpringSmoother sm, SparseMatrix A, double *node_weights, int dim, double *x);
/*------------------------------------------------------------------*/
-void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control ctrl, real *node_weights, real *x, int *flag);
+void post_process_smoothing(int dim, SparseMatrix A, spring_electrical_control ctrl, double *node_weights, double *x, int *flag);
/*-------------------- sparse stress majorizationp ------------------- */
typedef StressMajorizationSmoother SparseStressMajorizationSmoother;
void SparseStressMajorizationSmoother_delete(SparseStressMajorizationSmoother sm);
enum {WEIGHTING_SCHEME_NONE, WEIGHTING_SCHEME_INV_DIST, WEIGHTING_SCHEME_SQR_DIST};
-SparseStressMajorizationSmoother SparseStressMajorizationSmoother_new(SparseMatrix A, int dim, real lambda, real *x,
+SparseStressMajorizationSmoother SparseStressMajorizationSmoother_new(SparseMatrix A, int dim, double lambda, double *x,
int weighting_scheme, int scale_initial_coord);
-real SparseStressMajorizationSmoother_smooth(SparseStressMajorizationSmoother sm, int dim, real *x, int maxit_sm, real tol);
+double SparseStressMajorizationSmoother_smooth(SparseStressMajorizationSmoother sm, int dim, double *x, int maxit_sm, double tol);
-real get_stress(int m, int dim, int *iw, int *jw, real *w, real *d, real *x, real scaling, void *data, int weighted);
+double get_stress(int m, int dim, int *iw, int *jw, double *w, double *d, double *x, double scaling, void *data, int weighted);
/*--------------------------------------------------------------*/
/* getPos:
*/
-static real *getPos(Agraph_t * g)
+static double *getPos(Agraph_t * g)
{
Agnode_t *n;
- real *pos = N_NEW(Ndim * agnnodes(g), real);
+ double *pos = N_NEW(Ndim * agnnodes(g), double);
int ix, i;
if (agfindnodeattr(g, "pos") == NULL)
static void sfdpLayout(graph_t * g, spring_electrical_control ctrl,
int hops, pointf pad)
{
- real *sizes;
- real *pos;
+ double *sizes;
+ double *pos;
Agnode_t *n;
int flag, i;
int n_edge_label_nodes = 0, *edge_label_nodes = NULL;
break;
case METHOD_STRESS:{
int maxit = 200;
- real tol = 0.001;
+ double tol = 0.001;
int weighted = TRUE;
if (!D){
}
for (n = agfstnode(g); n; n = agnxtnode(g, n)) {
- real *npos = pos + (Ndim * ND_id(n));
+ double *npos = pos + (Ndim * ND_id(n));
for (i = 0; i < Ndim; i++) {
ND_pos(n)[i] = npos[i];
}
/* #define DEBUG_PRINT */
struct uniform_stress_matmul_data{
- real alpha;
+ double alpha;
SparseMatrix A;
};
-static real *Operator_uniform_stress_matmul_apply(Operator o, real *x, real *y){
+static double *Operator_uniform_stress_matmul_apply(Operator o, double *x, double *y){
struct uniform_stress_matmul_data *d = (struct uniform_stress_matmul_data*) (o->data);
SparseMatrix A = d->A;
- real alpha = d->alpha;
- real xsum = 0.;
+ double alpha = d->alpha;
+ double xsum = 0.;
int m = A->m, i;
SparseMatrix_multiply_vector(A, x, &y, FALSE);
-Operator Operator_uniform_stress_matmul(SparseMatrix A, real alpha){
+Operator Operator_uniform_stress_matmul(SparseMatrix A, double alpha){
Operator o;
struct uniform_stress_matmul_data *d;
}
-static real *Operator_matmul_apply(Operator o, real *x, real *y){
+static double *Operator_matmul_apply(Operator o, double *x, double *y){
SparseMatrix A = (SparseMatrix) o->data;
SparseMatrix_multiply_vector(A, x, &y, FALSE);
return y;
}
-static real* Operator_diag_precon_apply(Operator o, real *x, real *y){
+static double* Operator_diag_precon_apply(Operator o, double *x, double *y){
int i, m;
- real *diag = (real*) o->data;
+ double *diag = (double*) o->data;
m = (int) diag[0];
diag++;
for (i = 0; i < m; i++) y[i] = x[i]*diag[i];
}
-Operator Operator_uniform_stress_diag_precon_new(SparseMatrix A, real alpha){
+Operator Operator_uniform_stress_diag_precon_new(SparseMatrix A, double alpha){
Operator o;
- real *diag;
+ double *diag;
int i, j, m = A->m, *ia = A->ia, *ja = A->ja;
- real *a = (real*) A->a;
+ double *a = (double*) A->a;
assert(A->type == MATRIX_TYPE_REAL);
assert(a);
o = MALLOC(sizeof(struct Operator_struct));
- o->data = MALLOC(sizeof(real)*(m + 1));
- diag = (real*) o->data;
+ o->data = MALLOC(sizeof(double)*(m + 1));
+ diag = (double*) o->data;
diag[0] = m;
diag++;
static Operator Operator_diag_precon_new(SparseMatrix A){
Operator o;
- real *diag;
+ double *diag;
int i, j, m = A->m, *ia = A->ia, *ja = A->ja;
- real *a = (real*) A->a;
+ double *a = (double*) A->a;
assert(A->type == MATRIX_TYPE_REAL);
assert(a);
o = N_GNEW(1,struct Operator_struct);
- o->data = N_GNEW((A->m + 1),real);
- diag = (real*) o->data;
+ o->data = N_GNEW((A->m + 1),double);
+ diag = (double*) o->data;
diag[0] = m;
diag++;
free(o);
}
-static real conjugate_gradient(Operator A, Operator precon, int n, real *x, real *rhs, real tol, int maxit){
- real *z, *r, *p, *q, res = 10*tol, alpha;
- real rho = 1.0e20, rho_old = 1, res0, beta;
- real* (*Ax)(Operator o, real *in, real *out) = A->Operator_apply;
- real* (*Minvx)(Operator o, real *in, real *out) = precon->Operator_apply;
+static double conjugate_gradient(Operator A, Operator precon, int n, double *x, double *rhs, double tol, int maxit){
+ double *z, *r, *p, *q, res = 10*tol, alpha;
+ double rho = 1.0e20, rho_old = 1, res0, beta;
+ double* (*Ax)(Operator o, double *in, double *out) = A->Operator_apply;
+ double* (*Minvx)(Operator o, double *in, double *out) = precon->Operator_apply;
int iter = 0;
- z = N_GNEW(n,real);
- r = N_GNEW(n,real);
- p = N_GNEW(n,real);
- q = N_GNEW(n,real);
+ z = N_GNEW(n,double);
+ r = N_GNEW(n,double);
+ p = N_GNEW(n,double);
+ q = N_GNEW(n,double);
r = Ax(A, x, r);
r = vector_subtract_to(n, rhs, r);
beta = rho/rho_old;
p = vector_saxpy(n, z, p, beta);
} else {
- memcpy(p, z, sizeof(real)*n);
+ memcpy(p, z, sizeof(double)*n);
}
q = Ax(A, p, q);
return res;
}
-real cg(Operator Ax, Operator precond, int n, int dim, real *x0, real *rhs, real tol, int maxit){
- real *x, *b, res = 0;
+double cg(Operator Ax, Operator precond, int n, int dim, double *x0, double *rhs, double tol, int maxit){
+ double *x, *b, res = 0;
int k, i;
- x = N_GNEW(n, real);
- b = N_GNEW(n, real);
+ x = N_GNEW(n, double);
+ b = N_GNEW(n, double);
for (k = 0; k < dim; k++){
for (i = 0; i < n; i++) {
x[i] = x0[i*dim+k];
return res;
}
-static real* jacobi(SparseMatrix A, int dim, real *x0, real *rhs, int maxit, int *flag){
+static double* jacobi(SparseMatrix A, int dim, double *x0, double *rhs, int maxit, int *flag){
/* maxit iteration of jacobi */
- real *x, *y, *b, sum, diag, *a;
+ double *x, *y, *b, sum, diag, *a;
int k, i, j, n = A->n, *ia, *ja, iter;
- x = MALLOC(sizeof(real)*n);
- y = MALLOC(sizeof(real)*n);
- b = MALLOC(sizeof(real)*n);
+ x = MALLOC(sizeof(double)*n);
+ y = MALLOC(sizeof(double)*n);
+ b = MALLOC(sizeof(double)*n);
assert(A->type == MATRIX_TYPE_REAL);
- ia = A->ia; ja = A->ja; a = (real*) A->a;
+ ia = A->ia; ja = A->ja; a = (double*) A->a;
for (k = 0; k < dim; k++){
for (i = 0; i < n; i++) {
y[i] = (b[i] - sum)/diag;
}
- memcpy(x, y, sizeof(real)*n);
+ memcpy(x, y, sizeof(double)*n);
}
for (i = 0; i < n; i++) {
return rhs;
}
-real SparseMatrix_solve(SparseMatrix A, int dim, real *x0, real *rhs, real tol, int maxit, int method, int *flag){
+double SparseMatrix_solve(SparseMatrix A, int dim, double *x0, double *rhs, double tol, int maxit, int method, int *flag){
Operator Ax, precond;
int n = A->m;
- real res = 0;
+ double res = 0;
*flag = 0;
switch (method){
struct Operator_struct {
void *data;
- real* (*Operator_apply)(Operator o, real *in, real *out);
+ double* (*Operator_apply)(Operator o, double *in, double *out);
};
-real cg(Operator Ax, Operator precond, int n, int dim, real *x0, real *rhs, real tol, int maxit);
+double cg(Operator Ax, Operator precond, int n, int dim, double *x0, double *rhs, double tol, int maxit);
-real SparseMatrix_solve(SparseMatrix A, int dim, real *x0, real *rhs, real tol, int maxit, int method, int *flag);
+double SparseMatrix_solve(SparseMatrix A, int dim, double *x0, double *rhs, double tol, int maxit, int method, int *flag);
-Operator Operator_uniform_stress_matmul(SparseMatrix A, real alpha);
+Operator Operator_uniform_stress_matmul(SparseMatrix A, double alpha);
-Operator Operator_uniform_stress_diag_precon_new(SparseMatrix A, real alpha);
+Operator Operator_uniform_stress_diag_precon_new(SparseMatrix A, double alpha);
return opt;
}
-void oned_optimizer_train(oned_optimizer opt, real work){
+void oned_optimizer_train(oned_optimizer opt, double work){
int i = opt->i;
assert(i >= 0);
}
-real average_edge_length(SparseMatrix A, int dim, real *coord){
- real dist = 0, d;
+double average_edge_length(SparseMatrix A, int dim, double *coord){
+ double dist = 0, d;
int *ia = A->ia, *ja = A->ja, i, j, k;
assert(SparseMatrix_is_symmetric(A, TRUE));
}
#ifdef ENERGY
-static real spring_electrical_energy(int dim, SparseMatrix A, real *x, real p, real CRK, real KP){
+static double spring_electrical_energy(int dim, SparseMatrix A, double *x, double p, double CRK, double KP){
/* 1. Grad[||x-y||^k,x] = k||x-y||^(k-1)*0.5*(x-y)/||x-y|| = k/2*||x-y||^(k-2) (x-y)
which should equal to -force (force = -gradient),
hence energy for force ||x-y||^m (y-x) is ||x-y||^(m+2)*2/(m+2) where m != 2
*/
int i, j, k, *ia = A->ia, *ja = A->ja, n = A->m;
- real energy = 0, dist;
+ double energy = 0, dist;
for (i = 0; i < n; i++){
/* attractive force C^((2-p)/3) ||x_i-x_j||/K * (x_j - x_i) */
#endif
-void export_embedding(FILE *fp, int dim, SparseMatrix A, real *x, real *width){
+void export_embedding(FILE *fp, int dim, SparseMatrix A, double *x, double *width){
int i, j, k, *ia=A->ia, *ja = A->ja;
int ne = 0;
- real xsize, ysize, xmin, xmax, ymin, ymax;
+ double xsize, ysize, xmin, xmax, ymin, ymax;
xmax = xmin = x[0];
ymax = ymin = x[1];
}
-static real update_step(int adaptive_cooling, real step, real Fnorm, real Fnorm0, real cool){
+static double update_step(int adaptive_cooling, double step, double Fnorm, double Fnorm0, double cool){
if (!adaptive_cooling) {
return cool*step;
#define node_degree(i) (ia[(i)+1] - ia[(i)])
-static void check_real_array_size(real **a, int len, int *lenmax){
+static void check_real_array_size(double **a, int len, int *lenmax){
if (len >= *lenmax){
*lenmax = len + MAX((int) 0.2*len, 10);
- *a = REALLOC(*a, sizeof(real)*(*lenmax));
+ *a = REALLOC(*a, sizeof(double)*(*lenmax));
}
}
}
-static real get_angle(real *x, int dim, int i, int j){
+static double get_angle(double *x, int dim, int i, int j){
/* between [0, 2Pi)*/
int k;
- real y[2], res;
- real eps = 0.00001;
+ double y[2], res;
+ double eps = 0.00001;
for (k = 0; k < 2; k++){
y[k] = x[j*dim+k] - x[i*dim+k];
}
}
static int comp_real(const void *x, const void *y){
- const real *xx = (const real*) x;
- const real *yy = (const real*) y;
+ const double *xx = (const double*) x;
+ const double *yy = (const double*) y;
if (*xx > *yy){
return 1;
}
return 0;
}
-static void sort_real(int n, real *a){
- qsort(a, n, sizeof(real), comp_real);
+static void sort_real(int n, double *a){
+ qsort(a, n, sizeof(double), comp_real);
}
-static void set_leaves(real *x, int dim, real dist, real ang, int i, int j){
+static void set_leaves(double *x, int dim, double dist, double ang, int i, int j){
x[dim*j] = cos(ang)*dist + x[dim*i];
x[dim*j+1] = sin(ang)*dist + x[dim*i+1];
}
-static void beautify_leaves(int dim, SparseMatrix A, real *x){
+static void beautify_leaves(int dim, SparseMatrix A, double *x){
int m = A->m, i, j, *ia = A->ia, *ja = A->ja, k;
int p;
- real dist;
+ double dist;
int nleaves, nleaves_max = 10;
- real *angles, maxang, ang1 = 0, ang2 = 0, pad, step;
+ double *angles, maxang, ang1 = 0, ang2 = 0, pad, step;
int *leaves, nangles_max = 10, nangles;
assert(!SparseMatrix_has_diagonal(A));
bool *checked = gcalloc(sizeof(bool), m);
- angles = MALLOC(sizeof(real)*nangles_max);
+ angles = MALLOC(sizeof(double)*nangles_max);
leaves = MALLOC(sizeof(int)*nleaves_max);
free(leaves);
}
-void force_print(FILE *fp, int n, int dim, real *x, real *force){
+void force_print(FILE *fp, int n, int dim, double *x, double *force){
int i, k;
fprintf(fp,"Graphics[{");
}
-void spring_electrical_embedding_fast(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, int *flag){
+void spring_electrical_embedding_fast(int dim, SparseMatrix A0, spring_electrical_control ctrl, double *node_weights, double *x, int *flag){
/* x is a point to a 1D array, x[i*dim+j] gives the coordinate of the i-th node at dimension j. */
SparseMatrix A = A0;
int m, n;
int i, j, k;
- real p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP;
+ double p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP;
int *ia = NULL, *ja = NULL;
- real *xold = NULL;
- real *f = NULL, dist, F, Fnorm = 0, Fnorm0;
+ double *xold = NULL;
+ double *f = NULL, dist, F, Fnorm = 0, Fnorm0;
int iter = 0;
int adaptive_cooling = ctrl->adaptive_cooling;
QuadTree qt = NULL;
- real counts[4], *force = NULL;
+ double counts[4], *force = NULL;
#ifdef TIME
clock_t start, end, start0;
- real qtree_cpu = 0, qtree_cpu0 = 0, qtree_new_cpu = 0, qtree_new_cpu0 = 0;
- real total_cpu = 0;
+ double qtree_cpu = 0, qtree_cpu0 = 0, qtree_new_cpu = 0, qtree_new_cpu0 = 0;
+ double total_cpu = 0;
start0 = clock();
#endif
int max_qtree_level = ctrl->max_qtree_level;
KP = pow(K, 1 - p);
CRK = pow(C, (2.-p)/3.)/K;
- xold = MALLOC(sizeof(real)*dim*n);
- force = MALLOC(sizeof(real)*dim*n);
+ xold = MALLOC(sizeof(double)*dim*n);
+ force = MALLOC(sizeof(double)*dim*n);
do {
#ifdef TIME
#endif
iter++;
- memcpy(xold, x, sizeof(real)*dim*n);
+ memcpy(xold, x, sizeof(double)*dim*n);
Fnorm0 = Fnorm;
Fnorm = 0.;
}
#ifdef TIME
- qtree_new_cpu += ((real) (clock() - start))/CLOCKS_PER_SEC;
+ qtree_new_cpu += ((double) (clock() - start))/CLOCKS_PER_SEC;
#endif
/* repulsive force */
#ifdef TIME
end = clock();
- qtree_cpu += ((real) (end - start)) / CLOCKS_PER_SEC;
+ qtree_cpu += ((double) (end - start)) / CLOCKS_PER_SEC;
#endif
/* attractive force C^((2-p)/3) ||x_i-x_j||/K * (x_j - x_i) */
QuadTree_delete(qt);
#ifdef TIME
end = clock();
- qtree_new_cpu += ((real) (end - start)) / CLOCKS_PER_SEC;
+ qtree_new_cpu += ((double) (end - start)) / CLOCKS_PER_SEC;
#endif
#ifdef TIME
qtree_cpu0 = qtree_cpu - qtree_cpu0;
qtree_new_cpu0 = qtree_new_cpu - qtree_new_cpu0;
/* if (Verbose) fprintf(stderr, "\r iter=%d cpu=%.2f, quadtree=%.2f quad_force=%.2f other=%.2f counts={%.2f,%.2f,%.2f} step=%f Fnorm=%f nz=%d K=%f qtree_lev = %d",
- iter, ((real) (clock() - start2)) / CLOCKS_PER_SEC, qtree_new_cpu0,
- qtree_cpu0,((real) (clock() - start2))/CLOCKS_PER_SEC - qtree_cpu0 - qtree_new_cpu0,
+ iter, ((double) (clock() - start2)) / CLOCKS_PER_SEC, qtree_new_cpu0,
+ qtree_cpu0,((double) (clock() - start2))/CLOCKS_PER_SEC - qtree_cpu0 - qtree_new_cpu0,
counts[0], counts[1], counts[2],
step, Fnorm, A->nz,K,max_qtree_level);
*/
if (ctrl->beautify_leaves) beautify_leaves(dim, A, x);
#ifdef TIME
- total_cpu += ((real) (clock() - start0)) / CLOCKS_PER_SEC;
+ total_cpu += ((double) (clock() - start0)) / CLOCKS_PER_SEC;
if (Verbose) fprintf(stderr, "\n time for qtree = %f, qtree_force = %f, total cpu = %f\n",qtree_new_cpu, qtree_cpu, total_cpu);
#endif
free(force);
}
-static void spring_electrical_embedding_slow(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, int *flag){
+static void spring_electrical_embedding_slow(int dim, SparseMatrix A0, spring_electrical_control ctrl, double *node_weights, double *x, int *flag){
/* a version that does vertex moves in one go, instead of one at a time, use for debugging the fast version. Quadtree is not used. */
/* x is a point to a 1D array, x[i*dim+j] gives the coordinate of the i-th node at dimension j. */
SparseMatrix A = A0;
int m, n;
int i, j, k;
- real p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP;
+ double p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP;
int *ia = NULL, *ja = NULL;
- real *xold = NULL;
- real *f = NULL, dist, F, Fnorm = 0, Fnorm0;
+ double *xold = NULL;
+ double *f = NULL, dist, F, Fnorm = 0, Fnorm0;
int iter = 0;
int adaptive_cooling = ctrl->adaptive_cooling;
QuadTree qt = NULL;
int USE_QT = FALSE;
int nsuper = 0, nsupermax = 10;
- real *center = NULL, *supernode_wgts = NULL, *distances = NULL, nsuper_avg, counts = 0, counts_avg = 0;
- real *force;
+ double *center = NULL, *supernode_wgts = NULL, *distances = NULL, nsuper_avg, counts = 0, counts_avg = 0;
+ double *force;
#ifdef TIME
clock_t start, end, start0, start2;
- real qtree_cpu = 0, qtree_cpu0 = 0;
- real total_cpu = 0;
+ double qtree_cpu = 0, qtree_cpu0 = 0;
+ double total_cpu = 0;
start0 = clock();
#endif
int max_qtree_level = ctrl->max_qtree_level;
m = A->m, n = A->n;
if (n <= 0 || dim <= 0) return;
- force = MALLOC(sizeof(real)*n*dim);
+ force = MALLOC(sizeof(double)*n*dim);
if (n >= ctrl->quadtree_size) {
USE_QT = TRUE;
qtree_level_optimizer = oned_optimizer_new(max_qtree_level);
- center = MALLOC(sizeof(real)*nsupermax*dim);
- supernode_wgts = MALLOC(sizeof(real)*nsupermax);
- distances = MALLOC(sizeof(real)*nsupermax);
+ center = MALLOC(sizeof(double)*nsupermax*dim);
+ supernode_wgts = MALLOC(sizeof(double)*nsupermax);
+ distances = MALLOC(sizeof(double)*nsupermax);
}
USE_QT = FALSE;
*flag = 0;
}
#endif
- f = MALLOC(sizeof(real)*dim);
- xold = MALLOC(sizeof(real)*dim*n);
+ f = MALLOC(sizeof(double)*dim);
+ xold = MALLOC(sizeof(double)*dim*n);
do {
for (i = 0; i < dim*n; i++) force[i] = 0;
iter++;
- memcpy(xold, x, sizeof(real)*dim*n);
+ memcpy(xold, x, sizeof(double)*dim*n);
Fnorm0 = Fnorm;
Fnorm = 0.;
nsuper_avg = 0;
¢er, &supernode_wgts, &distances, &counts, flag);
#ifdef TIME
end = clock();
- qtree_cpu += ((real) (end - start)) / CLOCKS_PER_SEC;
+ qtree_cpu += ((double) (end - start)) / CLOCKS_PER_SEC;
#endif
counts_avg += counts;
nsuper_avg += nsuper;
counts_avg /= n;
#ifdef TIME
qtree_cpu0 = qtree_cpu - qtree_cpu0;
- if (Verbose && 0) fprintf(stderr, "\n cpu this outer iter = %f, quadtree time = %f other time = %f\n",((real) (clock() - start2)) / CLOCKS_PER_SEC, qtree_cpu0,((real) (clock() - start2))/CLOCKS_PER_SEC - qtree_cpu0);
+ if (Verbose && 0) fprintf(stderr, "\n cpu this outer iter = %f, quadtree time = %f other time = %f\n",((double) (clock() - start2)) / CLOCKS_PER_SEC, qtree_cpu0,((double) (clock() - start2))/CLOCKS_PER_SEC - qtree_cpu0);
qtree_cpu0 = qtree_cpu;
#endif
if (Verbose && 0) fprintf(stderr, "nsuper_avg=%f, counts_avg = %f 2*nsuper+counts=%f\n",nsuper_avg,counts_avg, 2*nsuper_avg+counts_avg);
if (ctrl->beautify_leaves) beautify_leaves(dim, A, x);
#ifdef TIME
- total_cpu += ((real) (clock() - start0)) / CLOCKS_PER_SEC;
+ total_cpu += ((double) (clock() - start0)) / CLOCKS_PER_SEC;
if (Verbose) fprintf(stderr, "time for supernode = %f, total cpu = %f\n",qtree_cpu, total_cpu);
#endif
-void spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, int *flag){
+void spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_control ctrl, double *node_weights, double *x, int *flag){
/* x is a point to a 1D array, x[i*dim+j] gives the coordinate of the i-th node at dimension j. */
SparseMatrix A = A0;
int m, n;
int i, j, k;
- real p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP;
+ double p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP;
int *ia = NULL, *ja = NULL;
- real *xold = NULL;
- real *f = NULL, dist, F, Fnorm = 0, Fnorm0;
+ double *xold = NULL;
+ double *f = NULL, dist, F, Fnorm = 0, Fnorm0;
int iter = 0;
int adaptive_cooling = ctrl->adaptive_cooling;
QuadTree qt = NULL;
int USE_QT = FALSE;
int nsuper = 0, nsupermax = 10;
- real *center = NULL, *supernode_wgts = NULL, *distances = NULL, nsuper_avg, counts = 0, counts_avg = 0;
+ double *center = NULL, *supernode_wgts = NULL, *distances = NULL, nsuper_avg, counts = 0, counts_avg = 0;
#ifdef TIME
clock_t start, end, start0, start2;
- real qtree_cpu = 0, qtree_cpu0 = 0;
- real total_cpu = 0;
+ double qtree_cpu = 0, qtree_cpu0 = 0;
+ double total_cpu = 0;
start0 = clock();
#endif
int max_qtree_level = ctrl->max_qtree_level;
if (n >= ctrl->quadtree_size) {
USE_QT = TRUE;
qtree_level_optimizer = oned_optimizer_new(max_qtree_level);
- center = MALLOC(sizeof(real)*nsupermax*dim);
- supernode_wgts = MALLOC(sizeof(real)*nsupermax);
- distances = MALLOC(sizeof(real)*nsupermax);
+ center = MALLOC(sizeof(double)*nsupermax*dim);
+ supernode_wgts = MALLOC(sizeof(double)*nsupermax);
+ distances = MALLOC(sizeof(double)*nsupermax);
}
*flag = 0;
if (m != n) {
}
#endif
- f = MALLOC(sizeof(real)*dim);
- xold = MALLOC(sizeof(real)*dim*n);
+ f = MALLOC(sizeof(double)*dim);
+ xold = MALLOC(sizeof(double)*dim*n);
do {
//#define VIS_MULTILEVEL
#endif
iter++;
- memcpy(xold, x, sizeof(real)*dim*n);
+ memcpy(xold, x, sizeof(double)*dim*n);
Fnorm0 = Fnorm;
Fnorm = 0.;
nsuper_avg = 0;
#ifdef TIME
end = clock();
- qtree_cpu += ((real) (end - start)) / CLOCKS_PER_SEC;
+ qtree_cpu += ((double) (end - start)) / CLOCKS_PER_SEC;
#endif
counts_avg += counts;
nsuper_avg += nsuper;
counts_avg /= n;
#ifdef TIME
qtree_cpu0 = qtree_cpu - qtree_cpu0;
- if (Verbose && 0) fprintf(stderr, "\n cpu this outer iter = %f, quadtree time = %f other time = %f\n",((real) (clock() - start2)) / CLOCKS_PER_SEC, qtree_cpu0,((real) (clock() - start2))/CLOCKS_PER_SEC - qtree_cpu0);
+ if (Verbose && 0) fprintf(stderr, "\n cpu this outer iter = %f, quadtree time = %f other time = %f\n",((double) (clock() - start2)) / CLOCKS_PER_SEC, qtree_cpu0,((double) (clock() - start2))/CLOCKS_PER_SEC - qtree_cpu0);
qtree_cpu0 = qtree_cpu;
#endif
if (Verbose & 0) fprintf(stderr, "nsuper_avg=%f, counts_avg = %f 2*nsuper+counts=%f\n",nsuper_avg,counts_avg, 2*nsuper_avg+counts_avg);
if (ctrl->beautify_leaves) beautify_leaves(dim, A, x);
#ifdef TIME
- total_cpu += ((real) (clock() - start0)) / CLOCKS_PER_SEC;
+ total_cpu += ((double) (clock() - start0)) / CLOCKS_PER_SEC;
if (Verbose) fprintf(stderr, "time for supernode = %f, total cpu = %f\n",qtree_cpu, total_cpu);
#endif
free(distances);
}
-static void scale_coord(int n, int dim, real *x, int *id, int *jd, real *d, real dj){
+static void scale_coord(int n, int dim, double *x, int *id, int *jd, double *d, double dj){
int i, j, k;
- real w_ij, dist, s = 0, stop = 0, sbot = 0., nz = 0;
+ double w_ij, dist, s = 0, stop = 0, sbot = 0., nz = 0;
if (dj == 0.) return;
for (i = 0; i < n; i++){
fprintf(stderr,"scaling factor = %f\n",s);
}
-static real dmean_get(int n, int *id, real* d){
- real dmean = 0;
+static double dmean_get(int n, int *id, double* d){
+ double dmean = 0;
int i, j;
if (!d) return 1.;
dmean += d[j];
}
}
- return dmean/((real) id[n]);
+ return dmean/((double) id[n]);
}
-static void spring_maxent_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *x, real rho, int *flag){
+static void spring_maxent_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control ctrl, double *node_weights, double *x, double rho, int *flag){
/* x is a point to a 1D array, x[i*dim+j] gives the coordinate of the i-th node at dimension j.
Minimize \Sum_{(i,j)\in E} w_ij (||x_i-x_j||-d_ij)^2 - \rho \Sum_{(i,j)\NotIn E} Log ||x_i-x_j||
SparseMatrix A = A0;
int m, n;
int i, j, k;
- real p = ctrl->p, C = ctrl->C, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, w_ij, dj = 1.;
+ double p = ctrl->p, C = ctrl->C, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, w_ij, dj = 1.;
int *ia = NULL, *ja = NULL;
int *id = NULL, *jd = NULL;
- real *d, dmean;
- real *xold = NULL;
- real *f = NULL, dist, F, Fnorm = 0, Fnorm0;
+ double *d, dmean;
+ double *xold = NULL;
+ double *f = NULL, dist, F, Fnorm = 0, Fnorm0;
int iter = 0;
int adaptive_cooling = ctrl->adaptive_cooling;
QuadTree qt = NULL;
int USE_QT = FALSE;
int nsuper = 0, nsupermax = 10;
- real *center = NULL, *supernode_wgts = NULL, *distances = NULL, nsuper_avg, counts = 0;
+ double *center = NULL, *supernode_wgts = NULL, *distances = NULL, nsuper_avg, counts = 0;
int max_qtree_level = 10;
#ifdef DEBUG
double stress = 0;
if (ctrl->tscheme != QUAD_TREE_NONE && n >= ctrl->quadtree_size) {
USE_QT = TRUE;
- center = MALLOC(sizeof(real)*nsupermax*dim);
- supernode_wgts = MALLOC(sizeof(real)*nsupermax);
- distances = MALLOC(sizeof(real)*nsupermax);
+ center = MALLOC(sizeof(double)*nsupermax*dim);
+ supernode_wgts = MALLOC(sizeof(double)*nsupermax);
+ distances = MALLOC(sizeof(double)*nsupermax);
}
*flag = 0;
if (D){
id = D->ia;
jd = D->ja;
- d = (real*) D->a;
+ d = (double*) D->a;
} else {
id = ia; jd = ja; d = NULL;
}
if (rho < 0) {
dmean = dmean_get(n, id, d);
- rho = rho*(id[n]/((((real) n)*((real) n)) - id[n]))/pow(dmean, p+1);
+ rho = rho*(id[n]/((((double) n)*((double) n)) - id[n]))/pow(dmean, p+1);
fprintf(stderr,"dmean = %f, rho = %f\n",dmean, rho);
}
}
#endif
- f = MALLOC(sizeof(real)*dim);
- xold = MALLOC(sizeof(real)*dim*n);
+ f = MALLOC(sizeof(double)*dim);
+ xold = MALLOC(sizeof(double)*dim*n);
do {
iter++;
- memcpy(xold, x, sizeof(real)*dim*n);
+ memcpy(xold, x, sizeof(double)*dim*n);
Fnorm0 = Fnorm;
Fnorm = 0.;
nsuper_avg = 0;
-void spring_electrical_spring_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *x, int *flag){
+void spring_electrical_spring_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control ctrl, double *node_weights, double *x, int *flag){
/* x is a point to a 1D array, x[i*dim+j] gives the coordinate of the i-th node at dimension j. Same as the spring-electrical except we also
introduce force due to spring length
*/
SparseMatrix A = A0;
int m, n;
int i, j, k;
- real p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP;
+ double p = ctrl->p, K = ctrl->K, C = ctrl->C, CRK, tol = ctrl->tol, maxiter = ctrl->maxiter, cool = ctrl->cool, step = ctrl->step, KP;
int *ia = NULL, *ja = NULL;
int *id = NULL, *jd = NULL;
- real *d;
- real *xold = NULL;
- real *f = NULL, dist, F, Fnorm = 0, Fnorm0;
+ double *d;
+ double *xold = NULL;
+ double *f = NULL, dist, F, Fnorm = 0, Fnorm0;
int iter = 0;
int adaptive_cooling = ctrl->adaptive_cooling;
QuadTree qt = NULL;
int USE_QT = FALSE;
int nsuper = 0, nsupermax = 10;
- real *center = NULL, *supernode_wgts = NULL, *distances = NULL, nsuper_avg, counts = 0;
+ double *center = NULL, *supernode_wgts = NULL, *distances = NULL, nsuper_avg, counts = 0;
int max_qtree_level = 10;
if (!A || maxiter <= 0) return;
if (n >= ctrl->quadtree_size) {
USE_QT = TRUE;
- center = MALLOC(sizeof(real)*nsupermax*dim);
- supernode_wgts = MALLOC(sizeof(real)*nsupermax);
- distances = MALLOC(sizeof(real)*nsupermax);
+ center = MALLOC(sizeof(double)*nsupermax*dim);
+ supernode_wgts = MALLOC(sizeof(double)*nsupermax);
+ distances = MALLOC(sizeof(double)*nsupermax);
}
*flag = 0;
if (m != n) {
ja = A->ja;
id = D->ia;
jd = D->ja;
- d = (real*) D->a;
+ d = (double*) D->a;
if (ctrl->random_start){
srand(ctrl->random_seed);
}
#endif
- f = MALLOC(sizeof(real)*dim);
- xold = MALLOC(sizeof(real)*dim*n);
+ f = MALLOC(sizeof(double)*dim);
+ xold = MALLOC(sizeof(double)*dim*n);
do {
iter++;
- memcpy(xold, x, sizeof(real)*dim*n);
+ memcpy(xold, x, sizeof(double)*dim*n);
Fnorm0 = Fnorm;
Fnorm = 0.;
nsuper_avg = 0;
-void print_matrix(real *x, int n, int dim){
+void print_matrix(double *x, int n, int dim){
int i, k;
printf("{");
for (i = 0; i < n; i++){
printf("}\n");
}
-void interpolate_coord(int dim, SparseMatrix A, real *x){
+void interpolate_coord(int dim, SparseMatrix A, double *x){
int i, j, k, *ia = A->ia, *ja = A->ja, nz;
- real alpha = 0.5, beta, *y;
+ double alpha = 0.5, beta, *y;
- y = MALLOC(sizeof(real)*dim);
+ y = MALLOC(sizeof(double)*dim);
for (i = 0; i < A->m; i++){
for (k = 0; k < dim; k++) y[k] = 0;
nz = 0;
free(y);
}
-static void prolongate(int dim, SparseMatrix A, SparseMatrix P, SparseMatrix R, real *x, real *y, int coarsen_scheme_used, real delta){
+static void prolongate(int dim, SparseMatrix A, SparseMatrix P, SparseMatrix R, double *x, double *y, int coarsen_scheme_used, double delta){
int nc, *ia, *ja, i, j, k;
SparseMatrix_multiply_dense(P, x, &y, dim);
return res;
}
-void pcp_rotate(int n, int dim, real *x){
+void pcp_rotate(int n, int dim, double *x){
int i, k,l;
- real y[4], axis[2], center[2], dist, x0, x1;
+ double y[4], axis[2], center[2], dist, x0, x1;
assert(dim == 2);
for (i = 0; i < dim*dim; i++) y[i] = 0;
}
-static void rotate(int n, int dim, real *x, real angle){
+static void rotate(int n, int dim, double *x, double angle){
int i, k;
- real axis[2], center[2], x0, x1;
- real radian = 3.14159/180;
+ double axis[2], center[2], x0, x1;
+ double radian = 3.14159/180;
assert(dim == 2);
for (i = 0; i < dim; i++) center[i] = 0;
}
-static void attach_edge_label_coordinates(int dim, SparseMatrix A, int n_edge_label_nodes, int *edge_label_nodes, real *x, real *x2){
+static void attach_edge_label_coordinates(int dim, SparseMatrix A, int n_edge_label_nodes, int *edge_label_nodes, double *x, double *x2){
int *mask;
int i, ii, j, k;
int nnodes = 0;
- real len;
+ double len;
mask = MALLOC(sizeof(int)*A->m);
}
}
- B = SparseMatrix_from_coordinate_arrays(nz, id, id, irn, jcn, NULL, MATRIX_TYPE_PATTERN, sizeof(real));
+ B = SparseMatrix_from_coordinate_arrays(nz, id, id, irn, jcn, NULL, MATRIX_TYPE_PATTERN, sizeof(double));
free(irn);
free(jcn);
}
-static void multilevel_spring_electrical_embedding_core(int dim, SparseMatrix A0, SparseMatrix D0, spring_electrical_control ctrl, real *node_weights, real *label_sizes,
- real *x, int n_edge_label_nodes, int *edge_label_nodes, int *flag){
+static void multilevel_spring_electrical_embedding_core(int dim, SparseMatrix A0, SparseMatrix D0, spring_electrical_control ctrl, double *node_weights, double *label_sizes,
+ double *x, int n_edge_label_nodes, int *edge_label_nodes, int *flag){
Multilevel_control mctrl = NULL;
int n, plg, coarsen_scheme_used;
SparseMatrix A = A0, D = D0, P = NULL;
Multilevel grid, grid0;
- real *xc = NULL, *xf = NULL;
+ double *xc = NULL, *xf = NULL;
struct spring_electrical_control_struct ctrl0;
#ifdef TIME
clock_t cpu;
&& n_edge_label_nodes > 0){
SparseMatrix A2;
- real *x2 = MALLOC(sizeof(real)*(A->m)*dim);
+ double *x2 = MALLOC(sizeof(double)*(A->m)*dim);
A2 = shorting_edge_label_nodes(A, n_edge_label_nodes, edge_label_nodes);
multilevel_spring_electrical_embedding(dim, A2, NULL, ctrl, NULL, NULL, x2, 0, NULL, flag);
if (Multilevel_is_finest(grid)){
xc = x;
} else {
- xc = MALLOC(sizeof(real)*grid->n*dim);
+ xc = MALLOC(sizeof(double)*grid->n*dim);
}
plg = power_law_graph(A);
if (Multilevel_is_finest(grid)){
xf = x;
} else {
- xf = MALLOC(sizeof(real)*grid->n*dim);
+ xf = MALLOC(sizeof(double)*grid->n*dim);
}
prolongate(dim, grid->A, P, grid->R, xc, xf, coarsen_scheme_used, (ctrl->K)*0.001);
free(xc);
#ifdef TIME
if (Verbose)
- fprintf(stderr, "layout time %f\n",((real) (clock() - cpu)) / CLOCKS_PER_SEC);
+ fprintf(stderr, "layout time %f\n",((double) (clock() - cpu)) / CLOCKS_PER_SEC);
cpu = clock();
#endif
Multilevel_delete(grid0);
}
-void multilevel_spring_electrical_embedding(int dim, SparseMatrix A, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *label_sizes,
- real *x, int n_edge_label_nodes, int *edge_label_nodes, int *flag){
+void multilevel_spring_electrical_embedding(int dim, SparseMatrix A, SparseMatrix D, spring_electrical_control ctrl, double *node_weights, double *label_sizes,
+ double *x, int n_edge_label_nodes, int *edge_label_nodes, int *flag){
multilevel_spring_electrical_embedding_core(dim, A, D, ctrl, node_weights, label_sizes, x, n_edge_label_nodes, edge_label_nodes, flag);
}
enum {METHOD_STA = -1, METHOD_SPRING_ELECTRICAL, METHOD_SPRING_MAXENT, METHOD_STRESS_MAXENT, METHOD_STRESS_APPROX, METHOD_STRESS, METHOD_UNIFORM_STRESS, METHOD_FULL_STRESS, METHOD_NONE, METHOD_STO};
struct spring_electrical_control_struct {
- real p;/*a negativve real number default to -1. repulsive force = dist^p */
- real q;/*a positive real number default to 2. attractive force = dist^q */
+ double p;/*a negativve real number default to -1. repulsive force = dist^p */
+ double q;/*a positive real number default to 2. attractive force = dist^q */
int random_start;/* whether to apply SE from a random layout, or from exisiting layout */
- real K;/* the natural distance. If K < 0, K will be set to the average distance of an edge */
- real C;/* another parameter. f_a(i,j) = C*dist(i,j)^2/K * d_ij, f_r(i,j) = K^(3-p)/dist(i,j)^(-p). By default C = 0.2. */
+ double K;/* the natural distance. If K < 0, K will be set to the average distance of an edge */
+ double C;/* another parameter. f_a(i,j) = C*dist(i,j)^2/K * d_ij, f_r(i,j) = K^(3-p)/dist(i,j)^(-p). By default C = 0.2. */
int multilevels;/* if <=1, single level */
int multilevel_coarsen_scheme;/* pass on to Multilevel_control->coarsen_scheme */
int multilevel_coarsen_mode;/* pass on to Multilevel_control->coarsen_mode */
int quadtree_size;/* cut off size above which quadtree approximation is used */
int max_qtree_level;/* max level of quadtree */
- real bh;/* Barnes-Hutt constant, if width(snode)/dist[i,snode] < bh, treat snode as a supernode. default 0.2*/
- real tol;/* minimum different between two subsequence config before terminating. ||x-xold|| < tol */
+ double bh;/* Barnes-Hutt constant, if width(snode)/dist[i,snode] < bh, treat snode as a supernode. default 0.2*/
+ double tol;/* minimum different between two subsequence config before terminating. ||x-xold|| < tol */
int maxiter;
- real cool;/* default 0.9 */
- real step;/* initial step size */
+ double cool;/* default 0.9 */
+ double step;/* initial step size */
int adaptive_cooling;
int random_seed;
int beautify_leaves;
int do_shrinking;
int tscheme; /* octree scheme. 0 (no octree), 1 (normal), 2 (fast) */
int method;/* spring_electical, spring_maxent */
- real initial_scaling;/* how to scale the layout of the graph before passing to overlap removal algorithm.
+ double initial_scaling;/* how to scale the layout of the graph before passing to overlap removal algorithm.
positive values are absolute in points, negative values are relative
to average label size.
*/
- real rotation;/* degree of rotation */
+ double rotation;/* degree of rotation */
int edge_labeling_scheme; /* specifying whether to treat node of the form |edgelabel|* as a special node representing an edge label.
0 (no action, default), 1 (penalty based method to make that kind of node close to the center of its neighbor),
1 (penalty based method to make that kind of node close to the old center of its neighbor),
spring_electrical_control spring_electrical_control_new(void);
void spring_electrical_control_print(spring_electrical_control ctrl);
-void spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, int *flag);
-void spring_electrical_embedding_fast(int dim, SparseMatrix A0, spring_electrical_control ctrl, real *node_weights, real *x, int *flag);
+void spring_electrical_embedding(int dim, SparseMatrix A0, spring_electrical_control ctrl, double *node_weights, double *x, int *flag);
+void spring_electrical_embedding_fast(int dim, SparseMatrix A0, spring_electrical_control ctrl, double *node_weights, double *x, int *flag);
-void multilevel_spring_electrical_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *label_sizes,
- real *x, int n_edge_label_nodes, int *edge_label_nodes, int *flag);
+void multilevel_spring_electrical_embedding(int dim, SparseMatrix A0, SparseMatrix D, spring_electrical_control ctrl, double *node_weights, double *label_sizes,
+ double *x, int n_edge_label_nodes, int *edge_label_nodes, int *flag);
-void export_embedding(FILE *fp, int dim, SparseMatrix A, real *x, real *width);
+void export_embedding(FILE *fp, int dim, SparseMatrix A, double *x, double *width);
void spring_electrical_control_delete(spring_electrical_control ctrl);
-void print_matrix(real *x, int n, int dim);
+void print_matrix(double *x, int n, int dim);
-real average_edge_length(SparseMatrix A, int dim, real *coord);
+double average_edge_length(SparseMatrix A, int dim, double *coord);
-void spring_electrical_spring_embedding(int dim, SparseMatrix A, SparseMatrix D, spring_electrical_control ctrl, real *node_weights, real *x, int *flag);
-void force_print(FILE *fp, int n, int dim, real *x, real *force);
+void spring_electrical_spring_embedding(int dim, SparseMatrix A, SparseMatrix D, spring_electrical_control ctrl, double *node_weights, double *x, int *flag);
+void force_print(FILE *fp, int n, int dim, double *x, double *force);
enum {MAX_I = 20, OPT_UP = 1, OPT_DOWN = -1, OPT_INIT = 0};
struct oned_optimizer_struct{
int i;
- real work[MAX_I+1];
+ double work[MAX_I+1];
int direction;
};
typedef struct oned_optimizer_struct *oned_optimizer;
void oned_optimizer_delete(oned_optimizer opt);
oned_optimizer oned_optimizer_new(int i);
-void oned_optimizer_train(oned_optimizer opt, real work);
+void oned_optimizer_train(oned_optimizer opt, double work);
int oned_optimizer_get(oned_optimizer opt);
-void interpolate_coord(int dim, SparseMatrix A, real *x);
+void interpolate_coord(int dim, SparseMatrix A, double *x);
int power_law_graph(SparseMatrix A);
-void pcp_rotate(int n, int dim, real *x);
+void pcp_rotate(int n, int dim, double *x);
#include <sfdpgen/post_process.h>
#include <sfdpgen/stress_model.h>
-static void stress_model_core(int dim, SparseMatrix B, real **x, int edge_len_weighted, int maxit_sm, real tol, int *flag){
+static void stress_model_core(int dim, SparseMatrix B, double **x, int edge_len_weighted, int maxit_sm, double tol, int *flag){
int m;
SparseStressMajorizationSmoother sm;
- real lambda = 0;
+ double lambda = 0;
int i;
SparseMatrix A = B;
*flag = 0;
m = A->m;
if (!x) {
- *x = MALLOC(sizeof(real)*m*dim);
+ *x = MALLOC(sizeof(double)*m*dim);
srand(123);
for (i = 0; i < dim*m; i++) (*x)[i] = drand();
}
if (A != B) SparseMatrix_delete(A);
}
-void stress_model(int dim, SparseMatrix A, SparseMatrix D, real **x, int edge_len_weighted, int maxit_sm, real tol, int *flag){
+void stress_model(int dim, SparseMatrix A, SparseMatrix D, double **x, int edge_len_weighted, int maxit_sm, double tol, int *flag){
stress_model_core(dim, D, x, edge_len_weighted, maxit_sm, tol, flag);
}
#pragma once
-void stress_model(int dim, SparseMatrix A, SparseMatrix D, real **x, int edge_len_weighted, int maxit, real tol, int *flag);
+void stress_model(int dim, SparseMatrix A, SparseMatrix D, double **x, int edge_len_weighted, int maxit, double tol, int *flag);
*/
-UniformStressSmoother UniformStressSmoother_new(int dim, SparseMatrix A, real *x, real alpha, real M, int *flag){
+UniformStressSmoother UniformStressSmoother_new(int dim, SparseMatrix A, double *x, double alpha, double M, int *flag){
UniformStressSmoother sm;
int i, j, k, m = A->m, *ia = A->ia, *ja = A->ja, *iw, *jw, *id, *jd;
int nz;
- real *d, *w, *a = (real*) A->a;
- real diag_d, diag_w, dist, epsilon = 0.01;
+ double *d, *w, *a = (double*) A->a;
+ double diag_d, diag_w, dist, epsilon = 0.01;
assert(SparseMatrix_is_symmetric(A, FALSE));
sm->data = NULL;
sm->scheme = SM_SCHEME_UNIFORM_STRESS;
sm->lambda = NULL;
- sm->data = MALLOC(sizeof(real)*2);
- ((real*) sm->data)[0] = alpha;
- ((real*) sm->data)[1] = M;
+ sm->data = MALLOC(sizeof(double)*2);
+ ((double*) sm->data)[0] = alpha;
+ ((double*) sm->data)[1] = M;
sm->data_deallocator = free;
sm->tol_cg = 0.01;
sm->maxit_cg = (int)sqrt((double) A->m);
iw = sm->Lw->ia; jw = sm->Lw->ja;
id = sm->Lwd->ia; jd = sm->Lwd->ja;
- w = (real*) sm->Lw->a; d = (real*) sm->Lwd->a;
+ w = (double*) sm->Lw->a; d = (double*) sm->Lwd->a;
iw[0] = id[0] = 0;
nz = 0;
}
-static real UniformStressSmoother_smooth(UniformStressSmoother sm, int dim, real *x, int maxit_sm) {
+static double UniformStressSmoother_smooth(UniformStressSmoother sm, int dim, double *x, int maxit_sm) {
return StressMajorizationSmoother_smooth(sm, dim, x, maxit_sm, 0.001);
}
-static SparseMatrix get_distance_matrix(SparseMatrix A, real scaling){
+static SparseMatrix get_distance_matrix(SparseMatrix A, double scaling){
/* get a distance matrix from a graph, at the moment we just symmetrize the matrix. At the moment if the matrix is not real,
we just assume distance of 1 among edges. Then we apply scaling to the entire matrix */
SparseMatrix B;
- real *val;
+ double *val;
int i;
if (A->type == MATRIX_TYPE_REAL){
} else {
B = SparseMatrix_get_real_adjacency_matrix_symmetrized(A);
}
- val = (real*) B->a;
+ val = (double*) B->a;
if (scaling != 1) for (i = 0; i < B->nz; i++) val[i] *= scaling;
return B;
}
-extern void scale_to_box(real xmin, real ymin, real xmax, real ymax, int n, int dim, real *x);
+extern void scale_to_box(double xmin, double ymin, double xmax, double ymax, int n, int dim, double *x);
-void uniform_stress(int dim, SparseMatrix A, real *x, int *flag){
+void uniform_stress(int dim, SparseMatrix A, double *x, int *flag){
UniformStressSmoother sm;
- real lambda0 = 10.1, M = 100, scaling = 1.;
+ double lambda0 = 10.1, M = 100, scaling = 1.;
int maxit = 300, samepoint = TRUE, i, k, n = A->m;
SparseMatrix B = NULL;
void UniformStressSmoother_delete(UniformStressSmoother sm);
-UniformStressSmoother UniformStressSmoother_new(int dim, SparseMatrix A, real *x, real alpha, real M, int *flag);
+UniformStressSmoother UniformStressSmoother_new(int dim, SparseMatrix A, double *x, double alpha, double M, int *flag);
-void uniform_stress(int dim, SparseMatrix A, real *x, int *flag);
+void uniform_stress(int dim, SparseMatrix A, double *x, int *flag);
#define ND_id(n) (((Agnodeinfo_t*)((n)->base.data))->id)
-static void color_string(int slen, char *buf, int dim, real *color){
+static void color_string(int slen, char *buf, int dim, double *color){
if (dim > 3 || dim < 1){
fprintf(stderr,"can only 1, 2 or 3 dimensional color space. with color value between 0 to 1\n");
assert(0);
}
}
-void attach_edge_colors(Agraph_t* g, int dim, real *colors){
+void attach_edge_colors(Agraph_t* g, int dim, double *colors){
/* colors is array of dim*nedges, with color for edge i at colors[dim*i, dim(i+1))
*/
Agsym_t* sym = agattr(g, AGEDGE, "color", 0);
* but not a->b and a->b
*/
SparseMatrix
-SparseMatrix_import_dot (Agraph_t* g, int dim, real **label_sizes, real **x, int *n_edge_label_nodes, int **edge_label_nodes, int format, SparseMatrix *D)
+SparseMatrix_import_dot (Agraph_t* g, int dim, double **label_sizes, double **x, int *n_edge_label_nodes, int **edge_label_nodes, int format, SparseMatrix *D)
{
SparseMatrix A = 0;
Agnode_t* n;
int i, row;
int* I;
int* J;
- real *val, *valD = NULL;
- real v;
+ double *val, *valD = NULL;
+ double v;
int type = MATRIX_TYPE_REAL;
- real padding = 10;
+ double padding = 10;
int nedge_nodes = 0;
A->nz = nedges;
I = A->ia;
J = A->ja;
- val = (real*) A->a;
+ val = (double*) A->a;
} else {
I = N_NEW(nedges, int);
J = N_NEW(nedges, int);
- val = N_NEW(nedges, real);
+ val = N_NEW(nedges, double);
}
sym = agattr(g, AGEDGE, "weight", NULL);
if (D) {
symD = agattr(g, AGEDGE, "len", NULL);
- valD = N_NEW(nedges, real);
+ valD = N_NEW(nedges, double);
}
i = 0;
for (n = agfstnode (g); n; n = agnxtnode (g, n)) {
}
- if (label_sizes) *label_sizes = MALLOC(sizeof(real)*2*nnodes);
+ if (label_sizes) *label_sizes = MALLOC(sizeof(double)*2*nnodes);
for (n = agfstnode (g); n; n = agnxtnode (g, n)) {
- real sz;
+ double sz;
i = ND_id(n);
if (edge_label_nodes && strncmp(agnameof(n), "|edgelabel|",11)==0) {
(*edge_label_nodes)[nedge_nodes++] = i;
bool has_positions = true;
char* pval;
if (!(*x)) {
- *x = MALLOC(sizeof(real)*dim*nnodes);
+ *x = MALLOC(sizeof(double)*dim*nnodes);
assert(*x);
}
for (n = agfstnode (g); n && has_positions; n = agnxtnode (g, n)) {
- real xx,yy, zz,ww;
+ double xx,yy, zz,ww;
int nitems;
i = ND_id(n);
if ((pval = agxget(n, psym)) && *pval) {
else if (x)
agerr (AGERR, "Error: graph %s has missing \"pos\" information", agnameof(g));
- size_t sz = sizeof(real);
+ size_t sz = sizeof(double);
if (format == FORMAT_CSR) A = SparseMatrix_from_coordinate_arrays(nedges, nnodes, nnodes, I, J, val, type, sz);
if (edge_label_nodes) *n_edge_label_nodes = nedge_nodes;
}
}
-SparseMatrix Import_coord_clusters_from_dot(Agraph_t* g, int maxcluster, int dim, int *nn, real **label_sizes, real **x, int **clusters, float **rgb_r, float **rgb_g, float **rgb_b, float **fsz, char ***labels, int default_color_scheme, int clustering_scheme, int useClusters){
+SparseMatrix Import_coord_clusters_from_dot(Agraph_t* g, int maxcluster, int dim, int *nn, double **label_sizes, double **x, int **clusters, float **rgb_r, float **rgb_g, float **rgb_b, float **fsz, char ***labels, int default_color_scheme, int clustering_scheme, int useClusters){
SparseMatrix A = 0;
Agnode_t* n;
Agedge_t* e;
int i, row, ic,nc, j;
int* I;
int* J;
- real* val;
- real v;
+ double* val;
+ double v;
int type = MATRIX_TYPE_REAL;
char scluster[100];
float ff;
/* form matrix */
I = N_NEW(nedges, int);
J = N_NEW(nedges, int);
- val = N_NEW(nedges, real);
+ val = N_NEW(nedges, double);
sym = agattr(g, AGEDGE, "weight", NULL);
clust_sym = agattr(g, AGNODE, "cluster", NULL);
}
}
A = SparseMatrix_from_coordinate_arrays(nedges, nnodes, nnodes, I, J, val,
- type, sizeof(real));
+ type, sizeof(double));
/* get clustering info */
*clusters = MALLOC(sizeof(int)*nnodes);
if (noclusterinfo) {
int use_value = TRUE, flag = 0;
- real modularity;
+ double modularity;
if (!clust_sym) clust_sym = agattr(g,AGNODE,"cluster","-1");
if (clustering_scheme == CLUSTERING_MQ){
}
}
- *label_sizes = MALLOC(sizeof(real)*dim*nnodes);
+ *label_sizes = MALLOC(sizeof(double)*dim*nnodes);
if (pal || (!noclusterinfo && clust_clr_sym)){
*rgb_r = MALLOC(sizeof(float)*(1+MAX_GRPS));
*rgb_g = MALLOC(sizeof(float)*(1+MAX_GRPS));
for (n = agfstnode (g); n; n = agnxtnode (g, n)) {
gvcolor_t color;
- real sz;
+ double sz;
i = ND_id(n);
if (agget(n, "width") && agget(n, "height")){
sscanf(agget(n, "width"), "%lf", &sz);
if (x){
bool has_position = false;
- *x = MALLOC(sizeof(real)*dim*nnodes);
+ *x = MALLOC(sizeof(double)*dim*nnodes);
for (n = agfstnode (g); n; n = agnxtnode (g, n)) {
- real xx,yy;
+ double xx,yy;
i = ND_id(n);
if (agget(n, "pos")){
has_position = true;
int i, row,nc;
int* I;
int* J;
- real* val;
- real v;
+ double* val;
+ double v;
int type = MATRIX_TYPE_REAL;
- size_t sz = sizeof(real);
+ size_t sz = sizeof(double);
int *clusters;
/* form matrix */
I = N_NEW(nedges, int);
J = N_NEW(nedges, int);
- val = N_NEW(nedges, real);
+ val = N_NEW(nedges, double);
sym = agattr(g, AGEDGE, "weight", NULL);
clust_sym = agattr(g, AGNODE, "cluster", NULL);
{
int use_value = TRUE, flag = 0;
- real modularity;
+ double modularity;
if (!clust_sym) clust_sym = agattr(g,AGNODE,"cluster","-1");
if (clustering_scheme == CLUSTERING_MQ){
extern void setDotNodeID (Agnode_t* n, int v);
extern int getDotNodeID (Agnode_t* n);
-extern void attach_edge_colors(Agraph_t* g, int dim, real *colors);
+extern void attach_edge_colors(Agraph_t* g, int dim, double *colors);
-extern SparseMatrix SparseMatrix_import_dot(Agraph_t* g, int dim, real **label_sizes, real **x, int *n_edge_label_nodes,
+extern SparseMatrix SparseMatrix_import_dot(Agraph_t* g, int dim, double **label_sizes, double **x, int *n_edge_label_nodes,
int **edge_label_nodes, int format, SparseMatrix *D);
-char * hue2rgb(real hue, char *color);
+char * hue2rgb(double hue, char *color);
-SparseMatrix Import_coord_clusters_from_dot(Agraph_t* g, int maxcluster, int dim, int *nn, real **label_sizes, real **x, int **clusters, float **rgb_r, float **rgb_g, float **rgb_b, float **fsz, char ***labels, int default_color_scheme, int clustering_scheme, int useClusters);
+SparseMatrix Import_coord_clusters_from_dot(Agraph_t* g, int maxcluster, int dim, int *nn, double **label_sizes, double **x, int **clusters, float **rgb_r, float **rgb_g, float **rgb_b, float **fsz, char ***labels, int default_color_scheme, int clustering_scheme, int useClusters);
void Dot_SetClusterColor(Agraph_t* g, float *rgb_r, float *rgb_g, float *rgb_b, int *clustering);
void attached_clustering(Agraph_t* g, int maxcluster, int clustering_scheme);
#include <sparse/LinkedList.h>
#include <sparse/QuadTree.h>
-extern real distance_cropped(real *x, int dim, int i, int j);
+extern double distance_cropped(double *x, int dim, int i, int j);
struct node_data_struct {
- real node_weight;
- real *coord;
- real id;
+ double node_weight;
+ double *coord;
+ double id;
void *data;
};
typedef struct node_data_struct *node_data;
-static node_data node_data_new(int dim, real weight, real *coord, int id){
+static node_data node_data_new(int dim, double weight, double *coord, int id){
node_data nd;
int i;
nd = MALLOC(sizeof(struct node_data_struct));
nd->node_weight = weight;
- nd->coord = MALLOC(sizeof(real)*dim);
+ nd->coord = MALLOC(sizeof(double)*dim);
nd->id = id;
for (i = 0; i < dim; i++) nd->coord[i] = coord[i];
nd->data = NULL;
free(nd);
}
-static real node_data_get_weight(void *d){
+static double node_data_get_weight(void *d){
node_data nd = (node_data) d;
return nd->node_weight;
}
-static real* node_data_get_coord(void *d){
+static double* node_data_get_coord(void *d){
node_data nd = (node_data) d;
return nd->coord;
}
#define node_data_get_data(d) (((node_data) (d))->data)
-static void check_or_realloc_arrays(int dim, int *nsuper, int *nsupermax, real **center, real **supernode_wgts, real **distances){
+static void check_or_realloc_arrays(int dim, int *nsuper, int *nsupermax, double **center, double **supernode_wgts, double **distances){
if (*nsuper >= *nsupermax) {
*nsupermax = *nsuper + MAX(10, (int) 0.2*(*nsuper));
- *center = REALLOC(*center, sizeof(real)*(*nsupermax)*dim);
- *supernode_wgts = REALLOC(*supernode_wgts, sizeof(real)*(*nsupermax));
- *distances = REALLOC(*distances, sizeof(real)*(*nsupermax));
+ *center = REALLOC(*center, sizeof(double)*(*nsupermax)*dim);
+ *supernode_wgts = REALLOC(*supernode_wgts, sizeof(double)*(*nsupermax));
+ *distances = REALLOC(*distances, sizeof(double)*(*nsupermax));
}
}
-static void QuadTree_get_supernodes_internal(QuadTree qt, real bh, real *point, int nodeid, int *nsuper, int *nsupermax, real **center, real **supernode_wgts, real **distances, real *counts, int *flag){
+static void QuadTree_get_supernodes_internal(QuadTree qt, double bh, double *point, int nodeid, int *nsuper, int *nsupermax, double **center, double **supernode_wgts, double **distances, double *counts, int *flag){
SingleLinkedList l;
- real *coord, dist;
+ double *coord, dist;
int dim, i;
(*counts)++;
}
-void QuadTree_get_supernodes(QuadTree qt, real bh, real *point, int nodeid, int *nsuper,
- int *nsupermax, real **center, real **supernode_wgts, real **distances, double *counts, int *flag){
+void QuadTree_get_supernodes(QuadTree qt, double bh, double *point, int nodeid, int *nsuper,
+ int *nsupermax, double **center, double **supernode_wgts, double **distances, double *counts, int *flag){
int dim = qt->dim;
(*counts) = 0;
*flag = 0;
*nsupermax = 10;
- if (!*center) *center = MALLOC(sizeof(real)*(*nsupermax)*dim);
- if (!*supernode_wgts) *supernode_wgts = MALLOC(sizeof(real)*(*nsupermax));
- if (!*distances) *distances = MALLOC(sizeof(real)*(*nsupermax));
+ if (!*center) *center = MALLOC(sizeof(double)*(*nsupermax)*dim);
+ if (!*supernode_wgts) *supernode_wgts = MALLOC(sizeof(double)*(*nsupermax));
+ if (!*distances) *distances = MALLOC(sizeof(double)*(*nsupermax));
QuadTree_get_supernodes_internal(qt, bh, point, nodeid, nsuper, nsupermax, center, supernode_wgts, distances, counts, flag);
}
-static real *get_or_assign_node_force(real *force, int i, SingleLinkedList l, int dim){
+static double *get_or_assign_node_force(double *force, int i, SingleLinkedList l, int dim){
- real *f = (real*) node_data_get_data(SingleLinkedList_get_data(l));
+ double *f = (double*) node_data_get_data(SingleLinkedList_get_data(l));
if (!f){
node_data_get_data(SingleLinkedList_get_data(l)) = &(force[i*dim]);
- f = (real*) node_data_get_data(SingleLinkedList_get_data(l));
+ f = (double*) node_data_get_data(SingleLinkedList_get_data(l));
}
return f;
}
-static real *get_or_alloc_force_qt(QuadTree qt, int dim){
+static double *get_or_alloc_force_qt(QuadTree qt, int dim){
int i;
- real *force = (real*) qt->data;
+ double *force = (double*) qt->data;
if (!force){
- qt->data = MALLOC(sizeof(real)*dim);
- force = (real*) qt->data;
+ qt->data = MALLOC(sizeof(double)*dim);
+ force = (double*) qt->data;
for (i = 0; i < dim; i++) force[i] = 0.;
}
return force;
}
-static void QuadTree_repulsive_force_interact(QuadTree qt1, QuadTree qt2, real *x, real *force, real bh, real p, real KP, real *counts){
+static void QuadTree_repulsive_force_interact(QuadTree qt1, QuadTree qt2, double *x, double *force, double bh, double p, double KP, double *counts){
/* calculate the all to all reopulsive force and accumulate on each node of the quadtree if an interaction is possible.
force[i*dim+j], j=1,...,dim is the force on node i
*/
SingleLinkedList l1, l2;
- real *x1, *x2, dist, wgt1, wgt2, f, *f1, *f2, w1, w2;
+ double *x1, *x2, dist, wgt1, wgt2, f, *f1, *f2, w1, w2;
int dim, i, j, i1, i2, k;
QuadTree qt11, qt12;
}
}
-static void QuadTree_repulsive_force_accumulate(QuadTree qt, real *force, real *counts){
+static void QuadTree_repulsive_force_accumulate(QuadTree qt, double *force, double *counts){
/* push down forces on cells into the node level */
- real wgt, wgt2;
- real *f, *f2;
+ double wgt, wgt2;
+ double *f, *f2;
SingleLinkedList l = qt->l;
int i, k, dim;
QuadTree qt2;
}
-void QuadTree_get_repulsive_force(QuadTree qt, real *force, real *x, real bh, real p, real KP, real *counts, int *flag){
+void QuadTree_get_repulsive_force(QuadTree qt, double *force, double *x, double bh, double p, double KP, double *counts, int *flag){
/* get repulsice force by a more efficient algortihm: we consider two cells, if they are well separated, we
calculate the overall repulsive force on the cell level, if not well separated, we divide one of the cell.
If both cells are at the leaf level, we calcuaulate repulsicve force among individual nodes. Finally
for (i = 0; i < 4; i++) counts[i] /= n;
}
-QuadTree QuadTree_new_from_point_list(int dim, int n, int max_level, real *coord, real *weight){
+QuadTree QuadTree_new_from_point_list(int dim, int n, int max_level, double *coord, double *weight){
/* form a new QuadTree data structure from a list of coordinates of n points
coord: of length n*dim, point i sits at [i*dim, i*dim+dim - 1]
weight: node weight of lentgth n. If NULL, unit weight assumed.
*/
- real *xmin, *xmax, *center, width;
+ double *xmin, *xmax, *center, width;
QuadTree qt = NULL;
int i, k;
- xmin = MALLOC(sizeof(real)*dim);
- xmax = MALLOC(sizeof(real)*dim);
- center = MALLOC(sizeof(real)*dim);
+ xmin = MALLOC(sizeof(double)*dim);
+ xmax = MALLOC(sizeof(double)*dim);
+ center = MALLOC(sizeof(double)*dim);
if (!xmin || !xmax || !center) {
free(xmin);
free(xmax);
return qt;
}
-QuadTree QuadTree_new(int dim, real *center, real width, int max_level){
+QuadTree QuadTree_new(int dim, double *center, double width, int max_level){
QuadTree q;
int i;
q = MALLOC(sizeof(struct QuadTree_struct));
q->dim = dim;
q->n = 0;
- q->center = MALLOC(sizeof(real)*dim);
+ q->center = MALLOC(sizeof(double)*dim);
for (i = 0; i < dim; i++) q->center[i] = center[i];
assert(width > 0);
q->width = width;
free(q);
}
-static int QuadTree_get_quadrant(int dim, real *center, real *coord){
+static int QuadTree_get_quadrant(int dim, double *center, double *coord){
/* find the quadrant that a point of coordinates coord is going into with reference to the center.
if coord - center == {+,-,+,+} = {1,0,1,1}, then it will sit in the i-quadrant where
i's binary representation is 1011 (that is, decimal 11).
return d;
}
-QuadTree QuadTree_new_in_quadrant(int dim, real *center, real width, int max_level, int i){
+QuadTree QuadTree_new_in_quadrant(int dim, double *center, double width, int max_level, int i){
/* a new quadtree in quadrant i of the original cell. The original cell is centered at 'center".
The new cell have width "width".
*/
return qt;
}
-static QuadTree QuadTree_add_internal(QuadTree q, real *coord, real weight, int id, int level){
+static QuadTree QuadTree_add_internal(QuadTree q, double *coord, double weight, int id, int level){
int i, dim = q->dim, ii;
node_data nd = NULL;
/* if this node is empty right now */
q->n = 1;
q->total_weight = weight;
- q->average = MALLOC(sizeof(real)*dim);
+ q->average = MALLOC(sizeof(double)*dim);
for (i = 0; i < q->dim; i++) q->average[i] = coord[i];
nd = node_data_new(q->dim, weight, coord, id);
assert(!(q->l));
}
-QuadTree QuadTree_add(QuadTree q, real *coord, real weight, int id){
+QuadTree QuadTree_add(QuadTree q, double *coord, double weight, int id){
if (!q) return q;
return QuadTree_add_internal(q, coord, weight, id, 0);
}
-static void draw_polygon(FILE *fp, int dim, real *center, real width){
+static void draw_polygon(FILE *fp, int dim, double *center, double width){
/* pliot the enclosing square */
if (dim < 2 || dim > 3) return;
fprintf(fp,"(*in c*){Line[{");
static void QuadTree_print_internal(FILE *fp, QuadTree q, int level){
/* dump a quad tree in Mathematica format. */
SingleLinkedList l, l0;
- real *coord;
+ double *coord;
int i, dim;
if (!q) return;
-static void QuadTree_get_nearest_internal(QuadTree qt, real *x, real *y, real *min, int *imin, int tentative, int *flag){
+static void QuadTree_get_nearest_internal(QuadTree qt, double *x, double *y, double *min, int *imin, int tentative, int *flag){
/* get the narest point years to {x[0], ..., x[dim]} and store in y.*/
SingleLinkedList l;
- real *coord, dist;
+ double *coord, dist;
int dim, i, iq = -1;
- real qmin;
- real *point = x;
+ double qmin;
+ double *point = x;
*flag = 0;
if (!qt) return;
if (qt->qts){
dist = point_distance(qt->center, point, dim);
- if (*min >= 0 && (dist - sqrt((real) dim) * qt->width > *min)){
+ if (*min >= 0 && (dist - sqrt((double) dim) * qt->width > *min)){
return;
} else {
if (tentative){/* quick first approximation*/
}
-void QuadTree_get_nearest(QuadTree qt, real *x, real *ymin, int *imin, real *min, int *flag){
+void QuadTree_get_nearest(QuadTree qt, double *x, double *ymin, int *imin, double *min, int *flag){
*flag = 0;
*min = -1;
2^dim QuadTrees qts. At the last level, all coordinates are stored in a single linked list l.
total_weight is the combined weights of the nodes */
int n;/* number of items */
- real total_weight;
+ double total_weight;
int dim;
- real *center;/* center of the bounding box, array of dimension dim. Allocated inside */
- real width;/* center +/- width gives the lower/upper bound, so really width is the
+ double *center;/* center of the bounding box, array of dimension dim. Allocated inside */
+ double width;/* center +/- width gives the lower/upper bound, so really width is the
"radius" */
- real *average;/* the average coordinates. Array of length dim. Allocated inside */
+ double *average;/* the average coordinates. Array of length dim. Allocated inside */
QuadTree *qts;/* subtree . If dim = 2, there are 4, dim = 3 gives 8 */
SingleLinkedList l;
int max_level;
};
-QuadTree QuadTree_new(int dim, real *center, real width, int max_level);
+QuadTree QuadTree_new(int dim, double *center, double width, int max_level);
void QuadTree_delete(QuadTree q);
-QuadTree QuadTree_add(QuadTree q, real *coord, real weight, int id);/* coord is copied in */
+QuadTree QuadTree_add(QuadTree q, double *coord, double weight, int id);/* coord is copied in */
void QuadTree_print(FILE *fp, QuadTree q);
-QuadTree QuadTree_new_from_point_list(int dim, int n, int max_level, real *coord, real *weight);
+QuadTree QuadTree_new_from_point_list(int dim, int n, int max_level, double *coord, double *weight);
-real point_distance(real *p1, real *p2, int dim);
+double point_distance(double *p1, double *p2, int dim);
-void QuadTree_get_supernodes(QuadTree qt, real bh, real *point, int nodeid, int *nsuper,
- int *nsupermax, real **center, real **supernode_wgts, real **distances, real *counts, int *flag);
+void QuadTree_get_supernodes(QuadTree qt, double bh, double *point, int nodeid, int *nsuper,
+ int *nsupermax, double **center, double **supernode_wgts, double **distances, double *counts, int *flag);
-void QuadTree_get_repulsive_force(QuadTree qt, real *force, real *x, real bh, real p, real KP, real *counts, int *flag);
+void QuadTree_get_repulsive_force(QuadTree qt, double *force, double *x, double bh, double p, double KP, double *counts, int *flag);
/* find the nearest point and put in ymin, index in imin and distance in min */
-void QuadTree_get_nearest(QuadTree qt, real *x, real *ymin, int *imin, real *min, int *flag);
+void QuadTree_get_nearest(QuadTree qt, double *x, double *ymin, int *imin, double *min, int *flag);
-QuadTree QuadTree_new_in_quadrant(int dim, real *center, real width, int max_level, int i);
+QuadTree QuadTree_new_in_quadrant(int dim, double *center, double width, int max_level, int i);
size_t size = 0;
switch (type){
case MATRIX_TYPE_REAL:
- size = sizeof(real);
+ size = sizeof(double);
break;
case MATRIX_TYPE_COMPLEX:
- size = 2*sizeof(real);
+ size = 2*sizeof(double);
break;
case MATRIX_TYPE_INTEGER:
size = sizeof(int);
switch (A->type){
case MATRIX_TYPE_REAL:{
- real *a = (real*) A->a;
- real *b = (real*) B->a;
+ double *a = (double*) A->a;
+ double *b = (double*) B->a;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
jb[ib[ja[j]]] = i;
break;
}
case MATRIX_TYPE_COMPLEX:{
- real *a = (real*) A->a;
- real *b = (real*) B->a;
+ double *a = (double*) A->a;
+ double *b = (double*) B->a;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
jb[ib[ja[j]]] = i;
switch (type){
case MATRIX_TYPE_REAL:{
- real *a = (real*) A->a;
- real *b = (real*) B->a;
+ double *a = (double*) A->a;
+ double *b = (double*) B->a;
for (i = 0; i <= m; i++) if (ia[i] != ib[i]) goto RETURN;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
break;
}
case MATRIX_TYPE_COMPLEX:{
- real *a = (real*) A->a;
- real *b = (real*) B->a;
+ double *a = (double*) A->a;
+ double *b = (double*) B->a;
for (i = 0; i <= m; i++) if (ia[i] != ib[i]) goto RETURN;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
}
static void SparseMatrix_print_csr(char *c, SparseMatrix A){
int *ia, *ja;
- real *a;
+ double *a;
int *ai;
int i, j, m = A->m;
a = A->a;
switch (A->type){
case MATRIX_TYPE_REAL:
- a = (real*) A->a;
+ a = (double*) A->a;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
printf("{%d, %d}->%f",i+1, ja[j]+1, a[j]);
}
break;
case MATRIX_TYPE_COMPLEX:
- a = (real*) A->a;
+ a = (double*) A->a;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
printf("{%d, %d}->%f + %f I",i+1, ja[j]+1, a[2*j], a[2*j+1]);
static void SparseMatrix_print_coord(char *c, SparseMatrix A){
int *ia, *ja;
- real *a;
+ double *a;
int *ai;
int i, m = A->m;
a = A->a;
switch (A->type){
case MATRIX_TYPE_REAL:
- a = (real*) A->a;
+ a = (double*) A->a;
for (i = 0; i < A->nz; i++){
printf("{%d, %d}->%f",ia[i]+1, ja[i]+1, a[i]);
if (i != A->nz - 1) printf(",");
printf("\n");
break;
case MATRIX_TYPE_COMPLEX:
- a = (real*) A->a;
+ a = (double*) A->a;
for (i = 0; i < A->nz; i++){
printf("{%d, %d}->%f + %f I",ia[i]+1, ja[i]+1, a[2*i], a[2*i+1]);
if (i != A->nz - 1) printf(",");
static void SparseMatrix_export_csr(FILE *f, SparseMatrix A){
int *ia, *ja;
- real *a;
+ double *a;
int *ai;
int i, j, m = A->m;
a = A->a;
switch (A->type){
case MATRIX_TYPE_REAL:
- a = (real*) A->a;
+ a = (double*) A->a;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
fprintf(f, "%d %d %16.8g\n",i+1, ja[j]+1, a[j]);
}
break;
case MATRIX_TYPE_COMPLEX:
- a = (real*) A->a;
+ a = (double*) A->a;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
fprintf(f, "%d %d %16.8g %16.8g\n",i+1, ja[j]+1, a[2*j], a[2*j+1]);
static void SparseMatrix_export_coord(FILE *f, SparseMatrix A){
int *ia, *ja;
- real *a;
+ double *a;
int *ai;
int i;
a = A->a;
switch (A->type){
case MATRIX_TYPE_REAL:
- a = (real*) A->a;
+ a = (double*) A->a;
for (i = 0; i < A->nz; i++){
fprintf(f, "%d %d %16.8g\n",ia[i]+1, ja[i]+1, a[i]);
}
break;
case MATRIX_TYPE_COMPLEX:
- a = (real*) A->a;
+ a = (double*) A->a;
for (i = 0; i < A->nz; i++){
fprintf(f, "%d %d %16.8g %16.8g\n",ia[i]+1, ja[i]+1, a[2*i], a[2*i+1]);
}
SparseMatrix A = NULL;
int *ia, *ja;
- real *a, *val;
+ double *a, *val;
int *ai, *vali;
int i;
switch (type){
case MATRIX_TYPE_REAL:
- val = (real*) val0;
- a = (real*) A->a;
+ val = (double*) val0;
+ a = (double*) A->a;
for (i = 0; i < nz; i++){
if (irn[i] < 0 || irn[i] >= m || jcn[i] < 0 || jcn[i] >= n) {
assert(0);
ia[0] = 0;
break;
case MATRIX_TYPE_COMPLEX:
- val = (real*) val0;
- a = (real*) A->a;
+ val = (double*) val0;
+ a = (double*) A->a;
for (i = 0; i < nz; i++){
if (irn[i] < 0 || irn[i] >= m || jcn[i] < 0 || jcn[i] >= n) {
assert(0);
ic[0] = 0;
switch (A->type){
case MATRIX_TYPE_REAL:{
- real *a = (real*) A->a;
- real *b = (real*) B->a;
- real *c = (real*) C->a;
+ double *a = (double*) A->a;
+ double *b = (double*) B->a;
+ double *c = (double*) C->a;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
mask[ja[j]] = nz;
break;
}
case MATRIX_TYPE_COMPLEX:{
- real *a = (real*) A->a;
- real *b = (real*) B->a;
- real *c = (real*) C->a;
+ double *a = (double*) A->a;
+ double *b = (double*) B->a;
+ double *c = (double*) C->a;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
mask[ja[j]] = nz;
return C;
}
-static void SparseMatrix_multiply_dense1(SparseMatrix A, real *v, real **res, int dim){
+static void SparseMatrix_multiply_dense1(SparseMatrix A, double *v, double **res, int dim){
/* A v where v a dense matrix of second dimension dim. Real only for now. */
int i, j, k, *ia, *ja, m;
- real *a, *u;
+ double *a, *u;
assert(A->format == FORMAT_CSR);
assert(A->type == MATRIX_TYPE_REAL);
- a = (real*) A->a;
+ a = (double*) A->a;
ia = A->ia;
ja = A->ja;
m = A->m;
u = *res;
- if (!u) u = MALLOC(sizeof(real)*((size_t) m)*((size_t) dim));
+ if (!u) u = MALLOC(sizeof(double)*((size_t) m)*((size_t) dim));
for (i = 0; i < m; i++){
for (k = 0; k < dim; k++) u[i*dim+k] = 0.;
for (j = ia[i]; j < ia[i+1]; j++){
*res = u;
}
-void SparseMatrix_multiply_dense(SparseMatrix A, real *v, real **res, int dim){
+void SparseMatrix_multiply_dense(SparseMatrix A, double *v, double **res, int dim){
/* A * V, with A dimension m x n, with V of dimension n x dim. v[i*dim+j] gives V[i,j]. Result of dimension m x dim
*/
SparseMatrix_multiply_dense1(A, v, res, dim);
}
-void SparseMatrix_multiply_vector(SparseMatrix A, real *v, real **res, int transposed){
+void SparseMatrix_multiply_vector(SparseMatrix A, double *v, double **res, int transposed){
/* A v or A^T v. Real only for now. */
int i, j, *ia, *ja, n, m;
- real *a, *u = NULL;
+ double *a, *u = NULL;
int *ai;
assert(A->format == FORMAT_CSR);
assert(A->type == MATRIX_TYPE_REAL || A->type == MATRIX_TYPE_INTEGER);
switch (A->type){
case MATRIX_TYPE_REAL:
- a = (real*) A->a;
+ a = (double*) A->a;
if (v){
if (!transposed){
- if (!u) u = MALLOC(sizeof(real)*((size_t)m));
+ if (!u) u = MALLOC(sizeof(double)*((size_t)m));
for (i = 0; i < m; i++){
u[i] = 0.;
for (j = ia[i]; j < ia[i+1]; j++){
}
}
} else {
- if (!u) u = MALLOC(sizeof(real)*((size_t)n));
+ if (!u) u = MALLOC(sizeof(double)*((size_t)n));
for (i = 0; i < n; i++) u[i] = 0.;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
} else {
/* v is assumed to be all 1's */
if (!transposed){
- if (!u) u = MALLOC(sizeof(real)*((size_t)m));
+ if (!u) u = MALLOC(sizeof(double)*((size_t)m));
for (i = 0; i < m; i++){
u[i] = 0.;
for (j = ia[i]; j < ia[i+1]; j++){
}
}
} else {
- if (!u) u = MALLOC(sizeof(real)*((size_t)n));
+ if (!u) u = MALLOC(sizeof(double)*((size_t)n));
for (i = 0; i < n; i++) u[i] = 0.;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
ai = (int*) A->a;
if (v){
if (!transposed){
- if (!u) u = MALLOC(sizeof(real)*((size_t)m));
+ if (!u) u = MALLOC(sizeof(double)*((size_t)m));
for (i = 0; i < m; i++){
u[i] = 0.;
for (j = ia[i]; j < ia[i+1]; j++){
}
}
} else {
- if (!u) u = MALLOC(sizeof(real)*((size_t)n));
+ if (!u) u = MALLOC(sizeof(double)*((size_t)n));
for (i = 0; i < n; i++) u[i] = 0.;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
} else {
/* v is assumed to be all 1's */
if (!transposed){
- if (!u) u = MALLOC(sizeof(real)*((size_t)m));
+ if (!u) u = MALLOC(sizeof(double)*((size_t)m));
for (i = 0; i < m; i++){
u[i] = 0.;
for (j = ia[i]; j < ia[i+1]; j++){
}
}
} else {
- if (!u) u = MALLOC(sizeof(real)*((size_t)n));
+ if (!u) u = MALLOC(sizeof(double)*((size_t)n));
for (i = 0; i < n; i++) u[i] = 0.;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
switch (type){
case MATRIX_TYPE_REAL:
{
- real *a = (real*) A->a;
- real *b = (real*) B->a;
- real *c = (real*) C->a;
+ double *a = (double*) A->a;
+ double *b = (double*) B->a;
+ double *c = (double*) C->a;
ic[0] = 0;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
break;
case MATRIX_TYPE_COMPLEX:
{
- real *a = (real*) A->a;
- real *b = (real*) B->a;
- real *c = (real*) C->a;
- a = (real*) A->a;
- b = (real*) B->a;
- c = (real*) C->a;
+ double *a = (double*) A->a;
+ double *b = (double*) B->a;
+ double *c = (double*) C->a;
+ a = (double*) A->a;
+ b = (double*) B->a;
+ c = (double*) C->a;
ic[0] = 0;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
switch (type){
case MATRIX_TYPE_REAL:
{
- real *a = (real*) A->a;
- real *b = (real*) B->a;
- real *c = (real*) C->a;
- real *d = (real*) D->a;
+ double *a = (double*) A->a;
+ double *b = (double*) B->a;
+ double *c = (double*) C->a;
+ double *d = (double*) D->a;
id[0] = 0;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
break;
case MATRIX_TYPE_COMPLEX:
{
- real *a = (real*) A->a;
- real *b = (real*) B->a;
- real *c = (real*) C->a;
- real *d = (real*) D->a;
+ double *a = (double*) A->a;
+ double *b = (double*) B->a;
+ double *c = (double*) C->a;
+ double *d = (double*) D->a;
id[0] = 0;
for (i = 0; i < m; i++){
for (j = ia[i]; j < ia[i+1]; j++){
switch (type){
case MATRIX_TYPE_REAL:
{
- real *a = (real*) A->a;
+ double *a = (double*) A->a;
nz = 0;
sta = ia[0];
for (i = 0; i < A->m; i++){
break;
case MATRIX_TYPE_COMPLEX:
{
- real *a = (real*) A->a;
+ double *a = (double*) A->a;
if (what_to_sum == SUM_REPEATED_ALL){
nz = 0;
sta = ia[0];
sta = ia[0];
switch (A->type){
case MATRIX_TYPE_REAL:{
- real *a = (real*) A->a;
+ double *a = (double*) A->a;
for (i = 0; i < A->m; i++){
for (j = sta; j < ia[i+1]; j++){
if (ja[j] != i){
break;
}
case MATRIX_TYPE_COMPLEX:{
- real *a = (real*) A->a;
+ double *a = (double*) A->a;
for (i = 0; i < A->m; i++){
for (j = sta; j < ia[i+1]; j++){
if (ja[j] != i){
sta = ia[0];
switch (A->type){
case MATRIX_TYPE_REAL:{
- real *a = (real*) A->a;
+ double *a = (double*) A->a;
for (i = 0; i < A->m; i++){
for (j = sta; j < ia[i+1]; j++){
if (ja[j] < i){
break;
}
case MATRIX_TYPE_COMPLEX:{
- real *a = (real*) A->a;
+ double *a = (double*) A->a;
for (i = 0; i < A->m; i++){
for (j = sta; j < ia[i+1]; j++){
if (ja[j] < i){
SparseMatrix SparseMatrix_divide_row_by_degree(SparseMatrix A){
int i, j, *ia, *ja;
- real deg;
+ double deg;
if (!A) return A;
ja = A->ja;
switch (A->type){
case MATRIX_TYPE_REAL:{
- real *a = (real*) A->a;
+ double *a = (double*) A->a;
for (i = 0; i < A->m; i++){
deg = ia[i+1] - ia[i];
for (j = ia[i]; j < ia[i+1]; j++){
break;
}
case MATRIX_TYPE_COMPLEX:{
- real *a = (real*) A->a;
+ double *a = (double*) A->a;
for (i = 0; i < A->m; i++){
deg = ia[i+1] - ia[i];
for (j = ia[i]; j < ia[i+1]; j++){
SparseMatrix SparseMatrix_get_real_adjacency_matrix_symmetrized(SparseMatrix A){
/* symmetric, all entries to 1, diaginal removed */
int i, *ia, *ja, nz, m, n;
- real *a;
+ double *a;
SparseMatrix B;
if (!A) return A;
A = SparseMatrix_symmetrize(B, TRUE);
SparseMatrix_delete(B);
A = SparseMatrix_remove_diagonal(A);
- A->a = MALLOC(sizeof(real)*((size_t)(A->nz)));
- a = (real*) A->a;
+ A->a = MALLOC(sizeof(double)*((size_t)(A->nz)));
+ a = (double*) A->a;
for (i = 0; i < A->nz; i++) a[i] = 1.;
A->type = MATRIX_TYPE_REAL;
- A->size = sizeof(real);
+ A->size = sizeof(double);
return A;
}
SparseMatrix SparseMatrix_apply_fun(SparseMatrix A, double (*fun)(double x)){
int i, j;
- real *a;
+ double *a;
if (!A) return A;
}
- a = (real*) A->a;
+ a = (double*) A->a;
for (i = 0; i < A->m; i++){
for (j = A->ia[i]; j < A->ia[i+1]; j++){
a[j] = fun(a[j]);
struct nodedata_struct {
- real dist;/* distance to root */
+ double dist;/* distance to root */
int id;/*node id */
};
typedef struct nodedata_struct* nodedata;
return 0;
}
-static int Dijkstra_internal(SparseMatrix A, int root, real *dist, int *nlist, int *list, real *dist_max, int *mask){
+static int Dijkstra_internal(SparseMatrix A, int root, double *dist, int *nlist, int *list, double *dist_max, int *mask){
/* Find the shortest path distance of all nodes to root. If khops >= 0, the shortest ath is of distance <= khops,
A: the nxn connectivity matrix. Entries are assumed to be nonnegative. Absolute value will be taken if
int m = A->m, i, j, jj, *ia = A->ia, *ja = A->ja, heap_id;
BinaryHeap h;
- real *a = NULL, *aa;
+ double *a = NULL, *aa;
int *ai;
nodedata ndata, ndata_min;
enum {UNVISITED = -2, FINISHED = -1};
switch (A->type){
case MATRIX_TYPE_COMPLEX:
- aa = (real*) A->a;
- a = MALLOC(sizeof(real)*((size_t)(A->nz)));
+ aa = (double*) A->a;
+ a = MALLOC(sizeof(double)*((size_t)(A->nz)));
for (i = 0; i < A->nz; i++) a[i] = aa[i*2];
break;
case MATRIX_TYPE_REAL:
- a = (real*) A->a;
+ a = (double*) A->a;
break;
case MATRIX_TYPE_INTEGER:
ai = (int*) A->a;
- a = MALLOC(sizeof(real)*((size_t)(A->nz)));
- for (i = 0; i < A->nz; i++) a[i] = (real) ai[i];
+ a = MALLOC(sizeof(double)*((size_t)(A->nz)));
+ for (i = 0; i < A->nz; i++) a[i] = (double) ai[i];
break;
case MATRIX_TYPE_PATTERN:
- a = MALLOC(sizeof(real)*((size_t)A->nz));
+ a = MALLOC(sizeof(double)*((size_t)A->nz));
for (i = 0; i < A->nz; i++) a[i] = 1.;
break;
default:
}
}
-static int Dijkstra(SparseMatrix A, int root, real *dist, int *nlist, int *list, real *dist_max){
+static int Dijkstra(SparseMatrix A, int root, double *dist, int *nlist, int *list, double *dist_max){
return Dijkstra_internal(A, root, dist, nlist, list, dist_max, NULL);
}
-static int Dijkstra_masked(SparseMatrix A, int root, real *dist, int *nlist, int *list, real *dist_max, int *mask){
+static int Dijkstra_masked(SparseMatrix A, int root, double *dist, int *nlist, int *list, double *dist_max, int *mask){
/* this makes the algorithm only consider nodes that are masked.
nodes are masked as 1, 2, ..., mask_max, which is (the number of hops from root)+1.
Only paths consists of nodes that are masked are allowed. */
switch (A->type){
case MATRIX_TYPE_REAL:{
- real *a = (real*) A->a;
- real *val;
+ double *a = (double*) A->a;
+ double *val;
irn = MALLOC(sizeof(int)*((size_t)nz));
jcn = MALLOC(sizeof(int)*((size_t)nz));
- val = MALLOC(sizeof(real)*((size_t)nz));
+ val = MALLOC(sizeof(double)*((size_t)nz));
nz = 0;
for (i = 0; i < m; i++){
break;
}
case MATRIX_TYPE_COMPLEX:{
- real *a = (real*) A->a;
- real *val;
+ double *a = (double*) A->a;
+ double *val;
irn = MALLOC(sizeof(int)*((size_t)nz));
jcn = MALLOC(sizeof(int)*((size_t)nz));
- val = MALLOC(sizeof(real)*2*((size_t)nz));
+ val = MALLOC(sizeof(double)*2*((size_t)nz));
nz = 0;
for (i = 0; i < m; i++){
}
SparseMatrix SparseMatrix_set_entries_to_real_one(SparseMatrix A){
- real *a;
+ double *a;
int i;
free(A->a);
- A->a = MALLOC(sizeof(real)*((size_t)A->nz));
- a = (real*) (A->a);
+ A->a = MALLOC(sizeof(double)*((size_t)A->nz));
+ a = (double*) (A->a);
for (i = 0; i < A->nz; i++) a[i] = 1.;
A->type = MATRIX_TYPE_REAL;
- A->size = sizeof(real);
+ A->size = sizeof(double);
return A;
}
-SparseMatrix SparseMatrix_from_dense(int m, int n, real *x){
+SparseMatrix SparseMatrix_from_dense(int m, int n, double *x){
/* wrap a mxn matrix into a sparse matrix. the {i,j} entry of the matrix is in x[i*n+j], 0<=i<m; 0<=j<n */
int i, j, *ja;
- real *a;
+ double *a;
SparseMatrix A = SparseMatrix_new(m, n, m*n, MATRIX_TYPE_REAL, FORMAT_CSR);
A->ia[0] = 0;
for (i = 1; i <= m; i++) (A->ia)[i] = (A->ia)[i-1] + n;
ja = A->ja;
- a = (real*) A->a;
+ a = (double*) A->a;
for (i = 0; i < m; i++){
for (j = 0; j < n; j++) {
ja[j] = j;
}
-int SparseMatrix_distance_matrix(SparseMatrix D0, int weighted, real **dist0){
+int SparseMatrix_distance_matrix(SparseMatrix D0, int weighted, double **dist0){
/*
Input:
D: the graph. If weighted, the entry values is used.
SparseMatrix D = D0;
int m = D->m, n = D->n;
int *levelset_ptr = NULL, *levelset = NULL, *mask = NULL;
- real *dist = NULL;
+ double *dist = NULL;
int nlist, *list = NULL;
int flag = 0, i, j, k, nlevel;
- real dmax;
+ double dmax;
if (!SparseMatrix_is_symmetric(D, FALSE)){
D = SparseMatrix_symmetrize(D, FALSE);
assert(m == n);
- if (!(*dist0)) *dist0 = MALLOC(sizeof(real)*n*n);
+ if (!(*dist0)) *dist0 = MALLOC(sizeof(double)*n*n);
for (i = 0; i < n*n; i++) (*dist0)[i] = -1;
if (!weighted){
SparseMatrix D = D0, B, C;
int m = D->m, n = D->n;
int *levelset_ptr = NULL, *levelset = NULL, *mask = NULL;
- real *dist = NULL;
+ double *dist = NULL;
int nlist, *list = NULL;
int flag = 0, i, j, k, itmp, nlevel;
- real dmax, dtmp;
+ double dmax, dtmp;
if (!SparseMatrix_is_symmetric(D, FALSE)){
D = SparseMatrix_symmetrize(D, FALSE);
}
} else {
list = MALLOC(sizeof(int)*n);
- dist = MALLOC(sizeof(real)*n);
+ dist = MALLOC(sizeof(double)*n);
/*
Dijkstra_khops(khops, D, 60, dist, &nlist, list, &dmax);
for (j = 0; j < nlist; j++){
SparseMatrix SparseMatrix_transpose(SparseMatrix A);
SparseMatrix SparseMatrix_symmetrize(SparseMatrix A, int pattern_symmetric_only);
SparseMatrix SparseMatrix_symmetrize_nodiag(SparseMatrix A);
-void SparseMatrix_multiply_vector(SparseMatrix A, real *v, real **res, int transposed);/* if v = NULL, v is assumed to be {1,1,...,1}*/
+void SparseMatrix_multiply_vector(SparseMatrix A, double *v, double **res, int transposed);/* if v = NULL, v is assumed to be {1,1,...,1}*/
SparseMatrix SparseMatrix_remove_diagonal(SparseMatrix A);
SparseMatrix SparseMatrix_remove_upper(SparseMatrix A);/* remove diag and upper diag */
SparseMatrix SparseMatrix_divide_row_by_degree(SparseMatrix A);
SparseMatrix SparseMatrix_get_real_adjacency_matrix_symmetrized(SparseMatrix A); /* symmetric, all entries to 1, diaginal removed */
-void SparseMatrix_multiply_dense(SparseMatrix A, real *v, real **res, int dim);
+void SparseMatrix_multiply_dense(SparseMatrix A, double *v, double **res, int dim);
SparseMatrix SparseMatrix_apply_fun(SparseMatrix A, double (*fun)(double x));/* for real only! */
SparseMatrix SparseMatrix_copy(SparseMatrix A);
int SparseMatrix_has_diagonal(SparseMatrix A);
SparseMatrix SparseMatrix_set_entries_to_real_one(SparseMatrix A);
-int SparseMatrix_distance_matrix(SparseMatrix A, int weighted, real **dist_matrix);
+int SparseMatrix_distance_matrix(SparseMatrix A, int weighted, double **dist_matrix);
SparseMatrix SparseMatrix_distance_matrix_khops(int khops, SparseMatrix A, int weighted);
-SparseMatrix SparseMatrix_from_dense(int m, int n, real *x);
+SparseMatrix SparseMatrix_from_dense(int m, int n, double *x);
#define SparseMatrix_set_undirected(A) set_flag((A)->property, MATRIX_UNDIRECTED)
#define SparseMatrix_set_symmetric(A) set_flag((A)->property, MATRIX_SYMMETRIC)
grid->next = NULL;
grid->prev = NULL;
grid->delete_top_level_A = FALSE;
- grid->matching = MALLOC(sizeof(real)*(n));
+ grid->matching = MALLOC(sizeof(double)*(n));
grid->deg = NULL;
grid->agglomerate_regardless = FALSE;
if (level == 0){
- real modularity = 0;
+ double modularity = 0;
int *ia = A->ia, *ja = A->ja, n = A->n;
- real deg_total = 0;
- real *deg, *a = (real*) (A->a);
- real *indeg;
+ double deg_total = 0;
+ double *deg, *a = (double*) (A->a);
+ double *indeg;
grid->deg_total = 0.;
- grid->deg = MALLOC(sizeof(real)*(n));
+ grid->deg = MALLOC(sizeof(double)*(n));
deg = grid->deg;
- indeg = MALLOC(sizeof(real)*n);
+ indeg = MALLOC(sizeof(double)*n);
for (i = 0; i < n; i++){
deg[i] = 0;
indeg[i] = 0.;
int *matching = grid->matching;
SparseMatrix A = grid->A;
int n = grid->n, level = grid->level, nc = 0;
- real modularity = 0;
+ double modularity = 0;
int *ia = A->ia, *ja = A->ja;
- real *a;
- real *deg = grid->deg;
- real *deg_new;
+ double *a;
+ double *deg = grid->deg;
+ double *deg_new;
int i, j, jj, jc, jmax;
- real inv_deg_total = 1./(grid->deg_total);
- real *deg_inter, gain;
+ double inv_deg_total = 1./(grid->deg_total);
+ double *deg_inter, gain;
int *mask;
- real maxgain;
- real total_gain = 0;
+ double maxgain;
+ double total_gain = 0;
modularity = grid->modularity;
- deg_new = MALLOC(sizeof(real)*n);
- deg_inter = MALLOC(sizeof(real)*n);
+ deg_new = MALLOC(sizeof(double)*n);
+ deg_inter = MALLOC(sizeof(double)*n);
mask = MALLOC(sizeof(int)*n);
for (i = 0; i < n; i++) mask[i] = -1;
/* gain in merging node i into cluster j is
deg(i,j)/deg_total - 2*deg(i)*deg(j)/deg_total^2
*/
- a = (real*) A->a;
+ a = (double*) A->a;
for (i = 0; i < n; i++){
if (matching[i] != UNMATCHED) continue;
/* accumulate connections between i and clusters */
if (nc >= 1 && (total_gain > 0 || nc < n)){
/* now set up restriction and prolongation operator */
SparseMatrix P, R, R0, B, cA;
- real one = 1.;
+ double one = 1.;
Multilevel_Modularity_Clustering cgrid;
R0 = SparseMatrix_new(nc, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD);
grid->R = R;
level++;
cgrid = Multilevel_Modularity_Clustering_init(cA, level);
- deg_new = REALLOC(deg_new, nc*sizeof(real));
+ deg_new = REALLOC(deg_new, nc*sizeof(double));
cgrid->deg = deg_new;
cgrid->modularity = grid->modularity + total_gain;
cgrid->deg_total = grid->deg_total;
static void hierachical_modularity_clustering(SparseMatrix A, int ncluster_target,
- int *nclusters, int **assignment, real *modularity, int *flag){
+ int *nclusters, int **assignment, double *modularity, int *flag){
/* find a clustering of vertices by maximize modularity
A: symmetric square matrix n x n. If real value, value will be used as edges weights, otherwise edge weights are considered as 1.
Multilevel_Modularity_Clustering grid, cgrid;
int *matching, i;
SparseMatrix P;
- real *u;
+ double *u;
assert(A->m == A->n);
*modularity = 0.;
}
/* project clustering up */
- u = MALLOC(sizeof(real)*cgrid->n);
- for (i = 0; i < cgrid->n; i++) u[i] = (real) (cgrid->matching)[i];
+ u = MALLOC(sizeof(double)*cgrid->n);
+ for (i = 0; i < cgrid->n; i++) u[i] = (double) (cgrid->matching)[i];
*nclusters = cgrid->n;
*modularity = cgrid->modularity;
while (cgrid->prev){
- real *v = NULL;
+ double *v = NULL;
P = cgrid->prev->P;
SparseMatrix_multiply_vector(P, u, &v, FALSE);
free(u);
void modularity_clustering(SparseMatrix A, int inplace, int ncluster_target, int use_value,
- int *nclusters, int **assignment, real *modularity, int *flag){
+ int *nclusters, int **assignment, double *modularity, int *flag){
/* find a clustering of vertices by maximize modularity
A: symmetric square matrix n x n. If real value, value will be used as edges weights, otherwise edge weights are considered as 1.
inplace: whether A can e modified. If true, A will be modified by removing diagonal.
Multilevel_Modularity_Clustering prev;
int delete_top_level_A;
int *matching; /* dimension n. matching[i] is the clustering assignment of node i */
- real modularity;
- real deg_total; /* total edge weights, including self-edges */
- real *deg;/* dimension n. deg[i] equal to the sum of edge weights connected to vertex i. I.e., sum of row i */
+ double modularity;
+ double deg_total; /* total edge weights, including self-edges */
+ double *deg;/* dimension n. deg[i] equal to the sum of edge weights connected to vertex i. I.e., sum of row i */
int agglomerate_regardless;/* whether to agglomerate nodes even if this causes modularity reduction. This is used if we want to force
agglomeration so as to get less clusters
*/
modularity: achieve modularity
*/
void modularity_clustering(SparseMatrix A, int inplace, int maxcluster, int use_value,
- int *nclusters, int **assignment, real *modularity, int *flag);
+ int *nclusters, int **assignment, double *modularity, int *flag);
}
}
-static real Hue2RGB(real v1, real v2, real H) {
+static double Hue2RGB(double v1, double v2, double H) {
if(H < 0.0) H += 1.0;
if(H > 1.0) H -= 1.0;
if((6.0*H) < 1.0) return (v1 + (v2 - v1) * 6.0 * H);
return v1;
}
-char * hue2rgb(real hue, char *color){
- real v1, v2, lightness = .5, saturation = 1;
+char * hue2rgb(double hue, char *color){
+ double v1, v2, lightness = .5, saturation = 1;
int red, blue, green;
if(lightness < 0.5)
void rgb2hex(float r, float g, float b, char *cstring, const char* opacity);
/* dimension of cstring must be >=7 */
-char* hue2rgb(real hue, char *color);
+char* hue2rgb(double hue, char *color);
double _statistics[10];
#endif
-real drand(){
- return rand()/(real) RAND_MAX;
+double drand(){
+ return rand()/(double) RAND_MAX;
}
int irand(int n){
}
-real* vector_subtract_to(int n, real *x, real *y){
+double* vector_subtract_to(int n, double *x, double *y){
/* y = x-y */
int i;
for (i = 0; i < n; i++) y[i] = x[i] - y[i];
return y;
}
-real vector_product(int n, real *x, real *y){
- real res = 0;
+double vector_product(int n, double *x, double *y){
+ double res = 0;
int i;
for (i = 0; i < n; i++) res += x[i]*y[i];
return res;
}
-real* vector_saxpy(int n, real *x, real *y, real beta){
+double* vector_saxpy(int n, double *x, double *y, double beta){
/* y = x+beta*y */
int i;
for (i = 0; i < n; i++) y[i] = x[i] + beta*y[i];
return y;
}
-real* vector_saxpy2(int n, real *x, real *y, real beta){
+double* vector_saxpy2(int n, double *x, double *y, double beta){
/* x = x+beta*y */
int i;
for (i = 0; i < n; i++) x[i] = x[i] + beta*y[i];
return x;
}
-void vector_print(char *s, int n, real *x){
+void vector_print(char *s, int n, double *x){
int i;
printf("%s{",s);
for (i = 0; i < n; i++) {
}
static int comp_ascend(const void *s1, const void *s2){
- const real *ss1 = s1;
- const real *ss2 = s2;
+ const double *ss1 = s1;
+ const double *ss2 = s2;
if ((ss1)[0] > (ss2)[0]){
return 1;
return 0;
}
-void vector_ordering(int n, real *v, int **p){
+void vector_ordering(int n, double *v, int **p){
/* give the position of the smallest, second smallest etc in vector v.
results in p. If *p == NULL, p is assigned.
*/
- real *u;
+ double *u;
int i;
if (!*p) *p = MALLOC(sizeof(int)*n);
- u = MALLOC(sizeof(real)*2*n);
+ u = MALLOC(sizeof(double)*2*n);
for (i = 0; i < n; i++) {
u[2*i+1] = i;
u[2*i] = v[i];
}
- qsort(u, n, sizeof(real)*2, comp_ascend);
+ qsort(u, n, sizeof(double)*2, comp_ascend);
for (i = 0; i < n; i++) (*p)[i] = (int) u[2*i+1];
free(u);
qsort(v, n, sizeof(int), comp_ascend_int);
}
-real distance_cropped(real *x, int dim, int i, int j){
+double distance_cropped(double *x, int dim, int i, int j){
int k;
- real dist = 0.;
+ double dist = 0.;
for (k = 0; k < dim; k++) dist += (x[i*dim+k] - x[j*dim + k])*(x[i*dim+k] - x[j*dim + k]);
dist = sqrt(dist);
return MAX(dist, MINDIST);
}
-real distance(real *x, int dim, int i, int j){
+double distance(double *x, int dim, int i, int j){
int k;
- real dist = 0.;
+ double dist = 0.;
for (k = 0; k < dim; k++) dist += (x[i*dim+k] - x[j*dim + k])*(x[i*dim+k] - x[j*dim + k]);
dist = sqrt(dist);
return dist;
}
-real point_distance(real *p1, real *p2, int dim){
+double point_distance(double *p1, double *p2, int dim){
int i;
- real dist;
+ double dist;
dist = 0;
for (i = 0; i < dim; i++) dist += (p1[i] - p2[i])*(p1[i] - p2[i]);
return sqrt(dist);
return s;
}
-void scale_to_box(real xmin, real ymin, real xmax, real ymax, int n, int dim, real *x){
- real min[3], max[3], min0[3], ratio = 1;
+void scale_to_box(double xmin, double ymin, double xmax, double ymax, int n, int dim, double *x){
+ double min[3], max[3], min0[3], ratio = 1;
int i, k;
for (i = 0; i < dim; i++) {
#include <common/memory.h>
#endif /* STANDALONE */
-#define real double
-
#define set_flag(a, flag) ((a)=((a)|(flag)))
#define test_flag(a, flag) ((a)&(flag))
#define clear_flag(a, flag) ((a) &=(~(flag)))
extern int irand(int n);
-extern real drand(void);
+extern double drand(void);
extern int *random_permutation(int n);/* random permutation of 0 to n-1 */
-real* vector_subtract_to(int n, real *x, real *y);/* y = x-y */
+double* vector_subtract_to(int n, double *x, double *y);/* y = x-y */
-real vector_product(int n, real *x, real *y);
+double vector_product(int n, double *x, double *y);
-real* vector_saxpy(int n, real *x, real *y, real beta); /* y = x+beta*y */
+double* vector_saxpy(int n, double *x, double *y, double beta); /* y = x+beta*y */
-real* vector_saxpy2(int n, real *x, real *y, real beta);/* x = x+beta*y */
+double* vector_saxpy2(int n, double *x, double *y, double beta);/* x = x+beta*y */
/* take m elements v[p[i]]],i=1,...,m and oput in u. u will be assigned if *u = NULL */
void vector_float_take(int n, float *v, int m, int *p, float **u);
/* give the position of the smallest, second smallest etc in vector v.
results in p. If *p == NULL, p is assigned.
*/
-void vector_ordering(int n, real *v, int **p);
+void vector_ordering(int n, double *v, int **p);
void vector_sort_int(int n, int *v);
-void vector_print(char *s, int n, real *x);
+void vector_print(char *s, int n, double *x);
#define MACHINEACC 1.0e-16
#define SQRT_MACHINEACC 1.0e-8
enum {UNMATCHED = -1};
-real distance(real *x, int dim, int i, int j);
-real distance_cropped(real *x, int dim, int i, int j);
+double distance(double *x, int dim, int i, int j);
+double distance_cropped(double *x, int dim, int i, int j);
-real point_distance(real *p1, real *p2, int dim);
+double point_distance(double *p1, double *p2, int dim);
char *strip_dir(char *s);
-void scale_to_box(real xmin, real ymin, real xmax, real ymax, int n, int dim, real *x);
+void scale_to_box(double xmin, double ymin, double xmax, double ymax, int n, int dim, double *x);
#include <sparse/LinkedList.h>
#include <string.h>
-static real get_mq(SparseMatrix A, int *assignment, int *ncluster0, real *mq_in0, real *mq_out0, real **dout0){
+static double get_mq(SparseMatrix A, int *assignment, int *ncluster0, double *mq_in0, double *mq_out0, double **dout0){
/* given a symmetric matrix representation of a graph and an assignment of nodes into clusters, calculate the modularity quality.
assignment: assignmenet[i] gives the cluster assignment of node i. 0 <= assignment[i] < ncluster.
ncluster: number of clusters
int n = A->m;
int test_pattern_symmetry_only = FALSE;
int *counts, *ia = A->ia, *ja = A->ja, k, i, j, jj;
- real mq_in = 0, mq_out = 0, *a = NULL, Vi, Vj;
+ double mq_in = 0, mq_out = 0, *a = NULL, Vi, Vj;
int c;
- real *dout;
+ double *dout;
assert(SparseMatrix_is_symmetric(A, test_pattern_symmetry_only));
assert(A->n == n);
- if (A->type == MATRIX_TYPE_REAL) a = (real*) A->a;
+ if (A->type == MATRIX_TYPE_REAL) a = (double*) A->a;
counts = CALLOC(n, sizeof(int));
}
/* calculate scaled out degree */
- dout = MALLOC(sizeof(real)*n);
+ dout = MALLOC(sizeof(double)*n);
for (i = 0; i < n; i++){
dout[i] = 0;
for (j = ia[i]; j < ia[i+1]; j++){
jj = ja[j];
if (jj == i) continue;
if (a){
- dout[i] += a[j]/(real) counts[assignment[jj]];
+ dout[i] += a[j]/(double) counts[assignment[jj]];
} else {
- dout[i] += 1./(real) counts[assignment[jj]];
+ dout[i] += 1./(double) counts[assignment[jj]];
}
}
}
grid->next = NULL;
grid->prev = NULL;
grid->delete_top_level_A = FALSE;
- matching = grid->matching = MALLOC(sizeof(real)*(n));
+ matching = grid->matching = MALLOC(sizeof(double)*(n));
grid->deg_intra = NULL;
grid->dout = NULL;
grid->wgt = NULL;
if (level == 0){
- real mq = 0, mq_in, mq_out;
+ double mq = 0, mq_in, mq_out;
int n = A->n, ncluster;
- real *deg_intra, *wgt, *dout;
+ double *deg_intra, *wgt, *dout;
- grid->deg_intra = MALLOC(sizeof(real)*(n));
+ grid->deg_intra = MALLOC(sizeof(double)*(n));
deg_intra = grid->deg_intra;
- grid->wgt = MALLOC(sizeof(real)*n);
+ grid->wgt = MALLOC(sizeof(double)*n);
wgt = grid->wgt;
for (i = 0; i < n; i++){
int *matching = grid->matching;
SparseMatrix A = grid->A;
int n = grid->n, level = grid->level, nc = 0, nclusters = n;
- real mq = 0, mq_in = 0, mq_out = 0, mq_new, mq_in_new, mq_out_new, mq_max = 0, mq_in_max = 0, mq_out_max = 0;
+ double mq = 0, mq_in = 0, mq_out = 0, mq_new, mq_in_new, mq_out_new, mq_max = 0, mq_in_max = 0, mq_out_max = 0;
int *ia = A->ia, *ja = A->ja;
- real *a, amax = 0;
- real *deg_intra = grid->deg_intra, *wgt = grid->wgt;
- real *deg_intra_new, *wgt_new = NULL;
+ double *a, amax = 0;
+ double *deg_intra = grid->deg_intra, *wgt = grid->wgt;
+ double *deg_intra_new, *wgt_new = NULL;
int i, j, k, jj, jc, jmax;
- real *deg_inter, gain = 0, *dout = grid->dout, *dout_new, deg_in_i, deg_in_j, wgt_i, wgt_j, a_ij, dout_i, dout_j, dout_max = 0, wgt_jmax = 0;
+ double *deg_inter, gain = 0, *dout = grid->dout, *dout_new, deg_in_i, deg_in_j, wgt_i, wgt_j, a_ij, dout_i, dout_j, dout_max = 0, wgt_jmax = 0;
int *mask;
- real maxgain = 0;
- real total_gain = 0;
+ double maxgain = 0;
+ double total_gain = 0;
SingleLinkedList *neighbors = NULL, lst;
mq_in = grid->mq_in;
mq_out = grid->mq_out;
- deg_intra_new = MALLOC(sizeof(real)*n);
- wgt_new = MALLOC(sizeof(real)*n);
- deg_inter = MALLOC(sizeof(real)*n);
+ deg_intra_new = MALLOC(sizeof(double)*n);
+ wgt_new = MALLOC(sizeof(double)*n);
+ deg_inter = MALLOC(sizeof(double)*n);
mask = MALLOC(sizeof(int)*n);
- dout_new = MALLOC(sizeof(real)*n);
+ dout_new = MALLOC(sizeof(double)*n);
for (i = 0; i < n; i++) mask[i] = -1;
assert(n == A->n);
mq_new = mq_in_new/(k-1) - mq_out_new/((k-1)*(k-2))
gain = mq_new - mq
*/
- a = (real*) A->a;
+ a = (double*) A->a;
for (i = 0; i < n; i++){
if (matching[i] != UNMATCHED) continue;
/* accumulate connections between i and clusters */
double mq2, mq_in2, mq_out2, *dout2;
int *matching2, nc2 = nc;
matching2 = MALLOC(sizeof(int)*A->m);
- memcpy(matching2, matching, sizeof(real)*A->m);
+ memcpy(matching2, matching, sizeof(double)*A->m);
if (jc != UNMATCHED) {
matching2[i] = jc;
} else {
if (nc >= 1 && (total_gain > 0 || nc < n)){
/* now set up restriction and prolongation operator */
SparseMatrix P, R, R0, B, cA;
- real one = 1.;
+ double one = 1.;
Multilevel_MQ_Clustering cgrid;
R0 = SparseMatrix_new(nc, n, 1, MATRIX_TYPE_REAL, FORMAT_COORD);
grid->R = R;
level++;
cgrid = Multilevel_MQ_Clustering_init(cA, level);
- deg_intra_new = REALLOC(deg_intra_new, nc*sizeof(real));
- wgt_new = REALLOC(wgt_new, nc*sizeof(real));
+ deg_intra_new = REALLOC(deg_intra_new, nc*sizeof(double));
+ wgt_new = REALLOC(wgt_new, nc*sizeof(double));
cgrid->deg_intra = deg_intra_new;
cgrid->mq = grid->mq + total_gain;
cgrid->wgt = wgt_new;
- dout_new = REALLOC(dout_new, nc*sizeof(real));
+ dout_new = REALLOC(dout_new, nc*sizeof(double));
cgrid->dout = dout_new;
cgrid = Multilevel_MQ_Clustering_establish(cgrid, maxcluster);
static void hierachical_mq_clustering(SparseMatrix A, int maxcluster,
- int *nclusters, int **assignment, real *mq, int *flag){
+ int *nclusters, int **assignment, double *mq, int *flag){
/* find a clustering of vertices by maximize mq
A: symmetric square matrix n x n. If real value, value will be used as edges weights, otherwise edge weights are considered as 1.
maxcluster: used to specify the maximum number of cluster desired, e.g., maxcluster=10 means that a maximum of 10 clusters
Multilevel_MQ_Clustering grid, cgrid;
int *matching, i;
SparseMatrix P;
- real *u;
+ double *u;
assert(A->m == A->n);
*mq = 0.;
}
/* project clustering up */
- u = MALLOC(sizeof(real)*cgrid->n);
- for (i = 0; i < cgrid->n; i++) u[i] = (real) (cgrid->matching)[i];
+ u = MALLOC(sizeof(double)*cgrid->n);
+ for (i = 0; i < cgrid->n; i++) u[i] = (double) (cgrid->matching)[i];
*nclusters = cgrid->n;
*mq = cgrid->mq;
while (cgrid->prev){
- real *v = NULL;
+ double *v = NULL;
P = cgrid->prev->P;
SparseMatrix_multiply_vector(P, u, &v, FALSE);
free(u);
void mq_clustering(SparseMatrix A, int inplace, int maxcluster, int use_value,
- int *nclusters, int **assignment, real *mq, int *flag){
+ int *nclusters, int **assignment, double *mq, int *flag){
/* find a clustering of vertices by maximize mq
A: symmetric square matrix n x n. If real value, value will be used as edges weights, otherwise edge weights are considered as 1.
inplace: whether A can e modified. If true, A will be modified by removing diagonal.
. = mq_in/k - mq_out/(k*(k-1))
*/
- real mq;
- real mq_in, mq_out;/* mqs(A) = deg_in(A)/|A|^2 - deg_out(A)/|A|/(|V|-|A|) */
+ double mq;
+ double mq_in, mq_out;/* mqs(A) = deg_in(A)/|A|^2 - deg_out(A)/|A|/(|V|-|A|) */
int ncluster; /* number of clusters */
- real *deg_intra;/* dimension n. deg[i] equal to the sum of edge weights within cluster i */
- real *dout;/* dimension n, dout[i] = \sum_{j -- i} a(i,j)/|j| is the scaled sum of outdegree */
- real *wgt; /* total vertex weight each coarse grid vertex represent */
+ double *deg_intra;/* dimension n. deg[i] equal to the sum of edge weights within cluster i */
+ double *dout;/* dimension n, dout[i] = \sum_{j -- i} a(i,j)/|j| is the scaled sum of outdegree */
+ double *wgt; /* total vertex weight each coarse grid vertex represent */
};
/* find a clustering of vertices by maximize modularity quality
mq: achieve modularity
*/
void mq_clustering(SparseMatrix A, int inplace, int maxcluster, int use_value,
- int *nclusters, int **assignment, real *mq, int *flag);
+ int *nclusters, int **assignment, double *mq, int *flag);
projects / graphviz / commitdiff