Category Archives: Julia

Building our own graph type in Julia

By: Julia on μβ

Re-posted from: https://matbesancon.github.io/post/2018-08-17-abstract_graph/


This is an adapted post on the talk we gave with James
at JuliaCon 2018 in London. You can see the
original slides,
the video still requires a bit of post-processing.

Last week JuliaCon in London was a great and very condensed experience.
The two talks on LightGraphs.jl
received a lot of positive feedback and more than that, we saw
how people are using the library for a variety of use cases which is a great
signal for the work on the JuliaGraphs ecosystem
(see the lightning talk).

I wanted to re-build the same graph for people who prefer a post version to
my clumsy live explanations on a laptop not handling dual-screen well
(those who prefer the latter are invited to see the live-stream of the talk).

Why abstractions?

The LightGraphs library is built to contain as few elements as possible to get
anyone going with graphs. This includes:

  • The interface a graph type has to comply with to be used
  • Essential algorithms implemented by any graph respecting that interface
  • A simple, battery-included implementation based on adjacency lists

The thing is, if you design an abstraction which in fact has just one
implementation, you’re doing abstraction wrong. This talks was also a
reality-check for LightGraphs, are we as composable, extensible as we promised?

The reason for abstraction is also that minimalism has its price.
The package was designed as the least amount of complexity required to get
graphs working. When people started to use it, obviously they needed more
features, some of which they could code themselves, some other required
extensions built within LightGraphs. By getting the core abstractions right,
you guarantee people will be able to use it and to build on top with minimal
friction, while keeping it simple to read and contribute to.

Our matrix graph type

Let’s recall that a graph is a collection of nodes and a collection of
edges between these nodes. To keep it simple, for a graph of $n$ edges,
we will consider they are numbered from 1 to n. An edge connects a node $i$
to a node $j$, therefore all the information of a graph can be kept as an
adjacency matrix $M_{ij}$ of size $n \times n$:

$$M_{ij} = \begin{cases} 1, & \mbox{if edge (i $\rightarrow$ j) exists} \\ 0 & \mbox{otherwise}\end{cases}$$

We don’t know what the use cases for our type will be, and therefore,
we will parametrize the graph type over the matrix type:

import LightGraphs; const lg = LightGraphs
mutable struct MatrixDiGraph{MT <: AbstractMatrix{Bool}} <: lg.AbstractGraph{Int}
  matrix::MT
end

The edges are simply mapping an entry (i,j) to a boolean (whether there is an
edge from i to j). Even though creating a graph type that can be directed
or undirected depending on the situation is possible, we are creating a type
that will be directed by default.

Implementing the core interface

We can now implement the core LightGraphs interface for this type, starting
with methods defined over the type itself, of the form function(g::MyType)

I’m not going to re-define each function here, their meaning can be found
by checking the help in a Julia REPL: ?LightGraphs.vertices or on the
documentation page.

lg.is_directed(::MatrixDiGraph) = true
lg.edgetype(::MatrixDiGraph) = lg.SimpleGraphs.SimpleEdge{Int}
lg.ne(g::MatrixDiGraph) = sum(g.m)
lg.nv(g::MatrixDiGraph) = size(g.m)[1]

lg.vertices(g::MatrixDiGraph) = 1:nv(g)

function lg.edges(g::MatrixDiGraph)
    n = lg.nv(g)
    return (lg.SimpleGraphs.SimpleEdge(i,j) for i in 1:n for j in 1:n if g.m[i,j])
end

Note the last function edges, for which the documentation specifies that we
need to return an iterator over edges. We don’t need to collect the comprehension
in a Vector, returning a lazy generator.

Some operations have to be defined on both the graph and a node, of the form
function(g::MyType, node).

lg.outneighbors(g::MatrixDiGraph, node) = [v for v in 1:lg.nv(g) if g.m[node, v]]
lg.inneighbors(g::MatrixDiGraph, node) = [v for v in 1:lg.nv(g) if g.m[v, node]]
lg.has_vertex(g::MatrixDiGraph, v::Integer) = v <= lg.nv(g) && v > 0

