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.1 by greg, Fri Oct 19 04:14:29 2012 UTC vs.
Revision 2.4 by greg, Thu Nov 8 23:11:41 2012 UTC

# Line 23 | Line 23 | int                    nprocs = 1;
23                                  /* number of children (-1 in child) */
24   static int              nchild = 0;
25  
26 < /* Compute (and allocate) migration price matrix for optimization */
27 < static float *
28 < price_routes(const RBFNODE *from_rbf, const RBFNODE *to_rbf)
29 < {
30 <        float   *pmtx = (float *)malloc(sizeof(float) *
31 <                                        from_rbf->nrbf * to_rbf->nrbf);
32 <        FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
33 <        int     i, j;
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 > } PRICEMAT;                     /* sorted pricing matrix */
31  
32 <        if ((pmtx == NULL) | (vto == NULL)) {
33 <                fprintf(stderr, "%s: Out of memory in migration_costs()\n",
37 <                                progname);
38 <                exit(1);
39 <        }
40 <        for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
41 <                ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
32 > #define pricerow(p,i)   ((p)->price + (i)*(p)->ncols)
33 > #define psortrow(p,i)   ((p)->sord + (i)*(p)->ncols)
34  
43 –        for (i = from_rbf->nrbf; i--; ) {
44 –            const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
45 –            FVECT               vfrom;
46 –            ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
47 –            for (j = to_rbf->nrbf; j--; )
48 –                pmtx[i*to_rbf->nrbf + j] = acos(DOT(vfrom, vto[j])) +
49 –                                fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
50 –        }
51 –        free(vto);
52 –        return(pmtx);
53 – }
54 –
55 – /* Comparison routine needed for sorting price row */
56 – static const float      *price_arr;
57 – static int
58 – msrt_cmp(const void *p1, const void *p2)
59 – {
60 –        float   c1 = price_arr[*(const int *)p1];
61 –        float   c2 = price_arr[*(const int *)p2];
62 –
63 –        if (c1 > c2) return(1);
64 –        if (c1 < c2) return(-1);
65 –        return(0);
66 – }
67 –
68 – /* Compute minimum (optimistic) cost for moving the given source material */
69 – static double
70 – min_cost(double amt2move, const double *avail, const float *price, int n)
71 – {
72 –        static int      *price_sort = NULL;
73 –        static int      n_alloc = 0;
74 –        double          total_cost = 0;
75 –        int             i;
76 –
77 –        if (amt2move <= FTINY)                  /* pre-emptive check */
78 –                return(0.);
79 –        if (n > n_alloc) {                      /* (re)allocate sort array */
80 –                if (n_alloc) free(price_sort);
81 –                price_sort = (int *)malloc(sizeof(int)*n);
82 –                if (price_sort == NULL) {
83 –                        fprintf(stderr, "%s: Out of memory in min_cost()\n",
84 –                                        progname);
85 –                        exit(1);
86 –                }
87 –                n_alloc = n;
88 –        }
89 –        for (i = n; i--; )
90 –                price_sort[i] = i;
91 –        price_arr = price;
92 –        qsort(price_sort, n, sizeof(int), &msrt_cmp);
93 –                                                /* move cheapest first */
94 –        for (i = 0; i < n && amt2move > FTINY; i++) {
95 –                int     d = price_sort[i];
96 –                double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
97 –
98 –                total_cost += amt * price[d];
99 –                amt2move -= amt;
100 –        }
101 –        return(total_cost);
102 – }
103 –
104 – /* Take a step in migration by choosing optimal bucket to transfer */
105 – static double
106 – migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const float *pmtx)
107 – {
108 –        const double    maxamt = .1;
109 –        const double    minamt = maxamt*.0001;
110 –        static double   *src_cost = NULL;
111 –        static int      n_alloc = 0;
112 –        struct {
113 –                int     s, d;   /* source and destination */
114 –                double  price;  /* price estimate per amount moved */
115 –                double  amt;    /* amount we can move */
116 –        } cur, best;
117 –        int             i;
118 –
119 –        if (mtx_nrows(mig) > n_alloc) {         /* allocate cost array */
120 –                if (n_alloc)
121 –                        free(src_cost);
122 –                src_cost = (double *)malloc(sizeof(double)*mtx_nrows(mig));
123 –                if (src_cost == NULL) {
124 –                        fprintf(stderr, "%s: Out of memory in migration_step()\n",
125 –                                        progname);
126 –                        exit(1);
127 –                }
128 –                n_alloc = mtx_nrows(mig);
129 –        }
130 –        for (i = mtx_nrows(mig); i--; )         /* starting costs for diff. */
131 –                src_cost[i] = min_cost(src_rem[i], dst_rem,
132 –                                        pmtx+i*mtx_ncols(mig), mtx_ncols(mig));
133 –
134 –                                                /* find best source & dest. */
135 –        best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
136 –        for (cur.s = mtx_nrows(mig); cur.s--; ) {
137 –            const float *price = pmtx + cur.s*mtx_ncols(mig);
138 –            double      cost_others = 0;
139 –            if (src_rem[cur.s] < minamt)
140 –                    continue;
141 –            cur.d = -1;                         /* examine cheapest dest. */
142 –            for (i = mtx_ncols(mig); i--; )
143 –                if (dst_rem[i] > minamt &&
144 –                                (cur.d < 0 || price[i] < price[cur.d]))
145 –                        cur.d = i;
146 –            if (cur.d < 0)
147 –                    return(.0);
148 –            if ((cur.price = price[cur.d]) >= best.price)
149 –                    continue;                   /* no point checking further */
150 –            cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
151 –                                src_rem[cur.s] : dst_rem[cur.d];
152 –            if (cur.amt > maxamt) cur.amt = maxamt;
153 –            dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
154 –            for (i = mtx_nrows(mig); i--; )
155 –                if (i != cur.s)
156 –                        cost_others += min_cost(src_rem[i], dst_rem,
157 –                                                price, mtx_ncols(mig))
158 –                                        - src_cost[i];
159 –            dst_rem[cur.d] += cur.amt;          /* undo trial move */
160 –            cur.price += cost_others/cur.amt;   /* adjust effective price */
161 –            if (cur.price < best.price)         /* are we better than best? */
162 –                    best = cur;
163 –        }
164 –        if ((best.s < 0) | (best.d < 0))
165 –                return(.0);
166 –                                                /* make the actual move */
167 –        mig->mtx[mtx_ndx(mig,best.s,best.d)] += best.amt;
168 –        src_rem[best.s] -= best.amt;
169 –        dst_rem[best.d] -= best.amt;
170 –        return(best.amt);
171 – }
172 –
173 – #ifdef DEBUG
174 – static char *
175 – thetaphi(const FVECT v)
176 – {
177 –        static char     buf[128];
178 –        double          theta, phi;
179 –
180 –        theta = 180./M_PI*acos(v[2]);
181 –        phi = 180./M_PI*atan2(v[1],v[0]);
182 –        sprintf(buf, "(%.0f,%.0f)", theta, phi);
183 –
184 –        return(buf);
185 – }
186 – #endif
187 –
35   /* Create a new migration holder (sharing memory for multiprocessing) */
36   static MIGRATION *
37   new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
# Line 286 | Line 133 | run_subprocess(void)
133  
134   #endif  /* ! _WIN32 */
135  
136 + /* Comparison routine needed for sorting price row */
137 + static int
138 + msrt_cmp(void *b, const void *p1, const void *p2)
139 + {
140 +        PRICEMAT        *pm = (PRICEMAT *)b;
141 +        int             ri = ((const short *)p1 - pm->sord) / pm->ncols;
142 +        float           c1 = pricerow(pm,ri)[*(const short *)p1];
143 +        float           c2 = pricerow(pm,ri)[*(const short *)p2];
144 +
145 +        if (c1 > c2) return(1);
146 +        if (c1 < c2) return(-1);
147 +        return(0);
148 + }
149 +
150 + /* Compute (and allocate) migration price matrix for optimization */
151 + static void
152 + price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
153 + {
154 +        FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
155 +        int     i, j;
156 +
157 +        pm->nrows = from_rbf->nrbf;
158 +        pm->ncols = to_rbf->nrbf;
159 +        pm->price = (float *)malloc(sizeof(float) * pm->nrows*pm->ncols);
160 +        pm->sord = (short *)malloc(sizeof(short) * pm->nrows*pm->ncols);
161 +        
162 +        if ((pm->price == NULL) | (pm->sord == NULL) | (vto == NULL)) {
163 +                fprintf(stderr, "%s: Out of memory in migration_costs()\n",
164 +                                progname);
165 +                exit(1);
166 +        }
167 +        for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
168 +                ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
169 +
170 +        for (i = from_rbf->nrbf; i--; ) {
171 +            const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
172 +            FVECT               vfrom;
173 +            ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
174 +            for (j = to_rbf->nrbf; j--; ) {
175 +                pricerow(pm,i)[j] = acos(DOT(vfrom, vto[j])) +
176 +                                fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
177 +                psortrow(pm,i)[j] = j;
178 +            }
179 +            qsort_r(psortrow(pm,i), pm->ncols, sizeof(short), pm, &msrt_cmp);
180 +        }
181 +        free(vto);
182 + }
183 +
184 + /* Free price matrix */
185 + static void
186 + free_routes(PRICEMAT *pm)
187 + {
188 +        free(pm->price); pm->price = NULL;
189 +        free(pm->sord); pm->sord = NULL;
190 + }
191 +
192 + /* Compute minimum (optimistic) cost for moving the given source material */
193 + static double
194 + min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
195 + {
196 +        double          total_cost = 0;
197 +        int             j;
198 +
199 +        if (amt2move <= FTINY)                  /* pre-emptive check */
200 +                return(0.);
201 +                                                /* move cheapest first */
202 +        for (j = 0; j < pm->ncols && amt2move > FTINY; j++) {
203 +                int     d = psortrow(pm,s)[j];
204 +                double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
205 +
206 +                total_cost += amt * pricerow(pm,s)[d];
207 +                amt2move -= amt;
208 +        }
209 +        return(total_cost);
210 + }
211 +
212 + /* Take a step in migration by choosing optimal bucket to transfer */
213 + static double
214 + migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const PRICEMAT *pm)
215 + {
216 +        const double    maxamt = 1./(double)pm->ncols;
217 +        const double    minamt = maxamt*5e-6;
218 +        static double   *src_cost = NULL;
219 +        static int      n_alloc = 0;
220 +        struct {
221 +                int     s, d;   /* source and destination */
222 +                double  price;  /* price estimate per amount moved */
223 +                double  amt;    /* amount we can move */
224 +        } cur, best;
225 +        int             i;
226 +
227 +        if (pm->nrows > n_alloc) {              /* allocate cost array */
228 +                if (n_alloc)
229 +                        free(src_cost);
230 +                src_cost = (double *)malloc(sizeof(double)*pm->nrows);
231 +                if (src_cost == NULL) {
232 +                        fprintf(stderr, "%s: Out of memory in migration_step()\n",
233 +                                        progname);
234 +                        exit(1);
235 +                }
236 +                n_alloc = pm->nrows;
237 +        }
238 +        for (i = pm->nrows; i--; )              /* starting costs for diff. */
239 +                src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
240 +
241 +                                                /* find best source & dest. */
242 +        best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
243 +        for (cur.s = pm->nrows; cur.s--; ) {
244 +            double      cost_others = 0;
245 +            if (src_rem[cur.s] <= minamt)
246 +                    continue;
247 +                                                /* examine cheapest dest. */
248 +            for (i = 0; i < pm->ncols; i++)
249 +                if (dst_rem[cur.d = psortrow(pm,cur.s)[i]] > minamt)
250 +                        break;
251 +            if (i >= pm->ncols)
252 +                    return(.0);
253 +            if ((cur.price = pricerow(pm,cur.s)[cur.d]) >= best.price)
254 +                    continue;                   /* no point checking further */
255 +            cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
256 +                                src_rem[cur.s] : dst_rem[cur.d];
257 +            if (cur.amt > maxamt) cur.amt = maxamt;
258 +            dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
259 +            for (i = pm->nrows; i--; )
260 +                if (i != cur.s)
261 +                        cost_others += min_cost(src_rem[i], dst_rem, pm, i)
262 +                                        - src_cost[i];
263 +            dst_rem[cur.d] += cur.amt;          /* undo trial move */
264 +            cur.price += cost_others/cur.amt;   /* adjust effective price */
265 +            if (cur.price < best.price)         /* are we better than best? */
266 +                    best = cur;
267 +        }
268 +        if ((best.s < 0) | (best.d < 0))
269 +                return(.0);
270 +                                                /* make the actual move */
271 +        mtx_coef(mig,best.s,best.d) += best.amt;
272 +        src_rem[best.s] -= best.amt;
273 +        dst_rem[best.d] -= best.amt;
274 +        return(best.amt);
275 + }
276 +
277 + #ifdef DEBUG
278 + static char *
279 + thetaphi(const FVECT v)
280 + {
281 +        static char     buf[128];
282 +        double          theta, phi;
283 +
284 +        theta = 180./M_PI*acos(v[2]);
285 +        phi = 180./M_PI*atan2(v[1],v[0]);
286 +        sprintf(buf, "(%.0f,%.0f)", theta, phi);
287 +
288 +        return(buf);
289 + }
290 + #endif
291 +
292   /* Compute and insert migration along directed edge (may fork child) */
293   static MIGRATION *
294   create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
295   {
296 <        const double    end_thresh = 0.1/(from_rbf->nrbf*to_rbf->nrbf);
297 <        const double    check_thresh = 0.01;
295 <        const double    rel_thresh = 5e-6;
296 <        float           *pmtx;
296 >        const double    end_thresh = 5e-6;
297 >        PRICEMAT        pmtx;
298          MIGRATION       *newmig;
299          double          *src_rem, *dst_rem;
300          double          total_rem = 1., move_amt;
# Line 307 | Line 308 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
308          newmig = new_migration(from_rbf, to_rbf);
309          if (run_subprocess())
310                  return(newmig);                 /* child continues */
311 <        pmtx = price_routes(from_rbf, to_rbf);
311 >        price_routes(&pmtx, from_rbf, to_rbf);
312          src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
313          dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
314          if ((src_rem == NULL) | (dst_rem == NULL)) {
# Line 318 | Line 319 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
319   #ifdef DEBUG
320          fprintf(stderr, "Building path from (theta,phi) %s ",
321                          thetaphi(from_rbf->invec));
322 <        fprintf(stderr, "to %s", thetaphi(to_rbf->invec));
323 <        /* if (nchild) */ fputc('\n', stderr);
322 >        fprintf(stderr, "to %s with %d x %d matrix\n",
323 >                        thetaphi(to_rbf->invec),
324 >                        from_rbf->nrbf, to_rbf->nrbf);
325   #endif
326                                                  /* starting quantities */
327          memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
# Line 328 | Line 330 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
330          for (i = to_rbf->nrbf; i--; )
331                  dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
332          do {                                    /* move a bit at a time */
333 <                move_amt = migration_step(newmig, src_rem, dst_rem, pmtx);
333 >                move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
334                  total_rem -= move_amt;
335   #ifdef DEBUG
336                  if (!nchild)
337 <                        /* fputc('.', stderr); */
336 <                        fprintf(stderr, "%.9f remaining...\r", total_rem);
337 >                        fprintf(stderr, "\r%.9f remaining...", total_rem);
338   #endif
339 <        } while (total_rem > end_thresh && (total_rem > check_thresh) |
339 <                                        (move_amt > rel_thresh*total_rem));
339 >        } while ((total_rem > end_thresh) & (move_amt > 0));
340   #ifdef DEBUG
341 <        if (!nchild) fputs("\ndone.\n", stderr);
341 >        if (!nchild) fputs("done.\n", stderr);
342          else fprintf(stderr, "finished with %.9f remaining\n", total_rem);
343   #endif
344          for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
# Line 347 | Line 347 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
347              if (nf <= FTINY) continue;
348              nf = from_rbf->vtotal / nf;
349              for (j = to_rbf->nrbf; j--; )
350 <                newmig->mtx[mtx_ndx(newmig,i,j)] *= nf;
350 >                mtx_coef(newmig,i,j) *= nf;
351          }
352          end_subprocess();                       /* exit here if subprocess */
353 <        free(pmtx);                             /* free working arrays */
353 >        free_routes(&pmtx);                     /* free working arrays */
354          free(src_rem);
355          free(dst_rem);
356          return(newmig);

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines