7_优化建模

一、基本模式分析（Geoda/PySAL

PySAL - Open Source Python Library fo https://https://www.360docs.net/doc/6915579358.html,/learning/tutorials https://www.360docs.net/doc/6915579358.html,/en/latest/· pysal.cg – Computational Geom · pysal.weights — Spatial Weigh · pysal.esda — Exploratory Spat · pysal.core — Core Data Structu · pysal.inequality — Spatial Ineq · pysal.region — Spatially Const · pysal.spatial_dynamics — Spat · pysal.spreg — Regression and · https://www.360docs.net/doc/6915579358.html,work — Network Con ·

pysal.contrib – Contributed M

1．空间权重问题（pysal.weights.W

1.1 Contiguity Based Weights

根据空间对象与周围邻居的邻接关系

0 1 2 3 4 5

第七章优化建模

SAL ）

ary for Spatial Analytical Functions

l Geometry Weights y Spatial Data Analysis Structures and IO al Inequality Analysis Constrained Clustering Spatial Dynamics n and Diagnostics rk Constrained Analysis ed Modules

hts.W ）接关系，计算权重。

10 11 12 13 14

15 16 17 18 19

20 21 22 23 24

import pysal

#强连接rook

w = https://www.360docs.net/doc/6915579358.html,t2W(5, 5)

w.weights[0]

w.neighbors[0]

w.histogram

#弱连接queen

wq = https://www.360docs.net/doc/6915579358.html,t2W(rook = False)

wq.weights[0]

wq.neighbors[0]

wq.histogram

>>> import pysal

>>> w = pysal.rook_from_shapefile("./columbus.shp")

>>> w.n #shape含49个多边形，多边形标示或引用，用顺序号表示！

>>> w.histogram

[(2, 7), (3, 10), (4, 17), (5, 8), (6, 3), (7, 3), (8, 0), (9, 1)]

>>> w.weights[0]

[1.0, 1.0] #高亮多边形有两个邻居，权重均设为1

>>>

>>> w.histogram

[(2, 5), (3, 9), (4, 12), (5, 5), (6, 9), (7, 3), (8, 4), (9, 1), (10, 1)]

1.2 Distance Based Weights (基于距离的权重)

（1）k-nearest neighbor weights

找出最邻近的k个邻居，计算权重

import numpy as np

x,y=np.indices((5,5))

x.shape=(25,1)

y.shape=(25,1)

data=np.hstack([x,y])

w3 = pysal.knnW(data, k = 3)

w3.neighbors[0]

w4 = pysal.knnW(data, k = 4)

w4.neighbors[0]

>>> wknn5 = pysal.knnW_from_shapefile("./columbus.shp")

>>> wknn5 = pysal.knnW_from_shapefile("./columbus.shp",3) #k=3

>>> wknn5.neighbors[0]

[2, 1]

（2）Distance band weights

#找出相邻距离小于1度的邻居（阈值1度，约100公里）

wthresh = pysal.threshold_binaryW_from_shapefile("./columbus.shp",1.0)

wthresh.neighbors[0]

[1, 2, 3]

1.3 Kernel Weights (核函数权重)

核函数有个参数bandwidth，起到距离阈值的作用，用于选取符合阈值范围内的邻居来计算权重。如：Gaussion核时bandwidth常指标准差值。

>>> points = [(10, 10), (20, 10), (40, 10), (15, 20), (30, 20), (30, 30)]

>>> kw = pysal.Kernel(points)

>>> kw.weights[0]

[1.0, 0.500000049999995, 0.4409830615267465]

>>> kw.neighbors[0]

[0, 1, 3]

>>> kw.bandwidth

array([[ 20.000002],

[ 20.000002],

[ 20.000002]])

#所有点采用用户指定的同一带宽

kw15 = pysal.Kernel(points,bandwidth = 15.0)

#所有点采用不同的带宽

bw = [25.0,15.0,25.0,16.0,14.5,25.0]

>>> kwa = pysal.Kernel(points,bandwidth = bw)

#指定采用的核函数，基于输入数据自适应计算带宽kweag = pysal.Kernel(points,fixed = False,function = 'gaussian')