Out MatrixDiGraph type is pretty straight-forward to work with and all
required methods are easy to relate to the way information is stored in the
adjacency matrix.

The last step is implementing methods on both the graph and an edge of the
form function(g::MatrixDiGraph,e). The only one we need here is:

lg.has_edge(g::MatrixDiGraph,i,j) = g.m[i,j]

Optional mutability

Mutating methods were removed from the core interace some time ago,
as they are not required to describe a graph-like behavior.
The general behavior for operations mutating a graph is to return whether
the operation succeded. They consist in adding or removing elements from
either the edges or nodes.

import LightGraphs: rem_edge!, rem_vertex!, add_edge!, add_vertex!

function add_edge!(g::MatrixDiGraph, e)
    has_edge(g,e) && return false
    n = nv(g)
    (src(e) > n || dst(e) > n) && return false
    g.m[src(e),dst(e)] = true
end

function rem_edge!(g::MatrixDiGraph,e)
    has_edge(g,e) || return false
    n = nv(g)
    (src(e) > n || dst(e) > n) && return false
    g.m[src(e),dst(e)] = false
    return true
end

function add_vertex!(g::MatrixDiGraph)
    n = nv(g)
    m = zeros(Bool,n+1,n+1)
    m[1:n,1:n] .= g.m
    g.m = m
    return true
end

Testing our graph type on real data

We will use the graph type to compute the PageRank of

import SNAPDatasets
data = SNAPDatasets.loadsnap(:ego_twitter_d)
twitter_graph = MatrixDiGraph(lg.adjacency_matrix(data)[1:10,1:10].==1);
ranks = lg.pagerank(twitter_graph)

Note the broadcast check .==1, adjacency_matrix is specified to yield a
matrix of Int, so we use this to cast the entries to boolean values.

I took only the first 10 nodes of the graph, but feel free to do the same with
500, 1000 or more nodes, depending on what your machine can stand ?

Overloading non-mandatory functions

Some methods are already implemented for free by implementing the core interface.
That does not mean it should be kept as-is in every case. Depending on your
graph type, some functions might have smarter implementations, let’s see one
example. What MatrixDiGraph is already an adjacency_matrix, so we know
there should be no computation required to return it (it’s almost a no-op).

using BenchmarkTools: @btime

@btime adjacency_matrix(bigger_twitter)
println("why did that take so long?")
lg.adjacency_matrix(g::MatrixDiGraph) = Int.(g.m)
@btime A = lg.adjacency_matrix(bigger_twitter)
println("that's better.")

This should yield roughly:

13.077 ms (5222 allocations: 682.03 KiB)
why did that take so long?
82.077 μs (6 allocations: 201.77 KiB)
that's better.

You can fall down to a no-op by storing the matrix entries as Int directly,
but the type ends up being a bit heavier in memory, your type, your trade-off.

Conclusion

We’ve implemented a graph type suited to our need in a couple lines of Julia,
guided by the LightGraphs interface specifying how to think about our
graph instead of getting in the way of what to store. A lighter version
of this post can be read as slides.

As usual, ping me on Twitter for any
question or comment.

Bonus

If you read this and want to try building your own graph type, here are two
implementations you can try out, put them out in a public repo and show them off
afterwards:
1. We created a type just for directed graphs, why bothering so much? You can create your own type which can be directed or not,
either by storing the information in the struct or by parametrizing the type
and getting the compiler to do the work for you.
2. We store the entries as an AbstractMatrix{Bool}, if your graph is dense
enough (how dense? No idea), it might be interesting to store entries as as
BitArray.


Image source: GraphPlot.jl

Building our own graph type in Julia

By: Julia on μβ

Re-posted from: https://matbesancon.xyz/post/2018-08-17-abstract_graph/


This is an adapted post on the talk we gave with James
at JuliaCon 2018 in London. You can see the
original slides,
the video still requires a bit of post-processing.

Last week JuliaCon in London was a great and very condensed experience.
The two talks on LightGraphs.jl
received a lot of positive feedback and more than that, we saw
how people are using the library for a variety of use cases which is a great
signal for the work on the JuliaGraphs ecosystem
(see the lightning talk).

