/* Last edited on 2008-01-08 22:31:53 by stolfi */ /* TRAVEL PROBABILITY MATRICES A /visit intensity matrix/ is a mobility matrix {mm} where the value of {mm[A,B]} is the expected number of times that site {B} will be encoutered during a random exploratory walk with certain parameters that starts at site {A}. We assume that, once the traveler gives up on a path, he returns home by retracing his steps; so we may as well count only the sites seen in the forward path. {pr(m,A,B)} for such an elementary move is a decreasing function of the relative difficulty of that move; which in turn depends on the mean altitude of the two sites (including whether the sites are on land or water), the slope along the move, and possibly other geographic factors (latitude, climate, etc.). In an elementary matrix, the probabilities {pr(m,A,A)} of staying put are zero. In a /transitive/ mobility matrix, the edges represent paths of zero or more elementary steps. The probability {pr(m,A,B)} may depend on the minimum difficulty of any path between them (the /min-diff/ criterion), or on the probability of reaching {B} from {A} by a random walk that starts at {A} with some stopping criterion (the /rand-walk/ criterion). which has a directed edge from site {A} to site {B} iff there is a non-negligible probability that a person or family in {A} will move to {B} in one random /travel step/. The edge is labeled with a weight directly proportional to the . */ /* ---------------------------------------------------------------------- */ /* For instance, the mobility matrix may be used to store the /a priori/ probability of that travel occurring in certain situations, ignoring any population-dependent factors. For a more especific example, the element {M[A,B]} of a mobility matrix {M} may represent the /a priori/ probability that a twin couple which was born at site {A} will consider settling down at site {B}, before checking whether {B} is occupied or not. In this example, to get the actual probabilities of the couple settling down at each site {B}, one would typically multiply the values {M[A,B]} by other factors that depend on the current occupants of {B}, and re-normalize the results to unit sum. */ /* ---------------------------------------------------------------------- */ double tr = 0.25; /* Transparency of language colors. */ double op = 1 - tr; /* Opacity of language colors. */ rgb[0] = tr*rgb[0] + op*rgbv[0]; rgb[1] = tr*rgb[1] + op*rgbv[1]; rgb[2] = tr*rgb[2] + op*rgbv[2]; /* ---------------------------------------------------------------------- */ #define langev_segregation_INFO \ "Segregation\n" \ "\n" \ " After reaching {pc}, the children move a short distance towards the" \ " centroid of adult people living near {pc} that speak languages similar to" \ " their own. This second stage of the migration takes place exclusively" \ " on dry land, landing at some dry land site {pd}.\n" \ "\n" \ /* ---------------------------------------------------------------------- */ void throw_in_disk(double radius, double *x, double *y); /* Generates a random point {(*x,*y)} in the dist of given {radius} centered at the origin. */ void throw_in_disk(double radius, double *rx, double *ry) { double x, y; do { x = 2*drandom() - 1; y = 2*drandom() - 1; } while (hypot(x,y) > 1.0); (*rx) = radius*x; (*ry) = radius*y; } /* ---------------------------------------------------------------------- */ /* Distribution of site of segregation: */ spm->segregate = new_distr_params ( /*water_rel_diff:*/ 1.0, /* Moving across across shallow water is OK. */ /*mean_trip_diff:*/ 2.0, /* Segregate by moving about this many cells. */ /*self:*/ 0.500, /* Home is not a bad place to stay after all. */ /*new_state:*/ TRUE, /* Segregation depends on new state. */ /*mean_lang_diff:*/ 0.20, /* Tolerable lang diff of peers. */ /*vacant:*/ 0.00 /* Vacant sites have no peers. */ ); /* ---------------------------------------------------------------------- */ #define PROB_LANG_IMITATE 0.90 /* Probability of a language trait being absorbed from a neighbor (apart from the language similarity and closeness factors). */ #define PROB_GENE_MUTATE 0.05 /* Probability of each language trait mutating spontaneously. */ #define PROB_GENE_LINK (0.75) /* Probability that a genetic trait is inherited from the same parent as the preceding adjacent trait. */ lang_t mix_languages(lang_t vc, qentry_t qt[], double wt[], int nw) { int c[3]; language_decode(vc, c); int i; for (i = 0; i < 3; i++) { if (drandom() < PROB_IMITATE) { int k = pick_elem(wt, nw); lang_t vt = qt[k].v; assert(vt != NULL_LANG); int t[3]; language_decode(vt, t); c[i] = t[i]; } } return language_encode(c); } lang_t mutate_language(lang_t v) { int c[3]; language_decode(v, c); c[0] = mutate_letter(c[0], N0); c[1] = mutate_letter(c[1], N1); c[2] = mutate_letter(c[2], N2); return language_encode(c); } int mutate_letter(int c, int N) { double r = drandom(); if (r < PROB_MUTATE/2) { c--; if (c < 0) { c = 1; } } else if (r > 1 - PROB_MUTATE/2) { c++; if (c >= N) { c = N - 2; } } return c; } void map_language_to_color(lang_t v, double rgb[]) { int c[3]; language_decode(v, c); /* Pick a zero-brightness direction {rgb}: */ double f0 = ((double)c[0])/((double)N0); double f1 = ((double)c[1] + f0)/((double)N1); double ang = (2*M_PI)*f1; double p = cos(ang); double q = sin(ang); rgb[0] = p*(-0.001164) + q*(+1.139704); rgb[1] = p*(-0.395735) + q*(-0.580624); rgb[2] = p*(+2.030875) + q*(-0.000606); /* Pick a brightness {y}: */ double f2 = ((double)c[2])/((double)N2 - 1); double y = 0.450 + 0.150*f2; /* Find the max scaling {s} such that {(y,y,y) + rgb} is valid: */ double s = +MAXDOUBLE; int i; for (i = 0; i < 3; i++) { if (rgb[i] > 0.0) { s = fmin(s, (1.0 - y)/rgb[i]); } if (rgb[i] < 0.0) { s = fmin(s, (0.0 - y)/rgb[i]); } } /* Combine {(y,y,y)} with {rgb} scaled by {s} and {sat}: */ double g0 = ((double)c[0])/((double)N0 - 1); double sat = 0.300 + 0.700*g0; for (i = 0; i < 3; i++) { rgb[i] = y + s*sat*rgb[i]; } } /* ---------------------------------------------------------------------- */ void gather_nearby_sites ( frame_t *frm, /* A frame. */ pnm_image_t *alt, /* Altitude map. */ pnm_sample_t level, /* Sea level. */ int xini, /* X coord of starting site. */ int yini, /* Y coord of starting site. */ step_diff_func_t st_diff, /* A function that defines the remoteness of a step. */ diff_weight_func_t diff_wt, /* A remoteness-related weight function. */ site_weight_func_t site_wt, /* A site-related weight function. */ qentry_vec_t *Q, /* (OUT) List of reachable sites. */ double_vec_t *W, /* (OUT) Weights (unnormalized) of those sites. */ int *nQ /* (OUT) Number of neighboring sites found. */ ) { int nx = frm->cols; int ny = frm->rows; assert((xini >= 0) & (xini < nx)); assert((yini >= 0) & (yini < ny)); /* Uses Dijkstra's algorithm. */ int nf = 0; /* Number of finalized sites. */ int nq = 0; /* Number sites seen so far. */ /* Sites {Q[0..nq-1]} are those seen so far, sorted by incr. remoteness. */ /* Sites {Q[0..nf-1]} are finalized and will not move. */ /* Sites {Q[nf..nq-1]} are not finalized and may move. */ auto void enqueue_and_mark(int xt, int yt, double dt, int pt); /* Appends the site {(xt,yt)} to the queue {Q} and sorts it in its proper place. Assumes that the easiest path to that site has remoteness {dt} and reaches it from site {Q[pt]}, which must be in the fixed part of {Q}. Also marks the site as seen (by ORing its language to {MARK_LANG}). */ auto bool_t is_marked(int xt, int yt); /* TRUE iff the site {(xt,yt)} has been seen (and added to {Q}). */ auto void clear_mark(int xt, int yt); /* Clears the mark of site {(xt,yt)}. */ /* Initialize queue with starting site {(x,y)}: */ enqueue_and_mark(xini, yini, 0.0, -1); double wtot = 0.0; /* Total weight of finalized sites. */ /* Dijkstra's loop: */ while (nf < nq) { /* Get the index {q} of the easiest site in queue: */ int q = nf; /* Current pivot site. */ /* Any unseen site has remoteness {Qd[q]} or more, so {Q[q]} is final: */ nf++; /* Get hold of data about {Q[q]}: */ qentry_t *Qq = &(Q->el[q]); int xq = Qq->x, yq = Qq->y; /* Compute weight factor {wdiff} due to trip remoteness: */ double wdiff = diff_wt(Qq->diff); demand(wdiff >= 0, "bad diff_wt()"); /* Compute weight factor {wsite} due to other factors: */ double wsite = site_wt(Qq->diff, xq, yq); demand((wsite >= 0) && (wsite <= 1), "bad site_wt()"); /* Set weight of {Q[q]} (which will not move any more): */ double wq = wdiff*wsite; double_vec_expand(W, q); W->el[q] = wq; /* Update total weight of fixed sites: */ wtot += wq; /* Estimate max total weight of unfixed sites: */ double wrem = (nx*ny - nq)*wdiff; /* If we got most of the available weight, stop: */ if (wrem < 0.001*wtot) { break; } /* Check sites adjacent to {(xq,yq)}: */ int sx, sy; for (sx = -1; sx <= +1; sx++) { for (sy = -1; sy <= +1; sy++) { /* Compute neighbor's X coord {xp}, in cylindrical topology: */ int xp = (xq + sx + nx) % nx; assert((xp >= 0) && (xp < nx)); /* Compute neighbor's Y coord {yp}: */ int yp = yq + sy; if ((xp != xq) && (yp != yq) && (yp >= 0) && (yp < ny) && (! is_marked(xp, yp))) { /* Site {(xp,yp)} exists and is not {(xq,yq)}. */ /* Compute the remoteness of the step from {(xq,yq)} to {(xp,yp)}: */ double dstep = st_diff(xq, yq, xp, yp); /* Compute the remoteness of the path from start to {(xp,yp)}: */ double dp = Qq->diff + dstep; /* Add {(xp,yp)} to queue: */ enqueue_and_mark(xp, yp, dp, q); } } } } /* Restore the frame to its original state, by clearing all marks: */ int k; for (k = 0; k < nq; k++) { qentry_t *Qk = &(Q->el[k]); clear_mark(Qk->x, Qk->y); } /* Return the sites that were fixed: */ (*nQ) = nf; return; bool_t is_marked(int xt, int yt) { lang_t vt = &(frm->site[yt*nx + xt].lang); return ((vt & MARK_LANG) != 0); } void enqueue_and_mark(int xt, int yt, double dt, int pt) { assert(pt < nf); /* Parent must be a fixed node. */ int k = nq; nq++; qentry_vec_expand(Q, k); qentry_t *Qk = &(Q->el[k]); Qk->x = xt; Qk->y = yt; Qk->diff = dt; Qk->prev = pt; /* Save language and mark site as seen: */ sibs_t *cpt = &(frm->site[yt*nx + xt]); assert((cpt->lang & MARK_LANG) == 0); cpt->lang |= MARK_LANG; /* Sorts the site into its proper place in the queue: */ while ((k > 0) && (Q->el[k].diff < Q->el[k-1].diff)) { assert(k > nf); qentry_t qt = Q->el[k-1]; Q->el[k-1] = Q->el[k]; Q->el[k] = qt; } } void clear_mark(int xt, int yt) { sibs_t *cpt = &(frm->site[yt*nx + xt]); assert((cpt->lang & MARK_LANG) != 0); cpt->lang &= (~ MARK_LANG); } } /* ---------------------------------------------------------------------- */ /* Wander off to a nearby dry land site, possibly over shallow water: */ while (try < maxtry) { try++; int xw = x + abrandom(-1, +1); while(xw < 0) { xw += nx; } while(xw >= nx) { xw -= nx; } int yw = y + abrandom(-1, +1); if (yw < 0) { yw = 0; } if (yw >= nx) { yw = nx - 1; } double hw = get_altitude(alt, level, xw, yw); if (debug) { fprintf(stderr, " checking (%d,%d) alt = %+6.3f\n", xw, yw, hw); } if (hw >= 0.0) { /* Dry land or shallow water, children can move there: */ x = xw; y = yw; if (hw > 0.0) { /* Reached dry land, try setting there: */ break; } } } /* ---------------------------------------------------------------------- */ double dtot = 0.0; while (dtot < d) { /* Try to take a step in the direction {(dx,dy)}: */ int nc = 0; /* Number of candidate steps. */ int xc[9], yc[9]; /* Candidate number {k} is {(xc[k],yc[k])}. */ double dc[9]; /* {dc[k]} is the difficulty of the step to cand. number {k}. */ double hc[9]; /* {hc[k]} is the altitude after step to cand. number {k}. */ double wc[9]; /* {wc[k]} is the prob. weight (unnorm.) of cand. number {k}. */ /* Gather candidates for the step (including `stay put'): */ int sx, sy; for (sx = -1; sx <= +1; sx++) { for (sy = -1; sy <= +1; sy++) { int xp = (xq + sx + nx) % nx; int yp = yq + sy; assert((xp >= 0) && (xp < nx)); if ((yp >= 0) && (yp < ny)) { double hp = get_altitude(alt, level, xp, yp); if (hp >= 0.0) { /* Site {(xp,yp)} is on land or shallow water. */ /* Get the difficulty of the step from {(x,y)} to {(xp,yp)}: */ double dstep = step_difficulty(alt, level, xq, yq, xp, yp); /* Compute the difficulty factor in the step's probability: */ double udif = dstep/MEAN_STEP_DIFF; double wdif = exp(-M_LN2*udif*udif); /* Compute the direction factor on the step's probability: */ double wdir = fmax(0.01, (sx*dx + sy*dy)/hypot(sx,sy)/hypot(dx,dy)); /* Save the step and its difficulty: */ xc[nc] = xp; yc[nc] = yp; dc[nc] = dstep; hc[nc] = hp; wc[nc] = wdif*wdir; nc++; } } } } /* Pick a step: */ int k = pick_elem(wc, nc); /* Take that step: */ xq = xc[k]; yq = yc[k]; /* Remember last dry land site: */ if (hc[k] > 0.0) { (*x) = xq; (*y) = yq; } /* Update the total difficulty of trip (including elapsed time): */ dtot += fmax(dc[k], 0.01); } /* ---------------------------------------------------------------------- */ #define MAXCHANNELS 3 #define N0 3 #define N1 6 #define N2 3 static int NC[3] = { N0, N1, N2 }; /* Num alternatives in each letter slot. */ #define NL (N0*N1*N2) /* Number of distinct languages. */ /* ---------------------------------------------------------------------- */ void write_pnm_image(FILE *wr, pnm_image_t *im) { int nx = im->cols; int ny = im->rows; int chns = im->cols; int format = 0; bool_t forceplain = FALSE; pnm_sample_t maxmaxval = 0; switch (chns) { case 1: maxmaxval = PGM_OVERALLMAXVAL; format = RPGM_FORMAT; break; case 3: maxmaxval = PPM_MAXMAXVAL; format = RPPM_FORMAT; break; default: pm_error("bad number of channels"); } if (im->maxval > maxmaxval) { pm_error("maxval must be at most %d", maxmaxval); } xel* otp; int y, x; xel xelrow[im->cols]; pnm_writepnminit(wr, nx, ny, im->maxval, format, forceplain); for (y = 0; y < ny; ++y) { pnm_sample_t *imrow = &(im->el[y*nx*chns]); for (x = 0, otp = &(xelrow[0]); x < nx; ++x, ++otp) { pnm_sample_t *imv = &(imrow[x*chns]); switch (PNM_FORMAT_TYPE(format)) { case PPM_TYPE: PPM_ASSIGN(*otp, imv[0], imv[1], imv[2]); break; case PGM_TYPE: PNM_ASSIGN1(*otp, imv[0]); break; } } pnm_writepnmrow(wr, xelrow, im->cols, im->maxval, format, forceplain); } fflush(wr); } pnm_image_t *read_pnm_image(FILE *rd) { int nx, ny; pnm_sample_t maxval; int format; pnm_readpnminit(rd, &nx, &ny, &maxval, &format); fprintf ( stderr, "input format = %c%c maxval = %d nx = %d ny = %d\n", (format >> 8)&255, format&255, maxval, ny, nx ); /* Decide number of channels: */ int chns = 0; switch (format) { case PPM_FORMAT: case RPPM_FORMAT: chns = 3; break; case PGM_FORMAT: case RPGM_FORMAT: case PBM_FORMAT: case RPBM_FORMAT: chns = 1; break; default: pm_error("bad input file format"); } pnm_image_t *im = allocate_image(chns, nx, ny, maxval); read_pnm_file_pixels(rd, format, im); return im; } pnm_image_t *allocate_image(int chns, int nx, int ny, pnm_sample_t maxval) { pnm_image_t *im; im = (pnm_image_t*)malloc(sizeof(pnm_image_t)); if (im == NULL) { pm_error("out of memory for image descriptor"); } im->sz[0] = chns; im->cols = nx; im->rows = ny; im->maxval = maxval; /* Allocate pixel arrays for floatized input image: */ im->el = (pnm_sample_t *)malloc(ny*nx*chns*sizeof(pnm_sample_t)); if (im->el == NULL) { pm_error("out of memory for image pixels"); } return im; } void read_pnm_file_pixels(FILE *f, int format, pnm_image_t *im) { int nx = im->cols; int ny = im->rows; int chns = im->sz[0]; xel xelrow[im->cols]; /* Scanline buffer. */ register xel *pp; /* Scans input pixels. */ int y, x; for (y = 0; y < ny; ++y) { pnm_readpnmrow(f, xelrow, im->cols, im->maxval, format); pnm_sample_t *imrow = &(im->el[y*nx*chns]); for (x = 0, pp = xelrow; x < im->cols; ++x, ++pp) { pnm_sample_t *imv = &(imrow[x*chns]); switch (format) { case PPM_FORMAT: case RPPM_FORMAT: imv[0] = PPM_GETR(*pp); imv[1] = PPM_GETG(*pp); imv[2] = PPM_GETB(*pp); break; case PGM_FORMAT: case RPGM_FORMAT: case PBM_FORMAT: case RPBM_FORMAT: imv[0] = PNM_GET1(*pp); break; default: pm_error("bad input file format"); } } } } /* ---------------------------------------------------------------------- */ int mutate_letter(int c, int N); /* Changes the letter {c} (assumed to range in {0..N-1}) to an adjecent letter, with probability {PROB_MUTATE}. */ void language_decode(lang_t v, int c[]); /* Breaks the language code {v} into its separate traits {c[0..2]}, where {c[i]} lies in {0..NC[i]-1}. */ lang_t language_encode(int c[]); /* Combines the traits {c[0..2]}, where {c[i]} is assumed to be in {0..NC[i]-1}, into a language code. */ void language_decode(lang_t v, int c[]) { assert(v != NULL_LANG); assert(v != MARKESE); int i; for (i = 2; i >= 0; i--) { c[i] = v % NC[i]; v /= NC[i]; } assert(v == 0); } lang_t language_encode(int c[]) { lang_t v = 0; int i; for (i = 2; i >= 0; i--) { assert((c[i] >= 0) && (c[i] < NC[i])); v = v*NC[i] + c[i]; } return v; } /* ---------------------------------------------------------------------- */ #define PROB_INHERIT 0.1 /* Probability of learning language from a parent instead of from peers. */ int mix_letters(int c, int m, int f, int N); /* Returns the letter {c} (assumed to vary in {0..N-1}), modified with probability {PROB_INHERIT} by admixture of letters {f} or {m}, and randomly mutated with probability {PROB_MUTATE}. */ int mix_letters(int c, int m, int f, int N) { double r = drandom(); if (r < PROB_MUTATE) { /* Random mutation: */ r = r/PROB_MUTATE; if (r < 0.50) { c--; if (c < 0) { c = 0; } } else { c++; if (c >= N) { c = N - 1; } } return c; } r = (r - PROB_MUTATE)/(1.0 - PROB_MUTATE); if (r < PROB_INHERIT) { /* Inherit something from parent: */ r = r/PROB_INHERIT; if (r < 0.50) { c = (int)(0.75*m + 0.25*c + 0.50); } else { c = (int)(0.75*f + 0.25*c + 0.50); } if (c < 0) { c = 0; } if (c >= N) { c = N - 1; } return c; } return c; }