提示：

queen_from_shapefile

rook_from_shapefile

knnW_from_shapefile

threshold_binaryW_from_shapefile

threshold_continuousW_from_shapefile

kernelW_from_shapefile

adaptive_kernelW_from_shapefile

2.空间自相关(Global autocorrelation)

#计算莫兰指数

import pysal

import numpy as np

f = pysal.open("./stl_hom.txt")

y = np.array(f.by_col['HR8893'])

w = pysal.open("./stl.gal").read()

mi = pysal.Moran(y, w, two_tailed=False)

print "%.3f"%mi.I

#expected value

print mi.EI

#p-value for one-side, so two_tailed is false

print "%.5f"%mi.p_norm

#计算全局G统计量，定义邻居的阈值0.6度（约60公里内）

import pysal

import numpy as np

dist_w = pysal.threshold_binaryW_from_shapefile('./stl_hom.shp',0.6)

#transform

#B: Binary R: Row-standardization (global sum=n) D: Double-standardization (global sum=1) #V: Variance stabilizing O: Restore original transformation (from instantiation)

dist_w.transform = "B"

f = pysal.open("./stl_hom.txt")

y = np.array(f.by_col['HR8893'])

from pysal.esda.getisord import G

g = G(y, dist_w)

print g.G

print g.EG

print g.z_norm

print g.p_norm

#计算局部G统计量

from pysal.esda.getisord import G_Local

np.random.seed(12345)

lg = G_Local(y, dist_w)

lg.n

#orginal G statistic in Getis & Ord (1992)

len(lg.Gs)

#expected value of Gs under normality assumption the values is scalar

len(lg.EGs)

#standardized Gs

len(lg.Zs)

#p-value under normality assumption (one-sided)

len(lg.p_norm)

lgstar = G_Local(y, dist_w, star=True)

len(lgstar.Gs)

二、优化建模（openopt/cvxopt）

1. 了解CVXOPT

（1）矩阵matrix

>>>from cvxopt import matrix

>>>A = matrix([1.0, 2.0, 3.0, 4.0, 5.0, 6.0], (2,3)) >>>print(A)

[ 1.00e+00 3.00e+00 5.00e+00]

[ 2.00e+00 4.00e+00 6.00e+00]

>>>A.size

(2, 3)

>>>B = matrix([ [1.0, 2.0], [3.0, 4.0] ])

>>>print(B)

[ 1.00e+00 3.00e+00]

[ 2.00e+00 4.00e+00]

（2）Numpy和CVXOPT

>>>from cvxopt import matrix

>>>from numpy import array

>>>A = matrix([[1,2,3],[4,5,6]])

>>>print(A)

[ 1 4]

[ 2 5]

[ 3 6]

>>>B = array(A)

>>>B

array([[1, 4],

[2, 5],

[3, 6]])

>>>C = matrix(B)

>>>print(C)

[ 1 4]

[ 2 5]

[ 3 6]

（3）线性规划的求解

例如：

>>> from cvxopt import matri >>> A = matrix([ [-1.0, -1.0>>> b = matrix([ 1.0, -2.0>>> c = matrix([ 2.0, 1.0 >>> sol=solvers.lp(c,A,b) pcost dcost 0: 2.6471e+00 -7.0588e-01 1: 3.0726e+00 2.8437e+00 2: 2.4891e+00 2.4808e+00 3: 2.4999e+00 2.4998e+00 4: 2.5000e+00 2.5000e+00 5: 2.5000e+00 2.5000e+00 >>> print(sol['x']) [ 5.00e-01] [ 1.50e+00]

（3）二次规划的求解

例如：

matrix, solvers 1.0, 0.0, 1.0], [1.0, -1.0, -1.0, -2.02.0, 0.0, 4.0 ]) ]) gap pres dres k/t 01 2e+01 8e-01 2e+00 1e+00 e+00 1e+00 1e-01 2e-01 3e-01 e+00 1e-01 1e-02 2e-02 5e-02 e+00 1e-03 1e-04 2e-04 5e-04 e+00 1e-05 1e-06 2e-06 5e-06 e+00 1e-07 1e-08 2e-08 5e-08

2.0] ])