I wanted to re-build the same graph for people who prefer a post version to
my clumsy live explanations on a laptop not handling dual-screen well
(those who prefer the latter are invited to see the live-stream of the talk).

Why abstractions?

The LightGraphs library is built to contain as few elements as possible to get
anyone going with graphs. This includes:

  • The interface a graph type has to comply with to be used
  • Essential algorithms implemented by any graph respecting that interface
  • A simple, battery-included implementation based on adjacency lists

The thing is, if you design an abstraction which in fact has just one
implementation, you’re doing abstraction wrong. This talks was also a
reality-check for LightGraphs, are we as composable, extensible as we promised?

The reason for abstraction is also that minimalism has its price.
The package was designed as the least amount of complexity required to get
graphs working. When people started to use it, obviously they needed more
features, some of which they could code themselves, some other required
extensions built within LightGraphs. By getting the core abstractions right,
you guarantee people will be able to use it and to build on top with minimal
friction, while keeping it simple to read and contribute to.

Our matrix graph type

Let’s recall that a graph is a collection of nodes and a collection of
edges between these nodes. To keep it simple, for a graph of $n$ edges,
we will consider they are numbered from 1 to n. An edge connects a node $i$
to a node $j$, therefore all the information of a graph can be kept as an
adjacency matrix $M_{ij}$ of size $n \times n$:

$$M_{ij} = \begin{cases} 1, & \mbox{if edge (i $\rightarrow$ j) exists} \\ 0 & \mbox{otherwise}\end{cases}$$

We don’t know what the use cases for our type will be, and therefore,
we will parametrize the graph type over the matrix type:

import LightGraphs; const lg = LightGraphs
mutable struct MatrixDiGraph{MT <: AbstractMatrix{Bool}} <: lg.AbstractGraph{Int}
  matrix::MT
end

The edges are simply mapping an entry (i,j) to a boolean (whether there is an
edge from i to j). Even though creating a graph type that can be directed
or undirected depending on the situation is possible, we are creating a type
that will be directed by default.

Implementing the core interface

We can now implement the core LightGraphs interface for this type, starting
with methods defined over the type itself, of the form function(g::MyType)

I’m not going to re-define each function here, their meaning can be found
by checking the help in a Julia REPL: ?LightGraphs.vertices or on the
documentation page.

lg.is_directed(::MatrixDiGraph) = true
lg.edgetype(::MatrixDiGraph) = lg.SimpleGraphs.SimpleEdge{Int}
lg.ne(g::MatrixDiGraph) = sum(g.m)
lg.nv(g::MatrixDiGraph) = size(g.m)[1]

lg.vertices(g::MatrixDiGraph) = 1:nv(g)

function lg.edges(g::MatrixDiGraph)
    n = lg.nv(g)
    return (lg.SimpleGraphs.SimpleEdge(i,j) for i in 1:n for j in 1:n if g.m[i,j])
end

Note the last function edges, for which the documentation specifies that we
need to return an iterator over edges. We don’t need to collect the comprehension
in a Vector, returning a lazy generator.

Some operations have to be defined on both the graph and a node, of the form
function(g::MyType, node).

lg.outneighbors(g::MatrixDiGraph, node) = [v for v in 1:lg.nv(g) if g.m[node, v]]
lg.inneighbors(g::MatrixDiGraph, node) = [v for v in 1:lg.nv(g) if g.m[v, node]]
lg.has_vertex(g::MatrixDiGraph, v::Integer) = v <= lg.nv(g) && v > 0

Out MatrixDiGraph type is pretty straight-forward to work with and all
required methods are easy to relate to the way information is stored in the
adjacency matrix.

The last step is implementing methods on both the graph and an edge of the
form function(g::MatrixDiGraph,e). The only one we need here is:

lg.has_edge(g::MatrixDiGraph,i,j) = g.m[i,j]

Optional mutability

Mutating methods were removed from the core interace some time ago,
as they are not required to describe a graph-like behavior.
The general behavior for operations mutating a graph is to return whether
the operation succeded. They consist in adding or removing elements from
either the edges or nodes.

import LightGraphs: rem_edge!, rem_vertex!, add_edge!, add_vertex!

