ViewVC Help

View File | Revision Log | Show Annotations | Download File | Root Listing

root/Development/ray/src/cv/bsdfmesh.c

Comparing ray/src/cv/bsdfmesh.c (file contents):
Revision 2.2 by greg, Sat Oct 20 07:02:00 2012 UTC vs.
Revision 2.11 by greg, Thu Oct 24 16:11:37 2013 UTC

#	Line 23 | Line 23 | int                    nprocs = 1;
23		/* number of children (-1 in child) */
24		static int              nchild = 0;
25
26	+	typedef struct {
27	+	int             nrows, ncols;   /* array size (matches migration) */
28	+	float           *price;         /* migration prices */
29	+	short           *sord;          /* sort for each row, low to high */
30	+	float           *prow;          /* current price row */
31	+	} PRICEMAT;                     /* sorted pricing matrix */
32	+
33	+	#define pricerow(p,i)   ((p)->price + (i)*(p)->ncols)
34	+	#define psortrow(p,i)   ((p)->sord + (i)*(p)->ncols)
35	+
36		/* Create a new migration holder (sharing memory for multiprocessing) */
37		static MIGRATION *
38		new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
#	Line 110 | Line 120 | run_subprocess(void)
120		if (pid < 0) {
121		fprintf(stderr, "%s: cannot fork subprocess\n",
122		progname);
123	+	await_children(nchild);
124		exit(1);
125		}
126		++nchild;                       /* subprocess started */
#	Line 124 | Line 135 | run_subprocess(void)
135
136		#endif  /* ! _WIN32 */
137
138	+	/* Comparison routine needed for sorting price row */
139	+	static int
140	+	msrt_cmp(void *b, const void *p1, const void *p2)
141	+	{
142	+	PRICEMAT        *pm = (PRICEMAT *)b;
143	+	float           c1 = pm->prow[*(const short *)p1];
144	+	float           c2 = pm->prow[*(const short *)p2];
145	+
146	+	if (c1 > c2) return(1);
147	+	if (c1 < c2) return(-1);
148	+	return(0);
149	+	}
150	+
151		/* Compute (and allocate) migration price matrix for optimization */
152	<	static float *
153	<	price_routes(const RBFNODE *from_rbf, const RBFNODE *to_rbf)
152	>	static void
153	>	price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
154		{
131	–	float   *pmtx = (float *)malloc(sizeof(float) *
132	–	from_rbf->nrbf * to_rbf->nrbf);
155		FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
156		int     i, j;
157
158	<	if ((pmtx == NULL) | (vto == NULL)) {
158	>	pm->nrows = from_rbf->nrbf;
159	>	pm->ncols = to_rbf->nrbf;
160	>	pm->price = (float *)malloc(sizeof(float) * pm->nrows*pm->ncols);
161	>	pm->sord = (short *)malloc(sizeof(short) * pm->nrows*pm->ncols);
162	>
163	>	if ((pm->price == NULL) | (pm->sord == NULL) | (vto == NULL)) {
164		fprintf(stderr, "%s: Out of memory in migration_costs()\n",
165		progname);
166		exit(1);
#	Line 144 | Line 171 | price_routes(const RBFNODE *from_rbf, const RBFNODE *t
171		for (i = from_rbf->nrbf; i--; ) {
172		const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
173		FVECT               vfrom;
174	+	short               *srow;
175		ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
176	<	for (j = to_rbf->nrbf; j--; )
177	<	pmtx[i*to_rbf->nrbf + j] = acos(DOT(vfrom, vto[j])) +
176	>	pm->prow = pricerow(pm,i);
177	>	srow = psortrow(pm,i);
178	>	for (j = to_rbf->nrbf; j--; ) {
179	>	double          dprod = DOT(vfrom, vto[j]);
180	>	pm->prow[j] = ((dprod >= 1.) ? .0 : acos(dprod)) +
181		fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
182	+	srow[j] = j;
183	+	}
184	+	qsort_r(srow, pm->ncols, sizeof(short), pm, &msrt_cmp);
185		}
186		free(vto);
153	–	return(pmtx);
187		}
188
189	<	/* Comparison routine needed for sorting price row */
190	<	static const float      *price_arr;
191	<	static int
159	<	msrt_cmp(const void *p1, const void *p2)
189	>	/* Free price matrix */
190	>	static void
191	>	free_routes(PRICEMAT *pm)
192		{
193	<	float   c1 = price_arr[*(const int *)p1];
194	<	float   c2 = price_arr[*(const int *)p2];
163	<
164	<	if (c1 > c2) return(1);
165	<	if (c1 < c2) return(-1);
166	<	return(0);
193	>	free(pm->price); pm->price = NULL;
194	>	free(pm->sord); pm->sord = NULL;
195		}
196
197		/* Compute minimum (optimistic) cost for moving the given source material */
198		static double
199	<	min_cost(double amt2move, const double *avail, const float *price, int n)
199	>	min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
200		{
201	<	static int      *price_sort = NULL;
202	<	static int      n_alloc = 0;
201	>	const short     *srow = psortrow(pm,s);
202	>	const float     *prow = pricerow(pm,s);
203		double          total_cost = 0;
204	<	int             i;
177	<
178	<	if (amt2move <= FTINY)                  /* pre-emptive check */
179	<	return(0.);
180	<	if (n > n_alloc) {                      /* (re)allocate sort array */
181	<	if (n_alloc) free(price_sort);
182	<	price_sort = (int *)malloc(sizeof(int)*n);
183	<	if (price_sort == NULL) {
184	<	fprintf(stderr, "%s: Out of memory in min_cost()\n",
185	<	progname);
186	<	exit(1);
187	<	}
188	<	n_alloc = n;
189	<	}
190	<	for (i = n; i--; )
191	<	price_sort[i] = i;
192	<	price_arr = price;
193	<	qsort(price_sort, n, sizeof(int), &msrt_cmp);
204	>	int             j;
205		/* move cheapest first */
206	<	for (i = 0; i < n && amt2move > FTINY; i++) {
207	<	int     d = price_sort[i];
206	>	for (j = 0; (j < pm->ncols) & (amt2move > FTINY); j++) {
207	>	int     d = srow[j];
208		double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
209
210	<	total_cost += amt * price[d];
210	>	total_cost += amt * prow[d];
211		amt2move -= amt;
212		}
213		return(total_cost);
214		}
215
216	<	/* Take a step in migration by choosing optimal bucket to transfer */
216	>	/* Compare entries by moving price */
217	>	static int
218	>	rmovcmp(void *b, const void *p1, const void *p2)
219	>	{
220	>	PRICEMAT        *pm = (PRICEMAT *)b;
221	>	const short     *ij1 = (const short *)p1;
222	>	const short     *ij2 = (const short *)p2;
223	>	float           price_diff;
224	>
225	>	if (ij1[1] < 0) return(ij2[1] >= 0);
226	>	if (ij2[1] < 0) return(-1);
227	>	price_diff = pricerow(pm,ij1[0])[ij1[1]] - pricerow(pm,ij2[0])[ij2[1]];
228	>	if (price_diff > 0) return(1);
229	>	if (price_diff < 0) return(-1);
230	>	return(0);
231	>	}
232	>
233	>	/* Take a step in migration by choosing reasonable bucket to transfer */
234		static double
235	<	migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const float *pmtx)
235	>	migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, PRICEMAT *pm)
236		{
237	<	const double    maxamt = .1;
237	>	const int       max2check = 100;
238	>	const double    maxamt = 1./(double)pm->ncols;
239		const double    minamt = maxamt*5e-6;
240	<	static double   *src_cost = NULL;
241	<	static int      n_alloc = 0;
240	>	double          *src_cost;
241	>	short           (*rord)[2];
242		struct {
243		int     s, d;   /* source and destination */
244		double  price;  /* price estimate per amount moved */
245		double  amt;    /* amount we can move */
246		} cur, best;
247	<	int             i;
248	<
249	<	if (mtx_nrows(mig) > n_alloc) {         /* allocate cost array */
250	<	if (n_alloc)
251	<	free(src_cost);
252	<	src_cost = (double *)malloc(sizeof(double)*mtx_nrows(mig));
253	<	if (src_cost == NULL) {
254	<	fprintf(stderr, "%s: Out of memory in migration_step()\n",
255	<	progname);
256	<	exit(1);
257	<	}
258	<	n_alloc = mtx_nrows(mig);
247	>	int             r2check, i, ri;
248	>	/*
249	>	* Check cheapest available routes only -- a higher adjusted
250	>	* destination price implies that another source is closer, so
251	>	* we can hold off considering more expensive options until
252	>	* some other (hopefully better) moves have been made.
253	>	*/
254	>	/* most promising row order */
255	>	rord = (short (*)[2])malloc(sizeof(short)*2*pm->nrows);
256	>	if (rord == NULL)
257	>	goto memerr;
258	>	for (ri = pm->nrows; ri--; ) {
259	>	rord[ri][0] = ri;
260	>	rord[ri][1] = -1;
261	>	if (src_rem[ri] <= minamt)          /* enough source material? */
262	>	continue;
263	>	for (i = 0; i < pm->ncols; i++)
264	>	if (dst_rem[ rord[ri][1] = psortrow(pm,ri)[i] ] > minamt)
265	>	break;
266	>	if (i >= pm->ncols) {               /* moved all we can? */
267	>	free(rord);
268	>	return(.0);
269	>	}
270		}
271	<	for (i = mtx_nrows(mig); i--; )         /* starting costs for diff. */
272	<	src_cost[i] = min_cost(src_rem[i], dst_rem,
273	<	pmtx+i*mtx_ncols(mig), mtx_ncols(mig));
274	<
271	>	if (pm->nrows > max2check)              /* sort if too many sources */
272	>	qsort_r(rord, pm->nrows, sizeof(short)*2, pm, &rmovcmp);
273	>	/* allocate cost array */
274	>	src_cost = (double *)malloc(sizeof(double)*pm->nrows);
275	>	if (src_cost == NULL)
276	>	goto memerr;
277	>	for (i = pm->nrows; i--; )              /* starting costs for diff. */
278	>	src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
279		/* find best source & dest. */
280		best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
281	<	for (cur.s = mtx_nrows(mig); cur.s--; ) {
282	<	const float *price = pmtx + cur.s*mtx_ncols(mig);
281	>	if ((r2check = pm->nrows) > max2check)
282	>	r2check = max2check;            /* put a limit on search */
283	>	for (ri = 0; ri < r2check; ri++) {      /* check each source row */
284		double      cost_others = 0;
285	<	if (src_rem[cur.s] <= minamt)
286	<	continue;
287	<	cur.d = -1;                         /* examine cheapest dest. */
288	<	for (i = mtx_ncols(mig); i--; )
289	<	if (dst_rem[i] > minamt &&
290	<	(cur.d < 0 || price[i] < price[cur.d]))
246	<	cur.d = i;
247	<	if (cur.d < 0)
248	<	return(.0);
249	<	if ((cur.price = price[cur.d]) >= best.price)
250	<	continue;                   /* no point checking further */
285	>	cur.s = rord[ri][0];
286	>	if ((cur.d = rord[ri][1]) < 0 ||
287	>	(cur.price = pricerow(pm,cur.s)[cur.d]) >= best.price) {
288	>	if (pm->nrows > max2check) break;       /* sorted end */
289	>	continue;                       /* else skip this one */
290	>	}
291		cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
292		src_rem[cur.s] : dst_rem[cur.d];
293	<	if (cur.amt > maxamt) cur.amt = maxamt;
294	<	dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
295	<	for (i = mtx_nrows(mig); i--; )
293	>	/* don't just leave smidgen */
294	>	if (cur.amt > maxamt*1.02) cur.amt = maxamt;
295	>	dst_rem[cur.d] -= cur.amt;          /* add up opportunity costs */
296	>	for (i = pm->nrows; i--; )
297		if (i != cur.s)
298	<	cost_others += min_cost(src_rem[i], dst_rem,
258	<	price, mtx_ncols(mig))
298	>	cost_others += min_cost(src_rem[i], dst_rem, pm, i)
299		- src_cost[i];
300		dst_rem[cur.d] += cur.amt;          /* undo trial move */
301		cur.price += cost_others/cur.amt;   /* adjust effective price */
302		if (cur.price < best.price)         /* are we better than best? */
303	<	best = cur;
303	>	best = cur;
304		}
305	<	if ((best.s < 0) | (best.d < 0))
305	>	free(src_cost);                         /* clean up */
306	>	free(rord);
307	>	if ((best.s < 0) | (best.d < 0))        /* nothing left to move? */
308		return(.0);
309	<	/* make the actual move */
309	>	/* else make the actual move */
310		mtx_coef(mig,best.s,best.d) += best.amt;
311		src_rem[best.s] -= best.amt;
312		dst_rem[best.d] -= best.amt;
313		return(best.amt);
314	+	memerr:
315	+	fprintf(stderr, "%s: Out of memory in migration_step()\n", progname);
316	+	exit(1);
317		}
318
274	–	#ifdef DEBUG
275	–	static char *
276	–	thetaphi(const FVECT v)
277	–	{
278	–	static char     buf[128];
279	–	double          theta, phi;
280	–
281	–	theta = 180./M_PI*acos(v[2]);
282	–	phi = 180./M_PI*atan2(v[1],v[0]);
283	–	sprintf(buf, "(%.0f,%.0f)", theta, phi);
284	–
285	–	return(buf);
286	–	}
287	–	#endif
288	–
319		/* Compute and insert migration along directed edge (may fork child) */
320		static MIGRATION *
321		create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
322		{
323		const double    end_thresh = 5e-6;
324	<	float           *pmtx;
324	>	PRICEMAT        pmtx;
325		MIGRATION       *newmig;
326		double          *src_rem, *dst_rem;
327		double          total_rem = 1., move_amt;
328	<	int             i;
328	>	int             i, j;
329		/* check if exists already */
330		for (newmig = from_rbf->ejl; newmig != NULL;
331		newmig = nextedge(from_rbf,newmig))
332		if (newmig->rbfv[1] == to_rbf)
333		return(NULL);
334		/* else allocate */
335	+	#ifdef DEBUG
336	+	fprintf(stderr, "Building path from (theta,phi) (%.0f,%.0f) ",
337	+	get_theta180(from_rbf->invec),
338	+	get_phi360(from_rbf->invec));
339	+	fprintf(stderr, "to (%.0f,%.0f) with %d x %d matrix\n",
340	+	get_theta180(to_rbf->invec),
341	+	get_phi360(to_rbf->invec),
342	+	from_rbf->nrbf, to_rbf->nrbf);
343	+	#endif
344		newmig = new_migration(from_rbf, to_rbf);
345		if (run_subprocess())
346		return(newmig);                 /* child continues */
347	<	pmtx = price_routes(from_rbf, to_rbf);
347	>	price_routes(&pmtx, from_rbf, to_rbf);
348		src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
349		dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
350		if ((src_rem == NULL) | (dst_rem == NULL)) {
#	Line 313 | Line 352 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
352		progname);
353		exit(1);
354		}
316	–	#ifdef DEBUG
317	–	fprintf(stderr, "Building path from (theta,phi) %s ",
318	–	thetaphi(from_rbf->invec));
319	–	fprintf(stderr, "to %s with %d x %d matrix\n",
320	–	thetaphi(to_rbf->invec),
321	–	from_rbf->nrbf, to_rbf->nrbf);
322	–	#endif
355		/* starting quantities */
356		memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
357		for (i = from_rbf->nrbf; i--; )
358		src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
359	<	for (i = to_rbf->nrbf; i--; )
360	<	dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
359	>	for (j = to_rbf->nrbf; j--; )
360	>	dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
361	>
362		do {                                    /* move a bit at a time */
363	<	move_amt = migration_step(newmig, src_rem, dst_rem, pmtx);
363	>	move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
364		total_rem -= move_amt;
332	–	#ifdef DEBUG
333	–	if (!nchild)
334	–	fprintf(stderr, "\r%.9f remaining...", total_rem);
335	–	#endif
365		} while ((total_rem > end_thresh) & (move_amt > 0));
366	<	#ifdef DEBUG
338	<	if (!nchild) fputs("done.\n", stderr);
339	<	else fprintf(stderr, "finished with %.9f remaining\n", total_rem);
340	<	#endif
366	>
367		for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
368	<	float       nf = rbf_volume(&from_rbf->rbfa[i]);
343	<	int         j;
368	>	double      nf = rbf_volume(&from_rbf->rbfa[i]);
369		if (nf <= FTINY) continue;
370		nf = from_rbf->vtotal / nf;
371		for (j = to_rbf->nrbf; j--; )
372	<	mtx_coef(newmig,i,j) *= nf;
372	>	mtx_coef(newmig,i,j) *= nf;     /* row now sums to 1.0 */
373		}
374		end_subprocess();                       /* exit here if subprocess */
375	<	free(pmtx);                             /* free working arrays */
375	>	free_routes(&pmtx);                     /* free working arrays */
376		free(src_rem);
377		free(dst_rem);
378		return(newmig);

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines (old)
>	Changed lines (new)