rm(list=ls())
library(igraph)
##
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
library(shp2graph)
library(geosphere)
library(rgdal)
## Loading required package: sp
## rgdal: version: 1.3-9, (SVN revision 794)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
## Linking to sp version: 1.3-1
library(leaflet)
library(scales)
library(lwgeom)
## Linking to liblwgeom 2.5.0dev r16016, GEOS 3.6.1, PROJ 4.9.3
library(shapefiles)
## Loading required package: foreign
##
## Attaching package: 'shapefiles'
## The following objects are masked from 'package:foreign':
##
## read.dbf, write.dbf
library(sp)
library(rgeos)
## rgeos version: 0.4-2, (SVN revision 581)
## GEOS runtime version: 3.6.1-CAPI-1.10.1
## Linking to sp version: 1.3-1
## Polygon checking: TRUE
library(osmdata)
## Data (c) OpenStreetMap contributors, ODbL 1.0. http://www.openstreetmap.org/copyright
library(sf)
## Linking to GEOS 3.6.1, GDAL 2.1.3, PROJ 4.9.3
library(ggmap)
## Loading required package: ggplot2
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library(maptools)
## Checking rgeos availability: TRUE
library(Matrix)
library(tripack)
##
## Attaching package: 'tripack'
## The following objects are masked from 'package:igraph':
##
## convex.hull, triangles
library(RColorBrewer)
library(units)
## udunits system database from /Library/Frameworks/R.framework/Versions/3.5/Resources/library/units/share/udunits
library(classInt)
dd <- data.frame(Id=c(1,1,1,2,2,2,3,3,4,4),
X=c(3,5,8,5,7,8,3,5,5,6),
Y=c(9,8,3,4,7,4,8,9,7,8))
dd
par(mfrow=c(1,2))
plot(dd$X,dd$Y,col="white",xlab="",ylab="")
for(i in unique(dd[,1])) lines(dd$X[dd$Id==i],dd$Y[dd$Id==i])
points(dd$X,dd$Y,pch=19,col="red")
library(RColorBrewer)
darkcols = brewer.pal(8, "Dark2")
plot(dd$X,dd$Y,xlab="",ylab="")
for(i in unique(dd[,"Id"])) lines(dd$X[dd$Id==i],dd$Y[dd$Id==i],col=darkcols[i],lwd=2)
extend_shp <- function(dd){
newdd1=NULL
newdd2=NULL
cpte=0
W=names(which(table(dd$Id)>1))
dd=dd[dd$Id %in% as.numeric(W),]
for(no in unique(dd$Id)){
cpte=cpte+1
G=NULL
DD=dd[dd$Id==no,2:3]
p1=c(min(DD[,1]),min(DD[,2]),max(DD[,1]),max(DD[,2]))
P1=cbind(c(p1[1],p1[1],p1[3],p1[3],p1[1]),c(p1[2],p1[4],p1[4],p1[2],p1[2]))
box_i=function(i){
b=dd[dd$Id==i,2:3]
p2=c(min(b[,1]),min(b[,2]),max(b[,1]),max(b[,2]))
P2=cbind(c(p2[1],p2[1],p2[3],p2[3]),c(p2[2],p2[4],p2[4],p2[2]))
ct=function(j) point.in.polygon(P2[j,1],P2[j,2],P1[,1],P2[,2])
sum(Vectorize(ct)(1:4))
}
vu = unique(dd$Id)[-which(unique(dd$Id)==no)]
for(t in 2:nrow(DD)){
P=DD[t-(1:0),]
L1=SpatialLines(list(Lines(list(Line(P)),"Line1")))
for(i in vu){
L2=SpatialLines(list(Lines(list(Line(dd[dd$Id==i,2:3])),"Line2")))
Inters=gIntersection(L1,L2)
if(length(Inters)>0){
if(class(Inters)=="SpatialPoints") G=rbind(G,c(Inters@coords,t))
if(class(Inters)=="SpatialLines") G=rbind(G,c(Inters@lines[[1]]@Lines[[1]]@coords,t))
}}
P2=DD[1,]
for(t in 2:nrow(DD)){
if(sum(G[,3]==t)==0) P2=rbind(P2,DD[t,])
if(sum(G[,3]==t)==1) P2=rbind(P2,G[G[,3]==t,1:2],DD[t,])
if(sum(G[,3]==t)>1){
sG=G[G[,3]==t,1:2]
if(DD[t-1,1]<DD[t,1]) sG=sG[order(sG[,1]),]
if(DD[t-1,1]>DD[t,1]) sG=sG[-order(sG[,1]),]
colnames(sG)=c("X","Y")
P2=rbind(P2,sG,DD[t,])
}
}
dd2=data.frame(Id=no,X=P2$X,Y=P2$Y)
dd2=dd2[!duplicated(dd2),]
newdd1=rbind(newdd1,dd2)
dd2=data.frame(Id=no,X1=P2$X[1:(nrow(dd2)-1)],Y1=P2$Y[1:(nrow(dd2)-1)],
X2=P2$X[2:(nrow(dd2))],Y2=P2$Y[2:(nrow(dd2))])
dd2=dd2[!duplicated(dd2),]
newdd2=rbind(newdd2,dd2)
}}
LIST_KNOTS = newdd1[!duplicated(newdd1[,2:3]),2:3]
LIST_KNOTS=data.frame(LIST_KNOTS,NO=1:nrow(LIST_KNOTS))
newdd2$Tr = 1:nrow(newdd2)
depart = data.frame(Tr=newdd2[,"Tr"],Id=newdd2[,"Id"],X=newdd2[,"X1"],Y=newdd2[,"Y1"])
arrivee = data.frame(Tr=newdd2[,"Tr"],Id=newdd2[,"Id"],X=newdd2[,"X2"],Y=newdd2[,"Y2"])
depart = merge(depart,LIST_KNOTS,all.x=TRUE)
arrivee = merge(arrivee,LIST_KNOTS,all.x=TRUE)
names(depart)[names(depart)=="NO"]="knot1"
names(arrivee)[names(arrivee)=="NO"]="knot2"
base=merge(newdd2,depart[,c("Tr","Id","knot1")],all.x=TRUE,by="Tr")
base=merge(base,arrivee[,c("Tr","Id","knot2")],all.x=TRUE,by="Tr")
base$distance = distHaversine( base[,c("X1","Y1")] , base[,c("X2","Y2")] , r =6378.137)
base2=base[,c("knot1","knot2","distance","Id")]
nb=base2[,c("knot2","knot1","distance","Id")];names(nb)=c("knot1","knot2","distance","Id")
base2=rbind(base2,nb)
base2=base2[!duplicated(base2),]
return(list(shp=newdd1,
knots=LIST_KNOTS,
connected=base2))}
dd2 <- extend_shp(dd)
par(mfrow=c(1,2))
plot(dd$X,dd$Y,col="white",xlab="",ylab="")
for(i in unique(dd[,1])) lines(dd$X[dd$Id==i],dd$Y[dd$Id==i])
points(dd2[["knots"]]$X,dd2[["knots"]]$Y,pch=19,col="blue")
points(dd$X,dd$Y,pch=19,col="red")
plot(dd2[["shp"]]$X,dd2[["shp"]]$Y,col="white",xlab="",ylab="",axes=FALSE,xlim=c(2.6,8.4))
text(dd2[["knots"]]$X,dd2[["knots"]]$Y,dd2[["knots"]]$NO)
points(dd2[["knots"]]$X,dd2[["knots"]]$Y,cex=3)
library(RColorBrewer)
darkcols = brewer.pal(8, "Dark2")
plot(dd2[["shp"]]$X,dd2[["shp"]]$Y,col="white",xlab="",ylab="",axes=FALSE)
for(i in 1:(nrow(dd2$connected)/2)) segments(dd2[["knots"]]$X[dd2[["connected"]]$knot1[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot1[i]],
dd2[["knots"]]$X[dd2[["connected"]]$knot2[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot2[i]],
col=darkcols[dd2[["connected"]]$Id[i]])
points(dd2[["knots"]]$X,dd2[["knots"]]$Y,pch=19,col="black",cex=3)
text(dd2[["knots"]]$X,dd2[["knots"]]$Y,dd2[["knots"]]$NO,col="white")
agents=data.frame(X=c(5.5,3.5,6.5,6.5,4.5),
Y=c(4,9,8,5.5,9),
Type=as.factor(c("Demand","Demand","Oferta","Oferta","Oferta")))
forme=c(15,19,0,1)
plot(dd2[["shp"]]$X,dd2[["shp"]]$Y,col="white",xlab="",ylab="",axes=FALSE)
for(i in 1:(nrow(dd2$connected)/2)) segments(dd2[["knots"]]$X[dd2[["connected"]]$knot1[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot1[i]],
dd2[["knots"]]$X[dd2[["connected"]]$knot2[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot2[i]])
points(dd2[["knots"]]$X,dd2[["knots"]]$Y,pch=19,col="white",cex=3)
points(dd2[["knots"]]$X,dd2[["knots"]]$Y,pch=1,col="black",cex=3)
text(dd2[["knots"]]$X,dd2[["knots"]]$Y,dd2[["knots"]]$NO,col="black")
points(agents$X,agents$Y,pch=forme[agents$Type],col=darkcols[agents$Type],cex=3)
plot(dd2[["shp"]]$X,dd2[["shp"]]$Y,col="white",xlab="",ylab="",axes=FALSE)
for(i in which(agents$Type=="Demand")) segments(agents[i,"X"],agents[i,"Y"],agents[agents$Type=="Oferta","X"],agents[agents$Type=="Oferta","Y"])
points(agents$X,agents$Y,pch=forme[agents$Type],col=darkcols[agents$Type],cex=3)
close_supply <- function(coord,knots=agents[agents$Type=="Oferta",]){
knots[which.min(distHaversine( coord[,c("X","Y")] , knots[,c("X","Y")] , r =6378.137)),]
}
list_demand <- which(agents$Type=="Demand")
close_supply(agents[list_demand[1],])
plot(dd2[["shp"]]$X,dd2[["shp"]]$Y,col="white",xlab="",ylab="",axes=FALSE)
for(i in which(agents$Type=="Demand")) segments(agents[i,"X"],agents[i,"Y"],agents[agents$Type=="Oferta","X"],agents[agents$Type=="Oferta","Y"])
for(i in list_demand){
supply=unlist(close_supply(agents[i,]))
segments(agents[i,"X"],agents[i,"Y"],supply["X"],supply["Y"],col=darkcols[3],lwd=6)
}
points(agents$X,agents$Y,pch=forme[agents$Type],col=darkcols[agents$Type],cex=3.2)
text(agents$X,agents$Y,1:nrow(agents),pch=forme[agents$Type],col="white")
close_knot=function(coord,knots=dd2$knots){
knots[which.min(distHaversine( coord[,c("X","Y")] , knots[,c("X","Y")] , r =6378.137)),]
}
shift_agents=agents
shift_agents$Xk=shift_agents$X
shift_agents$Yk=shift_agents$Y
shift_agents$NO=0
for(i in 1:nrow(shift_agents)) shift_agents[i,c("Xk","Yk","NO")]=close_knot(agents[i,])
shift_agents
plot(dd2[["shp"]]$X,dd2[["shp"]]$Y,col="white",xlab="",ylab="",axes=FALSE)
for(i in 1:(nrow(dd2$connected)/2)) segments(dd2[["knots"]]$X[dd2[["connected"]]$knot1[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot1[i]],
dd2[["knots"]]$X[dd2[["connected"]]$knot2[i]],dd2[["knots"]]$Y[dd2[["connected"]]$knot2[i]])
points(dd2[["knots"]]$X,dd2[["knots"]]$Y,pch=19,col="white",cex=3)
points(dd2[["knots"]]$X,dd2[["knots"]]$Y,pch=1,col="black",cex=3)
points(agents$Xk,agents$Yk,pch=forme[agents$Type],col=darkcols[agents$Type],cex=3)
text(dd2[["knots"]]$X,dd2[["knots"]]$Y,dd2[["knots"]]$NO,col="black")
points(shift_agents$Xk,shift_agents$Yk,pch=forme[shift_agents$Type],col=darkcols[shift_agents$Type],cex=3.2)
points(shift_agents$X,shift_agents$Y,pch=forme[2+as.numeric(shift_agents$Type)],col=darkcols[as.numeric(shift_agents$Type)],cex=3)
plot(0:1,0:1,col="white",xlab="",ylab="",axes=FALSE)
text(.125,.5,"DEMANDA",srt=90,col=darkcols[1])
text(.875,.5,"Oferta",srt=270,col=darkcols[2])
i=which(agents$Type=="Demand")
y1=rev(seq(1/(2*length(i)),1-1/(2*length(i)),by=1/(2*length(i)))[seq(1,2*length(i),by=2)])
i=which(agents$Type=="Oferta")
y2=rev(seq(1/(2*length(i)),1-1/(2*length(i)),by=1/(2*length(i)))[seq(1,2*length(i),by=2)])
id=which(agents$Type=="Demand")
for(i in 1:sum(agents$Type=="Demand")){
d=distHaversine( agents[id[i],c("X","Y")] , agents[agents$Type=="Oferta",c("X","Y")] , r =6378.137)
segments(.25,y1[i],.75,y2,col=c("black",darkcols[3])[1+(d==min(d))],lwd=c(1,6)[1+(d==min(d))])
text(.5,(y2+y1[i])/2,round(d),pos=1,col=c("black",darkcols[3])[1+(d==min(d))])
}
i=which(agents$Type=="Demand")
points(rep(.25,length(i)),y1,
pch=forme[1],col=darkcols[1],cex=3.2)
text(rep(.25,length(i)),y1,i,col="white")
i=which(agents$Type=="Oferta")
points(rep(.75,length(i)),y2,
pch=forme[2],col=darkcols[2],cex=3.2)
text(rep(.75,length(i)),y2,i,col="white")
Demand = shift_agents[shift_agents$Type=="Demand",c("NO","X","Y")]
Supply = shift_agents[shift_agents$Type=="Oferta",c("NO","X","Y")]
Nodes = dd2$knots[,c("NO","X","Y")]
Arcs = dd2$connected
library(igraph)
iarcs <- make_graph(edges = as.vector(rbind(Arcs[,1],Arcs[,2])),directed=TRUE)
plot(iarcs,layout=as.matrix(dd2$knots[,c("X","Y")]),edge.arrow.size=.2)
iarcs = set_edge_attr(iarcs,"distance", value=Arcs[,"distance"])
cat("The graph is connex :", is.connected(iarcs,mode="strong"),"\n")
## The graph is connex : TRUE
if(!is.connected(iarcs,mode="strong")){
nodestokeep=which(clusters(iarcs,mode="strong")[[1]]==which(clusters(iarcs,mode="strong")[[2]]==max(clusters(iarcs,mode="strong")[[2]])))
nodestokeep
Arcs = Arcs[(Arcs$knot1 %in% nodestokeep)&(Arcs$knot2 %in% nodestokeep),]
Nodes = Nodes[Nodes$NO %in% nodestokeep, ]
iarcs = graph.edgelist(as.matrix(Arcs[,c("knot1","knot2")]), directed=TRUE)
iarcs = set_edge_attr(iarcs,"distance", value=Arcs[,"distance"])
cat("The graph is connex :", is.connected(iarcs,mode="strong"),"\n")
plot(iarcs,layout=as.matrix(dd2$knots[dd2$knots$NO%in%nodestokeep,c("X","Y")]),edge.arrow.size=.2)}
AP=all_shortest_paths(iarcs,
from=Demand[1,"NO"],
to=Supply[3,"NO"])
L=AP$res[[1]]
n=length(L)
V(iarcs)$color="yellow"
V(iarcs)$color[L[2:(n-1)]]="light blue"
V(iarcs)$color[L[c(1,n)]]="blue"
plot(iarcs,layout=as.matrix(dd2$knots[,c("X","Y")]),edge.arrow.size=.2)
dist_graph = function(i,j) dd2$connected[(dd2$connected$knot1==i)&(dd2$connected$knot2==j),"distance"][1]
D=function(dem,sup){
AP=all_shortest_paths(iarcs,
from=Demand[dem,"NO"],
to=Supply[sup,"NO"],weights=dd2$connected$distance)
AP_path = AP$res[[1]]
d=0
for(u in 1:(length(AP_path)-1)) d=d+dist_graph(AP_path[u],AP_path[u+1])
d}
DF=expand.grid(dem=1:nrow(Demand),sup=1:nrow(Supply))
M=matrix(Vectorize(function(i) D(DF[i,"dem"],DF[i,"sup"]))(1:(nrow(Demand)*nrow(Supply))),nrow(Demand),nrow(Supply))
M
## [,1] [,2] [,3]
## [1,] 542.9799 253.2018 772.7032
## [2,] 425.6273 519.5013 246.6023
idx=apply(M,1,which.min)
for(i in 1:length(idx)){
AP=all_shortest_paths(iarcs,
from=Demand[i,"NO"],
to=Supply[idx[i],"NO"])
L=AP$res[[1]]
n=length(L)
V(iarcs)$color="yellow"
V(iarcs)$color[L[2:(n-1)]]="light blue"
V(iarcs)$color[L[c(1,n)]]="blue"
plot(iarcs,layout=as.matrix(dd2$knots[,c("X","Y")]),edge.arrow.size=.2,
main=paste("Menor caminho, demanda ",Demand[i,"NO"],", oferta ",Supply[idx[i],"NO"],sep=""))
}
arcs = read.csv(file = "arcs-metro.csv",sep=";", header=FALSE)
dim(arcs)
## [1] 868 3
nodes = read.csv("nodes-metro.csv",sep=";", header=FALSE) # loads the data
dim(nodes)
## [1] 366 6
par(mfrow=c(1,1))
subarcs=arcs[(arcs$V1<303)&(arcs$V2<303),]
plot(nodes[1:302,c("V3","V4")],pch=19,xlab="",ylab="",axes=FALSE)
for(i in 1:nrow(arcs)) segments(nodes[subarcs$V1[i],"V3"],nodes[subarcs$V1[i],"V4"],
nodes[subarcs$V2[i],"V3"],nodes[subarcs$V2[i],"V4"])
library(OpenStreetMap)
library(devtools)
## Loading required package: usethis
#install_github("ropensci/osmdata")
library(osmdata)
library(sf)
library(units)
library(scales)
#map=openmap(c(lat=48.9,lon=2.26),c(lat=48.8,lon=2.45))
#map=openproj(map,projection="+init=epsg:4326")
#plot(map)
# see https://rcarto.github.io/caRtosm/index.html
q0 = opq(bbox = c(2.2247, 48.8188, 2.4611, 48.9019))
q1 = add_osm_feature(opq = q0, key = 'name', value = "Paris")
res1 = osmdata_sf(q1)
paris = st_geometry(res1$osm_multipolygons[1,])
q2 <- add_osm_feature(opq = q0, key = 'name', value = "La Seine")
res2 <- osmdata_sf(q2)
seine1 <- st_geometry(res2$osm_multilines)
q2b <- add_osm_feature(opq = q0, key = 'name',
value = "La Seine - Bras de la Monnaie")
res2b <- osmdata_sf(q2b)
seine2 <- st_geometry(res2b$osm_lines)
q3 <- add_osm_feature(opq = q0, key = 'leisure', value = "park")
res3 <- osmdata_sf(q3)
parc1 <- st_geometry(res3$osm_polygons)
parc2 <- st_geometry(res3$osm_multipolygons)
q4 <- add_osm_feature(opq = q0, key = 'landuse', value = "cemetery")
res4 <- osmdata_sf(q4)
parc3 <- st_geometry(res4$osm_polygons)
q5 <- add_osm_feature(opq = q0, key = 'admin_level', value = "10")
res5 <- osmdata_sf(q5)
quartier <- res5$osm_multipolygons
parc <- do.call(c, list(parc1, parc2, parc3))
parc <- st_union(x = st_buffer(parc,0), by_feature = F)
## Warning in st_buffer.sfc(parc, 0): st_buffer does not correctly buffer
## longitude/latitude data
## dist is assumed to be in decimal degrees (arc_degrees).
parc <- st_cast(parc, "POLYGON")
parc <- parc[st_area(parc)>=set_units(10000, "m^2")]
parc <- st_intersection(x = parc, y = paris)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
seine <- st_intersection(x = seine1, y = paris)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
seine <- c(st_cast(seine[1])[2:5], seine[2])
seine <-c(seine, seine2)
back_paris=function(a=1,background=FALSE){
plot(paris, col=alpha("#D9D0C9",a), border = NA, add=background)
plot(parc, col=alpha("#CDEBB2",a), border = NA, add=TRUE)
plot(seine, col=alpha("#AAD3DF",a), add=TRUE, lwd = 4)
plot(st_geometry(quartier), col = NA,lty = 2, lwd = 0.2, add=TRUE)
}
back_paris()
points(nodes[1:302,c("V3","V4")],pch=19)
for(i in 1:nrow(arcs)) segments(nodes[subarcs$V1[i],"V3"],nodes[subarcs$V1[i],"V4"],
nodes[subarcs$V2[i],"V3"],nodes[subarcs$V2[i],"V4"])
library(ggmap)
adresses <- c("47 boulevard de l'hopital Paris","20 rue leblanc Paris","2 rue ambroise pare paris","1 place parvis notre dame paris","149 rue de sevres paris","1 avenue claude vellafaux paris","184 faubourg saint antoine paris","46 rue henri huchard paris")
# hospital <- lapply(adresses,geocode)
hospital = list(c(lon=2.3610258,lat=48.8401821),
c(lon=2.2719068,lat=48.838736),
c(lon=2.3507367,lat=48.8824963),
c(lon=2.346142,lat=48.8539541),
c(lon=2.3125916,lat=48.8455934),
c(lon=2.3676277,lat=48.873488),
c(lon=2.3806986,lat=48.8500029),
c(lon=2.3300414,lat=48.8982232))
adresses <- c("Place Denfert-Rochereau, Paris, France")
# sick <- lapply(adresses,geocode)
sick = list(c(lat=48.8342694,lon=2.3300964))
back_paris()
points(nodes[1:302,c("V3","V4")],pch=19)
for(i in 1:nrow(arcs)) segments(nodes[subarcs$V1[i],"V3"],nodes[subarcs$V1[i],"V4"],
nodes[subarcs$V2[i],"V3"],nodes[subarcs$V2[i],"V4"],col="grey")
points(unlist(hospital)[names(unlist(hospital))=="lon"],unlist(hospital)[names(unlist(hospital))=="lat"],pch=forme[2],col=darkcols[2],cex=2)
points(unlist(sick)[names(unlist(sick))=="lon"],unlist(sick)[names(unlist(sick))=="lat"],
pch=forme[1],col=darkcols[1],cex=2)
close_knot=function(coord,knots=data.frame(X=nodes$V3,Y=nodes$V4,NO=1:nrow(nodes))){
knots[which.min(distHaversine( coord[c("lon","lat")] , knots[,c("X","Y")] , r =6378.137)),]
}
newhospital=data.frame(X=rep(NA,length(hospital)),Y=rep(NA,length(hospital)),NO=rep(NA,length(hospital)))
for(i in 1:length(hospital)) newhospital[i,]=as.numeric(close_knot(hospital[[i]]))
newsick=data.frame(X=rep(NA,length(sick)),Y=rep(NA,length(sick)),NO=rep(NA,length(sick)))
for(i in 1:length(sick)) newsick[i,]=as.numeric(close_knot(sick[[i]]))
back_paris()
points(nodes[1:302,c("V3","V4")],pch=19)
for(i in 1:nrow(arcs)) segments(nodes[subarcs$V1[i],"V3"],nodes[subarcs$V1[i],"V4"],
nodes[subarcs$V2[i],"V3"],nodes[subarcs$V2[i],"V4"],col="grey")
points(newhospital[,"X"],newhospital[,"Y"],pch=forme[2],col=darkcols[2],cex=2)
points(newsick[,"X"],newsick[,"Y"],pch=forme[1],col=darkcols[1],cex=2)
library(igraph)
g <- make_graph(edges = as.vector(rbind(arcs[,1],arcs[,2])),directed=FALSE)
sp <- shortest_paths(g, newsick[1,"NO"], newhospital[1,"NO"], weights= arcs[,3])
pth <- as.numeric(sp$vpath[[1]])
back_paris()
points(nodes[,c("V3","V4")],pch=19)
for(i in 1:nrow(arcs)) segments(nodes[arcs$V1[i],"V3"],nodes[arcs$V1[i],"V4"],
nodes[arcs$V2[i],"V3"],nodes[arcs$V2[i],"V4"],col="grey")
points(nodes[pth,c("V3","V4")],col=darkcols[3],pch=19)
for(i in 1:(length(pth)-1)) segments(nodes[pth[i],"V3"],nodes[pth[i],"V4"],
nodes[pth[i+1],"V3"],nodes[pth[i+1],"V4"],col=darkcols[3],lwd=5)
PCH=rep(1,nrow(newhospital))
PCH[1]=forme[2]
points(newhospital[,"X"],newhospital[,"Y"],pch=PCH,col=darkcols[2],cex=2,lwd=2)
points(newsick[,"X"],newsick[,"Y"],pch=forme[1],col=darkcols[1],cex=2)
library(igraph)
g <- make_graph(edges = as.vector(rbind(arcs[,1],arcs[,2])),directed=FALSE)
sp <- shortest_paths(g, newsick[1,"NO"], newhospital[3,"NO"], weights= arcs[,3])
pth <- as.numeric(sp$vpath[[1]])
back_paris()
points(nodes[,c("V3","V4")],pch=19)
for(i in 1:nrow(arcs)) segments(nodes[arcs$V1[i],"V3"],nodes[arcs$V1[i],"V4"],
nodes[arcs$V2[i],"V3"],nodes[arcs$V2[i],"V4"],col="grey")
points(nodes[pth,c("V3","V4")],col=darkcols[3],pch=19)
for(i in 1:(length(pth)-1)) segments(nodes[pth[i],"V3"],nodes[pth[i],"V4"],
nodes[pth[i+1],"V3"],nodes[pth[i+1],"V4"],col=darkcols[3],lwd=5)
PCH=rep(1,nrow(newhospital))
PCH[3]=forme[2]
points(newhospital[,"X"],newhospital[,"Y"],pch=PCH,col=darkcols[2],cex=2,lwd=2)
points(newsick[,"X"],newsick[,"Y"],pch=forme[1],col=darkcols[1],cex=2)
distances(g, newsick[1,"NO"], newhospital[1,"NO"], weights= arcs[,3])
## [,1]
## [1,] 2872.404
distances(g, newsick[1,"NO"], newhospital[3,"NO"], weights= arcs[,3])
## [,1]
## [1,] 6029.571
plot(0:1,0:1,col="white",xlab="",ylab="",axes=FALSE)
text(.125,.5,"DEMANDA",srt=90,col=darkcols[1])
text(.875,.5,"OFERTA",srt=270,col=darkcols[2])
agents=data.frame(X=c(newhospital$X,newsick$X),Y=c(newhospital$Y,newsick$Y),
Type=c(rep("Oferta",nrow(newhospital)),c(rep("Demanda",nrow(newsick)))))
i=which(agents$Type=="Demanda")
y1=rev(seq(1/(2*length(i)),1-1/(2*length(i)),by=1/(2*length(i)))[seq(1,2*length(i),by=2)])
i=which(agents$Type=="Oferta")
y2=rev(seq(1/(2*length(i)),1-1/(2*length(i)),by=1/(2*length(i)))[seq(1,2*length(i),by=2)])
id=which(agents$Type=="Demanda")
for(i in 1:sum(agents$Type=="Demanda")){
d=rep(NA,sum(agents$Type=="Oferta"))
for(j in 1:sum(agents$Type=="Oferta")){
d[j]=distances(g, newsick[i,"NO"], newhospital[j,"NO"], weights= arcs[,3])
}
segments(.25,y1[i],.75,y2,col=c("black",darkcols[3])[1+(d==min(d))],lwd=c(1,6)[1+(d==min(d))])
text(.5,(y2+y1[i])/2,round(d),pos=1,col=c("black",darkcols[3])[1+(d==min(d))])
}
i=which(agents$Type=="Demanda")
points(rep(.25,length(i)),y1,
pch=forme[1],col=darkcols[1],cex=3.2)
text(rep(.25,length(i)),y1,i,col="white")
i=which(agents$Type=="Oferta")
points(rep(.75,length(i)),y2,
pch=forme[2],col=darkcols[2],cex=3.2)
text(rep(.75,length(i)),y2,i,col="white")
which(d==min(d))
## [1] 5
library(igraph)
g <- make_graph(edges = as.vector(rbind(arcs[,1],arcs[,2])),directed=FALSE)
sp <- shortest_paths(g, newsick[1,"NO"], newhospital[which.min(d),"NO"], weights= arcs[,3])
pth <- as.numeric(sp$vpath[[1]])
back_paris()
points(nodes[,c("V3","V4")],pch=19)
for(i in 1:nrow(arcs)) segments(nodes[arcs$V1[i],"V3"],nodes[arcs$V1[i],"V4"],
nodes[arcs$V2[i],"V3"],nodes[arcs$V2[i],"V4"],col="grey")
points(nodes[pth,c("V3","V4")],col=darkcols[3],pch=19)
for(i in 1:(length(pth)-1)) segments(nodes[pth[i],"V3"],nodes[pth[i],"V4"],
nodes[pth[i+1],"V3"],nodes[pth[i+1],"V4"],col=darkcols[3],lwd=5)
PCH=rep(1,nrow(newhospital))
PCH[which(d==min(d))]=forme[2]
points(newhospital[,"X"],newhospital[,"Y"],pch=PCH,col=darkcols[2],cex=2,lwd=2)
points(newsick[,"X"],newsick[,"Y"],pch=forme[1],col=darkcols[1],cex=2)
