--- ray/src/cv/bsdfmesh.c	2012/10/19 04:14:29	2.1
+++ ray/src/cv/bsdfmesh.c	2012/10/20 07:02:00	2.2
@@ -1,5 +1,5 @@
 #ifndef lint
-static const char RCSid[] = "$Id: bsdfmesh.c,v 2.1 2012/10/19 04:14:29 greg Exp $";
+static const char RCSid[] = "$Id: bsdfmesh.c,v 2.2 2012/10/20 07:02:00 greg Exp $";
 #endif
 /*
  * Create BSDF advection mesh from radial basis functions.
@@ -23,6 +23,107 @@ int			nprocs = 1;
 				/* number of children (-1 in child) */
 static int		nchild = 0;
 
+/* Create a new migration holder (sharing memory for multiprocessing) */
+static MIGRATION *
+new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
+{
+	size_t		memlen = sizeof(MIGRATION) +
+				sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1);
+	MIGRATION	*newmig;
+#ifdef _WIN32
+	if (nprocs > 1)
+		fprintf(stderr, "%s: warning - multiprocessing not supported\n",
+				progname);
+	nprocs = 1;
+	newmig = (MIGRATION *)malloc(memlen);
+#else
+	if (nprocs <= 1) {			/* single process? */
+		newmig = (MIGRATION *)malloc(memlen);
+	} else {				/* else need to share memory */
+		newmig = (MIGRATION *)mmap(NULL, memlen, PROT_READ|PROT_WRITE,
+						MAP_ANON|MAP_SHARED, -1, 0);
+		if ((void *)newmig == MAP_FAILED)
+			newmig = NULL;
+	}
+#endif
+	if (newmig == NULL) {
+		fprintf(stderr, "%s: cannot allocate new migration\n", progname);
+		exit(1);
+	}
+	newmig->rbfv[0] = from_rbf;
+	newmig->rbfv[1] = to_rbf;
+						/* insert in edge lists */
+	newmig->enxt[0] = from_rbf->ejl;
+	from_rbf->ejl = newmig;
+	newmig->enxt[1] = to_rbf->ejl;
+	to_rbf->ejl = newmig;
+	newmig->next = mig_list;		/* push onto global list */
+	return(mig_list = newmig);
+}
+
+#ifdef _WIN32
+#define await_children(n)	(void)(n)
+#define run_subprocess()	0
+#define end_subprocess()	(void)0
+#else
+
+/* Wait for the specified number of child processes to complete */
+static void
+await_children(int n)
+{
+	int	exit_status = 0;
+
+	if (n > nchild)
+		n = nchild;
+	while (n-- > 0) {
+		int	status;
+		if (wait(&status) < 0) {
+			fprintf(stderr, "%s: missing child(ren)!\n", progname);
+			nchild = 0;
+			break;
+		}
+		--nchild;
+		if (status) {			/* something wrong */
+			if ((status = WEXITSTATUS(status)))
+				exit_status = status;
+			else
+				exit_status += !exit_status;
+			fprintf(stderr, "%s: subprocess died\n", progname);
+			n = nchild;		/* wait for the rest */
+		}
+	}
+	if (exit_status)
+		exit(exit_status);
+}
+
+/* Start child process if multiprocessing selected */
+static pid_t
+run_subprocess(void)
+{
+	int	status;
+	pid_t	pid;
+
+	if (nprocs <= 1)			/* any children requested? */
+		return(0);
+	await_children(nchild + 1 - nprocs);	/* free up child process */
+	if ((pid = fork())) {
+		if (pid < 0) {
+			fprintf(stderr, "%s: cannot fork subprocess\n",
+					progname);
+			exit(1);
+		}
+		++nchild;			/* subprocess started */
+		return(pid);
+	}
+	nchild = -1;
+	return(0);				/* put child to work */
+}
+
+/* If we are in subprocess, call exit */
+#define	end_subprocess()	if (nchild < 0) _exit(0); else
+
+#endif	/* ! _WIN32 */
+
 /* Compute (and allocate) migration price matrix for optimization */
 static float *
 price_routes(const RBFNODE *from_rbf, const RBFNODE *to_rbf)