function add_edge!(g::MatrixDiGraph, e)
    has_edge(g,e) && return false
    n = nv(g)
    (src(e) > n || dst(e) > n) && return false
    g.m[src(e),dst(e)] = true
end

function rem_edge!(g::MatrixDiGraph,e)
    has_edge(g,e) || return false
    n = nv(g)
    (src(e) > n || dst(e) > n) && return false
    g.m[src(e),dst(e)] = false
    return true
end

function add_vertex!(g::MatrixDiGraph)
    n = nv(g)
    m = zeros(Bool,n+1,n+1)
    m[1:n,1:n] .= g.m
    g.m = m
    return true
end

Testing our graph type on real data

We will use the graph type to compute the PageRank of

import SNAPDatasets
data = SNAPDatasets.loadsnap(:ego_twitter_d)
twitter_graph = MatrixDiGraph(lg.adjacency_matrix(data)[1:10,1:10].==1);
ranks = lg.pagerank(twitter_graph)

Note the broadcast check .==1, adjacency_matrix is specified to yield a
matrix of Int, so we use this to cast the entries to boolean values.

I took only the first 10 nodes of the graph, but feel free to do the same with
500, 1000 or more nodes, depending on what your machine can stand ?

Overloading non-mandatory functions

Some methods are already implemented for free by implementing the core interface.
That does not mean it should be kept as-is in every case. Depending on your
graph type, some functions might have smarter implementations, let’s see one
example. What MatrixDiGraph is already an adjacency_matrix, so we know
there should be no computation required to return it (it’s almost a no-op).

using BenchmarkTools: @btime

@btime adjacency_matrix(bigger_twitter)
println("why did that take so long?")
lg.adjacency_matrix(g::MatrixDiGraph) = Int.(g.m)
@btime A = lg.adjacency_matrix(bigger_twitter)
println("that's better.")

This should yield roughly:

13.077 ms (5222 allocations: 682.03 KiB)
why did that take so long?
82.077 μs (6 allocations: 201.77 KiB)
that's better.

You can fall down to a no-op by storing the matrix entries as Int directly,
but the type ends up being a bit heavier in memory, your type, your trade-off.

Conclusion

We’ve implemented a graph type suited to our need in a couple lines of Julia,
guided by the LightGraphs interface specifying how to think about our
graph instead of getting in the way of what to store. A lighter version
of this post can be read as slides.

As usual, ping me on Twitter for any
question or comment.

Bonus

If you read this and want to try building your own graph type, here are two
implementations you can try out, put them out in a public repo and show them off
afterwards:

  1. We created a type just for directed graphs, why bothering so much? You can create your own type which can be directed or not,
    either by storing the information in the struct or by parametrizing the type
    and getting the compiler to do the work for you.
  2. We store the entries as an AbstractMatrix{Bool}, if your graph is dense
    enough (how dense? No idea), it might be interesting to store entries as as
    BitArray.

Image source: GraphPlot.jl

Validation of the linear regressor model

By: Sören Dobberschütz

Re-posted from: https://tensorflowjulia.blogspot.com/2018/08/validation-of-linear-regressor-model.html

The third part of the Machine Learning Crash Course deals with validation of the model.

The Jupyter notebook can be downloaded here. For the version displayed below, I needed to remove some scatter plots, which are contained in the original file.


This notebook is based on the file Validation programming exercise, which is part of Google’s Machine Learning Crash Course.
In [0]:
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# https://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

Validation

Learning Objectives:
  • Use multiple features, instead of a single feature, to further improve the effectiveness of a model
  • Debug issues in model input data
  • Use a test data set to check if a model is overfitting the validation data
As in the prior exercises, we’re working with the California housing data set, to try and predict median_house_value at the city block level from 1990 census data.

Setup

First off, let’s load up and prepare our data. This time, we’re going to work with multiple features, so we’ll modularize the logic for preprocessing the features a bit:
In [14]:
# Load packages
using Plots
gr()
using DataFrames
using TensorFlow
import CSV

