<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 
Step1:  Automating Multiple Single-Objective Spatial Optimization Models for Efficiency and Reproducibility
James D. Gaboardi &nbsp;&nbsp; |  &nbsp;&nbsp; Association of American Geographers  2016
Florida State University  &nbsp;&nbsp; |  &nbsp;&nbsp; Department of Geography
Outline
Background & Information
Models
PMP
PCP
CentDian
PMCP Method
Data & Processing
Solutions
Visualizations
Future Work
COIN-OR
Set Path
Step2:  Background & Information
$\Rightarrow$ Automating solutions for p-median, p-center, and p-centdian problems with p={1, 2, ..., n} facilities
$\Rightarrow$ Compare coverage and costs numerically, graphically, and visually 
PySAL 1.11.0  $\Rightarrow$ PySAL.Network
Python Spatial Analysis Library
[https
Step3:  Data & Processing
Process Imports
Step4:  Data preparation
Step5:  What we have at this point
Step6:  Call Client to Service Matrix Function
Step7:  Solutions
Solve all
Step8:  Visualizations
Draw PMP figure [p=1] large $ \rightarrow $ small [p=15]
Step9:  Pandas PMP Data Frame
Step11:  Bokeh PMP [p vs. cost] trade off
Step12:  Draw PCP figure [p=1] large $ \rightarrow $ small [p=15]
Step13:  Pandas PCP Data Frame
Step15:  Bokeh PCP [p vs. cost] trade off
Step16:  Draw CentDian figure [p=1] large $ \rightarrow $ small [p=15]
Step17:  Pandas CentDian Data Frame
Step19:  Bokeh CentDian [p vs. cost] trade off
Step20:  Draw PMCP figure
Step21:  Pandas PMCP Data Frame
Step22:  Bokeh PMP & PCP [p vs. cost] comparision
Step23:  Convert Service Facilities Back to Longitude/Latitude for Google Maps Plots
Step24:  Create Lists of Selected Locations for Google Maps Plot
Step26:  Google Maps Plot
Step27:  Future Work & Vision
$\Longrightarrow$ Provide more accessible alternatives to spatial optimization & GIS
$\Longrightarrow$ Develop a python library to bring together spatial analysis & spatial optimization in one package [spanoptpy], potentially incorporating
Step28:  COIN warning that will be debugged in the next release
System Specs
Step29:  Thanks for your time! Please share with your students and local government GIS departments, etc.!
email $\Longrightarrow$ jgaboardi@fsu.edu
GitHub $\Longrightarrow$ https
<ASSISTANT_TASK:>
Python Code:
#    GNU LESSER GENERAL PUBLIC LICENSE
#    Version 3, 29 June 2007
#    Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
#    Everyone is permitted to copy and distribute verbatim copies
#    of this license document, but changing it is not allowed.
#    James Gaboardi, 2016
Explanation: 
End of explanation
import IPython.display as IPd
# Local path on user's machine
path = '/Users/jgaboardi/AAG_16/Data/'
Explanation: Automating Multiple Single-Objective Spatial Optimization Models for Efficiency and Reproducibility
James D. Gaboardi &nbsp;&nbsp; |  &nbsp;&nbsp; Association of American Geographers  2016
Florida State University  &nbsp;&nbsp; |  &nbsp;&nbsp; Department of Geography
Outline
Background & Information
Models
PMP
PCP
CentDian
PMCP Method
Data & Processing
Solutions
Visualizations
Future Work
COIN-OR
Set Path
End of explanation
# Conceptual Model Workflow
workflow = IPd.Image(path+'/AAG_16.png')
workflow
Explanation: Background & Information
$\Rightarrow$ Automating solutions for p-median, p-center, and p-centdian problems with p={1, 2, ..., n} facilities
$\Rightarrow$ Compare coverage and costs numerically, graphically, and visually 
PySAL 1.11.0  $\Rightarrow$ PySAL.Network
Python Spatial Analysis Library
[https://www.pysal.readthedocs.org]
Sergio Rey at Arizona State University leads the PySAL project.
[https://geoplan.asu.edu/people/sergio-j-rey]. 
PySAL.Network was principally developed by Jay Laura at Arizona State Universty and the United States Geological Suvery.  [https://geoplan.asu.edu/people/jay-laura]
"PySAL is an open source library of spatial analysis functions written in Python intended to support the development of high level applications. PySAL is open source under the BSD License." [https://pysal.readthedocs.org/en/latest/]
Gurobi 6.5.0 $\Rightarrow$ gurobipy
Relatively new company founded by optimization experts formerly at key positions with CPLEX.
[http://www.gurobi.com] [http://www.gurobi.com/company/about-gurobi]
Models
The p-Median Problem  $\Rightarrow$ "minimum total cost" $\Rightarrow$ efficiency
The objective of the p-median problem, also know as the minimsum problem and the PMP, is to minimize the total weighted cost while siting [p] facilities to serve all demand/client nodes.  It was originally proposed by Hakimi (1964) and is well-studied in Geography, Operations Research, Mathematics, etc.  In this particular project the network-based vertice PMP is used meaning the cost will be calculated on a road network and solutions will be determined based on discrete locations.  Cost is generally defined as either travel time or distance and it is the latter in the project.  Population (demand) utilized as a weight at each client node.  The average cost can be calculated by dividing the minimized total cost by the total demand.
For more information refer to references section.
Minimize
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\displaystyle {Z} = {\sum_{i \in 1}^n\sum_{j\in 1}^m a_i c_{ij} x_{ij}}$
Subject to
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\displaystyle\sum_{j\in m} x_{ij}  = 1 ,$&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;  $\forall i \in n$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\displaystyle\sum_{j \in m} y_j = p$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $x_{ij} - y_j \geq 0,$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\forall i \in n, j \in m$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $x_{ij}, y_j \in {0,1}$ &nbsp;&nbsp; $\forall i \in n , j \in m$
where
- $i$ = a specific origin
- $j$ = a specific destination
- $n$ = the set of origins
- $m$ = the set of destinations
- $a_i$ = weight at each node
- $c_{ij}$ = travel costs between nodes
- $x_{ij}$ = the decision variable at each node in the matrix
- $y_j$ = nodes chosen as service facilities
- $p$ = the number of facilities to be sited
Adapted from:
- Daskin, M. S. 1995. Network and Discrete Location: Models, Algorithms, and Applications. Hoboken, NJ, USA: John Wiley & Sons, Inc.
The p-Center Problem $\Rightarrow$ "minimum worst case" $\Rightarrow$ equity
The objective of the p-center problem, also know as the minimax problem and the PCP, is to 
minimize the worst case cost (W) scenario while siting [p] facilities to serve all demand/client nodes.  It was originally proposed by Minieka (1970) and, like the PMP, is well-studied in Geography, Operations Research, Mathematics, etc.  In this particular project the network-based vertice PCP is used meaning the cost will be calculated on a road network and solutions will be determined based on discrete locations.  Cost is generally defined as either travel time or distance and it is the latter in the project.
For more information refer to references section.
Minimize
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $W$
Subject to
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\displaystyle\sum_{j\in m} x_{ij} = 1,$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;   $\forall i \in n$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\displaystyle\sum_{j \in m} y_j = p$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $x_{ij} - y_j \geq 0,$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;  $\forall i\in n, j \in m$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\displaystyle W \geq \sum_{j \in m} c_{ij} x_{ij}$&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\forall i \in n$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $x_{ij}, y_j \in {0,1}$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\forall i \in n, j \in m$
where
$W$ = the worst case cost between a client and a service node
$i$ = a specific origin
$j$ = a specific destination
$n$ = the set of origins
$m$ = the set of destinations
$a_i$ = weight at each node
$c_{ij}$ = travel costs between nodes
$x_{ij}$ = the decision variable at each node in the matrix
$y_j$ = nodes chosen as service facilities
$p$ = the number of facilities to be sited
Adapted from:
- Daskin, M. S. 1995. Network and Discrete Location: Models, Algorithms, and Applications. Hoboken, NJ, USA: John Wiley & Sons, Inc.
The p-CentDian Problem
The p-CentDian Problem was first descibed by Halpern (1976).  It is a combination of the p-median problem and the p-center problem with a dual objective of minimizing both the worst case scenario and the total travel distance. The objective used for the model in this demonstration is the average of (1) the p-center objective function and (2) the p-median objective function divided by the total demand.  An alternative formulation is the p-$\lambda$-CentDian Problem, where ( $\lambda$ ) represents the weight attributed to the p-center objective function and (1 - $\lambda$) represents the weight attributed to the p-median objective function which was was proposed by Pérez-Brito, et al (1997). 
For more information refer to references section.
Minimize
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $ \displaystyle {W + {Z \over \sum_{i=1}a_i} \over 2}$
Subject to
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\displaystyle\sum_{j\in m} x_{ij} = 1,$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;   $\forall i \in n$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\displaystyle\sum_{j \in m} y_j = p$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $x_{ij} - y_j \geq 0,$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;  $\forall i\in n, j \in m$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\displaystyle W \geq \sum_{j \in m} c_{ij} x_{ij}$&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\forall i \in n$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $x_{ij}, y_j \in {0,1}$ &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $\forall i \in n, j \in m$
where
$W$ = the maximum travel cost between client and service nodes
$Z$ = the minimized total travel cost $\big({\sum_{i \in 1}^n\sum_{j\in 1}^m a_i c_{ij} x_{ij}}\big)$
$i$ = a specific origin
$j$ = a specific destination
$n$ = the set of origins
$m$ = the set of destinations
$a_i$ = weight at each node
$c_{ij}$ = travel costs between nodes
$x_{ij}$ = the decision variable at each node in the matrix
$y_j$ = nodes chosen as service facilities
$p$ = the number of facilities to be sited
Adapted from:
Halpern, J. 1976. The Location of a Center-Median Convex Combination on an Undirected Tree*. Journal of Regional Science 16 (2):237–245
The PMCP Method
$\Rightarrow$ solve the PMP and PCP to determine with optimal locations with equivalent [p]
$\Rightarrow$ "poor man's" p-CentDian Problem?
automated & efficient decision making for those who don't have access to multiple-objective capable solvers or the know-how to formulate a mutliple objectives into a single objective funtion
what it is:
a comparision to determine equivalent site selection of single objective solutions
probably best used with low cost sites
a opportunity for finding optimal solutions without sacrificing either efficiency or equity
what it is not:
an optimization solution with multiple objective functions
capable of a true 'best solution' trade-off between efficiency and equity
guaranteed to find identical solutions
Workflow for PMCP
End of explanation
import pysal as ps
import geopandas as gpd
import numpy as np
import networkx as nx
import shapefile as shp
from shapely.geometry import Point
import shapely
from collections import OrderedDict
import pandas as pd
import qgrid
import gurobipy as gbp
import time
import bokeh
from bokeh.plotting import figure, show, ColumnDataSource
from bokeh.io import output_notebook, output_file, show
from bokeh.models import (HoverTool, BoxAnnotation, GeoJSONDataSource, 
                          GMapPlot, GMapOptions, ColumnDataSource, Circle, 
                          DataRange1d, PanTool, WheelZoomTool, BoxSelectTool,
                          ResetTool, MultiLine)
import utm
from cylp.cy import CyCbcModel, CyClpSimplex
import matplotlib.pyplot as plt
import matplotlib as mpl
%matplotlib inline
mpl.rcParams['figure.figsize']=8,8
output_notebook()
Explanation: Data & Processing
Process Imports
End of explanation
import spanoptpy as sp
print dir(sp)
from spanoptpy import AAG16_DataPrep 
AAG16_DataPrep.data_prep()
Explanation: Data preparation
End of explanation
Buff.plot()
STREETS.plot()
CLIENT.plot()
SERVICE.plot(colormap=True)
Explanation: What we have at this point
End of explanation
# Call Client to Service Matrix Function
c_s_matrix() 
Explanation: Call Client to Service Matrix Function
End of explanation
Gurobi_PMCP(len(SER), Ai, AiSum, All_Dist_MILES)
Explanation: Solutions
Solve all
End of explanation
# Draw Network Actual Roads and Nodes
for e in ntw.edges:
    g.add_edge(*e)
nx.draw(g, ntw.node_coords, node_size=5, alpha=0.25, edge_color='r', width=2)
#PMP
size = 3000
for i,j in P_Med_Graphs.iteritems():
    size=size-120
    # p-Median
    P_Med = ps.open(path+'Results/Selected_Locations_Pmedian'+str(i)+'.shp')
    points_median = {}
    for idx, coords in enumerate(P_Med):
        P_Med_Graphs[i].add_node(idx)
        points_median[idx] = coords
        P_Med_Graphs[i].node[idx] = coords
    nx.draw(P_Med_Graphs[i], 
                points_median, 
                node_size=size, 
                alpha=.1, 
                node_color='k')
# Legend (Ordered Dictionary)
LEGEND = OrderedDict()
LEGEND['Network Nodes']=g
LEGEND['Roads']=g
for i in P_Med_Graphs:
    LEGEND['Optimal PMP '+str(i)]=P_Med_Graphs[i]
plt.legend(LEGEND, 
       loc='upper left', 
       fancybox=True, 
       framealpha=0.5, 
       scatterpoints=1)
# Title
plt.title('Waverly Hills\n Tallahassee, Florida', family='Times New Roman', 
      size=40, color='k', backgroundcolor='w', weight='bold')
# North Arrow and 'N' --> Must be changed for different spatial resolutions, etc.
plt.arrow(624000, 164050, 0.0, 500, width=50, head_width=125, 
          head_length=75, fc='k', ec='k',alpha=0.75,)
plt.annotate('N', xy=(623900, 164700), fontstyle='italic', fontsize='xx-large',
            fontweight='heavy', alpha=0.75)
Explanation: Visualizations
Draw PMP figure [p=1] large $ \rightarrow $ small [p=15]
End of explanation
qgrid.show_grid(pydf_M, show_toolbar=True)
Explanation: Pandas PMP Data Frame
End of explanation
source_m = ColumnDataSource(
        data=dict(
            x=range(1, len(SER)+1),
            y=AVG_PMP,
            avg=AVG_PMP,
            desc=p_list,
            change=PMP_Avg_Diff))
TOOLS = 'wheel_zoom, pan, reset, crosshair, save'
hover = HoverTool(line_policy="nearest", mode="hline", tooltips=
        <div>
            <div>
                
            </div>
            <div>
                <span style="font-size: 17px; font-weight: bold;">@desc</span> 
            </div>
            <div>
                <span style="font-size: 15px;">Average Minimized Cost</span>
                <span style="font-size: 15px; font-weight: bold; color: #ff4d4d;">[@avg]</span>
            </div>
            <div>
                <span style="font-size: 15px;">Variation</span>
                <span style="font-size: 15px; font-weight: bold; color: #ff4d4d;">[@change]</span>
            </div>
        </div>)
# Instantiate Plot
pmp_plot = figure(plot_width=600, plot_height=600, tools=[TOOLS,hover],
           title="Average Distance vs. p-Facilities", y_range=(0,2))
# Create plot points and set source
pmp_plot.circle('x', 'y', size=15, color='red',source=source_m, 
                legend='Total Minimized Cost / Total Demand')
pmp_plot.line('x', 'y', line_width=2, color='red', alpha=.5, source=source_m, 
              legend='Total Minimized Cost / Total Demand')
pmp_plot.xaxis.axis_label = '[p = n]'
pmp_plot.yaxis.axis_label = 'Miles'
one_quarter = BoxAnnotation(plot=pmp_plot, top=.35, 
                            fill_alpha=0.1, fill_color='green')
half = BoxAnnotation(plot=pmp_plot, bottom=.35, top=.7, 
                            fill_alpha=0.1, fill_color='blue')
three_quarter = BoxAnnotation(plot=pmp_plot, bottom=.7, top=1.05,
                            fill_alpha=0.1, fill_color='gray')
pmp_plot.renderers.extend([one_quarter, half, three_quarter])
# Display the figure
show(pmp_plot)
Explanation: Bokeh PMP [p vs. cost] trade off
End of explanation
# Draw Network Actual Roads and Nodes
for e in ntw.edges:
    g.add_edge(*e)
nx.draw(g, ntw.node_coords, node_size=5, alpha=0.25, edge_color='r', width=2)
#PCP
size = 3000
for i,j in P_Cent_Graphs.iteritems():
    size=size-150
    # p-Center
    P_Cent = ps.open(path+'Results/Selected_Locations_Pcenter'+str(i)+'.shp')
    points_center = {}
    for idx, coords in enumerate(P_Cent):
        P_Cent_Graphs[i].add_node(idx)
        points_center[idx] = coords
        P_Cent_Graphs[i].node[idx] = coords
    nx.draw(P_Cent_Graphs[i], 
                points_center, 
                node_size=size, 
                alpha=.1, 
                node_color='k')
# Legend (Ordered Dictionary)
LEGEND = OrderedDict()
LEGEND['Network Nodes']=g
LEGEND['Roads']=g
for i in P_Cent_Graphs:
    LEGEND['Optimal PCP '+str(i)]=P_Cent_Graphs[i]
plt.legend(LEGEND, 
       loc='upper left', 
       fancybox=True, 
       framealpha=0.5, 
       scatterpoints=1)
# Title
plt.title('Waverly Hills\n Tallahassee, Florida', family='Times New Roman', 
      size=40, color='k', backgroundcolor='w', weight='bold')
# North Arrow and 'N' --> Must be changed for different spatial resolutions, etc.
plt.arrow(624000, 164050, 0.0, 500, width=50, head_width=125, 
          head_length=75, fc='k', ec='k',alpha=0.75,)
plt.annotate('N', xy=(623900, 164700), fontstyle='italic', fontsize='xx-large',
            fontweight='heavy', alpha=0.75)
Explanation: Draw PCP figure [p=1] large $ \rightarrow $ small [p=15]
End of explanation
qgrid.show_grid(pydf_C, show_toolbar=True)
Explanation: Pandas PCP Data Frame
End of explanation
source_c = ColumnDataSource(
        data=dict(
            x=range(1, len(SER)+1),
            y=VAL_PCP,
            obj=VAL_PCP,
            desc=p_list,
            change=PCP_Diff))
    
TOOLS = 'wheel_zoom, pan, reset, crosshair, save'
    
hover = HoverTool(line_policy="nearest", mode="vline", tooltips=
        <div>
            <div>
                
            </div>
            <div>
                <span style="font-size: 17px; font-weight: bold;">@desc</span> 
            </div>
            <div>
                <span style="font-size: 15px;">Worst Case Cost</span>
                <span style="font-size: 15px; font-weight: bold; color: #00b300;">[@obj]</span>
            </div>
            <div>
                <span style="font-size: 15px;">Variation</span>
                <span style="font-size: 15px; font-weight: bold; color: #00b300;">[@change]</span>
            </div>
        </div>)
# Instantiate Plot
pcp_plot = figure(plot_width=600, plot_height=600, tools=[TOOLS,hover],
           title="Worst Case Distance vs. p-Facilities", y_range=(0,2))
# Create plot points and set source
pcp_plot.circle('x', 'y', size=15, color='green', source=source_c,
                   legend='Minimized Worst Case')
pcp_plot.line('x', 'y', line_width=2, color='green', alpha=.5, source=source_c,
                   legend='Minimized Worst Case')
pcp_plot.xaxis.axis_label = '[p = n]'
pcp_plot.yaxis.axis_label = 'Miles'
one_quarter = BoxAnnotation(plot=pcp_plot, top=.35, 
                            fill_alpha=0.1, fill_color='green')
half = BoxAnnotation(plot=pcp_plot, bottom=.35, top=.7, 
                            fill_alpha=0.1, fill_color='blue')
three_quarter = BoxAnnotation(plot=pcp_plot, bottom=.7, top=1.05,
                            fill_alpha=0.1, fill_color='gray')
pcp_plot.renderers.extend([one_quarter, half, three_quarter])
# Display the figure
show(pcp_plot)
Explanation: Bokeh PCP [p vs. cost] trade off
End of explanation
# Draw Network Actual Roads and Nodes
for e in ntw.edges:
    g.add_edge(*e)
nx.draw(g, ntw.node_coords, node_size=5, alpha=0.25, edge_color='r', width=2)
#CentDian
size = 3000
for i,j in P_CentDian_Graphs.iteritems():
    size=size-150
    P_CentDian = ps.open(path+'Results/Selected_Locations_CentDian'+str(i)+'.shp')
    points_centdian = {}
    for idx, coords in enumerate(P_CentDian):
        P_CentDian_Graphs[i].add_node(idx)
        points_centdian[idx] = coords
        P_CentDian_Graphs[i].node[idx] = coords
    nx.draw(P_CentDian_Graphs[i], 
                points_centdian, 
                node_size=size, 
                alpha=.1, 
                node_color='k')
# Legend (Ordered Dictionary)
LEGEND = OrderedDict()
LEGEND['Network Nodes']=g
LEGEND['Roads']=g
for i in P_CentDian_Graphs:
    LEGEND['Optimal CentDian '+str(i)]=P_CentDian_Graphs[i]
plt.legend(LEGEND, 
       loc='upper left', 
       fancybox=True, 
       framealpha=0.5, 
       scatterpoints=1)
# Title
plt.title('Waverly Hills\n Tallahassee, Florida', family='Times New Roman', 
      size=40, color='k', backgroundcolor='w', weight='bold')
# North Arrow and 'N' --> Must be changed for different spatial resolutions, etc.
plt.arrow(624000, 164050, 0.0, 500, width=50, head_width=125, 
          head_length=75, fc='k', ec='k',alpha=0.75,)
plt.annotate('N', xy=(623900, 164700), fontstyle='italic', fontsize='xx-large',
            fontweight='heavy', alpha=0.75)
Explanation: Draw CentDian figure [p=1] large $ \rightarrow $ small [p=15]
End of explanation
qgrid.show_grid(pydf_CentDian, show_toolbar=True)
Explanation: Pandas CentDian Data Frame
End of explanation
source_centdian = ColumnDataSource(
        data=dict(
            x=range(1, len(SER)+1),
            y=VAL_CentDian,
            obj=VAL_CentDian,
            desc=p_list,
            change=CentDian_Diff))
    
TOOLS = 'wheel_zoom, pan, reset, crosshair, save'
    
hover = HoverTool(line_policy="nearest", mode="vline" ,tooltips=
        <div>
            <div>
                
            </div>
            <div>
                <span style="font-size: 17px; font-weight: bold;">@desc</span> 
            </div>
            <div>
                <span style="font-size: 15px;">Center Median Cost</span>
                <span style="font-size: 15px; font-weight: bold; color: #3385ff;">[@obj]</span>
            </div>
            <div>
                <span style="font-size: 15px;">Variation</span>
                <span style="font-size: 15px; font-weight: bold; color: #3385ff;">[@change]</span>
            </div>
        </div>)
# Instantiate Plot
centdian_plot = figure(plot_width=600, plot_height=600, tools=[TOOLS,hover],
           title="Center Median Distance vs. p-Facilities", y_range=(0,2))
# Create plot points and set source
centdian_plot.circle('x', 'y', size=15, color='blue', source=source_centdian,
                   legend='Center Median')
centdian_plot.line('x', 'y', line_width=2, color='blue', alpha=.5, source=source_centdian,
                   legend='Center Median')
centdian_plot.xaxis.axis_label = '[p = n]'
centdian_plot.yaxis.axis_label = 'Miles'
one_quarter = BoxAnnotation(plot=pcp_plot, top=.35, 
                            fill_alpha=0.1, fill_color='green')
half = BoxAnnotation(plot=pcp_plot, bottom=.35, top=.7, 
                            fill_alpha=0.1, fill_color='blue')
three_quarter = BoxAnnotation(plot=pcp_plot, bottom=.7, top=1.05,
                            fill_alpha=0.1, fill_color='gray')
centdian_plot.renderers.extend([one_quarter, half, three_quarter])
# Display the figure
show(centdian_plot)
Explanation: Bokeh CentDian [p vs. cost] trade off
End of explanation
# Draw Network Actual Roads and Nodes
for e in ntw.edges:
    g.add_edge(*e)
nx.draw(g, ntw.node_coords, node_size=5, alpha=0.25, edge_color='r', width=2)
size = 500
shape = 'sdh^Vp<8>'
counter = -1
for i,j in PMCP_Graphs.iteritems():
    if int(i) <= len(SER)-1:
        counter = counter+1
        pmcp = ps.open(path+'Results/Selected_Locations_PMCP'+str(i)+'.shp')
        points_pmcp = {}
        for idx, coords in enumerate(pmcp):
            PMCP_Graphs[i].add_node(idx)
            points_pmcp[idx] = coords
            PMCP_Graphs[i].node[idx] = coords
        nx.draw(PMCP_Graphs[i], 
                    points_pmcp, 
                    node_size=size,
                    node_shape=shape[counter],
                    alpha=.5, 
                    node_color='k')
    else:
        pass
    
# Legend (Ordered Dictionary)
LEGEND = OrderedDict()
LEGEND['Network Nodes']=g
LEGEND['Roads']=g
for i in PMCP_Graphs:
    if int(i) <= len(SER)-1:
        LEGEND['PMP/PCP == '+str(i)]=PMCP_Graphs[i]
plt.legend(LEGEND, 
       loc='upper left', 
       fancybox=True, 
       framealpha=0.5, 
       scatterpoints=1)
# Title
plt.title('Waverly Hills\n Tallahassee, Florida', family='Times New Roman', 
      size=40, color='k', backgroundcolor='w', weight='bold')
# North Arrow and 'N' --> Must be changed for different spatial resolutions, etc.
plt.arrow(624000, 164050, 0.0, 500, width=50, head_width=125, 
          head_length=75, fc='k', ec='k',alpha=0.75,)
plt.annotate('N', xy=(623900, 164700), fontstyle='italic', fontsize='xx-large',
            fontweight='heavy', alpha=0.75)
Explanation: Draw PMCP figure
End of explanation
qgrid.show_grid(pydf_MC, show_toolbar=True)
Explanation: Pandas PMCP Data Frame
End of explanation
TOOLS = 'wheel_zoom, pan, reset, crosshair, save, hover'
bokeh_df_PMCP = pd.DataFrame()
bokeh_df_PMCP['p'] = [int(i[2:]) for i in p_dens]
bokeh_df_PMCP['Total Obj. Value'] = [VAL_PMCP[x][0] for x in range(len(VAL_PMCP))]
bokeh_df_PMCP['Avg. Obj. Value'] = [VAL_PMCP[x][1] for x in range(len(VAL_PMCP))]
bokeh_df_PMCP['Worst Case Obj. Value'] = [VAL_PMCP[x][2] for x in range(len(VAL_PMCP))]
bokeh_df_PMCP['CentDian Obj. Value'] = [VAL_PMCP[x][3] for x in range(len(VAL_PMCP))]
plot_PMCP = figure(title="Optimal PMP & PCP Selections without Sacrifice", 
                        plot_width=800, plot_height=600, tools=[TOOLS], y_range=(0,2))
plot_PMCP.circle('x', 'y', size=5, color='red', source=source_m, legend='PMP')
plot_PMCP.line('x', 'y', 
                 color="#ff4d4d", alpha=0.2, line_width=2, source=source_m, legend='PMP')
plot_PMCP.circle('x', 'y', size=5, color='green', source=source_c, legend='PCP')
plot_PMCP.line('x', 'y', 
                 color='#00b300', alpha=0.2, line_width=2, source=source_c, legend='PCP')
plot_PMCP.circle('x', 'y', size=5, color='blue', source=source_centdian, legend='CentDian')
plot_PMCP.line('x', 'y', 
                 color='#3385ff', alpha=0.2, line_width=2, source=source_centdian, legend='CentDian')
plot_PMCP.circle_x(bokeh_df_PMCP['p'], 
                 bokeh_df_PMCP['Avg. Obj. Value'], 
                 legend="Location PMP=PCP for PM+CP", 
                 color="#ff4d4d",
                 fill_alpha=0.2,
                 size=15)
plot_PMCP.circle_x(bokeh_df_PMCP['p'], 
             bokeh_df_PMCP['Worst Case Obj. Value'], 
             legend="Location PCP=PMP for PM+CP", 
             color='#00b300',
             fill_alpha=0.2,
             size=15)
plot_PMCP.circle_x(bokeh_df_PMCP['p'], 
             bokeh_df_PMCP['CentDian Obj. Value'], 
             legend="Location CentDian = PMCP", 
             color='#3385ff',
             fill_alpha=0.2,
             size=15)
plot_PMCP.xaxis.axis_label = '[p = n]'
plot_PMCP.yaxis.axis_label = 'Miles'
one_quarter = BoxAnnotation(plot=plot_PMCP, top=.35, 
                            fill_alpha=0.1, fill_color='green')
half = BoxAnnotation(plot=plot_PMCP, bottom=.35, top=.7, 
                            fill_alpha=0.1, fill_color='blue')
three_quarter = BoxAnnotation(plot=plot_PMCP, bottom=.7, top=1.05,
                            fill_alpha=0.1, fill_color='gray')
plot_PMCP.renderers.extend([one_quarter, half, three_quarter])
show(plot_PMCP)
Explanation: Bokeh PMP & PCP [p vs. cost] comparision
End of explanation
points = SERVICE
points.to_crs(epsg=32616, inplace=True) # UTM 16N
LonLat_Dict = OrderedDict()
LonLat_List = []
for i,j in points['geometry'].iteritems():
    LonLat_Dict[y_list[i]] = utm.to_latlon(j.xy[0][-1], j.xy[1][-1], 16, 'N')  
    LonLat_List.append((utm.to_latlon(j.xy[0][-1], j.xy[1][-1], 16, 'N')))
Service_Lat_List = []
Service_Lon_List = []
    
for i in LonLat_List:
    Service_Lat_List.append(i[0])
for i in LonLat_List:
    Service_Lon_List.append(i[1])
Explanation: Convert Service Facilities Back to Longitude/Latitude for Google Maps Plots
End of explanation
# p-Median Selected Sites
list_of_p_MEDIAN = []
for y in range(len(y_list)):
    list_of_p_MEDIAN.append([])
    for p in range(len(p_list)):
        if pydf_M[y_list[y]][p_list[p]] == u'\u2588':
            list_of_p_MEDIAN[y].append([p_list[p]])
            
# p-Center Selected Sites
list_of_p_CENTER = []
for y in range(len(y_list)):
    list_of_p_CENTER.append([])
    for p in range(len(p_list)):
        if pydf_C[y_list[y]][p_list[p]] == u'\u2588':
            list_of_p_CENTER[y].append([p_list[p]])
            
# p-CentDian Selected Sites
list_of_p_CentDian = []
for y in range(len(y_list)):
    list_of_p_CentDian.append([])
    for p in range(len(p_list)):
        if pydf_CentDian[y_list[y]][p_list[p]] == u'\u2588':
            list_of_p_CentDian[y].append([p_list[p]])
# PMCP Selected Sites
list_of_PMCP = []
for y in range(len(y_list)):
    list_of_PMCP.append([])
    for p in p_dens:
        if pydf_MC[y_list[y]][p] == u'\u2588':
            list_of_PMCP[y].append(p)
Explanation: Create Lists of Selected Locations for Google Maps Plot
End of explanation
map_options = GMapOptions(lat=30.4855, lng=-84.265, map_type="hybrid", zoom=14)
plot = GMapPlot(
    x_range=DataRange1d(), y_range=DataRange1d(), map_options=map_options, title="Waverly Hills")
hover = HoverTool(tooltips=
        <div>
            <div>
                
            </div>
            <div>
                <span style="font-size: 30px; font-weight: bold;">Site @desc</span> 
            </div>
            <div>
                <span> \b </span>
            </div>
            <div>
                <span style="font-size: 17px; font-weight: bold;">p-Median: </span> 
            </div>
                <div>
                <span style="font-size: 15px; font-weight: bold; color: #ff4d4d;">@p_select_median</span>
            </div>
            <div>
                <span> \b </span>
            </div>
            <div>
                <span style="font-size: 17px; font-weight: bold;">p-Center</span> 
            
            <div>
                <span style="font-size: 14px; font-weight: bold; color: #00b300;">@p_select_center</span>
            </div>
            <div>
                <span> \b </span>
            </div>
            <div>
                <span style="font-size: 17px; font-weight: bold;">p-CentDian</span> 
            </div>
            <div>
                <span style="font-size: 14px; font-weight: bold; color: #3385ff;">@p_select_centdian</span>
            </div>
            <div>
                <span> \b </span>
            </div>
            <span style="font-size: 17px; font-weight: bold;">PMCP Method</span> 
            </div>
            <div>
                <span style="font-size: 14px; font-weight: bold; color: 'gray';">@p_select_pmcp</span>
            </div>
        </div>)
source_1 = ColumnDataSource(
        data=dict(
        lat=Service_Lat_List,
        lon=Service_Lon_List,
        desc=y_list,
        p_select_center=list_of_p_CENTER,
        p_select_median=list_of_p_MEDIAN,
        p_select_centdian= list_of_p_CentDian,
        p_select_pmcp=list_of_PMCP))
#source_2 = ColumnDataSource(
#    data=dict(
#        xs=line1x,
#        ys=line1y))
facilties = Circle(x="lon", y="lat", size=10, fill_color="yellow", fill_alpha=0.6, line_color=None)
#streets = MultiLine(xs="xs", ys="ys", line_width=20, line_color='red')
#plot.title = "Waverly"
plot.add_glyph(source_1, facilties)
#plot.add_glyph(source_2, streets)
plot.add_tools(PanTool(), WheelZoomTool(), BoxSelectTool(), ResetTool(), hover)
output_file("gmap_plot.html")
show(plot)
Explanation: Google Maps Plot
End of explanation
# Gurobi
m = gbp.Model()
m.setParam('OutputFlag', False) 
x = m.addVar(vtype=gbp.GRB.CONTINUOUS, name='x')
y = m.addVar(vtype=gbp.GRB.CONTINUOUS, name='y')
m.update()
m.setObjective(3*x + 2*y, gbp.GRB.MINIMIZE)
m.addConstr(x >= 3)
m.addConstr(y >= 5)
m.addConstr(x - y <= 20)
m.optimize()
#m.write('path_m.lp')
print m.objVal
print m.getVars()
# CyLP
s = CyClpSimplex()
x = s.addVariable('x', 1)
y = s.addVariable('y', 1)
s.addConstraint(x >= 3)
s.addConstraint(y >= 5)
s.addConstraint(x - y <= 20)
s.objective = 3 * x + 2 * y
s.primal()
#s.writeLp('path_s')
print s.objectiveValue
print s.primalVariableSolution
print 'Gurobi & CLP Objective Values match? --> ', m.objVal == s.objectiveValue
Explanation: Future Work & Vision
$\Longrightarrow$ Provide more accessible alternatives to spatial optimization & GIS
$\Longrightarrow$ Develop a python library to bring together spatial analysis & spatial optimization in one package [spanoptpy], potentially incorporating:
|QGIS|  PySAL  |  NetworkX  |  Pandas & QGrid |  GeoPandas  |  NumPy / SciPy  |  Shapely  |  Bokeh | CyLP | PostgreSQL & PostGIS
|----|---------|-------------------------------------------------------------------------
|GIS|network analysis|network analysis|data frames|geo dataframes|scientific computing|geometric objects|visualizations|optimization | DBMS
$\Longrightarrow$ Need PySAL.Network to be able to handle larger networks
$\Longrightarrow$ Intergration with PostgreSQL & PostGIS for data storage
$\Longrightarrow$ Develop functionality within a Linux environment
$\Longrightarrow$  scipy.spatial.cKDTree(dist_matrix)
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$\ast$  query_ball_point() for close neighbors of the selected sites
$\Longrightarrow$ Master CyLP from COIN-OR or develop an open-source optimization suite
interface with CLP, CBC, CGL
[ http://mpy.github.io/CyLPdoc/ ]
relatively steep learning curve 
Computational Optimization Infrastructure for Operations Research
[ http://www.coin-or.org ]
$\ast$ CyLP example
Minimize
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $ 3x + 2y $
Subject To
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $ x \geq 3$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $ y \geq 5$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; $ x - y \leq 20$
End of explanation
import datetime as dt
import os
import platform
import sys
import bokeh
import cylp
names = ['OSX', 'Processor ', 'Machine ', 'Python ','PySAL ','Gurobi ','Pandas ','GeoPandas ',
         'Shapely ', 'NumPy ', 'Bokeh ', 'CyLP', 'Date & Time']
versions = [platform.mac_ver()[0], platform.processor(), platform.machine(), platform.python_version(),
            ps.version, gbp.gurobi.version(), pd.__version__, gpd.__version__, 
            str(shapely.__version__), np.__version__, 
            bokeh.__version__, '0.7.1', dt.datetime.now()]
specs = pd.DataFrame(index=names, columns=['Version'])
specs.columns.name = 'Platform & Software Specs'
specs['Version'] = versions
specs # Pandas DF of specifications
Explanation: COIN warning that will be debugged in the next release
System Specs
End of explanation
IPd.HTML('https://github.com/jGaboardi')
Explanation: Thanks for your time! Please share with your students and local government GIS departments, etc.!
email $\Longrightarrow$ jgaboardi@fsu.edu
GitHub $\Longrightarrow$ https://github.com/jGaboardi/AAG_16
End of explanation
<END_TASK>