@@ -106,7 +207,7 @@ static double
 migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const float *pmtx)
 {
 	const double	maxamt = .1;
-	const double	minamt = maxamt*.0001;
+	const double	minamt = maxamt*5e-6;
 	static double	*src_cost = NULL;
 	static int	n_alloc = 0;
 	struct {
@@ -136,7 +237,7 @@ migration_step(MIGRATION *mig, double *src_rem, double
 	for (cur.s = mtx_nrows(mig); cur.s--; ) {
 	    const float	*price = pmtx + cur.s*mtx_ncols(mig);
 	    double	cost_others = 0;
-	    if (src_rem[cur.s] < minamt)
+	    if (src_rem[cur.s] <= minamt)
 		    continue;
 	    cur.d = -1;				/* examine cheapest dest. */
 	    for (i = mtx_ncols(mig); i--; )
@@ -164,7 +265,7 @@ migration_step(MIGRATION *mig, double *src_rem, double
 	if ((best.s < 0) | (best.d < 0))
 		return(.0);
 						/* make the actual move */
-	mig->mtx[mtx_ndx(mig,best.s,best.d)] += best.amt;
+	mtx_coef(mig,best.s,best.d) += best.amt;
 	src_rem[best.s] -= best.amt;
 	dst_rem[best.d] -= best.amt;
 	return(best.amt);
@@ -185,114 +286,11 @@ thetaphi(const FVECT v)
 }
 #endif
 
-/* Create a new migration holder (sharing memory for multiprocessing) */
-static MIGRATION *
-new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
-{
-	size_t		memlen = sizeof(MIGRATION) +
-				sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1);
-	MIGRATION	*newmig;
-#ifdef _WIN32
-	if (nprocs > 1)
-		fprintf(stderr, "%s: warning - multiprocessing not supported\n",
-				progname);
-	nprocs = 1;
-	newmig = (MIGRATION *)malloc(memlen);
-#else
-	if (nprocs <= 1) {			/* single process? */
-		newmig = (MIGRATION *)malloc(memlen);
-	} else {				/* else need to share memory */
-		newmig = (MIGRATION *)mmap(NULL, memlen, PROT_READ|PROT_WRITE,
-						MAP_ANON|MAP_SHARED, -1, 0);
-		if ((void *)newmig == MAP_FAILED)
-			newmig = NULL;
-	}
-#endif
-	if (newmig == NULL) {
-		fprintf(stderr, "%s: cannot allocate new migration\n", progname);
-		exit(1);
-	}
-	newmig->rbfv[0] = from_rbf;
-	newmig->rbfv[1] = to_rbf;
-						/* insert in edge lists */
-	newmig->enxt[0] = from_rbf->ejl;
-	from_rbf->ejl = newmig;
-	newmig->enxt[1] = to_rbf->ejl;
-	to_rbf->ejl = newmig;
-	newmig->next = mig_list;		/* push onto global list */
-	return(mig_list = newmig);
-}
-
-#ifdef _WIN32
-#define await_children(n)	(void)(n)
-#define run_subprocess()	0
-#define end_subprocess()	(void)0
-#else
-
-/* Wait for the specified number of child processes to complete */
-static void
-await_children(int n)
-{
-	int	exit_status = 0;
-
-	if (n > nchild)
-		n = nchild;
-	while (n-- > 0) {
-		int	status;
-		if (wait(&status) < 0) {
-			fprintf(stderr, "%s: missing child(ren)!\n", progname);
-			nchild = 0;
-			break;
-		}
-		--nchild;
-		if (status) {			/* something wrong */
-			if ((status = WEXITSTATUS(status)))
-				exit_status = status;
-			else
-				exit_status += !exit_status;
-			fprintf(stderr, "%s: subprocess died\n", progname);
-			n = nchild;		/* wait for the rest */
-		}
-	}
-	if (exit_status)
-		exit(exit_status);
-}
-
-/* Start child process if multiprocessing selected */
-static pid_t
-run_subprocess(void)
-{
-	int	status;
-	pid_t	pid;
-
-	if (nprocs <= 1)			/* any children requested? */
-		return(0);
-	await_children(nchild + 1 - nprocs);	/* free up child process */
-	if ((pid = fork())) {
-		if (pid < 0) {
-			fprintf(stderr, "%s: cannot fork subprocess\n",
-					progname);
-			exit(1);
-		}
-		++nchild;			/* subprocess started */
-		return(pid);
-	}
-	nchild = -1;
-	return(0);				/* put child to work */
-}
-
-/* If we are in subprocess, call exit */
-#define	end_subprocess()	if (nchild < 0) _exit(0); else
-
-#endif	/* ! _WIN32 */
-
 /* Compute and insert migration along directed edge (may fork child) */
 static MIGRATION *
 create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
 {
-	const double	end_thresh = 0.1/(from_rbf->nrbf*to_rbf->nrbf);
-	const double	check_thresh = 0.01;
-	const double	rel_thresh = 5e-6;
+	const double	end_thresh = 5e-6;
 	float		*pmtx;
 	MIGRATION	*newmig;
 	double		*src_rem, *dst_rem;
@@ -318,8 +316,9 @@ create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
 #ifdef DEBUG
 	fprintf(stderr, "Building path from (theta,phi) %s ",
 			thetaphi(from_rbf->invec));
-	fprintf(stderr, "to %s", thetaphi(to_rbf->invec));
-	/* if (nchild) */ fputc('\n', stderr);
+	fprintf(stderr, "to %s with %d x %d matrix\n",
+			thetaphi(to_rbf->invec), 
+			from_rbf->nrbf, to_rbf->nrbf);
 #endif
 						/* starting quantities */
 	memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
@@ -332,13 +331,11 @@ create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
 		total_rem -= move_amt;
 #ifdef DEBUG
 		if (!nchild)
-			/* fputc('.', stderr); */
-			fprintf(stderr, "%.9f remaining...\r", total_rem);
+			fprintf(stderr, "\r%.9f remaining...", total_rem);
 #endif
-	} while (total_rem > end_thresh && (total_rem > check_thresh) |
-					(move_amt > rel_thresh*total_rem));
+	} while ((total_rem > end_thresh) & (move_amt > 0));
 #ifdef DEBUG
-	if (!nchild) fputs("\ndone.\n", stderr);
+	if (!nchild) fputs("done.\n", stderr);
 	else fprintf(stderr, "finished with %.9f remaining\n", total_rem);
 #endif
 	for (i = from_rbf->nrbf; i--; ) {	/* normalize final matrix */
@@ -347,7 +344,7 @@ create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
 	    if (nf <= FTINY) continue;
 	    nf = from_rbf->vtotal / nf;
 	    for (j = to_rbf->nrbf; j--; )
-		newmig->mtx[mtx_ndx(newmig,i,j)] *= nf;
+		mtx_coef(newmig,i,j) *= nf;
 	}
 	end_subprocess();			/* exit here if subprocess */
 	free(pmtx);				/* free working arrays */