# Start a TensorFlow session and load the data
sess=Session()
california_housing_dataframe = CSV.read("california_housing_train.csv", delim=",");
#california_housing_dataframe = california_housing_dataframe[shuffle(1:size(california_housing_dataframe, 1)),:];
In [2]:
function preprocess_features(california_housing_dataframe)
"""Prepares input features from California housing data set.

Args:
california_housing_dataframe: A DataFrame expected to contain data
from the California housing data set.
Returns:
A DataFrame that contains the features to be used for the model, including
synthetic features.
"""
selected_features = california_housing_dataframe[
[:latitude,
:longitude,
:housing_median_age,
:total_rooms,
:total_bedrooms,
:population,
:households,
:median_income]]
processed_features = selected_features
# Create a synthetic feature.
processed_features[:rooms_per_person] = (
california_housing_dataframe[:total_rooms] ./
california_housing_dataframe[:population])
return processed_features
end

function preprocess_targets(california_housing_dataframe)
"""Prepares target features (i.e., labels) from California housing data set.

Args:
california_housing_dataframe: A DataFrame expected to contain data
from the California housing data set.
Returns:
A DataFrame that contains the target feature.
"""
output_targets = DataFrame()
# Scale the target to be in units of thousands of dollars.
output_targets[:median_house_value] = (
california_housing_dataframe[:median_house_value] ./ 1000.0)
return output_targets
end
Out[2]:
preprocess_targets (generic function with 1 method)
2018-08-13 20:33:55.100558: I tensorflow/core/platform/cpu_feature_guard.cc:140] Your CPU supports instructions that this TensorFlow binary was not compiled to use: SSE4.2 AVX AVX2 FMA
For the training set, we’ll choose the first 12000 examples, out of the total of 17000.
In [15]:
training_examples = preprocess_features(head(california_housing_dataframe,12000))
describe(training_examples)
Out[15]:
variable mean min median max nunique nmissing eltype
1 latitude 35.6415 32.54 34.255 41.95 0 Float64
2 longitude -119.583 -124.35 -118.52 -114.31 0 Float64
3 housing_median_age 28.6681 1.0 29.0 52.0 0 Float64
4 total_rooms 2644.53 11.0 2139.5 28258.0 0 Float64
5 total_bedrooms 540.689 3.0 436.0 4819.0 0 Float64
6 population 1427.05 3.0 1166.0 35682.0 0 Float64
7 households 501.714 2.0 410.0 4769.0 0 Float64
8 median_income 3.8858 0.4999 3.5494 15.0001 0 Float64
9 rooms_per_person 1.98433 0.0616054 1.94325 34.2143 Float64
In [16]:
training_targets = preprocess_targets(head(california_housing_dataframe,12000))
describe(training_targets)
Out[16]:
variable mean min median max nunique nmissing eltype
1 median_house_value 208.244 14.999 181.3 500.001 Float64
For the validation set, we’ll choose the last 5000 examples, out of the total of 17000.
In [17]:
validation_examples = preprocess_features(tail(california_housing_dataframe,5000))
describe(validation_examples)
Out[17]:
variable mean min median max nunique nmissing eltype
1 latitude 35.5861 32.55 34.235 41.95 0 Float64
2 longitude -119.511 -124.27 -118.45 -114.58 0 Float64
3 housing_median_age 28.4004 2.0 29.0 52.0 0 Float64
4 total_rooms 2641.58 2.0 2110.0 37937.0 0 Float64
5 total_bedrooms 536.343 1.0 430.0 6445.0 0 Float64
6 population 1435.63 6.0 1168.0 28566.0 0 Float64
7 households 500.04 1.0 407.0 6082.0 0 Float64
8 median_income 3.87824 0.4999 3.5318 15.0001 0 Float64
9 rooms_per_person 1.97284 0.0180649 1.93763 55.2222 Float64
In [18]:
validation_targets = preprocess_targets(tail(california_housing_dataframe,5000))
describe(validation_targets)
Out[18]:
variable mean min median max nunique nmissing eltype
1 median_house_value 205.038 14.999 177.85 500.001 Float64

Task 1: Examine the Data

