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.26 by greg, Wed Mar 26 00:11:30 2014 UTC vs.
Revision 2.27 by greg, Wed Mar 26 02:52:31 2014 UTC

#	Line 27 | Line 27 | int                    nprocs = 1;
27		/* number of children (-1 in child) */
28		static int              nchild = 0;
29
30	–	typedef struct {
31	–	int             nrows, ncols;   /* array size (matches migration) */
32	–	float           *price;         /* migration prices */
33	–	short           *sord;          /* sort for each row, low to high */
34	–	float           *prow;          /* current price row */
35	–	} PRICEMAT;                     /* sorted pricing matrix */
36	–
37	–	#define pricerow(p,i)   ((p)->price + (i)*(p)->ncols)
38	–	#define psortrow(p,i)   ((p)->sord + (i)*(p)->ncols)
39	–
30		/* Create a new migration holder (sharing memory for multiprocessing) */
31		static MIGRATION *
32		new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
#	Line 180 | Line 170 | neighborhood_dist2(int x0, int y0, int x1, int y1)
170		return(sum2 / (4*rad*(rad+1) + 1));
171		}
172
173	<	/* Comparison routine needed for sorting price row */
174	<	static int
175	<	msrt_cmp(void *b, const void *p1, const void *p2)
173	>	/* Compute distance between two RBF lobes */
174	>	double
175	>	lobe_distance(RBFVAL *rbf1, RBFVAL *rbf2)
176		{
177	<	PRICEMAT        *pm = (PRICEMAT *)b;
178	<	float           c1 = pm->prow[*(const short *)p1];
179	<	float           c2 = pm->prow[*(const short *)p2];
180	<
181	<	if (c1 > c2) return(1);
182	<	if (c1 < c2) return(-1);
183	<	return(0);
177	>	FVECT   vfrom, vto;
178	>	double  d, res;
179	>	/* quadratic cost function */
180	>	ovec_from_pos(vfrom, rbf1->gx, rbf1->gy);
181	>	ovec_from_pos(vto, rbf2->gx, rbf2->gy);
182	>	d = Acos(DOT(vfrom, vto));
183	>	res = d*d;
184	>	d = R2ANG(rbf2->crad) - R2ANG(rbf1->crad);
185	>	res += d*d;
186	>	/* neighborhood difference */
187	>	res += NEIGH_FACT2 * neighborhood_dist2( rbf1->gx, rbf1->gy,
188	>	rbf2->gx, rbf2->gy );
189	>	return(res);
190		}
191
196	–	/* Compute (and allocate) migration price matrix for optimization */
197	–	static void
198	–	price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
199	–	{
200	–	FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
201	–	int     i, j;
192
203	–	compute_nDSFs(from_rbf, to_rbf);
204	–	pm->nrows = from_rbf->nrbf;
205	–	pm->ncols = to_rbf->nrbf;
206	–	pm->price = (float *)malloc(sizeof(float) * pm->nrows*pm->ncols);
207	–	pm->sord = (short *)malloc(sizeof(short) * pm->nrows*pm->ncols);
208	–
209	–	if ((pm->price == NULL) | (pm->sord == NULL) | (vto == NULL)) {
210	–	fprintf(stderr, "%s: Out of memory in migration_costs()\n",
211	–	progname);
212	–	exit(1);
213	–	}
214	–	for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
215	–	ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
216	–
217	–	for (i = from_rbf->nrbf; i--; ) {
218	–	const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
219	–	FVECT               vfrom;
220	–	short               *srow;
221	–	ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
222	–	pm->prow = pricerow(pm,i);
223	–	srow = psortrow(pm,i);
224	–	for (j = to_rbf->nrbf; j--; ) {
225	–	double  d;                      /* quadratic cost function */
226	–	d = Acos(DOT(vfrom, vto[j]));
227	–	pm->prow[j] = d*d;
228	–	d = R2ANG(to_rbf->rbfa[j].crad) - from_ang;
229	–	pm->prow[j] += d*d;
230	–	/* neighborhood difference */
231	–	pm->prow[j] += NEIGH_FACT2 * neighborhood_dist2(
232	–	from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy,
233	–	to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy );
234	–	srow[j] = j;
235	–	}
236	–	qsort_r(srow, pm->ncols, sizeof(short), pm, &msrt_cmp);
237	–	}
238	–	free(vto);
239	–	}
240	–
241	–	/* Free price matrix */
242	–	static void
243	–	free_routes(PRICEMAT *pm)
244	–	{
245	–	free(pm->price); pm->price = NULL;
246	–	free(pm->sord); pm->sord = NULL;
247	–	}
248	–
249	–	/* Compute minimum (optimistic) cost for moving the given source material */
250	–	static double
251	–	min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
252	–	{
253	–	const short     *srow = psortrow(pm,s);
254	–	const float     *prow = pricerow(pm,s);
255	–	double          total_cost = 0;
256	–	int             j;
257	–	/* move cheapest first */
258	–	for (j = 0; (j < pm->ncols) & (amt2move > FTINY); j++) {
259	–	int     d = srow[j];
260	–	double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
261	–
262	–	total_cost += amt * prow[d];
263	–	amt2move -= amt;
264	–	}
265	–	return(total_cost);
266	–	}
267	–
268	–	typedef struct {
269	–	short   s, d;           /* source and destination */
270	–	float   dc;             /* discount to push inventory */
271	–	} ROWSENT;              /* row sort entry */
272	–
273	–	/* Compare entries by discounted moving price */
274	–	static int
275	–	rmovcmp(void *b, const void *p1, const void *p2)
276	–	{
277	–	PRICEMAT        *pm = (PRICEMAT *)b;
278	–	const ROWSENT   *re1 = (const ROWSENT *)p1;
279	–	const ROWSENT   *re2 = (const ROWSENT *)p2;
280	–	double          price_diff;
281	–
282	–	if (re1->d < 0) return(re2->d >= 0);
283	–	if (re2->d < 0) return(-1);
284	–	price_diff = re1->dc*pricerow(pm,re1->s)[re1->d] -
285	–	re2->dc*pricerow(pm,re2->s)[re2->d];
286	–	if (price_diff > 0) return(1);
287	–	if (price_diff < 0) return(-1);
288	–	return(0);
289	–	}
290	–
291	–	/* Take a step in migration by choosing reasonable bucket to transfer */
292	–	static double
293	–	migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, PRICEMAT *pm)
294	–	{
295	–	const int       max2check = 100;
296	–	const double    maxamt = 1./(double)pm->ncols;
297	–	const double    minamt = maxamt*1e-4;
298	–	double          *src_cost;
299	–	ROWSENT         *rord;
300	–	struct {
301	–	int     s, d;   /* source and destination */
302	–	double  price;  /* cost per amount moved */
303	–	double  amt;    /* amount we can move */
304	–	} cur, best;
305	–	int             r2check, i, ri;
306	–	/*
307	–	* Check cheapest available routes only -- a higher adjusted
308	–	* destination price implies that another source is closer, so
309	–	* we can hold off considering more expensive options until
310	–	* some other (hopefully better) moves have been made.
311	–	* A discount based on source remaining is supposed to prioritize
312	–	* movement from large lobes, but it doesn't seem to do much,
313	–	* so we have it set to 1.0 at the moment.
314	–	*/
315	–	#define discount(qr)    1.0
316	–	/* most promising row order */
317	–	rord = (ROWSENT *)malloc(sizeof(ROWSENT)*pm->nrows);
318	–	if (rord == NULL)
319	–	goto memerr;
320	–	for (ri = pm->nrows; ri--; ) {
321	–	rord[ri].s = ri;
322	–	rord[ri].d = -1;
323	–	rord[ri].dc = 1.f;
324	–	if (src_rem[ri] <= minamt)          /* enough source material? */
325	–	continue;
326	–	for (i = 0; i < pm->ncols; i++)
327	–	if (dst_rem[ rord[ri].d = psortrow(pm,ri)[i] ] > minamt)
328	–	break;
329	–	if (i >= pm->ncols) {               /* moved all we can? */
330	–	free(rord);
331	–	return(.0);
332	–	}
333	–	rord[ri].dc = discount(src_rem[ri]);
334	–	}
335	–	if (pm->nrows > max2check)              /* sort if too many sources */
336	–	qsort_r(rord, pm->nrows, sizeof(ROWSENT), pm, &rmovcmp);
337	–	/* allocate cost array */
338	–	src_cost = (double *)malloc(sizeof(double)*pm->nrows);
339	–	if (src_cost == NULL)
340	–	goto memerr;
341	–	for (i = pm->nrows; i--; )              /* starting costs for diff. */
342	–	src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
343	–	/* find best source & dest. */
344	–	best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
345	–	if ((r2check = pm->nrows) > max2check)
346	–	r2check = max2check;            /* put a limit on search */
347	–	for (ri = 0; ri < r2check; ri++) {      /* check each source row */
348	–	double      cost_others = 0;
349	–	cur.s = rord[ri].s;
350	–	if ((cur.d = rord[ri].d) < 0 ||
351	–	rord[ri].dc*pricerow(pm,cur.s)[cur.d] >= best.price) {
352	–	if (pm->nrows > max2check) break;       /* sorted end */
353	–	continue;                       /* else skip this one */
354	–	}
355	–	cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
356	–	src_rem[cur.s] : dst_rem[cur.d];
357	–	/* don't just leave smidgen */
358	–	if (cur.amt > maxamt*1.02) cur.amt = maxamt;
359	–	dst_rem[cur.d] -= cur.amt;          /* add up opportunity costs */
360	–	for (i = pm->nrows; i--; )
361	–	if (i != cur.s)
362	–	cost_others += min_cost(src_rem[i], dst_rem, pm, i)
363	–	- src_cost[i];
364	–	dst_rem[cur.d] += cur.amt;          /* undo trial move */
365	–	/* discount effective price */
366	–	cur.price = ( pricerow(pm,cur.s)[cur.d] + cost_others/cur.amt ) *
367	–	rord[ri].dc;
368	–	if (cur.price < best.price)         /* are we better than best? */
369	–	best = cur;
370	–	}
371	–	free(src_cost);                         /* clean up */
372	–	free(rord);
373	–	if ((best.s < 0) | (best.d < 0))        /* nothing left to move? */
374	–	return(.0);
375	–	/* else make the actual move */
376	–	mtx_coef(mig,best.s,best.d) += best.amt;
377	–	src_rem[best.s] -= best.amt;
378	–	dst_rem[best.d] -= best.amt;
379	–	return(best.amt);
380	–	memerr:
381	–	fprintf(stderr, "%s: Out of memory in migration_step()\n", progname);
382	–	exit(1);
383	–	#undef discount
384	–	}
385	–
386	–	#ifdef DUMP_MATRIX
387	–	/* Dump transport plan and corresponding price matrix to a text file */
388	–	static void
389	–	dump_matrix(const MIGRATION *me, const PRICEMAT *pm)
390	–	{
391	–	char    fname[256];
392	–	FILE    *fp;
393	–	int     i, j;
394	–
395	–	sprintf(fname, "edge_%d-%d.txt", me->rbfv[0]->ord, me->rbfv[1]->ord);
396	–	if ((fp = fopen(fname, "w")) == NULL)
397	–	return;
398	–	for (j = 0; j < 2; j++) {
399	–	fprintf(fp, "Available from %d source RBF lobes in node %d:\n",
400	–	me->rbfv[j]->nrbf, me->rbfv[j]->ord);
401	–	for (i = 0; i < me->rbfv[j]->nrbf; i++)
402	–	fprintf(fp, " %.4e", rbf_volume(&me->rbfv[j]->rbfa[i]) /
403	–	me->rbfv[j]->vtotal);
404	–	fputc('\n', fp);
405	–	}
406	–	fprintf(fp, "Price (quadratic distance metric) matrix:\n");
407	–	for (i = 0; i < pm->nrows; i++) {
408	–	for (j = 0; j < pm->ncols; j++)
409	–	fprintf(fp, " %.4e", pricerow(pm,i)[j]);
410	–	fputc('\n', fp);
411	–	}
412	–	fprintf(fp, "Solution matrix (transport plan):\n");
413	–	for (i = 0; i < mtx_nrows(me); i++) {
414	–	for (j = 0; j < mtx_ncols(me); j++)
415	–	fprintf(fp, " %.4e", mtx_coef(me,i,j));
416	–	fputc('\n', fp);
417	–	}
418	–	fclose(fp);
419	–	}
420	–	#endif
421	–
193		/* Compute and insert migration along directed edge (may fork child) */
194		static MIGRATION *
195		create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
196		{
426	–	const double    end_thresh = 5e-6;
427	–	PRICEMAT        pmtx;
197		MIGRATION       *newmig;
429	–	double          *src_rem, *dst_rem;
430	–	double          total_rem = 1., move_amt;
198		int             i, j;
199		/* check if exists already */
200		for (newmig = from_rbf->ejl; newmig != NULL;
#	Line 447 | Line 214 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
214		newmig = new_migration(from_rbf, to_rbf);
215		if (run_subprocess())
216		return(newmig);                 /* child continues */
450	–	price_routes(&pmtx, from_rbf, to_rbf);
451	–	src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
452	–	dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
453	–	if ((src_rem == NULL) | (dst_rem == NULL)) {
454	–	fprintf(stderr, "%s: Out of memory in create_migration()\n",
455	–	progname);
456	–	exit(1);
457	–	}
458	–	/* starting quantities */
459	–	memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
460	–	for (i = from_rbf->nrbf; i--; )
461	–	src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
462	–	for (j = to_rbf->nrbf; j--; )
463	–	dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
217
218	<	do {                                    /* move a bit at a time */
219	<	move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
220	<	total_rem -= move_amt;
468	<	} while ((total_rem > end_thresh) & (move_amt > 0));
218	>	/* compute transport plan */
219	>	compute_nDSFs(from_rbf, to_rbf);
220	>	plan_transport(newmig);
221
222		for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
223		double      nf = rbf_volume(&from_rbf->rbfa[i]);
#	Line 474 | Line 226 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
226		for (j = to_rbf->nrbf; j--; )
227		mtx_coef(newmig,i,j) *= nf;     /* row now sums to 1.0 */
228		}
477	–	#ifdef DUMP_MATRIX
478	–	dump_matrix(newmig, &pmtx);
479	–	#endif
229		end_subprocess();                       /* exit here if subprocess */
481	–	free_routes(&pmtx);                     /* free working arrays */
482	–	free(src_rem);
483	–	free(dst_rem);
230		return(newmig);
231		}
232

Diff Legend

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