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.7 by greg, Sat Nov 10 19:47:42 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 272 | Line 119 | run_subprocess(void)
119                  if (pid < 0) {
120                          fprintf(stderr, "%s: cannot fork subprocess\n",
121                                          progname);
122 +                        await_children(nchild);
123                          exit(1);
124                  }
125                  ++nchild;                       /* subprocess started */
# Line 286 | Line 134 | run_subprocess(void)
134  
135   #endif  /* ! _WIN32 */
136  
137 + /* Comparison routine needed for sorting price row */
138 + static int
139 + msrt_cmp(void *b, const void *p1, const void *p2)
140 + {
141 +        PRICEMAT        *pm = (PRICEMAT *)b;
142 +        int             ri = ((const short *)p1 - pm->sord) / pm->ncols;
143 +        float           c1 = pricerow(pm,ri)[*(const short *)p1];
144 +        float           c2 = pricerow(pm,ri)[*(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 void
153 + price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
154 + {
155 +        FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
156 +        int     i, j;
157 +
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);
167 +        }
168 +        for (j = to_rbf->nrbf; j--; )           /* save repetitive ops. */
169 +                ovec_from_pos(vto[j], to_rbf->rbfa[j].gx, to_rbf->rbfa[j].gy);
170 +
171 +        for (i = from_rbf->nrbf; i--; ) {
172 +            const double        from_ang = R2ANG(from_rbf->rbfa[i].crad);
173 +            FVECT               vfrom;
174 +            ovec_from_pos(vfrom, from_rbf->rbfa[i].gx, from_rbf->rbfa[i].gy);
175 +            for (j = to_rbf->nrbf; j--; ) {
176 +                double          dprod = DOT(vfrom, vto[j]);
177 +                pricerow(pm,i)[j] = ((dprod >= 1.) ? .0 : acos(dprod)) +
178 +                                fabs(R2ANG(to_rbf->rbfa[j].crad) - from_ang);
179 +                psortrow(pm,i)[j] = j;
180 +            }
181 +            qsort_r(psortrow(pm,i), pm->ncols, sizeof(short), pm, &msrt_cmp);
182 +        }
183 +        free(vto);
184 + }
185 +
186 + /* Free price matrix */
187 + static void
188 + free_routes(PRICEMAT *pm)
189 + {
190 +        free(pm->price); pm->price = NULL;
191 +        free(pm->sord); pm->sord = NULL;
192 + }
193 +
194 + /* Compute minimum (optimistic) cost for moving the given source material */
195 + static double
196 + min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
197 + {
198 +        double          total_cost = 0;
199 +        int             j;
200 +
201 +        if (amt2move <= FTINY)                  /* pre-emptive check */
202 +                return(.0);
203 +                                                /* move cheapest first */
204 +        for (j = 0; j < pm->ncols && amt2move > FTINY; j++) {
205 +                int     d = psortrow(pm,s)[j];
206 +                double  amt = (amt2move < avail[d]) ? amt2move : avail[d];
207 +
208 +                total_cost += amt * pricerow(pm,s)[d];
209 +                amt2move -= amt;
210 +        }
211 +        return(total_cost);
212 + }
213 +
214 + /* Take a step in migration by choosing optimal bucket to transfer */
215 + static double
216 + migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const PRICEMAT *pm)
217 + {
218 +        const double    maxamt = 1./(double)pm->ncols;
219 +        const double    minamt = maxamt*5e-6;
220 +        double          *src_cost;
221 +        struct {
222 +                int     s, d;   /* source and destination */
223 +                double  price;  /* price estimate per amount moved */
224 +                double  amt;    /* amount we can move */
225 +        } cur, best;
226 +        int             i;
227 +                                                /* allocate cost array */
228 +        src_cost = (double *)malloc(sizeof(double)*pm->nrows);
229 +        if (src_cost == NULL) {
230 +                fprintf(stderr, "%s: Out of memory in migration_step()\n",
231 +                                progname);
232 +                exit(1);
233 +        }
234 +        for (i = pm->nrows; i--; )              /* starting costs for diff. */
235 +                src_cost[i] = min_cost(src_rem[i], dst_rem, pm, i);
236 +
237 +                                                /* find best source & dest. */
238 +        best.s = best.d = -1; best.price = FHUGE; best.amt = 0;
239 +        for (cur.s = pm->nrows; cur.s--; ) {
240 +            double      cost_others = 0;
241 +
242 +            if (src_rem[cur.s] <= minamt)
243 +                    continue;
244 +                                                /* examine cheapest dest. */
245 +            for (i = 0; i < pm->ncols; i++)
246 +                if (dst_rem[ cur.d = psortrow(pm,cur.s)[i] ] > minamt)
247 +                        break;
248 +            if (i >= pm->ncols)
249 +                break;
250 +            if ((cur.price = pricerow(pm,cur.s)[cur.d]) >= best.price)
251 +                continue;                       /* no point checking further */
252 +            cur.amt = (src_rem[cur.s] < dst_rem[cur.d]) ?
253 +                                src_rem[cur.s] : dst_rem[cur.d];
254 +            if (cur.amt > maxamt) cur.amt = maxamt;
255 +            dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
256 +            for (i = pm->nrows; i--; )
257 +                if (i != cur.s)
258 +                        cost_others += min_cost(src_rem[i], dst_rem, pm, i)
259 +                                        - src_cost[i];
260 +            dst_rem[cur.d] += cur.amt;          /* undo trial move */
261 +            cur.price += cost_others/cur.amt;   /* adjust effective price */
262 +            if (cur.price < best.price)         /* are we better than best? */
263 +                    best = cur;
264 +        }
265 +        free(src_cost);                         /* finish up */
266 +
267 +        if ((best.s < 0) | (best.d < 0))        /* nothing left to move? */
268 +                return(.0);
269 +                                                /* else make the actual move */
270 +        mtx_coef(mig,best.s,best.d) += best.amt;
271 +        src_rem[best.s] -= best.amt;
272 +        dst_rem[best.d] -= best.amt;
273 +        return(best.amt);
274 + }
275 +
276   /* Compute and insert migration along directed edge (may fork child) */
277   static MIGRATION *
278   create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
279   {
280 <        const double    end_thresh = 0.1/(from_rbf->nrbf*to_rbf->nrbf);
281 <        const double    check_thresh = 0.01;
295 <        const double    rel_thresh = 5e-6;
296 <        float           *pmtx;
280 >        const double    end_thresh = 5e-6;
281 >        PRICEMAT        pmtx;
282          MIGRATION       *newmig;
283          double          *src_rem, *dst_rem;
284          double          total_rem = 1., move_amt;
285 <        int             i;
285 >        int             i, j;
286                                                  /* check if exists already */
287          for (newmig = from_rbf->ejl; newmig != NULL;
288                          newmig = nextedge(from_rbf,newmig))
289                  if (newmig->rbfv[1] == to_rbf)
290                          return(NULL);
291                                                  /* else allocate */
292 + #ifdef DEBUG
293 +        fprintf(stderr, "Building path from (theta,phi) (%.0f,%.0f) ",
294 +                        get_theta180(from_rbf->invec),
295 +                        get_phi360(from_rbf->invec));
296 +        fprintf(stderr, "to (%.0f,%.0f) with %d x %d matrix\n",
297 +                        get_theta180(to_rbf->invec),
298 +                        get_phi360(to_rbf->invec),
299 +                        from_rbf->nrbf, to_rbf->nrbf);
300 + #endif
301          newmig = new_migration(from_rbf, to_rbf);
302          if (run_subprocess())
303                  return(newmig);                 /* child continues */
304 <        pmtx = price_routes(from_rbf, to_rbf);
304 >        price_routes(&pmtx, from_rbf, to_rbf);
305          src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
306          dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
307          if ((src_rem == NULL) | (dst_rem == NULL)) {
# Line 315 | Line 309 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
309                                  progname);
310                  exit(1);
311          }
318 – #ifdef DEBUG
319 –        fprintf(stderr, "Building path from (theta,phi) %s ",
320 –                        thetaphi(from_rbf->invec));
321 –        fprintf(stderr, "to %s", thetaphi(to_rbf->invec));
322 –        /* if (nchild) */ fputc('\n', stderr);
323 – #endif
312                                                  /* starting quantities */
313          memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
314          for (i = from_rbf->nrbf; i--; )
315                  src_rem[i] = rbf_volume(&from_rbf->rbfa[i]) / from_rbf->vtotal;
316 <        for (i = to_rbf->nrbf; i--; )
317 <                dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
316 >        for (j = to_rbf->nrbf; j--; )
317 >                dst_rem[j] = rbf_volume(&to_rbf->rbfa[j]) / to_rbf->vtotal;
318 >
319          do {                                    /* move a bit at a time */
320 <                move_amt = migration_step(newmig, src_rem, dst_rem, pmtx);
320 >                move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
321                  total_rem -= move_amt;
322 < #ifdef DEBUG
323 <                if (!nchild)
335 <                        /* fputc('.', stderr); */
336 <                        fprintf(stderr, "%.9f remaining...\r", total_rem);
337 < #endif
338 <        } while (total_rem > end_thresh && (total_rem > check_thresh) |
339 <                                        (move_amt > rel_thresh*total_rem));
340 < #ifdef DEBUG
341 <        if (!nchild) fputs("\ndone.\n", stderr);
342 <        else fprintf(stderr, "finished with %.9f remaining\n", total_rem);
343 < #endif
322 >        } while ((total_rem > end_thresh) & (move_amt > 0));
323 >
324          for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
325 <            float       nf = rbf_volume(&from_rbf->rbfa[i]);
346 <            int         j;
325 >            double      nf = rbf_volume(&from_rbf->rbfa[i]);
326              if (nf <= FTINY) continue;
327              nf = from_rbf->vtotal / nf;
328              for (j = to_rbf->nrbf; j--; )
329 <                newmig->mtx[mtx_ndx(newmig,i,j)] *= nf;
329 >                mtx_coef(newmig,i,j) *= nf;     /* row now sums to 1.0 */
330          }
331          end_subprocess();                       /* exit here if subprocess */
332 <        free(pmtx);                             /* free working arrays */
332 >        free_routes(&pmtx);                     /* free working arrays */
333          free(src_rem);
334          free(dst_rem);
335          return(newmig);

Diff Legend

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