Okay, let’s look at the data above. We have 9 input features that we can use.
Take a quick skim over the table of values. Everything look okay? See how many issues you can spot. Don’t worry if you don’t have a background in statistics; common sense will get you far.
After you’ve had a chance to look over the data yourself, check the solution for some additional thoughts on how to verify data.

Solution

Let’s check our data against some baseline expectations:
  • For some values, like median_house_value, we can check to see if these values fall within reasonable ranges (keeping in mind this was 1990 data — not today!).
  • For other values, like latitude and longitude, we can do a quick check to see if these line up with expected values from a quick Google search.
If you look closely, you may see some oddities:
  • median_income is on a scale from about 3 to 15. It’s not at all clear what this scale refers to—looks like maybe some log scale? It’s not documented anywhere; all we can assume is that higher values correspond to higher income.
  • The maximum median_house_value is 500,001. This looks like an artificial cap of some kind.
  • Our rooms_per_person feature is generally on a sane scale, with a 75th percentile value of about 2. But there are some very large values, like 18 or 55, which may show some amount of corruption in the data.
We’ll use these features as given for now. But hopefully these kinds of examples can help to build a little intuition about how to check data that comes to you from an unknown source.

Task 2: Plot Latitude/Longitude vs. Median House Value

Let’s take a close look at two features in particular: latitude and longitude. These are geographical coordinates of the city block in question.
This might make a nice visualization — let’s plot latitude and longitude, and use color to show the median_house_value.
In [30]:
ax1=scatter(validation_examples[:longitude],
validation_examples[:latitude],
color=:coolwarm,
zcolor=validation_targets[:median_house_value] ./ maximum(validation_targets[:median_house_value]),
ms=5,
markerstrokecolor=false,
title="Validation Data",
ylim=[32,43],
xlim=[-126,-112])

ax2=scatter(training_examples[:longitude],
training_examples[:latitude],
color=:coolwarm,
zcolor=training_targets[:median_house_value] ./ maximum(training_targets[:median_house_value]),
markerstrokecolor=false,
ms=5,
title="Training Data",
ylim=[32,43],
xlim=[-126,-112]);

#plot(ax1, ax2, legend=false, colorbar=false, layout=(1,2))
Wait a second…this should have given us a nice map of the state of California, with red showing up in expensive areas like the San Francisco and Los Angeles.
The training set sort of does, compared to a real map, but the validation set clearly doesn’t.
Go back up and look at the data from Task 1 again.
Do you see any other differences in the distributions of features or targets between the training and validation data?

Solution

Looking at the tables of summary stats above, it’s easy to wonder how anyone would do a useful data check. What’s the right 75th percentile value for total_rooms per city block?
The key thing to notice is that for any given feature or column, the distribution of values between the train and validation splits should be roughly equal.
The fact that this is not the case is a real worry, and shows that we likely have a fault in the way that our train and validation split was created.

Task 3: Return to the Data Importing and Pre-Processing Code, and See if You Spot Any Bugs

If you do, go ahead and fix the bug. Don’t spend more than a minute or two looking. If you can’t find the bug, check the solution.
When you’ve found and fixed the issue, re-run latitude / longitude plotting cell above and confirm that our sanity checks look better.
By the way, there’s an important lesson here.
Debugging in ML is often data debugging rather than code debugging.
If the data is wrong, even the most advanced ML code can’t save things.

Solution

Take a look at how the data is randomized when it’s read in.
If we don’t randomize the data properly before creating training and validation splits, then we may be in trouble if the data is given to us in some sorted order, which appears to be the case here.

Task 4: Train and Evaluate a Model

Spend 5 minutes or so trying different hyperparameter settings. Try to get the best validation performance you can.
Next, we’ll train a linear regressor using all the features in the data set, and see how well we do.
Let’s define the same input function we’ve used previously for loading the data into a TensorFlow model.
In [19]:
function create_batches(features, targets, steps, batch_size=5, num_epochs=0)

if(num_epochs==0)
num_epochs=ceil(batch_size*steps/size(features,1))
end

names_features=names(features);
names_targets=names(targets);

features_batches=copy(features)
target_batches=copy(targets)

for i=1:num_epochs

select=shuffle(1:size(features,1))

if i==1
features_batches=(features[select,:])
target_batches=(targets[select,:])
else

