source: subversion/applications/rendering/gosmore/density6.cpp @ 28848

Last change on this file since 28848 was 19893, checked in by nic, 10 years ago

Rebuild: Filter out more bulk imports
Density / bboxes: Smaller boxes, better algorithm. Replace png imagemap with OpenLayers?
Windows: Improve installer script.
Better default background

  • Property svn:executable set to *
File size: 8.1 KB
Line 
1// Finds rectangles that cover the planet. Requires density.csv, which is created by
2// make CFLAGS="-DMKDENSITY_CSV -O2" && bzcat planet.osm.bz2 | ./gosmore sortRelations
3// g++ -O2 density6.cpp && time nice -n 19 ./a.out >density.txt
4// Run time is approximately 5 minutes and will require 5 GB of RAM
5
6#include <string.h>
7#include <stdlib.h>
8#include <stdio.h>
9#include <vector>
10#include <queue>
11#include <assert.h>
12using namespace::std;
13
14int b[4000000][8], cnt = 0;
15int dp[1025][1025], target, maks = 0, col;
16
17int BCmp (int *a, int *b)
18{
19  return *a - *b;
20}
21
22struct state {
23  long long nCovered;
24  int **s, **best, c, mActive;
25  state () {}
26};
27
28bool operator < (const state &a, const state &b)
29{
30//  return a.c == 0 ? false : b.c == 0 ? true : a.nCovered / a.c < b.nCovered / b.c;
31  int cadj = col > 800 ? 112165 : col > 620 ? 3165 : col < 140 ? 15 :
32    (col > 340 && col < 410) ? 1516 : 5565;
33  return (a.nCovered + 60000000) / (a.c + cadj) <
34         (b.nCovered + 60000000) / (b.c + cadj);
35}
36
37int IPtrCmp (int **a, int **b)
38{
39  return **b - **a;
40}
41
42vector<state> states;
43priority_queue<state> que;
44vector<int*> active;
45
46
47#define I(a,b,e,d) (dp[a+1][b+1] - dp[e][b+1] - dp[a+1][d] + dp[e][d])
48#define M 11000000
49
50void Try (void)
51{
52  while (states.size () < 2000 && !que.empty ()) {
53    //printf ("%lld %d\n", que.top ().nCovered, states.size ());
54    int rh[1024], i, j, **s = que.top ().s, l, c = que.top ().c;
55    int **best = que.top ().best, nCovered = que.top ().nCovered;
56    memset (rh, 0, sizeof (rh));
57    for (i = 0, l = 0; s[i]; i++) {
58      for (j = s[i][1]; j < s[i][3]; j++) {
59        if (rh[j] < s[i][2]) rh[j] = s[i][2];
60      }
61      if (s[i][2] > col) s[l++] = s[i]; // Remove rectangles we passed.
62    }
63    s[l] = NULL; // Remove rectangles we passed.
64    for (i = 0; i < 1024 && rh[i] > col; i++) {}
65    if (i == 1024) { // The state covers this column
66      states.push_back (state ()); //que.top ());
67      memcpy (&states.back (), &que.top (), sizeof (que.top ()));
68    }
69    //printf ("%d\n", i);
70//          assert (states.empty () || states[0].best[0]);
71    j = que.top ().mActive;
72    que.pop ();
73//          assert (states.empty () || states[0].best[0]);
74    for (; i < 1024 && j < active.size (); j++) { // 'i' is the first row not covered.
75      for (l = 0x3fffff; l >= active.size (); l/=2) {}
76      int *a = active[j < l ? j ^ (/*nCovered &*/(col*1717) & l) : j];
77      if (a[1] <= i && i < a[3]) {
78        int newNcov = nCovered;
79        for (l = a[1]; l < a[3]; l++) {
80          if (rh[l] < a[2]) newNcov += I (l, a[2] - 1, l,
81            (rh[l] < col ? col : rh[l]));
82        }
83       //assert (newNcov < nCovered + 14000000);
84        //assert (newNcov >= I (1023, col - 1, 0, 0));
85        //assert (a[2] > 800 ||newNcov <= I (1023, a[2] + 168, 0, 0));
86        //if (newNcov > nCovered + (col < 200 ? 2000000: col < 600 ? 10000*0:
87        //      col < 650 ? 200000:  col < 700 ? 1200000:col<780 ? 500000: col < 850 ? 100000: 10)) { // (c + 1) * 1000000) { //target) {
88        if (newNcov > nCovered + que.size () / (col < 512 ? 20 :10)) {
89          state nxt;
90          for (l = 0; s[l]; l++) {}
91          nxt.s = new int *[l+2];//(int**) malloc ((l + 2) * sizeof (*s));
92          nxt.s[0] = a;
93          memcpy (nxt.s + 1, s, (l + 1) * sizeof (*s));
94
95          for (l = 0; best[l]; l++) {}
96          assert (l ==c);
97          nxt.best = new int *[l+2]; //(int**) malloc ((l + 2) * sizeof (*s));
98          nxt.best[0] = a;
99          memcpy (nxt.best + 1, best, (l + 1) * sizeof (*s));
100          nxt.c = c + 1;
101          nxt.nCovered = newNcov;
102          nxt.mActive = j + 1;
103          if (0&&newNcov < maks) {
104            fprintf (stderr, "%d %d %d %d\n", newNcov, c, que.size (), i);
105            maks = newNcov;
106          }
107          que.push (nxt);
108//          assert (states.empty () || states[0].best[0]);
109        }
110      }
111    } // Look for a bbox that can cover the gap.
112    if (i < 1024) delete s; //free (s);
113//          assert (states.empty () || (states[0].best[0] && best != states[0].best));
114    if (i < 1024) delete best; //free (best);
115  } // While creating new states
116  for (; !que.empty (); que.pop ()) {
117    delete que.top().s; //free (que.top ().s);
118    delete que.top().best;//free (que.top ().best);
119  }
120}
121
122main ()
123{
124  FILE *in = fopen ("density.csv", "r");
125  FILE *bf = fopen ("b.bin", "r+");
126  if (!bf) bf = fopen ("b.bin", "w");
127//  memset (c, 0, sizeof (c));
128//  memset (dc, 0, sizeof (dc));
129  int d, i, j, k, l;
130  #if 1
131  for (i = 0; i < 1025; i++) dp[i][0] = dp[0][i] = 0;
132  for (i = 0; i < 1024; i++) {
133    for (j = 0; j < 1024; j++) {
134      fscanf (in, j == 1023 ? "%d\n" : "%d ", &d);
135      dp[i+1][j+1] = d + dp[i][j+1] + dp[i+1][j] - dp[i][j];
136    }
137  }
138  #endif
139  //printf ("%d\n", dp[1024][1024]);
140#if 1
141//#define O 2
142  for (i = 0; i < 1024; i++) {
143    for (j = 0; j < 1024; j++) {
144      int maxk = i < 415 && j > 412 ? 462 : i < 666 && j > 430 ? 712 : 1024, l2;
145      // Reject rectangles that cross the South Atlantic or the Indian Ocean
146      for (k = i, l2 = 1023; k < maxk; k++) {
147        int mini = j > 430 && k > 712 ? 666 : j > 412 && k > 462 ? 415 : 0;
148        maxk = min (maxk, i + 3 * (l2 - j));
149        for (; l2 >= j && I (l2, k, j, i) > M; l2--) {}
150        l = min (l2, j + 3 * (k - i));
151        if (l >= j && // Area > 0
152            (i == mini || I (l, k, j, i - 1) > M) && // Not expandable Westwards
153            (j == 0 || I (l, k, j - 1, i) > M) && // Not expandable Northwards
154            (k == maxk - 1 || I (l, k + 1, j, i) > M)) { // Not expandable Eastwards
155          int O = 1; //i > 150 && k < 320 ? 0 : 2;
156          // Some parts of US so crowded that 4x4 is too large
157          if ((l >= 1023 || j == 0 || l >= j + O*1) && 
158              (k >= 1023 || i == 0 || k >= i + O*1)) {
159            b[cnt][0] = i == 0 ? i : i + (k - i + O*4)/9;
160            b[cnt][1] = j == 0 ? j : j + (l - j + O*4)/9;
161            b[cnt][2] = k == 1023 ? k + 1 : k + 1 - (k - i + O*4)/9;
162            b[cnt][3] = l == 1023 ? l + 1 : l + 1 - (l - j + O*4)/9;
163            b[cnt][4] = i;
164            b[cnt][5] = j;
165            b[cnt][6] = k + 1;
166            b[cnt++][7] = l + 1;
167          }
168          //printf ("%4d %4d %4d %4d %d\n", i, j, k, l, I (l, k, j, i));
169        }
170      }
171    }
172  }
173  qsort (b, cnt, sizeof (b[0]), (int (*)(const void *, const void *)) BCmp);
174  fwrite (b, cnt, sizeof (b[0]), bf);
175#else 
176  cnt = fread (b, sizeof (b[0]), sizeof (b) / sizeof (b[0]), bf);
177#endif
178  fprintf (stderr, "%d bboxes\n", cnt);
179  states.push_back (state ());
180  states.back ().s = new int*[1];//(int**) calloc (sizeof (*states.back ().s), 1);
181  states.back ().s[0] = NULL;
182  states.back ().best = new int*[1]; //(int**) calloc (sizeof (*states.back ().best), 1);
183  states.back ().best[0] = NULL;
184  states.back ().nCovered = 0;
185  states.back ().c = 0;
186  for (col = 0, j = 0; col < 1024 && !states.empty (); col++) {
187    for (;j < cnt && b[j][0] <= col; j++) active.push_back (&b[j][0]);
188    target = 0; //10000000 / (state.size () + 3);
189    for (k = 0; k < states.size (); k++) {
190      target += (states[k].nCovered + 12000000 + states.size()*5000) / (states[k].c + 1) /
191        states.size ();
192      //printf ("%d %d\n", states[k].nCovered,
193    }
194    fprintf (stderr, "%4d %d active, %d states %d %d\n", col, active.size (), states.size (),
195      target, states[0].c);
196    for (k = states.size () - 1; k >= 0; k--) {
197      for (l = 0; states[k].s[l] && states[k].s[l][2] > col; l++) {}
198      if (states[k].s[l] || l == 0 /* bootstrap */) {
199        states[k].mActive = 0;
200        que.push (states[k]);
201        memcpy (&states[k], &states.back (), sizeof (states[k]));
202        states.pop_back ();
203      }
204    }
205    Try ();
206/*    for (k = states.size () - 1; k >= 0; k--) {
207          for (l = 0; states[k].best[l]; l++) {}
208          assert (l ==states[k].c);
209    }*/
210    for (k = active.size () - 1; k >= 0; k--) {
211      if (active[k][2] <= col + 1) {
212        active[k] = active.back();
213        active.pop_back ();
214      }
215    }
216  }
217  for (i = 1, j = 0; i < states.size (); i++) {
218    if (states[i].c < states[j].c) j = i;
219  }
220  fprintf (stderr, "%d rectangles\n", states[j].c);
221  for (i = 0; states[j].best[i]; i++) {
222    printf ("%d %d %d %d\n", states[j].best[i][4], states[j].best[i][5],
223      states[j].best[i][6], states[j].best[i][7]);
224  }
225 
226}
Note: See TracBrowser for help on using the repository browser.