The optimize-surface script attempts to determine unobserved allocation parameters by optimization.


In [1]:
include("../src/optimize-surface.jl")


Loading from saved region network...
WARNING: module Main should explicitly import + from Base
WARNING: imported binding for edges overwritten in module Main
Loading from saved water network...
Optimize a model with 25668 rows, 24707 columns and 909122 nonzeros
Coefficient statistics:
  Matrix range    [1e+00, 1e+00]
  Objective range [1e-02, 1e-02]
  Bounds range    [0e+00, 0e+00]
  RHS range       [1e-02, 3e+08]

Concurrent LP optimizer: dual simplex and barrier
Showing barrier log only...

Presolve removed 12231 rows and 5630 columns
Presolve time: 1.54s
Presolved: 13437 rows, 19077 columns, 690318 nonzeros

Ordering time: 0.40s

Barrier performed 0 iterations in 2.40 seconds
Barrier solve interrupted - model solved by another algorithm


Solved with dual simplex
Solved in 4670 iterations and 2.54 seconds
Optimal objective  1.049115029e+06
  5.621560 seconds (1.96 M allocations: 108.189 MB, 7.13% gc time)
waterfromsupersource
[21232.142923059582,0.0,0.0,0.0,0.0,0.0,0.0,9360.43849363865,0.0,0.0,65.00170948720825,56110.81879542234,0.0,2002.0857514751724,0.0,0.0,1.7209988912970556e6,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,90515.48652000807,0.0,0.0,0.0,358392.72527034086,0.0,3981.86496,125727.37583285896,1.4348887809746969e6,0.0,0.0,0.0,68793.59049035238,0.0,3.020221249168929e6,0.0,6171.165947854569,60053.31784274879,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,21028.257741429174,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,10499.886867604593,0.0,0.0,0.0,0.0,0.0,341631.8258520396,0.0,0.0,247466.52290471178,1382.592,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
Sum: 1.0491150290368535e8
withdrawals
[30158.80171694042,0.0,0.0,0.0,25221.447657311724,0.0,0.0,0.0,0.0,6000.44928,0.0,0.0,2143.0176,0.0,0.0,38483.10270516199,0.0,34849.576974838,0.0,0.0,0.0,0.0,0.0,594.51456,41.187106361348384,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1755.8815297707333,1520.8615102292665,0.0,0.0,0.2647625296391568,2658.4397679831527,359.08924345869696,0.0,0.0,36.62700111896295,11971.989305567411,17408.09069443259,0.0,0.0,1216.68096,5101.76448,0.0,0.0,0.0,0.0,134785.22070294435,0.0,0.0,0.0,276.51840000000004,0.0,0.0,0.0,0.0,122.93852729135982,1.4947527086401746,0.0,246.89193806441537,1063.8348753185605,113.09159294032209,0.0,0.0,3138.7351936767013,0.0,0.0,884.85888,47837.6832,0.0,0.0,0.0,2225.97312,0.0,0.0,0.0,0.0,0.0,0.0,0.0,50713.47456,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,10922.4768,0.0,0.0,0.0]
Sum: 6.407456937892703e8

In [2]:
recorded = readtable("../data/extraction/USGS-2010.csv");

In [3]:
using RCall
R"library(ggplot2)"
R"library(PBSmapping)"
R"shapes <- importShapefile('../data/mapping/US_county_2000-simple')"
R"polydata <- attributes(shapes)$PolyData"
R"polydata$STATE <- as.numeric(levels(polydata$STATE))[polydata$STATE]"
R"polydata$COUNTY <- as.numeric(levels(polydata$COUNTY))[polydata$COUNTY]"
R"shapes$id <- polydata$STATE[shapes$PID] * 100 + polydata$COUNTY[shapes$PID] / 10";
R"names(shapes) <- tolower(names(shapes))";


WARNING: RCall.jl Warning: package ‘ggplot2’ was built under R version 3.2.3
WARNING: RCall.jl 
-----------------------------------------------------------
PBS Mapping 2.69.76 -- Copyright (C) 2003-2016 Fisheries and Oceans Canada

PBS Mapping comes with ABSOLUTELY NO WARRANTY;
for details see the file COPYING.
This is free software, and you are welcome to redistribute
it under certain conditions, as outlined in the above file.

A complete user guide 'PBSmapping-UG.pdf' is located at 
/Library/Frameworks/R.framework/Versions/3.2/Resources/library/PBSmapping/doc/PBSmapping-UG.pdf

Packaged on 2015-04-23
Pacific Biological Station, Nanaimo

All available PBS packages can be found at
http://code.google.com/p/pbs-software/

To see demos, type '.PBSfigs()'.
-----------------------------------------------------------


WARNING: RCall.jl Loading required package: maptools
Loading required package: sp
Checking rgeos availability: FALSE
 	Note: when rgeos is not available, polygon geometry 	computations in maptools depend on gpclib,
 	which has a restricted licence. It is disabled by default;
 	to enable gpclib, type gpclibPermit()
Loading required package: foreign

In [24]:
recorded[:supersource] = sol.sol[1:varlens[1]]
R"ggplot($recorded, aes(fill=supersource)) + geom_map(aes(map_id=FIPS), map=shapes) +
expand_limits(x=c(-2500000, 2500000), y=c(-1.4e6, 1.6e6)) +
scale_fill_gradient(name='Supersource', trans='log', low='white', high='#132B43', breaks=c(1, 10, 1e3, 1e4, 1e5, 1e6)) +
theme_bw() + theme(legend.justification=c(0,0), legend.position=c(0,0)) + xlab('') + ylab('')"


Out[24]:
RCall.RObject{RCall.VecSxp}

In [ ]: