ViewVC Help
View File | Revision Log | Show Annotations | Download File | Root Listing
root/radiance/ray/src/cv/bsdfmesh.c
(Generate patch)

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
> Changed lines