append!(features_batches, features[select,:])
append!(target_batches, targets[select,:])
end
end

return features_batches, target_batches
end
Out[19]:
create_batches (generic function with 3 methods)
In [20]:
function next_batch(features_batches, targets_batches, batch_size, iter)

select=mod((iter-1)*batch_size+1, size(features_batches,1)):mod(iter*batch_size, size(features_batches,1));

ds=features_batches[select,:];
target=targets_batches[select,:];

return ds, target
end
Out[20]:
next_batch (generic function with 1 method)
In [21]:
function my_input_fn(features_batches, targets_batches, iter, batch_size=5, shuffle_flag=1):
"""Trains a linear regression model of one feature.

Args:
features: DataFrame of features
targets: DataFrame of targets
batch_size: Size of batches to be passed to the model
shuffle: True or False. Whether to shuffle the data.
num_epochs: Number of epochs for which data should be repeated. None = repeat indefinitely
Returns:
Tuple of (features, labels) for next data batch
"""

# Convert pandas data into a dict of np arrays.
#features = {key:np.array(value) for key,value in dict(features).items()}

# Construct a dataset, and configure batching/repeating.
#ds = Dataset.from_tensor_slices((features,targets)) # warning: 2GB limit
ds, target = next_batch(features_batches, targets_batches, batch_size, iter)

# Shuffle the data, if specified.
if shuffle_flag==1
select=shuffle(1:size(ds, 1));
ds = ds[select,:]
target = target[select, :]
end

# Return the next batch of data.
# features, labels = ds.make_one_shot_iterator().get_next()
return ds, target
end
Out[21]:
my_input_fn (generic function with 3 methods)
Because we’re now working with multiple input features, let’s modularize our code for configuring feature columns into a separate function. (For now, this code is fairly simple, as all our features are numeric, but we’ll build on this code as we use other types of features in future exercises.)
In [23]:
function construct_columns(input_features)
"""Construct the TensorFlow Feature Columns.

Args:
input_features: A dataframe of numerical input features to use.
Returns:
A set of feature columns
"""
out=convert(Array, input_features[:,:])
return convert.(Float64,out)

end
Out[23]:
construct_columns (generic function with 1 method)
Next, we use the train_model() code below to set up the input functions and calculate predictions.
Compare the losses on training data and validation data. With a single raw feature, our best root mean squared error (RMSE) was of about 180.
See how much better you can do now that we can use multiple features.
Check the data using some of the methods we’ve looked at before. These might include:
  • Comparing distributions of predictions and actual target values
  • Creating a scatter plot of predictions vs. target values
  • Creating two scatter plots of validation data using latitude and longitude:
    • One plot mapping color to actual target median_house_value
    • A second plot mapping color to predicted median_house_value for side-by-side comparison.
In [24]:
function train_model(learning_rate,
steps,
batch_size,
training_examples,
training_targets,
validation_examples,
validation_targets)
"""Trains a linear regression model of one feature.

Args:
learning_rate: A `float`, the learning rate.
steps: A non-zero `int`, the total number of training steps. A training step
consists of a forward and backward pass using a single batch.
batch_size: A non-zero `int`, the batch size.
training_examples: A dataframe of training examples.
training_targets: A column of training targets.
validation_examples: A dataframe of validation examples.
validation_targets: A column of validation targets.
"""

periods = 10
steps_per_period = steps / periods

# Create feature columns.
feature_columns = placeholder(Float32)
target_columns = placeholder(Float32)

# Create a linear regressor object.
m=Variable(zeros(length(training_examples),1))
b=Variable(0.0)
y=(feature_columns*m) .+ b
loss=reduce_sum((target_columns - y).^2)
run(sess, global_variables_initializer())
features_batches, targets_batches = create_batches(training_examples, training_targets, steps, batch_size)

# Advanced gradient decent with gradient clipping
my_optimizer=(train.GradientDescentOptimizer(learning_rate))
gvs = train.compute_gradients(my_optimizer, loss)
capped_gvs = [(clip_by_norm(grad, 5.), var) for (grad, var) in gvs]
my_optimizer = train.apply_gradients(my_optimizer,capped_gvs)


