/* Last edited on 2005-03-16 14:46:38 by stolfi */ /* Sweep algorithm for rasterization. */ auto int cmpy (int i, int j); /* Compares points {p[i]} and {p[j]} in Y-coordinate, then X: */ int cmpy (int i, int j) { return cmpy(&(p.el[i]), &(p.el[j])); } /* Sort vertices by increasing {y}. */ int ix[np]; int i; for (i = 0; i < np; i++) { ix[i] = i; } isrt_heapsort(ix, np, cmpy, +1); /* List of active edges: */ int ne = 0; /* Num of active edges, indexed {0..ne-1}. */ int v[n]; /* Lower vertex of edge {k} is {p[v[k]]}. */ int u[n]; /* Upper vertex of edge {k} is {p[u[k]]}. */ /* Agenda of future events (a heap sorted by Y): */ int_vec_t ag = int_vect_new(2*np); /* In theory {np^2}, but... */ int na = 0; /* Insert vertex events in agenda: */ /* Sweep with line in {+Y} direction: */ double ymin = Y(p.el[ix[0]]); /* Minimum Y of polygon. */ double ymax = p.el[ix[np-1]].c[1]; /* Maximum Y of polygon. */ int iy = (int)floor((ymin - yorg)/dy); /* Index of current scan line. */ int ia = 0; /* Vertices {ix[ia..np-1]} haven't been swept over yet. */ int ns = 0; /* next vertex to sweep over is {p[ix[ns]]} */ while (iy < ny) { /* Compute Y-coord of row {iy}: */ double y = Y(s->org) + iy*dy; /* Process events with ordinate {\leq y}: */ while ((ia < np) && (Y(p.el[ix[ia]]) <= y)) { /* Sweep over vertex {p[ix[ia]]}: */ cpk_update_active_list(u, v, &ne, p, ix[ia]); ia++; } /* Rasterize row: */ typedef hxg_color_t hxg_color_fn_t(hxg_pixel_t v); /* A procedure that selects a color for a canvas pixel based on its value {v}. */ typedef double hxg_radius_fn_t(hxg_pixel_t v); /* A procedure that selects the radius of a canvas pixel based on its value {v}. */ ---------------------------------------------------------------------- void LL_to_EUTM_orig(double Lat, double Lon, double *dx, double *dy, double refLon); void EUTM_to_LL_orig (double dx, double dy, double *Lat, double *Lon, double refLon); /* ROTINAS ORIGINAIS (QUASE) */ void LL_to_EUTM_orig(double Lat, double Lon, double *dx, double *dy, double refLon) { double LatRad; double e2; double e2P; double k0; double a; double N; double T; double C; double A; double M; a = EquatorialRadius; e2 = EccentricitySquared; k0 = 0.9996; LatRad = Lat*deg2rad; e2P = (e2)/(1-e2); N = a/sqrt(1-e2*sin(LatRad)*sin(LatRad)); T = tan(LatRad)*tan(LatRad); C = e2P*cos(LatRad)*cos(LatRad); A = cos(LatRad)*(Lon-refLon)*deg2rad; M = a * ( ( 1 - e2/4 - 3*e2*e2/64 - 5*e2*e2*e2/256 )*LatRad - (3*e2/8 + 3*e2*e2/32 + 45*e2*e2*e2/1024)*sin(2*LatRad) + (15*e2*e2/256 + 45*e2*e2*e2/1024)*sin(4*LatRad) - (35*e2*e2*e2/3072)*sin(6*LatRad) ); (*dx) = k0*N*(A+(1-T+C)*A*A*A/6 + (5-18*T+T*T+72*C-58*e2P)*A*A*A*A*A/120); (*dy) = k0*(M+N*tan(LatRad)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24 + (61-58*T+T*T+600*C-330*e2P)*A*A*A*A*A*A/720)); } void EUTM_to_LL_orig (double dx, double dy, double *Lat, double *Lon, double refLon) { double k0; double a; double e2; double e2P; double e1; double N1; double T1; double C1; double R1; double D; double M; double mu; double phi1; double phi1Rad; k0 = 0.9996; a = EquatorialRadius; e2 = EccentricitySquared; e1 = (1-sqrt(1-e2))/(1+sqrt(1-e2)); e2P = (e2)/(1-e2); M = dy / k0; mu = M/(a*(1-e2/4-3*e2*e2/64-5*e2*e2*e2/256)); phi1Rad = mu + (3*e1/2-27*e1*e1*e1/32)*sin(2*mu) + (21*e1*e1/16-55*e1*e1*e1*e1/32)*sin(4*mu) + (151*e1*e1*e1/96)*sin(6*mu); phi1 = phi1Rad*rad2deg; N1 = a/sqrt(1-e2*sin(phi1Rad)*sin(phi1Rad)); T1 = tan(phi1Rad)*tan(phi1Rad); C1 = e2P*cos(phi1Rad)*cos(phi1Rad); R1 = a*(1-e2)/pow(1-e2*sin(phi1Rad)*sin(phi1Rad), 1.5); D = dx/(N1*k0); (*Lat) = phi1Rad - (N1*tan(phi1Rad)/R1)*(D*D/2-(5+3*T1+10*C1-4*C1*C1-9*e2P)*D*D*D*D/24 +(61+90*T1+298*C1+45*T1*T1-252*e2P-3*C1*C1)*D*D*D*D*D*D/720); (*Lat) = (*Lat) * rad2deg; (*Lon) = (D-(1+2*T1+C1)*D*D*D/6+(5-2*C1+28*T1-3*C1*C1+8*e2P+24*T1*T1) *D*D*D*D*D/120)/cos(phi1Rad); (*Lon) = refLon + (*Lon)*rad2deg; } ---------------------------------------------------------------------- /* Checking: */ double xOld, yOld; LL_to_EUTM_orig(Lat, Lon, &xOld, &yOld, *refLon); affirm(fabs(x - xOld) < 0.001, "bug in x formula"); affirm(fabs(y - yOld) < 0.001, "bug in y formula"); ---------------------------------------------------------------------- St->tamS = 2; St->S[0] = v1; St->S[1] = v2; /* Constrói o conjunto admissivel {St->Livre} inicial */ { int t3; for(t3 = 0; t3Livre[t3] = TRUE; } } St->Livre[St->S[0]] = FALSE; St->Livre[St->S[1]] = FALSE; /* bloqueia os vértices adjacentes a St->S[0] e St->S[1] */ for(t3 = St->G.fnb[St->S[0]];t3G.fnb[St->S[0]+1];t3++) St->Livre[St->G.nbr[t3]] = FALSE; for(t3 = St->G.fnb[St->S[1]];t3G.fnb[St->S[1]+1];t3++) St->Livre[St->G.nbr[t3]] = FALSE; /* inicializa o vetor admissiveis contendo os graus dinâmicos dos vértices. Este vetor terá como elementos estruturas do tipo Tgrau. Talvez pudesse ser implementado como uma estrutura do tipo heap mas, por simplicidade do código, optou-se por um vetor simples. */ tamAdmissiveis = 0; for(t3 = 0;t3Livre[t3]){ admissiveis[tamAdmissiveis].vertice = t3; i = 0; for(t2 = St->G.fnb[t3];t2G.fnb[t3+1];t2++) if(St->Livre[St->G.nbr[t2]]) i++; admissiveis[tamAdmissiveis].grau = i; tamAdmissiveis++; } } ---------------------------------------------------------------------- void grasp_mis_greedy_build(grasp_mis_state_t *St, int v1, int v2) { int nV = St->G.nV; if (St->verbose) { cpk_trace(); } /* Pega vetores de trabalho: */ int *S = St->S; /* Vértices selecionados. */ bool_t *pertenceS = St->pertenceS; /* Indicadores de pertinência em {S}. */ int *L = St->L; /* Vértices admissíveis. */ bool_t *pertenceL = St->pertenceL; /* Indicadores de pertinência em {L}. */ int *grauL = St->grauL; /* Graus no subgrafo {G[L]}. */ /* No decorrer da malha, a solução corrente (ainda não maximal) é {S[0..tamS-1]}. Os vértices admissíveis para acrescentar a {S} -- isto é, que não estão em {S} e não são vizinhos a ninguém de {S} -- são {L[0..tamL]}. O booleano {pertenceL[v]} é TRUE se e somente se {v} está em {L}. Para todo vértice admissível {v}, {grauL[v]} é o grau de {v} no grafo induzido {G[L]}; ou seja, é o número de vizinhos de {v} que também são admissíveis. No começo de cada iteração, a lista {L} está ordenada por {grauL} crescente. Quando um vértice é retirado de {L} e colocado em {S}. */ /* Inicializa a solução corrente {S} com o conjunto vazio: */ int tamS = 0; double pesoS = 0; { int v; for (v = 0; v < nV; v++) { pertenceS[v] = FALSE; pertenceL[v] = TRUE; grauL[v] = St->G.deg[v]; L[v] = v; mudou[v] = FALSE; } } int tamL = nV; int tamMud = 0; auto void elimina(int v); /* Se {v} ainda é admissível, elimina-o do conjunto de admissíveis. Isso implica decrementar {grauL[w]} para todo vizinho {w} de {v}. Os vértices {w} que mudaram são marcados com {mudou[w] = TRUE} e acrescentados ao conjunto {Mud[0..tamMud-1]}. */ /* Malha do algoritmo semi-guloso: */ while (pqueue_count(&Q) > 0) { /* Neste ponto a fila deve estar atualizada: */ assert(tamMud == 0); /* Escolhe o próximo candidato: */ { int u = pqueue_smallest(&Q); assert(pertenceL[u]); elimina(u); /* Add {u} to {J}: */ J.el[nJ] = u; nJ++; /* Kill all neighbors of {u}: */ int fu = G->fnb[u], du = G->deg[u]; int iv; for (iv = 0; iv < du; iv++) { int v = G->nbr[fu + iv]; elimina(v); } } /* Recompute values of cands that had their {(deg - grauL)} changed: */ while (tamMud > 0) { tamMud--; int v = tofix[tamMud]; if (pertenceL[v]) { int Kv = G->deg[v] - (deg - grauL)[v]; /* Num of neighbors of {u} that are cands. */ assert(Kv >= 0); double Rv = N->el[v]/((double)1+Kv); /* New goodness of {v}. */ pqueue_set_value(&Q, v, -Rv); } mudou[v] = FALSE; } } int_vec_trim(&J, nJ); /* Return set: */ return J; void elimina(int v) { if (pertenceL[v]) { /* Kill candidate {v} */ pertenceL[v] = FALSE; pqueue_delete(&Q, v); /* Increment {(deg - grauL)[w]} for every neighbor {w} of {v}: */ int fv = G->fnb[v], dv = G->deg[v]; int iw; for (iw = 0; iw < dv; iw++) { int w = G->nbr[fv + iw]; (deg - grauL)[w]++; /* Rememeber that we must recompute the goodness of {w}: */ if (! mudou[w]) { tofix[tamMud] = w; tamMud++; mudou[w]=TRUE; } } } } void acrescentaS(int v); /* Acrescenta o vértice {v} ao conjunto {S}, atualizando {pertenceS,pertenceL,grauL}.*/ /* Inicializa a solução com {v1,v2}: */ acrescentaS(v1); acrescentaS(v2); /* enquanto o conjunto de L nao for nulo, aumenta o conjunto independente corrente St->S com o vértice de menor grau (dinâmico) presente neste conjunto */ while (tamL > 0) { /* ordena a lista de L em ordem de {grauL} crescente: */ qsort(L, tamL, sizeof(int), &ComparaGraus); if (St->verbose) { cpk_print_vertex_set(stderr, "vértices admissíveis", tamL, L, 25); } /* Escolhe um dos {JANELA} primeiros vértices de {L}: */ int j = abrandom(0, min(tamL, JANELA)-1); /* Retira-o da lista {L} /* adiciona o j-esimo elemento de L ao IS corrente */ u = L[j].vertice; St->S[St->tamS] = u; St->tamS++; /* L deve ser atualizado. u e seus vizinhos deixarao de ser L. Os demais vértices terao seus graus mudoudos se estiverem a distância 2 de u e forem livres. O vetor mudou indicara quais vértices precisam ter seu grau mudoudo. */ for(t2 = 0;t2pertenceL[u] = FALSE; for(t2 = St->G.fnb[u];t2G.fnb[u+1];t2++) { i = St->G.nbr[t2]; if (St->pertenceL[i]) { St->pertenceL[i] = FALSE; for(t3 = St->G.fnb[i];t3G.fnb[i+1];t3++){ j = St->G.nbr[t3]; if (St->pertenceL[j]) mudou[j] = TRUE; } } } t2 = 0; while (t2pertenceL[i]) { /* mudou grau dos vértices livres se necessario*/ if (mudou[i]) { j = 0; for(t3 = St->G.fnb[i];t3G.fnb[i+1];t3++) if (St->pertenceL[St->G.nbr[t3]]) j++; L[t2].grau = j; } t2++; /* avanca */ } else{ /* quando vértice nao for mais livre, ele sera removido de L (substituído pelo ultimo elemento) */ L[t2].vertice = L[tamL-1].vertice; L[t2].grau = L[tamL-1].grau; /* diminui o conjunto L mas ainda precisa testar a posição t2 e, por isso, nao avanca. */ tamL--; } } /* while */ /* === TESTE */ /* TSTimprimeS(); */ /* === FIM TESTE */ } /* repete ate chegar a um IS corrente (St->S) que seja maximal */ } ---------------------------------------------------------------------- { int u; for (u = 0; u < nV; u++) { degQ[u] = G->deg[u]; needfix[u] = FALSE; int Ku = [u]; /* Num of neighbors of {u} that are cands. */ double Ru = W->el[u]/((double)1+Ku); /* Goodness of {u}. */ pqueue_insert(&Q, u, -Ru); } } auto void kill_cand(int v); /* If {v} is still a candidate, removes it from the candidate set. Then, for each neighbor {w} od {f}, increments {NK[w]}; if {needfix[w]} was false, sets it TRUE and saves {w} in {tofix}. */ /* Greedy loop: */ while (pqueue_count(&Q) > 0) { /* At this point, no cand should need fixing: */ assert(ntofix == 0); /* Pop best candidate: */ { int u = pqueue_smallest(&Q); assert(ok[u]); kill_cand(u); /* Add {u} to {S}: */ S.el[nS] = u; nS++; /* Kill all neighbors of {u}: */ int fu = G->fnb[u], du = G->deg[u]; int iv; for (iv = 0; iv < du; iv++) { int v = G->nbr[fu + iv]; kill_cand(v); } } /* Recompute values of cands that had their {NK} changed: */ while (ntofix > 0) { ntofix--; int v = tofix[ntofix]; if (ok[v]) { int Kv = G->deg[v] - NK[v]; /* Num of neighbors of {u} that are cands. */ assert(Kv >= 0); double Rv = W->el[v]/((double)1+Kv); /* New goodness of {v}. */ pqueue_set_value(&Q, v, -Rv); } needfix[v] = FALSE; } } void kill_cand(int v) { if (ok[v]) { /* Kill candidate {v} */ ok[v] = FALSE; pqueue_delete(&Q, v); /* Increment {NK[w]} for every neighbor {w} of {v}: */ int fv = G->fnb[v], dv = G->deg[v]; int iw; for (iw = 0; iw < dv; iw++) { int w = G->nbr[fv + iw]; NK[w]++; /* Rememeber that we must recompute the goodness of {w}: */ if (! needfix[w]) { tofix[ntofix] = w; ntofix++; needfix[w]=TRUE; } } } } ---------------------------------------------------------------------- int_vec_t cpk_greedy_find_indep_set_A ( cpk_graph_t *G, /* The incompatibilty graph. */ r2_vec_t *V, /* Vertex voordinates. */ bool_t verbose /* TRUE prints debugging diagnostics. */ ); /* Uses a greedy heuristic based on vertex positions. Each vertex of {V} is examined, in order of increasing distance from from an arbitrary point, and is added to {S} if it is compatible with all vertices already in {S}. */ int_vec_t cpk_greedy_find_indep_set_A ( cpk_graph_t *G, /* The incompatibilty graph. */ r2_vec_t *V, /* Vertex voordinates. */ double_vec_t *W, /* Weight of each vertex. */ unsigned seed, /* Seed for random choices. */ bool_t verbose /* TRUE prints debugging diagnostics. */ ) { int nV = G->nV; assert(V->nel == nV); assert(W->nel == nV); /* Pick a reference point: */ r2_t ctr = cpk_greedy_pick_origin(V); /* Compute vertex distances {dst[0..nV-1]} from {ctr}: */ double dst[nV]; { int u; for (u = 0; u < nV; u++) { dst[u] = r2_dist(&(V->el[u]), &ctr); } } /* Sort the graph vertices by distance from {ctr}: */ int ix[nV]; { int u; for (u = 0; u < nV; u++) { ix[u] = u; } } auto int dcmp(int u, int v); /* Compares {dst[u]} with {dst[v]}, returns -1, 0, or +1. */ int dcmp(int u, int v) { if (dst[u] < dst[v]) { return -1; } else if (dst[u] > dst[v] ) { return +1; } else { return 0; } } isrt_heapsort(ix, nV, dcmp, +1); /* Greedy packing: */ int nS = 0; int_vec_t S = int_vec_new(nV); /* Selected vertices are {S[0..nS-1]}. */ bool_t ok[nV]; /* {ok[u]} is true iff vertex {u} is compatible with {S[0..nS-1]}. */ { int u; for (u = 0; u < nV; u++) { ok[u] = TRUE; } } { int i; for (i = 0; i < nV; i++) { int u = ix[i]; { if (ok[u]) { /* Add {u} to {S}: */ S.el[nS] = u; nS++; /* Kill all neighbors of {u}: */ { int fu = G->fnb[u], du = G->deg[u]; int j; for (j = 0; j < du; j++) { int v = G->nbr[fu + j]; ok[v] = FALSE; } } } } } } int_vec_trim(&S, nS); /* Return set: */ return S; } ---------------------------------------------------------------------- r2_t cpk_greedy_pick_origin(r2_vec_t *V); /* Given the candidates {V} for auction points, chooses the reference point {opk} for the geometric greedy packing heuristic. Candidates will be ordered by distance from {opk}. */ r2_t cpk_greedy_pick_origin(r2_vec_t *V) { int nV = V->nel; /* Compute vertex barycenter {b}: */ r2_t b; { double xb = 0, yb = 0; int i; for (i = 0; i < nV; i++) { r2_t *Vi = &(V->el[i]); xb += X(*Vi); yb += Y(*Vi); } b = (r2_t){{ xb/nV, yb/nV }}; } /* Compute radius {dMax}: */ double dMax = 0; { int i; for (i = 0; i < nV; i++) { r2_t *Vi = &(V->el[i]); double di = r2_dist(Vi, &b); if (di > dMax) { dMax = di; } } } /* Pick an arbitrary direction {dir}: */ r2_t dir = (r2_t) {{ 3, 4 }}; /* Pick a point far away in direction {dir}: */ r2_t opk; r2_mix(1.0, &b, -50*dMax/r2_norm(&dir), &dir, &opk); return opk; } ---------------------------------------------------------------------- int_vec_t cpk_greedy_find_indep_set_Z ( cpk_graph_t *G, /* The incompatibilty graph. */ r2_vec_t *V, /* Vertex voordinates. */ double_vec_t *W, /* Weight of each vertex. */ bool_t verbose /* TRUE prints debugging diagnostics. */ ); /* Uses the canonical greedy heuristic. Each vertex of {V} is examined, in order of decreasing weight {W}, and is added to {S} if it is compatible with all vertices already in {S}. */ ---------------------------------------------------------------------- /* calcula peso da solução St->S corrente */ St->pesoS = 0; for(t2 = 0;t2tamS;t2++) St->pesoS = St->pesoS+St->peso[St->S[t2]]; ---------------------------------------------------------------------- /* INFINITE HEXAGONAL GRIDS */ typedef struct cpk_grid_spec_t { r2_t o; /* Grid's reference point. */ r2_t dir; /* Grid's main period vector. */ } cpk_grid_spec_t; /* A {cpk_grid_spec_t} describes an infinite regular grid that includes the point {o} and has six fundamental period vectors of equal length, spaced 60 degrees apart, which include the vector {dir}. */ ---------------------------------------------------------------------- grasp_mis_state_t *grasp_mis_state_init ( cpk_graph_t *G, r2_vec_t *V, double_vec_t *W, double cutoff, bool_t verbose, bool_t validate ) { int nV = G->nV; /* Cria o registro de estado: */ grasp_mis_state_t *St; St = (grasp_mis_state_t *)notnull(malloc(sizeof(grasp_mis_state_t)), "no mem"); /* Armazena os parâmetros no registro de estado: */ St->G = *G; St->V = V->el; St->cutoff = cutoff; St->verbose = verbose; St->validate = validate; /* Pré-aloca vetores de trabalho: */ St->PertenceS = (bool_t *)notnull(malloc(nV*sizeof(bool_t)), "no mem"); St->Si = (int *)notnull(malloc(nV*sizeof(int)), "no mem"); St->Sij = (int *)notnull(malloc(nV*sizeof(int)), "no mem"); St->grauAux = (int *)notnull(malloc(nV*sizeof(int)), "no mem"); St->Sred = (int *)notnull(malloc(nV*sizeof(int)), "no mem"); St->Pares = (int *)notnull(malloc(2*nV*sizeof(int)), "no mem"); St->pertenceL = (bool_t *)notnull(malloc(nV*sizeof(bool_t)), "no mem"); /* inicializa o gerador de numeros pseudo-aleatórios */ /* srand(time(NULL)); */ } void grasp_mis_state_free(grasp_mis_state_t *St) { /* libera memoria local */ free(St->recalcula); /* liberação da memoria usada pelos vetores das buscas locais */ free(St->Si); free(St->Sij); free(St->PertenceS); free(St->grauAux); free(St->Sred); free(St->Pares); /* liberação da memoria usada aqui e em algumas buscas */ free(St->pertenceL); /* retorna o conjunto J */ } ---------------------------------------------------------------------- typedef struct grasp_mis_state_t { /* 1. variaveis de descricao da instancia: */ cpk_graph_t G; /* Grafo de incompatibilidades. */ r2_t *V; /* Vetor com as coordenadas dos vértices. */ double *peso; /* Vetor com os pesos dos vértices. */ /* 2. Parâmetros para a otimização: */ double cutoff; /* A busca local só é efetuada quando a fase de construção conseguiu um peso maior ou igual a {wtbest - cutoff} onde {wtbest} é o peso da melhor solução encontrada até o momento. */ bool_t verbose; /* Quando TRUE, causa a impressão de mensagens de diagnóstico em {stderr} no decorrer da otimização. */ bool_t validate; /* Quand TRUE, efetua vários testes (demorados!) de auto-consistência no decorrer da otimização. */ /* Áreas de trabalho de várias rotinas: */ bool_t *pertenceL; /* Vetor de tamanho {G->nV}, usado para indicar se um vértice eh independente da solucão corrente. */ bool_t *PertenceS; /* indica se um vértice esta na solucão corrente S ou não */ /* Conjunto dos vizinhos de um vértice {i}: */ int tamSi; /* Tamanho do conjunto (número de vértices). */ int *Si; /* Lista dos vértices. */ /* Conjunto dos vizinhos de dois vértices {i,j}: */ int tamSij; /* Tamanho do conjunto (número de vértices). */ int *Sij; /* Lista dos vértices. */ /* Conjunto independente dentro de un conjunto restrito de vértices: */ int tamSred; /* Tamanho do conjunto (número de vértices). */ int *Sred; /* Lista dos vértices. */ int *grauAux; /* vetor auxiliar para armazenar o grau dinamico dos vértices. O grau dinamico eh o grau original do vértice, menos o numero dos seus adjacentes que sao vizinhos a alguem que eh vizinho da solucão corrente S. */ int *Pares; /* este vetor sera usado para armazenar um ou dois vizinhos de vértices adjacentes a cada vértice i do grafo que estejam na solucão corrente (IS). */ } grasp_mis_state_t; ---------------------------------------------------------------------- void cpk_greedy_add_adm(cpk_greedy_state_t *SQ, int u); /* Adds a vertex {u} to the set {Q}. Automatically deletes {u} and its neighbors from {S}, if they are there. May change the order vertices in {SQ->S}. */ void cpk_greedy_add_adm(cpk_greedy_state_t *SQ, int u) { if (! pqueue_has(&(SQ->Q), u)) { /* Remove {u} from {S}, if it is there: */ cpk_greedy_simply_delete_indep(SQ, u); /* Get neighbors of {u}: */ int du = SQ->G->deg[u]; int *nbu = &(SQ->G->nbr[SQ->G->fnb[u]]); /* Delete from {S} any neighbor {w} of {u} that is in {S}: */ int i; for (i = 0; i < du; i++) { int w = nbu[i]; cpk_greedy_simply_delete_indep(SQ, w); } /* Add {u} to {Q}: */ cpk_greedy_simply_add_adm(SQ, u); } } ---------------------------------------------------------------------- /* In either case, the score will include small fractional biases to break ties towards auctions that are (1) as close as possible to the declared demand points, and (2) whose service areas include a larger fraction of the urbanized area. In any case, every proposed auction point will lie within the urbanized area {C->Urb}, and the corresponding auction center will be at the specified distances from the other obstacles: namely, {dMun} away from municipality borders {C->Mun}; {dNat} away from international borders {C->Nat}; {dMin+rAuc} away from other stations or predefined auction disks. */ ---------------------------------------------------------------------- /* The integer part of a {cpk_wcoeff_t} is interpreted as a priority {pr}: solutions are ranked by the terms with {pr = 0}, then ties are broken by terms with {pr = 1}, etc. The fractional part is interpreted as a multiplicative coefficient {alpha}, used to combine terms with same priority; except that a zero fraction is interpreted as {alpha = 1}. The {alpha}s of terms with same {pr} are normalized to unit sum. */ ---------------------------------------------------------------------- /* The weight of each vertex {v} is set to {W[i] = DNum(v) + Dist(v)*MDist + Area(v,rS)*MArea} for each vertex {v}. The terms are {DNum(v)}: the number of demand points {DI[k]} such that {dist(v,DI[k]) \leq rA}, where {rA = P->rAuc}. This term favors candidate auction centers that will be interesting to as many demand points as possible. {Dist(v)}: sum of {floor(100*(1-(dist(v,DI[k])/rA)^2))} for all of the above demand points {DI[k]}. Roughly this term favors candidate auction centers that are as close as possible to the demand points. {Area(v)}: approximately proportional to the area of the intersection of the urban polygon {C->R} and the circle with center {v} and radius {rS = P->dMin/2 + P->rAuc}. This term favors candidates whose service area is expected to cover a larger chunk of the urbanized area. The multipliers {MDist,MArea} are such that maximizing the total weight {W} means primarily maximizing the number of demand points that can compete for the auctions. In case of ties, auctions closer to the demand points are favored. The area criterion is hardly ever used. */ ---------------------------------------------------------------------- /* The weight is computed by the same formula, except that {DNum(v)} is replaced by {DNum(v)+1}. Maximizing the total weight then means essentially maximizing the number of auctions plus the number of declared demand points that can compete for them. */ ---------------------------------------------------------------------- void cpk_elim_no_demand_points(hxg_canvas_t *cvs); /* Eliminates all candidates to auction points that to not meet any declared demand. Assumes that a pixel has value {\geq 2} if and only if it meets some demand, that is, iff its auction disk contains at least one declared demand point. Thus the procedure merely sets all pixels less than 2 to 0, and subtracts 1 from the rest. */ ---------------------------------------------------------------------- #define MDist (0.0001) /* Multiplier for {Dist} term in candidate weight formula. */ #define MArea (0.000001) /* Multiplier for {Area} term in candidate weight formula. */ ---------------------------------------------------------------------- void cpk_elim_no_demand_points(hxg_canvas_t *cvs) { auto void op_elim_bg(int ix, int iy, hxg_pixel_t *b, hxg_pixel_t v); /* Subtracts {v} from {*b}, clipping at zero. */ hxg_paint_canvas(cvs, op_elim_bg, 1.0); /* Painting op implementations: */ void op_elim_bg(int ix, int iy, hxg_pixel_t *b, hxg_pixel_t v) { if ((*b) < 2) { (*b) = 0; } else { (*b) -= 1; } } } ---------------------------------------------------------------------- hxg_canvas_t hxg_canvas_new(r2_t org, r2_t step, i2_t size, hxg_pixel_t v) { int nx = X(size), ny = Y(size); demand(nx <= MAX_CANVAS_SIZE, "canvas too big"); demand(ny <= MAX_CANVAS_SIZE, "canvas too big"); int nt = nx*ny; hxg_pixel_t *pix = (hxg_pixel_t *)notnull(malloc(nt*sizeof(hxg_pixel_t)), "no mem"); int i; for (i = 0; i < nt; i++) { pix[i] = v; } return (hxg_canvas_t){ pix, org, step, size }; } ---------------------------------------------------------------------- double xKey = lo(B[0]) + 0.02*(hi(B[0]) - lo(B[0])); double yKey = hi(B[1]) - 0.02*(hi(B[1]) - lo(B[1])); double ySep = (ptSize/72.0)*((hi(B[0]) - lo(B[0]))/6.5); /* Guess. */ ---------------------------------------------------------------------- #define cpk_CUTOFF_WD 3.0 /* GRASP local search cutoff when DIs are specified (DI-driven mode). */ #define cpk_CUTOFF_ND 0.0 /* GRASP local search cutoff when running without DIs (area fill mode). */ ---------------------------------------------------------------------- double cutoff = (P->for_demand > 0 ? cpk_CUTOFF_WD : cpk_CUTOFF_ND ); ---------------------------------------------------------------------- double WR[nV]; /* Pesos perturbados para greedy semi-aleatório: */ double *WE;/* Pesos usados por {escore}; {WR} ou {W}, conforme fase. */ ---------------------------------------------------------------------- /* Calcula pesos perturbados {WR[v]} para greedy semi-aletório: */ { double R = cpk_grasp_REL_WT_NOISE; int v; for (v = 0; v < nV; v++) { WR[v] = (1.0 + R*drandom())*W[v]; } } ---------------------------------------------------------------------- WE = WR; ---------------------------------------------------------------------- WE = W; ---------------------------------------------------------------------- r2_vec_t cpk_read_points(char *inDir, char *tag, double *refLon, double *refY); /* Reads a list of points from the file "{inDir}/{tag}.txt". Ignores separator lines (whose first non-blank character is other than [0-9.+-]). */ /* Both procedures convert the {Lon,Lat} coordinates obtained from the file to {X,Y} coordinates (in meters) by EUTM projection. The parameter {refLon} should be the mean longitude (which defines the projection cylinder). The parameter {refY} should be a mean Y coordinate, that will be subtracted from the EUTM Y. If either of these parameters is infinite, they are both set to the appropriate averages of the input points. */ ---------------------------------------------------------------------- void cpk_write_solution ( char *outDir, char *solTag, char *demTag, r2_vec_t *V, double refLon, double refY, double_vec_t *W, int_vec_t *J ); /* Writes a proposed solution {J} (list of candidate vertices proposed as auction centers) to disk, as file "{outDir}/{solTag}/{demTag}.sol". The other parameters are: the EUTM X/Y coordinates {V} of the candidate auction points (grid points), the reference longitude {refLon} and reference EUTM Y {refY}, and the candidate weights {W}. */ void cpk_write_solution ( char *outDir, char *solTag, char *demTag, r2_vec_t *V, double refLon, double refY, double_vec_t *W, int_vec_t *J ) { /* Compose output file name: */ char *path = NULL; asprintf(&path, "%s/%s-%s.sol", outDir, solTag, demTag); FILE *wr = open_write(path); int i; double Lon = 777, Lat = 888; for (i = 0; i < J->nel; i++) { int u = J->el[i]; r2_t *up = &(V->el[u]); double xu = X(*up), yu = Y(*up); EUTM_to_LL(xu, yu + refY, &Lat, &Lon, refLon); fprintf(wr, LL_FMT " " LL_FMT, Lon, Lat); fprintf(wr, " "); fprintf(wr, XY_FMT " " XY_FMT, xu, yu); fprintf(wr, " "); fprintf(wr, "%12.8f", W->el[u]); fprintf(wr, "\n"); } fclose(wr); free(path); } ---------------------------------------------------------------------- double_vec_t cpk_copy_doubles(double_vec_t x); /* Makes a copy of a {double_vec_t}. */ double_vec_t cpk_copy_doubles(double_vec_t x) { double_vec_t y = double_vec_new(x.nel); int i; for (i = 0; i < x.nel; i++) { y.el[i] = x.el[i]; } return y; } ---------------------------------------------------------------------- typedef double cpk_mis_grader_t(int u, int d); /* A client-provided routine that is used internally by the routines below to sort the vertices of {Q}. When this procedure is called, the vertex {u} is in the set {Q}, and {d} is the number neighbors of {u} which are in {Q}. The procedure may be called several times for the same {u}, namely whenever {u} is added to {Q} and whenever {d} changes. */ cpk_mis_grader_t *grader; /* Defines the score of an elem of {Q}, or NULL. */ SQ->grader = NULL; void cpk_mis_set_grader(cpk_mis_state_t *SQ, cpk_mis_grader_t *grader); /* Specifies the scoring formula to be used to sort the elements of {Q}. The set {Q} must be empty when this procedure is called. Time: O(1). */ void cpk_mis_set_grader(cpk_mis_state_t *SQ, cpk_mis_grader_t *grader) { demand(pqueue_count(&(SQ->Q)) == 0, "Q not empty"); SQ->grader = grader; } cpk_mis_set_grader(SQ, &escore); cpk_mis_set_grader(SQ, &ratio_score); auto double escore(int u, int d); /* Função usada para selecionar o próximo vértice nas fases de contrução e busca local. Definida como sendo o peso corrente {W[u]} dividido por {d+1}, o número de candidatos que seriam eliminados se {u} fosse selecionado. */ double escore(int u, int d) { /* Calcula razão peso/consumo do vértice {u}: */ return W[u]/(d + 1); } ---------------------------------------------------------------------- void cpk_lopt_delete_indep(cpk_mis_state_t *SQ, double *WS, int v, double Wv, bool_t verbose) { assert(SQ->locS[v] >= 0); /* Retire {v} do conjunto {S}, acrescentando os novos admissíveis a {Q}: */ cpk_mis_delete_indep(SQ, v); /* Atualize o peso de {S}: */ (*WS) -= Wv; } void cpk_lopt_add_indep(cpk_mis_state_t *SQ, double *WS, int v, double Wv, bool_t verbose) { assert(SQ->locS[v] < 0); /* Coloque {v} no conjunto {S}, eliminando de {Q} os que ficarem inadmissíveis: */ cpk_mis_add_indep(SQ, v); /* Atualize o peso de {S}: */ (*WS) += Wv; } ---------------------------------------------------------------------- bool_t cpk_local_opt_multi_pass ( cpk_mis_state_t *SQ, /* Solução corrente {S} e estruturas auxiliares. */ int nRem, /* Número de vértice a trocar em cada tentativa. */ bool_t verbose, /* TRUE para imprimir mensagens de diagnóstico. */ /* Areas de trabalho: */ bool_t piv[], int QPiv[], int T[], bool_t temp[] ); /* Faz uma passada da otimização local com remoção de vértices múltiplos. Enumera todos os subconjuntos {T} de {S} com exatamente {nRem} elementos, tais que existe pelo menos um vértice {x} de {G-S} cujos vizinhos em {S} são exatamente {T}. Para cada {T} encontrado, tenta retirar {T} de {S}, e recobir o buraco usando greedy guloso com lookahead. Retorna TRUE assim que uma dessas tentativas melhora a solução; se nenhuma tentativa der certo, retorna FALSE. Os parâmetros {piv}, {Qpiv}, {T} e {temp} são áreas de trabalho que devem ser alocadas pelo cliente, com pelo menos {nV} elementos. O vetor {temp} deve ser todo FALSE na entrada, e será todo FALSE na saída. */ bool_t cpk_local_opt_multi_pass ( cpk_mis_state_t *SQ, /* Solução corrente {S} e estruturas auxiliares. */ int nRem, /* Número de vértices a trocar em cada tentativa. */ bool_t verbose, /* TRUE para imprimir mensagens de diagnóstico. */ /* Areas de trabalho: */ bool_t piv[], int QPiv[], int T[], bool_t temp[] ) { cpk_graph_t *G = SQ->G; int nV = G->nV; if (verbose) { fprintf(stderr, "\n"); } if (verbose) { fprintf(stderr, " Procurando grupos de tamanho %d...\n", nRem); } /* O conjunto {S} deve ser maximal: */ assert(pqueue_count(&(SQ->Q)) == 0); /* Constrói a lista de candidatos para pivô {QPiv[0..nQPiv-1]}. */ /* Também inicializa indicador de {piv[0..nV-1]} dos candidatos elegíveis. */ int nQPiv = 0; { int v; for (v = 0; v < nV; v++) { /* Pega vértices de {Q-S} com {2..maxRem} vizinhos em {S}: */ piv[v] = (SQ->locS[v] < 0) & (SQ->degS[v] == nRem); if (piv[v]) { QPiv[nQPiv] = v; nQPiv++; } } } if (nQPiv == 0) { if (verbose) { fprintf(stderr, " Não encontrou nenhum grupo\n"); } return FALSE; } /* Ordena os candidatos a pivôs, na ordem em que serão considerados: */ { auto int compara_grau_S_cre(const void *a, const void *b); /* Para ordenar em ordem CRESCENTE de {|T|}. */ auto int compara_peso_dec(const void *a, const void *b); /* Para ordenar em ordem DECRESCENTE de {W}. */ qsort(QPiv, nQPiv, sizeof(int), &compara_grau_S_cre); /* qsort(QPiv, nQPiv, sizeof(int), &compara_peso_dec); */ int compara_grau_S_cre(const void *a, const void *b) { int u = *((int*)a), dSu = SQ->degS[u]; int v = *((int*)b), dSv = SQ->degS[v]; return intcmp(dSu, dSv); } int compara_peso_dec(const void *a, const void *b) { int u = *((int*)a); int v = *((int*)b); return -dblcmp(SQ->W[u],SQ->W[v]); } } /* Malha sobre os pivôs: */ int iQPiv = 0; while (iQPiv < nQPiv) { /* Pega um candidato a pivô {x}: */ int x = QPiv[iQPiv]; if (piv[x]) { /* O vértice {x} ainda não entrou em nenhum {Q}. */ if (cpk_local_opt_at_pivot(SQ, x, verbose, piv, T, temp)) { if (verbose) { TRACE_SOL(" Melhorou solução corrente", SQ->nS, SQ->WS); } return TRUE; } } /* Próximo candidato: */ iQPiv++; } /* Não conseguimos melhorar: */ return FALSE; } ----------------------------------------------------------------------