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.3 by greg, Thu Nov 8 22:05:04 2012 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 + } PRICEMAT;                     /* sorted pricing matrix */
31 +
32 + #define pricerow(p,i)   ((p)->price + (i)*(p)->ncols)
33 + #define psortrow(p,i)   ((p)->sord + (i)*(p)->ncols)
34 +
35 + /* Create a new migration holder (sharing memory for multiprocessing) */
36 + static MIGRATION *
37 + new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
38 + {
39 +        size_t          memlen = sizeof(MIGRATION) +
40 +                                sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1);
41 +        MIGRATION       *newmig;
42 + #ifdef _WIN32
43 +        if (nprocs > 1)
44 +                fprintf(stderr, "%s: warning - multiprocessing not supported\n",
45 +                                progname);
46 +        nprocs = 1;
47 +        newmig = (MIGRATION *)malloc(memlen);
48 + #else
49 +        if (nprocs <= 1) {                      /* single process? */
50 +                newmig = (MIGRATION *)malloc(memlen);
51 +        } else {                                /* else need to share memory */
52 +                newmig = (MIGRATION *)mmap(NULL, memlen, PROT_READ|PROT_WRITE,
53 +                                                MAP_ANON|MAP_SHARED, -1, 0);
54 +                if ((void *)newmig == MAP_FAILED)
55 +                        newmig = NULL;
56 +        }
57 + #endif
58 +        if (newmig == NULL) {
59 +                fprintf(stderr, "%s: cannot allocate new migration\n", progname);
60 +                exit(1);
61 +        }
62 +        newmig->rbfv[0] = from_rbf;
63 +        newmig->rbfv[1] = to_rbf;
64 +                                                /* insert in edge lists */
65 +        newmig->enxt[0] = from_rbf->ejl;
66 +        from_rbf->ejl = newmig;
67 +        newmig->enxt[1] = to_rbf->ejl;
68 +        to_rbf->ejl = newmig;
69 +        newmig->next = mig_list;                /* push onto global list */
70 +        return(mig_list = newmig);
71 + }
72 +
73 + #ifdef _WIN32
74 + #define await_children(n)       (void)(n)
75 + #define run_subprocess()        0
76 + #define end_subprocess()        (void)0
77 + #else
78 +
79 + /* Wait for the specified number of child processes to complete */
80 + static void
81 + await_children(int n)
82 + {
83 +        int     exit_status = 0;
84 +
85 +        if (n > nchild)
86 +                n = nchild;
87 +        while (n-- > 0) {
88 +                int     status;
89 +                if (wait(&status) < 0) {
90 +                        fprintf(stderr, "%s: missing child(ren)!\n", progname);
91 +                        nchild = 0;
92 +                        break;
93 +                }
94 +                --nchild;
95 +                if (status) {                   /* something wrong */
96 +                        if ((status = WEXITSTATUS(status)))
97 +                                exit_status = status;
98 +                        else
99 +                                exit_status += !exit_status;
100 +                        fprintf(stderr, "%s: subprocess died\n", progname);
101 +                        n = nchild;             /* wait for the rest */
102 +                }
103 +        }
104 +        if (exit_status)
105 +                exit(exit_status);
106 + }
107 +
108 + /* Start child process if multiprocessing selected */
109 + static pid_t
110 + run_subprocess(void)
111 + {
112 +        int     status;
113 +        pid_t   pid;
114 +
115 +        if (nprocs <= 1)                        /* any children requested? */
116 +                return(0);
117 +        await_children(nchild + 1 - nprocs);    /* free up child process */
118 +        if ((pid = fork())) {
119 +                if (pid < 0) {
120 +                        fprintf(stderr, "%s: cannot fork subprocess\n",
121 +                                        progname);
122 +                        exit(1);
123 +                }
124 +                ++nchild;                       /* subprocess started */
125 +                return(pid);
126 +        }
127 +        nchild = -1;
128 +        return(0);                              /* put child to work */
129 + }
130 +
131 + /* If we are in subprocess, call exit */
132 + #define end_subprocess()        if (nchild < 0) _exit(0); else
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 float *
152 < price_routes(const RBFNODE *from_rbf, const RBFNODE *to_rbf)
151 > static void
152 > price_routes(PRICEMAT *pm, const RBFNODE *from_rbf, const RBFNODE *to_rbf)
153   {
30 –        float   *pmtx = (float *)malloc(sizeof(float) *
31 –                                        from_rbf->nrbf * to_rbf->nrbf);
154          FVECT   *vto = (FVECT *)malloc(sizeof(FVECT) * to_rbf->nrbf);
155          int     i, j;
156  
157 <        if ((pmtx == NULL) | (vto == NULL)) {
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);
# Line 44 | Line 171 | price_routes(const RBFNODE *from_rbf, const RBFNODE *t
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 <                pmtx[i*to_rbf->nrbf + j] = acos(DOT(vfrom, vto[j])) +
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);
52 –        return(pmtx);
182   }
183  
184 < /* Comparison routine needed for sorting price row */
185 < static const float      *price_arr;
186 < static int
58 < msrt_cmp(const void *p1, const void *p2)
184 > /* Free price matrix */
185 > static void
186 > free_routes(PRICEMAT *pm)
187   {
188 <        float   c1 = price_arr[*(const int *)p1];
189 <        float   c2 = price_arr[*(const int *)p2];
62 <
63 <        if (c1 > c2) return(1);
64 <        if (c1 < c2) return(-1);
65 <        return(0);
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 float *price, int n)
194 > min_cost(double amt2move, const double *avail, const PRICEMAT *pm, int s)
195   {
72 –        static int      *price_sort = NULL;
73 –        static int      n_alloc = 0;
196          double          total_cost = 0;
197 <        int             i;
197 >        int             j;
198  
199          if (amt2move <= FTINY)                  /* pre-emptive check */
200                  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);
201                                                  /* move cheapest first */
202 <        for (i = 0; i < n && amt2move > FTINY; i++) {
203 <                int     d = price_sort[i];
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 * price[d];
206 >                total_cost += amt * pricerow(pm,s)[d];
207                  amt2move -= amt;
208          }
209          return(total_cost);
# Line 103 | Line 211 | min_cost(double amt2move, const double *avail, const f
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 float *pmtx)
214 > migration_step(MIGRATION *mig, double *src_rem, double *dst_rem, const PRICEMAT *pm)
215   {
216          const double    maxamt = .1;
217 <        const double    minamt = maxamt*.0001;
217 >        const double    minamt = maxamt*5e-6;
218          static double   *src_cost = NULL;
219          static int      n_alloc = 0;
220          struct {
# Line 116 | Line 224 | migration_step(MIGRATION *mig, double *src_rem, double
224          } cur, best;
225          int             i;
226  
227 <        if (mtx_nrows(mig) > n_alloc) {         /* allocate cost array */
227 >        if (pm->nrows > n_alloc) {              /* allocate cost array */
228                  if (n_alloc)
229                          free(src_cost);
230 <                src_cost = (double *)malloc(sizeof(double)*mtx_nrows(mig));
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 = mtx_nrows(mig);
236 >                n_alloc = pm->nrows;
237          }
238 <        for (i = mtx_nrows(mig); i--; )         /* starting costs for diff. */
239 <                src_cost[i] = min_cost(src_rem[i], dst_rem,
132 <                                        pmtx+i*mtx_ncols(mig), mtx_ncols(mig));
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 = mtx_nrows(mig); cur.s--; ) {
244 <            const float *price = pmtx + cur.s*mtx_ncols(mig);
243 >        for (cur.s = pm->nrows; cur.s--; ) {
244 >            const float *price = pricerow(pm,cur.s);
245              double      cost_others = 0;
246 <            if (src_rem[cur.s] < minamt)
246 >            if (src_rem[cur.s] <= minamt)
247                      continue;
248              cur.d = -1;                         /* examine cheapest dest. */
249 <            for (i = mtx_ncols(mig); i--; )
249 >            for (i = pm->ncols; i--; )
250                  if (dst_rem[i] > minamt &&
251                                  (cur.d < 0 || price[i] < price[cur.d]))
252                          cur.d = i;
# Line 151 | Line 258 | migration_step(MIGRATION *mig, double *src_rem, double
258                                  src_rem[cur.s] : dst_rem[cur.d];
259              if (cur.amt > maxamt) cur.amt = maxamt;
260              dst_rem[cur.d] -= cur.amt;          /* add up differential costs */
261 <            for (i = mtx_nrows(mig); i--; )
261 >            for (i = pm->nrows; i--; )
262                  if (i != cur.s)
263 <                        cost_others += min_cost(src_rem[i], dst_rem,
157 <                                                price, mtx_ncols(mig))
263 >                        cost_others += min_cost(src_rem[i], dst_rem, pm, i)
264                                          - src_cost[i];
265              dst_rem[cur.d] += cur.amt;          /* undo trial move */
266              cur.price += cost_others/cur.amt;   /* adjust effective price */
# Line 164 | Line 270 | migration_step(MIGRATION *mig, double *src_rem, double
270          if ((best.s < 0) | (best.d < 0))
271                  return(.0);
272                                                  /* make the actual move */
273 <        mig->mtx[mtx_ndx(mig,best.s,best.d)] += best.amt;
273 >        mtx_coef(mig,best.s,best.d) += best.amt;
274          src_rem[best.s] -= best.amt;
275          dst_rem[best.d] -= best.amt;
276          return(best.amt);
# Line 185 | Line 291 | thetaphi(const FVECT v)
291   }
292   #endif
293  
188 – /* Create a new migration holder (sharing memory for multiprocessing) */
189 – static MIGRATION *
190 – new_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
191 – {
192 –        size_t          memlen = sizeof(MIGRATION) +
193 –                                sizeof(float)*(from_rbf->nrbf*to_rbf->nrbf - 1);
194 –        MIGRATION       *newmig;
195 – #ifdef _WIN32
196 –        if (nprocs > 1)
197 –                fprintf(stderr, "%s: warning - multiprocessing not supported\n",
198 –                                progname);
199 –        nprocs = 1;
200 –        newmig = (MIGRATION *)malloc(memlen);
201 – #else
202 –        if (nprocs <= 1) {                      /* single process? */
203 –                newmig = (MIGRATION *)malloc(memlen);
204 –        } else {                                /* else need to share memory */
205 –                newmig = (MIGRATION *)mmap(NULL, memlen, PROT_READ|PROT_WRITE,
206 –                                                MAP_ANON|MAP_SHARED, -1, 0);
207 –                if ((void *)newmig == MAP_FAILED)
208 –                        newmig = NULL;
209 –        }
210 – #endif
211 –        if (newmig == NULL) {
212 –                fprintf(stderr, "%s: cannot allocate new migration\n", progname);
213 –                exit(1);
214 –        }
215 –        newmig->rbfv[0] = from_rbf;
216 –        newmig->rbfv[1] = to_rbf;
217 –                                                /* insert in edge lists */
218 –        newmig->enxt[0] = from_rbf->ejl;
219 –        from_rbf->ejl = newmig;
220 –        newmig->enxt[1] = to_rbf->ejl;
221 –        to_rbf->ejl = newmig;
222 –        newmig->next = mig_list;                /* push onto global list */
223 –        return(mig_list = newmig);
224 – }
225 –
226 – #ifdef _WIN32
227 – #define await_children(n)       (void)(n)
228 – #define run_subprocess()        0
229 – #define end_subprocess()        (void)0
230 – #else
231 –
232 – /* Wait for the specified number of child processes to complete */
233 – static void
234 – await_children(int n)
235 – {
236 –        int     exit_status = 0;
237 –
238 –        if (n > nchild)
239 –                n = nchild;
240 –        while (n-- > 0) {
241 –                int     status;
242 –                if (wait(&status) < 0) {
243 –                        fprintf(stderr, "%s: missing child(ren)!\n", progname);
244 –                        nchild = 0;
245 –                        break;
246 –                }
247 –                --nchild;
248 –                if (status) {                   /* something wrong */
249 –                        if ((status = WEXITSTATUS(status)))
250 –                                exit_status = status;
251 –                        else
252 –                                exit_status += !exit_status;
253 –                        fprintf(stderr, "%s: subprocess died\n", progname);
254 –                        n = nchild;             /* wait for the rest */
255 –                }
256 –        }
257 –        if (exit_status)
258 –                exit(exit_status);
259 – }
260 –
261 – /* Start child process if multiprocessing selected */
262 – static pid_t
263 – run_subprocess(void)
264 – {
265 –        int     status;
266 –        pid_t   pid;
267 –
268 –        if (nprocs <= 1)                        /* any children requested? */
269 –                return(0);
270 –        await_children(nchild + 1 - nprocs);    /* free up child process */
271 –        if ((pid = fork())) {
272 –                if (pid < 0) {
273 –                        fprintf(stderr, "%s: cannot fork subprocess\n",
274 –                                        progname);
275 –                        exit(1);
276 –                }
277 –                ++nchild;                       /* subprocess started */
278 –                return(pid);
279 –        }
280 –        nchild = -1;
281 –        return(0);                              /* put child to work */
282 – }
283 –
284 – /* If we are in subprocess, call exit */
285 – #define end_subprocess()        if (nchild < 0) _exit(0); else
286 –
287 – #endif  /* ! _WIN32 */
288 –
294   /* Compute and insert migration along directed edge (may fork child) */
295   static MIGRATION *
296   create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
297   {
298 <        const double    end_thresh = 0.1/(from_rbf->nrbf*to_rbf->nrbf);
299 <        const double    check_thresh = 0.01;
295 <        const double    rel_thresh = 5e-6;
296 <        float           *pmtx;
298 >        const double    end_thresh = 5e-6;
299 >        PRICEMAT        pmtx;
300          MIGRATION       *newmig;
301          double          *src_rem, *dst_rem;
302          double          total_rem = 1., move_amt;
# Line 307 | Line 310 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
310          newmig = new_migration(from_rbf, to_rbf);
311          if (run_subprocess())
312                  return(newmig);                 /* child continues */
313 <        pmtx = price_routes(from_rbf, to_rbf);
313 >        price_routes(&pmtx, from_rbf, to_rbf);
314          src_rem = (double *)malloc(sizeof(double)*from_rbf->nrbf);
315          dst_rem = (double *)malloc(sizeof(double)*to_rbf->nrbf);
316          if ((src_rem == NULL) | (dst_rem == NULL)) {
# Line 318 | Line 321 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
321   #ifdef DEBUG
322          fprintf(stderr, "Building path from (theta,phi) %s ",
323                          thetaphi(from_rbf->invec));
324 <        fprintf(stderr, "to %s", thetaphi(to_rbf->invec));
325 <        /* if (nchild) */ fputc('\n', stderr);
324 >        fprintf(stderr, "to %s with %d x %d matrix\n",
325 >                        thetaphi(to_rbf->invec),
326 >                        from_rbf->nrbf, to_rbf->nrbf);
327   #endif
328                                                  /* starting quantities */
329          memset(newmig->mtx, 0, sizeof(float)*from_rbf->nrbf*to_rbf->nrbf);
# Line 328 | Line 332 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
332          for (i = to_rbf->nrbf; i--; )
333                  dst_rem[i] = rbf_volume(&to_rbf->rbfa[i]) / to_rbf->vtotal;
334          do {                                    /* move a bit at a time */
335 <                move_amt = migration_step(newmig, src_rem, dst_rem, pmtx);
335 >                move_amt = migration_step(newmig, src_rem, dst_rem, &pmtx);
336                  total_rem -= move_amt;
337   #ifdef DEBUG
338                  if (!nchild)
339 <                        /* fputc('.', stderr); */
336 <                        fprintf(stderr, "%.9f remaining...\r", total_rem);
339 >                        fprintf(stderr, "\r%.9f remaining...", total_rem);
340   #endif
341 <        } while (total_rem > end_thresh && (total_rem > check_thresh) |
339 <                                        (move_amt > rel_thresh*total_rem));
341 >        } while ((total_rem > end_thresh) & (move_amt > 0));
342   #ifdef DEBUG
343 <        if (!nchild) fputs("\ndone.\n", stderr);
343 >        if (!nchild) fputs("done.\n", stderr);
344          else fprintf(stderr, "finished with %.9f remaining\n", total_rem);
345   #endif
346          for (i = from_rbf->nrbf; i--; ) {       /* normalize final matrix */
# Line 347 | Line 349 | create_migration(RBFNODE *from_rbf, RBFNODE *to_rbf)
349              if (nf <= FTINY) continue;
350              nf = from_rbf->vtotal / nf;
351              for (j = to_rbf->nrbf; j--; )
352 <                newmig->mtx[mtx_ndx(newmig,i,j)] *= nf;
352 >                mtx_coef(newmig,i,j) *= nf;
353          }
354          end_subprocess();                       /* exit here if subprocess */
355 <        free(pmtx);                             /* free working arrays */
355 >        free_routes(&pmtx);                     /* free working arrays */
356          free(src_rem);
357          free(dst_rem);
358          return(newmig);

Diff Legend

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