# Train the model, but do so inside a loop so that we can periodically assess
# loss metrics.
println("Training model...")
println("RMSE (on training data):")
training_rmse = []
validation_rmse=[]
for period in 1:periods
# Train the model, starting from the prior state.
for i=1:steps_per_period
features, labels = my_input_fn(features_batches, targets_batches, convert(Int,(period-1)*steps_per_period+i), batch_size)
run(sess, my_optimizer, Dict(feature_columns=>construct_columns(features), target_columns=>construct_columns(labels)))
end
# Take a break and compute predictions.
training_predictions = run(sess, y, Dict(feature_columns=> construct_columns(training_examples)));
validation_predictions = run(sess, y, Dict(feature_columns=> construct_columns(validation_examples)));

# Compute loss.
training_mean_squared_error = mean((training_predictions- construct_columns(training_targets)).^2)
training_root_mean_squared_error = sqrt(training_mean_squared_error)
validation_mean_squared_error = mean((validation_predictions- construct_columns(validation_targets)).^2)
validation_root_mean_squared_error = sqrt(validation_mean_squared_error)
# Occasionally print the current loss.
println(" period ", period, ": ", training_root_mean_squared_error)
# Add the loss metrics from this period to our list.
push!(training_rmse, training_root_mean_squared_error)
push!(validation_rmse, validation_root_mean_squared_error)
end

weight = run(sess,m)
bias = run(sess,b)
println("Model training finished.")

# Output a graph of loss metrics over periods.
p1=plot(training_rmse, label="training", title="Root Mean Squared Error vs. Periods", ylabel="RMSE", xlabel="Periods")
p1=plot!(validation_rmse, label="validation")

println("Final RMSE (on training data): ", training_rmse[end])
println("Final Weight (on training data): ", weight)
println("Final Bias (on training data): ", bias)

return weight, bias, p1
end
Out[24]:
train_model (generic function with 1 method)
In [25]:
weight, bias, p1 = train_model(
# TWEAK THESE VALUES TO SEE HOW MUCH YOU CAN IMPROVE THE RMSE
0.00003, #learning rate
500, #steps
5, #batch_size
training_examples,
training_targets,
validation_examples,
validation_targets)
Training model...
RMSE (on training data):
period 1: 218.21101557986623
period 2: 200.39219050211705
period 3: 187.48228248649704
period 4: 177.86646056587998
period 5: 171.31757059486895
period 6: 167.42319001197586
period 7: 166.09887670830182
period 8: 165.48684651754442
period 9: 165.77122987589004
period 10: 166.47520437942347
Model training finished.
Final RMSE (on training data): 166.47520437942347
Final Weight (on training data):
Out[25]:
([0.00133516; -0.0045199; … ; 0.000193309; 8.00184e-5], 0.06642270821680482, Plot{Plots.GRBackend() n=2})
[0.00133516; -0.0045199; 0.0012989; 0.0423281; 0.00791081; 0.0207483; 0.0079124; 0.000193309; 8.00184e-5]
Final Bias (on training data): 0.06642270821680482
In [26]:
plot(p1)
Out[26]:
246810170180190200210Root Mean Squared Error vs. PeriodsPeriodsRMSEtrainingvalidation

Task 5: Evaluate on Test Data

In the cell below, load in the test data set and evaluate your model on it.
We’ve done a lot of iteration on our validation data. Let’s make sure we haven’t overfit to the pecularities of that particular sample. The test data set is located here.
How does your test performance compare to the validation performance? What does this say about the generalization performance of your model?
In [27]:
california_housing_test_data  = CSV.read("california_housing_test.csv", delim=",");

test_examples = preprocess_features(california_housing_test_data)
test_targets = preprocess_targets(california_housing_test_data)

test_predictions = construct_columns(test_examples)*weight .+ bias

test_mean_squared_error = mean((test_predictions- construct_columns(test_targets)).^2)
test_root_mean_squared_error = sqrt(test_mean_squared_error)

print("Final RMSE (on test data): ", test_root_mean_squared_error)
Final RMSE (on test data): 161.49519916004172
In [28]:
# end of file