>>> from cvxopt import matri >>> Q = 2*matrix([ [2, .5],>>> p = matrix([1.0, 1.0])>>> G = matrix([[-1.0,0.0],[>>> h = matrix([0.0,0.0]) >>> A = matrix([1.0, 1.0],>>> b = matrix(1.0)

>>> sol=solvers.qp(Q, p, G pcost dcost 0: 0.0000e+00 0.0000e+00 1: 9.9743e-01 1.4372e+00 2: 1.8062e+00 1.8319e+00 3: 1.8704e+00 1.8693e+00 4: 1.8749e+00 1.8748e+00 5: 1.8750e+00 1.8750e+00 6: 1.8750e+00 1.8750e+00 >>> print(sol['x']) [ 2.50e-01] [ 7.50e-01]

例子：混合像元分解

matrix, solvers

], [.5, 1] ]) ]) ],[0.0,-1.0]]) ], (1,2)) G, h, A, b) gap pres dres e+00 3e+00 1e+00 0e+00 e+00 5e-01 4e-01 3e-16 e+00 5e-02 4e-02 5e-16 e+00 6e-03 2e-03 1e-15 e+00 2e-04 6e-05 6e-16 e+00 2e-06 6e-07 7e-16

e+00 2e-08 6e-09 1e-15

GF2图像–多光谱数据，4米地面分辨率

环境星HJ-1A，高光谱106个波段（共115波段），地面分辨率100米

import numpy as np

from cvxopt import matrix, solvers

from gwp_image import IMAGE

solvers.options['show_progress'] = False

#member - rows(band number), cols(endmember number), mixed - band number def unmixture(member,mixed):

bandnum,endnum = member.shape

mat = np.zeros((endnum,endnum))

for i in range(bandnum):

mat += member[i,:][:, np.newaxis]*member[i,:]

Q = 2*mat

p = -2*np.sum((mixed[:,np.newaxis]*member),axis=0)

G_01 = -np.diag(np.ones(endnum,dtype=np.int))

h_01 = np.zeros(endnum)

G_02 = np.diag(np.ones(endnum,dtype=np.int))

h_02 = np.ones(endnum)

G = np.concatenate([G_01,G_02],axis=0)

h = np.concatenate([h_01,h_02],axis=0)

Q_ =matrix(Q*1.0)

p_= matrix(p*1.0)

G_= matrix(G*1.0)

h_= matrix(h*1.0)

sol = solvers.qp(Q_, p_, G_, h_)

return np.array(sol["x"]).ravel()

#total endmember is 4, total band number is 106

fh = open('spectrum1.txt')

lns = fh.readlines()

fh.close()

lns = [ln.strip() for ln in lns]

data = []

for ln in lns:

data.append([float(v) for v in ln.split('\t')])

spectrum = np.array(data)

#read tif

drv = IMAGE()

im_proj,im_geotrans,im_data = drv.read_img("sample.tif")

im_data = np.swapaxes(np.swapaxes(im_data,0,1),1,2)

row,col,bnd = im_data.shape

bandnum,endnum = spectrum.shape

percent = np.zeros((row,col,endnum))

for i in range(row):

for j in range(col):

mixed = im_data[i,j,:]

r = unmixture(spectrum,mixed)

percent[i,j,:] = r

print percent.shape

percent = np.swapaxes(np.swapaxes(percent,1,2),0,1)

for i in range(endnum):

drv.write_img("endmember_%02d" % (i),im_proj,im_geotrans,percent[i])

练习

（1）改进计算局部G统计量的例子，将结果存入SHAPE文件。

（2）最优路径问题，试着用线性规划求解从点0到点5的最短路径。

建模提示：

21Z = 60;

x02 = 0; x04 = 1; x05 = 0; x12 = 0; x23 = 0; x35 = 1; x43 = 1; x45 = 0

Edge (i, j)0, 20,40,51,22,33,54,34,5Length c(i,j)1030100550102060Label x(i,j)

x02

x04

x05

x12

x23

x35

x43

x45

Nodes/Label x02x04x05x12x23x35x43x45

112-1

13-1

-14-1

-1

0-1

0000-1