diff --git a/.travis.yml b/.travis.yml index 53574a9..8a39511 100644 --- a/.travis.yml +++ b/.travis.yml @@ -5,9 +5,12 @@ notifications: email: false before_install: - sudo add-apt-repository ppa:staticfloat/julia-deps -y - - sudo add-apt-repository ppa:staticfloat/julianightlies -y + - sudo add-apt-repository ppa:staticfloat/juliareleases -y - sudo apt-get update -qq -y - sudo apt-get install libpcre3-dev julia -y script: - julia -e 'Pkg.init(); run(`ln -s $(pwd()) $(Pkg.dir("LinearLeastSquares"))`); Pkg.pin("LinearLeastSquares"); Pkg.resolve()' - julia -e 'using LinearLeastSquares; @assert isdefined(:LinearLeastSquares); @assert typeof(LinearLeastSquares) === Module' + - julia --code-coverage test/runtests.jl +after_success: + - julia -e 'cd(Pkg.dir("LinearLeastSquares")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())' diff --git a/REQUIRE b/REQUIRE index 93db550..e3edc9e 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,2 +1,2 @@ -julia 0.2+ +julia 0.3 diff --git a/examples/regression.jl b/examples/regression.jl deleted file mode 100644 index aeb26a1..0000000 --- a/examples/regression.jl +++ /dev/null @@ -1,45 +0,0 @@ -# Regression problem -# Originally written by Keegan Go for lsqpy -# Translated into LinearLeastSquares.jl by Karanveer Mohan and David Zeng - -""" -Fit a large polynomial to some given data -Play with regularization to understand -Question: Why does increasing regularization not help very much -near the end of the polynomial (farther from zero)? -""" - -using LinearLeastSquares -import PyPlot.plt - -# Set the random seed to get consistent data -srand(1); - -# Number of examples to use -n = 50 - -# Generate data -x_data = rand(n, 1) * 18; -y_data = sin(x_data * 2) + cos(x_data)- 4 * cos(x_data) / 4 + 3 * sin(x_data / 3) + 0.2 * rand(n, 1); - -num_powers = 20; - -t_vals = 0 : 0.1 : 18; -t = reshape(t_vals, length(t_vals), 1); -T = hcat([t .^ i for i in 1 : num_powers]...); - -# We will regress using different powers of x -X = hcat([x_data .^ i for i in 1 : num_powers]...); - -# Solve the problem -mu = 0; -a = Variable(num_powers); -optval = minimize!(sum_squares(X * a - y_data) + mu * sum_squares(a)); - -# Plot our regressed function -plt.plot(x_data, y_data, "bo"); -plt.plot(t, T * a.value, "r"); -plt.xlabel("x"); -plt.ylabel("y"); -plt.ylim([-5, 5]); -plt.show(); diff --git a/examples/regression/data.jl b/examples/regression/data.jl new file mode 100644 index 0000000..b9b73de --- /dev/null +++ b/examples/regression/data.jl @@ -0,0 +1,22 @@ +using Gadfly + +# Set the random seed to get consistent data +srand(1) + +# Number of examples to use +n = 100 + +# Specify the true value of the variable +true_coeffs = [2; -2; 0.5]; + +# Generate data +x_data = rand(n, 1) * 5; +x_data_expanded = hcat([x_data .^ i for i in 1 : 3]...) +y_data = x_data_expanded * true_coeffs + 0.5 * rand(n, 1) + + +p = plot( + x=x_data, y=y_data, Geom.point, + Theme(panel_fill=color("white")) +) +draw(PNG("data.png", 16cm, 12cm), p) diff --git a/examples/regression/data.png b/examples/regression/data.png new file mode 100644 index 0000000..7ae7d09 Binary files /dev/null and b/examples/regression/data.png differ diff --git a/examples/regression/linear_regression.jl b/examples/regression/linear_regression.jl new file mode 100644 index 0000000..7a26a5e --- /dev/null +++ b/examples/regression/linear_regression.jl @@ -0,0 +1,35 @@ +# Simple Linear Regression +# Originally written by Keegan Go for lsqpy +# Translated into LinearLeastSquares.jl by Karanveer Mohan and David Zeng + +using LinearLeastSquares +using Gadfly +# Set the random seed to get consistent data +srand(1) + +# Number of examples to use +n = 100 + +# Specify the true value of the variable +true_coeffs = [2; -2; 0.5] + +# Generate data +x_data = rand(n) * 5 +x_data_expanded = hcat([x_data .^ i for i in 1 : 3]...) +y_data = x_data_expanded * true_coeffs + 0.5 * rand(n) + +slope = Variable() +offset = Variable() +line = offset + x_data * slope +residuals = line - y_data +error = sum_squares(residuals) +optval = minimize!(error) + +# plot the data and the line +t = [0; 5; 0.1] +p = plot( + layer(x=x_data, y=y_data, Geom.point), + layer(x=t, y=evaluate(slope) * t + evaluate(offset), Geom.line), + Theme(panel_fill=color("white")) +) +# draw(PNG("linear_regression.png", 16cm, 12cm), p) diff --git a/examples/regression/linear_regression.png b/examples/regression/linear_regression.png new file mode 100644 index 0000000..f20a214 Binary files /dev/null and b/examples/regression/linear_regression.png differ diff --git a/examples/simple_lin_and_quad_reg/quadratic_regression.jl b/examples/regression/quadratic_regression.jl similarity index 54% rename from examples/simple_lin_and_quad_reg/quadratic_regression.jl rename to examples/regression/quadratic_regression.jl index a5ea0e5..0dbe8f1 100644 --- a/examples/simple_lin_and_quad_reg/quadratic_regression.jl +++ b/examples/regression/quadratic_regression.jl @@ -3,6 +3,7 @@ # Translated into LinearLeastSquares.jl by Karanveer Mohan and David Zeng using LinearLeastSquares +using Gadfly # Set the random seed to get consistent data srand(1) @@ -10,28 +11,29 @@ srand(1) n = 100 # Specify the true value of the variable -true_coeffs = [2; -2; 0.5]; +true_coeffs = [2; -2; 0.5] # Generate data x_data = rand(n, 1) * 5 x_data_expanded = hcat([x_data .^ i for i in 1 : 3]...) y_data = x_data_expanded * true_coeffs + 0.5 * rand(n, 1) -quadratic_coeff = Variable(); -slope = Variable(); -offset = Variable(); -quadratic = offset + x_data * slope + quadratic_coeff * x_data .^ 2; -residuals = quadratic - y_data; -error = sum_squares(residuals); -optval = minimize!(error); +quadratic_coeff = Variable() +slope = Variable() +offset = Variable() +quadratic = offset + x_data * slope + quadratic_coeff * x_data .^ 2 +residuals = quadratic - y_data +error = sum_squares(residuals) +optval = minimize!(error) # Create some evenly spaced points for plotting, again replicate powers -t = reshape([0 : 0.1 : 5], length([0 : 0.1 : 5]), 1); -t_squared = t .^ 2; +t = reshape([0 : 0.1 : 5], length([0 : 0.1 : 5]), 1) +t_squared = t .^ 2 # Plot our regressed function -plt.plot(x_data, y_data, "ro") -plt.plot(t, evaluate(offset) + t * evaluate(slope) + t_squared * evaluate(quadratic_coeff), "b") -plt.xlabel("x") -plt.ylabel("y") -plt.show() +p = plot( + layer(x=x_data, y=y_data, Geom.point), + layer(x=t, y=evaluate(offset) + t * evaluate(slope) + t_squared * evaluate(quadratic_coeff), Geom.line), + Theme(panel_fill=color("white")) +) +# draw(PNG("quadratic_regression.png", 16cm, 12cm), p) diff --git a/examples/regression/quadratic_regression.png b/examples/regression/quadratic_regression.png new file mode 100644 index 0000000..5f8a561 Binary files /dev/null and b/examples/regression/quadratic_regression.png differ diff --git a/examples/simple_classifier.jl b/examples/simple_classifier.jl deleted file mode 100644 index 501b9cd..0000000 --- a/examples/simple_classifier.jl +++ /dev/null @@ -1,34 +0,0 @@ -# Simple Classifier -# Originally written by Keegan Go for lsqpy -# Translated into LinearLeastSquares.jl by Karanveer Mohan and David Zeng - -using LinearLeastSquares -import PyPlot.plt - -# Set the random seed to get consistent data -srand(1); - -# Number of examples to use -n = 200; - -# Specify the true value of the variable -true_vect = [-1; 1]; - -# Create data and labels -X = rand(n, 2) * 2; -y = sign(X * true_vect + rand(n, 1)); - -a = Variable(2); -optval = minimize!(sum_squares(X * a - y)); -println(a.value); - -# Plot the line we found -plt.plot(X[y[:, 1] .>= 0, 1], X[y[:, 1] .>= 0, 2], "ro"); -plt.plot(X[y[:, 1] .< 0, 1], X[y[:, 1] .< 0, 2], "bo"); - -t = [-5 : 0.1 : 5]; -plt.plot(t, -a.value[1, 1] * t / a.value[2, 1]); -plt.xlabel("x"); -plt.ylabel("y"); -plt.xlim([0, 2]); -plt.ylim([0, 2]); diff --git a/examples/simple_lin_and_quad_reg/data.jl b/examples/simple_lin_and_quad_reg/data.jl deleted file mode 100644 index 174f83c..0000000 --- a/examples/simple_lin_and_quad_reg/data.jl +++ /dev/null @@ -1,18 +0,0 @@ -import PyPlot.plt - -# Set the random seed to get consistent data -srand(1) - -# Number of examples to use -n = 100 - -# Specify the true value of the variable -true_coeffs = [2; -2; 0.5]; - -# Generate data -x_data = rand(n, 1) * 5; -x_data_expanded = hcat([x_data .^ i for i in 1 : 3]...); -y_data = x_data_expanded * true_coeffs + 0.5 * rand(n, 1); - -plt.plot(x_data, y_data, "ro"); -plt.show() diff --git a/examples/simple_lin_and_quad_reg/data.png b/examples/simple_lin_and_quad_reg/data.png deleted file mode 100644 index 1f124f0..0000000 Binary files a/examples/simple_lin_and_quad_reg/data.png and /dev/null differ diff --git a/examples/simple_lin_and_quad_reg/linear.png b/examples/simple_lin_and_quad_reg/linear.png deleted file mode 100644 index d81f398..0000000 Binary files a/examples/simple_lin_and_quad_reg/linear.png and /dev/null differ diff --git a/examples/simple_lin_and_quad_reg/linear_regression.jl b/examples/simple_lin_and_quad_reg/linear_regression.jl deleted file mode 100644 index 2507014..0000000 --- a/examples/simple_lin_and_quad_reg/linear_regression.jl +++ /dev/null @@ -1,34 +0,0 @@ -# Simple Linear Regression -# Originally written by Keegan Go for lsqpy -# Translated into LinearLeastSquares.jl by Karanveer Mohan and David Zeng - -using LinearLeastSquares -import PyPlot.plt -# Set the random seed to get consistent data -srand(1) - -# Number of examples to use -n = 100 - -# Specify the true value of the variable -true_coeffs = [2; -2; 0.5]; - -# Generate data -x_data = rand(n, 1) * 5; -x_data_expanded = hcat([x_data .^ i for i in 1 : 3]...); -y_data = x_data_expanded * true_coeffs + 0.5 * rand(n, 1); - -slope = Variable(); -offset = Variable(); -line = offset + x_data * slope; -residuals = line - y_data; -error = sum_squares(residuals); -optval = minimize!(error); - -# plot the data and the line -t = [0; 5; 0.1]; -plt.plot(x_data, y_data, "ro"); -plt.plot(t, evaluate(slope) .* t .+ evaluate(offset)); -plt.xlabel("x"); -plt.ylabel("y"); -plt.show(); diff --git a/examples/simple_lin_and_quad_reg/quadratic.png b/examples/simple_lin_and_quad_reg/quadratic.png deleted file mode 100644 index 01eb799..0000000 Binary files a/examples/simple_lin_and_quad_reg/quadratic.png and /dev/null differ diff --git a/examples/temperature/CRUTEM4v-gl.dat b/examples/temperature/CRUTEM4v-gl.dat deleted file mode 100644 index 6c2b99e..0000000 --- a/examples/temperature/CRUTEM4v-gl.dat +++ /dev/null @@ -1,328 +0,0 @@ -1851 0.757 0.333 -0.511 0.464 0.056 0.079 0.342 0.277 0.196 0.643 -0.109 -0.181 0.196 -1851 3 3 3 3 3 3 3 3 3 3 3 3 -1852 0.166 0.186 -0.042 -0.813 0.409 0.459 0.517 0.095 0.166 -0.201 0.126 1.686 0.230 -1852 3 3 3 3 3 3 3 3 3 3 3 3 -1853 0.675 -0.596 -1.249 -0.460 -0.194 -0.174 0.362 0.282 -0.317 0.062 -0.500 -1.150 -0.272 -1853 3 3 3 3 3 3 3 3 3 3 3 3 -1854 -0.427 0.207 0.175 -0.330 0.383 0.082 0.552 0.334 0.131 0.617 -0.487 0.310 0.129 -1854 3 3 3 3 3 3 3 3 3 3 3 3 -1855 0.020 -1.295 -0.523 0.379 0.180 0.142 0.221 0.083 -0.301 0.287 -0.431 -1.540 -0.231 -1855 3 3 3 3 3 3 3 3 3 3 3 3 -1856 -0.277 -0.713 -1.209 -0.061 -0.574 0.183 -0.395 -0.300 -0.575 -0.572 -1.397 -0.153 -0.504 -1856 4 4 3 3 3 3 3 3 3 3 3 3 -1857 -1.081 0.203 -0.685 -1.363 -0.964 -0.320 -0.036 0.101 -0.329 -0.398 -0.717 0.560 -0.419 -1857 3 4 3 4 4 4 4 4 4 4 4 4 -1858 -0.003 -1.310 -0.723 -0.395 -0.377 0.111 -0.168 -0.348 -0.351 -0.080 -1.719 -0.371 -0.478 -1858 4 4 4 4 4 4 4 4 4 4 4 4 -1859 -0.156 0.148 -0.034 -0.176 -0.021 -0.187 0.018 0.008 -0.535 -0.180 -0.316 -1.036 -0.206 -1859 4 4 4 4 4 4 4 4 4 4 4 4 -1860 0.106 -1.067 -1.486 -0.533 -0.296 0.001 -0.280 -0.068 -0.282 -0.348 -1.218 -1.536 -0.584 -1860 4 4 4 4 4 4 4 4 4 4 4 4 -1861 -1.899 0.202 0.048 -0.774 -0.980 0.113 -0.141 0.029 -0.323 -0.277 -0.691 -0.438 -0.428 -1861 4 4 4 4 4 4 4 4 4 4 4 4 -1862 -1.221 -1.272 -0.444 -0.318 -0.180 -0.395 -0.457 -0.702 -0.367 -0.517 -0.972 -1.115 -0.663 -1862 4 4 4 3 3 3 4 4 4 4 4 4 -1863 1.068 0.418 -0.369 -0.308 -0.230 -0.381 -0.471 -0.255 -0.470 -0.578 -0.367 -0.363 -0.192 -1863 3 4 4 4 4 4 4 3 3 3 3 3 -1864 -1.350 -0.583 -0.494 -0.738 -0.621 0.097 -0.077 -0.477 -0.304 -1.297 -1.263 -1.225 -0.694 -1864 4 4 4 4 4 4 4 4 4 4 4 4 -1865 -0.131 -1.306 -1.385 -0.179 -0.168 -0.211 0.106 -0.351 0.236 -0.595 0.023 -0.542 -0.375 -1865 4 4 4 4 4 4 4 4 4 4 4 4 -1866 0.547 -0.509 -0.897 -0.226 -0.729 0.034 -0.037 -0.506 -0.042 -0.818 -0.628 -0.374 -0.349 -1866 4 4 4 4 4 4 4 4 4 4 4 4 -1867 -0.711 0.560 -1.827 -0.492 -1.138 -0.270 -0.341 -0.184 -0.104 -0.109 -0.631 -0.962 -0.517 -1867 4 4 4 4 4 4 4 4 4 4 4 4 -1868 -1.348 -0.833 -0.208 -0.664 0.159 0.089 0.423 0.098 -0.354 -0.153 -0.980 -0.203 -0.331 -1868 4 4 4 5 5 4 5 5 5 4 5 5 -1869 -0.166 0.827 -1.125 -0.352 -0.444 -0.345 -0.167 0.192 -0.071 -0.972 -0.613 -0.415 -0.304 -1869 5 5 5 5 5 5 4 4 5 5 4 4 -1870 0.103 -0.959 -1.112 -0.446 -0.018 0.185 0.123 -0.704 -0.318 -0.482 -0.149 -1.758 -0.461 -1870 5 5 5 5 5 5 5 5 5 5 5 5 -1871 -0.997 -1.155 -0.007 -0.335 -0.556 -0.396 0.008 0.110 -0.593 -0.604 -1.268 -1.203 -0.583 -1871 5 5 5 5 5 5 5 5 5 5 5 5 -1872 -0.124 -0.498 -0.624 -0.124 -0.044 -0.089 -0.132 -0.422 -0.194 -0.177 -0.504 -0.789 -0.310 -1872 5 5 6 5 6 6 5 5 5 5 6 5 -1873 -0.031 -0.629 -0.530 -0.988 -0.587 0.203 -0.155 0.004 -0.398 -0.422 -0.918 -0.020 -0.373 -1873 6 6 6 6 6 6 6 6 6 6 6 6 -1874 0.256 -0.469 -1.071 -0.670 -0.840 -0.404 -0.248 -0.388 -0.201 -0.250 -0.806 -0.388 -0.457 -1874 6 6 6 6 6 6 6 6 6 6 6 6 -1875 -1.212 -1.581 -1.624 -1.047 -0.254 -0.195 -0.324 -0.298 -0.640 -0.814 -1.263 -0.888 -0.845 -1875 6 7 7 6 7 7 7 7 7 7 7 7 -1876 -0.270 -0.228 -0.458 -0.374 -0.745 -0.082 0.090 -0.172 -0.216 -0.446 -1.265 -1.450 -0.468 -1876 7 7 7 7 7 7 7 7 7 7 7 7 -1877 -0.317 -0.046 -0.637 -0.739 -0.885 -0.069 0.146 0.185 -0.384 -0.561 0.058 0.098 -0.263 -1877 7 7 7 7 7 7 7 7 7 7 7 7 -1878 -0.342 0.404 0.598 0.369 -0.435 -0.114 -0.095 -0.026 0.240 0.107 -0.155 -0.965 -0.034 -1878 7 7 7 7 7 7 7 7 7 7 7 7 -1879 -0.483 -0.482 -0.455 -0.625 -0.527 -0.587 -0.377 -0.368 -0.346 -0.124 -0.876 -1.181 -0.536 -1879 8 8 8 8 8 8 8 8 8 8 8 8 -1880 0.086 -0.266 -0.499 -0.343 -0.110 -0.128 -0.201 0.121 -0.181 -0.975 -1.041 -0.616 -0.346 -1880 8 8 8 8 8 8 8 8 8 8 8 8 -1881 -1.244 -0.620 -0.721 -0.477 -0.054 -0.602 -0.294 -0.163 -0.148 -0.666 -0.689 0.092 -0.465 -1881 8 8 8 9 9 9 9 8 8 9 8 8 -1882 0.598 0.080 0.099 -0.693 -0.720 -0.664 -0.538 -0.118 -0.122 -0.391 -0.758 -0.894 -0.343 -1882 9 9 9 9 9 9 9 9 9 9 9 9 -1883 -1.072 -0.913 -1.105 -0.710 -0.566 0.011 -0.279 -0.449 -0.548 -0.561 -0.573 -0.362 -0.594 -1883 9 9 9 9 9 9 9 9 9 9 9 9 -1884 -0.601 -0.415 -0.808 -1.056 -0.993 -0.684 -0.717 -0.269 -0.421 -0.515 -0.996 -0.637 -0.676 -1884 9 9 9 9 9 9 9 9 9 9 9 9 -1885 -1.315 -0.794 -1.034 -0.883 -0.882 -0.703 -0.367 -0.633 -0.373 -0.526 -0.462 -0.012 -0.665 -1885 9 9 9 9 9 9 9 9 9 9 10 10 -1886 -0.785 -1.042 -0.964 -0.311 -0.413 -0.647 -0.337 -0.399 -0.225 -0.567 -0.494 -0.302 -0.541 -1886 10 10 10 10 10 10 10 10 10 10 10 10 -1887 -0.712 -0.779 -0.577 -0.658 -0.434 -0.425 -0.299 -0.464 -0.315 -0.879 -0.669 -0.610 -0.569 -1887 10 10 10 10 10 10 10 10 10 10 10 10 -1888 -1.170 -1.003 -1.434 -0.344 -0.747 -0.469 -0.417 -0.452 -0.241 -0.416 -0.498 -0.318 -0.626 -1888 10 11 10 11 11 10 10 10 10 10 11 11 -1889 -0.446 -0.458 -0.304 -0.012 0.108 -0.171 -0.368 -0.400 -0.566 -0.353 -0.693 0.024 -0.303 -1889 11 11 11 11 11 11 11 11 11 11 11 11 -1890 -0.152 -0.372 -0.509 -0.270 -0.617 -0.299 -0.433 -0.440 -0.263 -0.531 -0.801 -0.705 -0.449 -1890 11 11 11 11 11 11 11 11 11 11 11 11 -1891 -0.935 -0.834 -0.691 -0.721 -0.474 -0.502 -0.519 -0.529 -0.232 -0.513 -0.922 0.160 -0.559 -1891 12 12 12 12 12 12 12 12 12 12 12 11 -1892 -0.726 -0.082 -0.884 -0.865 -0.805 -0.436 -0.494 -0.361 -0.280 -0.478 -0.784 -1.325 -0.627 -1892 12 12 12 12 12 12 12 12 12 12 12 12 -1893 -2.046 -1.506 -0.603 -0.677 -0.759 -0.487 -0.205 -0.390 -0.496 -0.265 -0.535 -0.387 -0.696 -1893 12 12 12 12 12 12 12 12 12 12 12 12 -1894 -0.552 -0.457 -0.212 -0.419 -0.538 -0.591 -0.425 -0.391 -0.741 -0.439 -0.443 -0.411 -0.468 -1894 12 12 12 12 12 12 12 12 12 12 12 12 -1895 -1.076 -1.525 -0.761 -0.458 -0.488 -0.377 -0.652 -0.356 -0.183 -0.345 -0.459 -0.277 -0.580 -1895 12 12 12 12 12 12 12 13 12 13 13 13 -1896 -0.348 -0.353 -0.842 -0.797 -0.371 -0.256 -0.244 -0.342 -0.313 -0.055 -0.950 -0.142 -0.418 -1896 13 13 13 13 13 13 13 13 13 13 13 13 -1897 -0.599 -0.375 -0.802 -0.089 -0.059 -0.172 -0.131 -0.273 -0.063 -0.294 -0.502 -0.580 -0.328 -1897 13 13 13 13 13 13 13 13 13 13 13 13 -1898 0.324 -0.394 -1.290 -0.659 -0.613 -0.213 -0.330 -0.134 -0.134 -0.570 -0.416 -0.135 -0.380 -1898 14 14 14 14 14 14 14 14 14 14 14 14 -1899 -0.124 -0.730 -0.809 -0.245 -0.466 -0.554 -0.266 -0.224 -0.005 -0.199 0.318 -0.610 -0.326 -1899 14 14 14 14 14 14 14 14 14 14 14 14 -1900 -0.397 -0.426 -0.470 -0.409 -0.279 -0.030 -0.284 -0.157 -0.236 0.288 -0.193 0.139 -0.204 -1900 14 14 14 14 14 14 14 14 14 14 14 14 -1901 -0.179 -0.545 -0.178 -0.070 -0.144 0.057 -0.020 -0.090 -0.152 -0.178 -0.329 -0.457 -0.190 -1901 14 14 14 14 14 14 14 14 14 14 14 14 -1902 0.045 -0.102 -0.311 -0.480 -0.455 -0.371 -0.445 -0.497 -0.465 -0.542 -0.604 -0.842 -0.422 -1902 15 15 15 15 15 15 14 15 15 15 15 15 -1903 -0.115 0.110 -0.120 -0.523 -0.556 -0.583 -0.527 -0.488 -0.400 -0.556 -0.591 -0.732 -0.424 -1903 15 15 15 15 15 15 15 15 15 15 15 15 -1904 -0.783 -0.696 -0.883 -0.535 -0.523 -0.505 -0.538 -0.467 -0.572 -0.389 -0.283 -0.484 -0.555 -1904 15 15 15 15 15 15 15 15 15 15 15 15 -1905 -0.568 -1.229 -0.518 -0.696 -0.397 -0.231 -0.336 -0.345 -0.250 -0.544 -0.067 -0.029 -0.434 -1905 15 15 15 15 16 15 15 15 15 15 15 15 -1906 0.100 -0.471 -0.604 0.122 -0.023 -0.134 -0.194 -0.216 -0.333 -0.252 -0.621 -0.154 -0.232 -1906 15 16 16 16 16 15 15 16 15 15 16 16 -1907 -0.700 -0.785 -0.528 -0.917 -0.967 -0.706 -0.496 -0.602 -0.420 -0.147 -0.782 -0.546 -0.633 -1907 16 16 16 16 16 16 16 16 16 16 16 16 -1908 -0.207 -0.393 -0.867 -0.567 -0.457 -0.398 -0.418 -0.576 -0.314 -0.595 -0.726 -0.410 -0.494 -1908 16 16 16 16 16 16 16 16 16 16 16 16 -1909 -0.631 -0.661 -0.872 -0.838 -0.797 -0.408 -0.437 -0.078 -0.152 -0.139 -0.091 -0.686 -0.482 -1909 16 16 16 16 16 16 16 16 17 16 16 16 -1910 -0.094 -0.377 -0.244 -0.224 -0.347 -0.312 -0.222 -0.259 -0.224 -0.413 -0.786 -0.662 -0.347 -1910 17 17 17 17 17 17 17 17 17 17 17 17 -1911 -0.587 -0.914 -0.799 -0.715 -0.410 -0.434 -0.282 -0.264 -0.296 -0.388 -0.202 -0.075 -0.447 -1911 17 18 18 18 18 17 17 17 17 18 17 18 -1912 -0.485 -0.099 -0.588 -0.442 -0.331 -0.144 -0.603 -0.693 -0.686 -0.810 -0.572 -0.275 -0.477 -1912 18 18 18 18 18 18 18 18 18 18 18 18 -1913 -0.257 -0.723 -0.638 -0.244 -0.828 -0.669 -0.379 -0.211 -0.339 -0.372 0.187 0.271 -0.350 -1913 18 18 18 18 18 18 18 18 18 18 18 18 -1914 0.454 0.138 -0.149 -0.309 -0.127 -0.115 -0.195 -0.155 -0.211 -0.000 -0.336 -0.344 -0.112 -1914 18 18 18 18 18 18 17 18 18 18 18 18 -1915 -0.088 0.302 -0.660 0.063 -0.340 -0.242 -0.034 -0.215 -0.075 -0.250 -0.034 -0.106 -0.140 -1915 17 17 17 17 17 17 17 17 17 17 17 17 -1916 0.017 -0.115 -0.726 -0.357 -0.381 -0.525 -0.227 -0.259 -0.330 -0.381 -0.466 -0.930 -0.390 -1916 18 18 18 18 18 18 17 17 18 18 18 18 -1917 -0.629 -1.340 -1.214 -0.765 -1.183 -0.440 -0.107 -0.330 -0.173 -0.609 -0.333 -1.267 -0.699 -1917 18 18 18 18 18 18 18 18 18 18 18 18 -1918 -0.788 -0.565 -0.510 -0.614 -0.736 -0.415 -0.466 -0.371 -0.306 0.121 -0.302 -0.319 -0.439 -1918 18 18 18 18 18 18 18 18 18 18 18 17 -1919 -0.080 -0.317 -0.599 -0.115 -0.467 -0.145 -0.234 -0.200 -0.072 -0.294 -0.835 -0.559 -0.326 -1919 17 17 17 17 17 17 17 17 17 17 17 17 -1920 -0.130 -0.417 -0.112 -0.346 -0.210 -0.280 -0.157 -0.225 -0.106 -0.374 -0.546 -0.495 -0.283 -1920 17 17 17 17 17 18 17 17 17 17 17 17 -1921 0.357 -0.118 -0.005 -0.033 0.053 0.064 -0.024 -0.369 -0.199 -0.211 -0.609 -0.145 -0.103 -1921 18 18 18 18 18 18 18 18 18 19 19 19 -1922 -0.631 -0.513 -0.367 -0.117 -0.165 -0.163 -0.281 -0.274 -0.202 -0.286 -0.161 -0.191 -0.279 -1922 19 19 19 19 19 19 19 19 19 19 19 19 -1923 0.071 -0.602 -0.376 -0.709 -0.377 -0.298 -0.433 -0.467 -0.189 -0.060 0.137 0.200 -0.259 -1923 19 19 19 19 19 19 19 19 19 19 20 19 -1924 -0.488 -0.492 -0.497 -0.630 -0.331 -0.207 -0.212 -0.237 -0.091 -0.166 -0.146 -0.583 -0.340 -1924 20 20 20 20 20 19 20 20 20 20 20 20 -1925 -0.327 -0.055 -0.287 -0.125 -0.395 -0.314 -0.395 -0.120 -0.245 -0.541 -0.084 0.170 -0.226 -1925 20 20 20 20 20 20 20 20 20 20 20 20 -1926 0.388 0.475 0.017 -0.372 -0.509 -0.211 -0.250 -0.077 -0.042 -0.089 0.074 -0.447 -0.087 -1926 20 20 20 20 20 20 20 20 20 20 20 20 -1927 -0.267 -0.174 -0.392 -0.429 -0.501 -0.259 -0.079 -0.171 0.020 0.215 -0.067 -0.860 -0.247 -1927 20 20 20 20 20 20 20 20 20 20 20 20 -1928 0.099 -0.113 -0.501 -0.305 -0.341 -0.569 -0.107 -0.130 -0.042 -0.069 0.104 -0.007 -0.165 -1928 20 20 20 20 21 20 20 20 20 20 20 20 -1929 -0.700 -1.366 -0.481 -0.695 -0.443 -0.431 -0.417 -0.164 -0.277 0.023 0.009 -0.699 -0.470 -1929 21 21 21 21 21 21 21 21 21 21 21 21 -1930 -0.347 -0.122 -0.057 -0.151 -0.315 -0.053 -0.007 0.037 -0.162 -0.219 0.104 -0.115 -0.117 -1930 21 21 21 21 21 21 21 21 21 21 21 21 -1931 0.024 -0.414 -0.330 -0.402 -0.342 -0.001 0.093 0.046 -0.053 0.154 -0.116 0.121 -0.102 -1931 21 21 21 22 22 22 22 22 22 22 21 21 -1932 0.878 -0.422 -0.676 0.030 -0.221 -0.195 -0.084 -0.041 0.080 0.057 -0.272 -0.147 -0.084 -1932 22 22 22 22 22 22 22 22 22 22 22 22 -1933 -0.450 -0.639 -0.541 -0.297 -0.257 -0.230 -0.055 -0.080 -0.072 -0.034 -0.365 -0.874 -0.325 -1933 22 22 23 23 23 23 23 23 23 23 23 23 -1934 -0.005 0.201 -0.399 -0.350 0.217 -0.036 0.014 -0.014 -0.141 0.035 0.219 0.122 -0.011 -1934 23 23 23 23 23 23 23 23 23 23 23 23 -1935 -0.435 0.723 -0.071 -0.563 -0.483 -0.200 -0.064 -0.040 -0.144 0.124 -0.508 -0.305 -0.164 -1935 23 23 23 23 23 23 23 23 23 23 23 23 -1936 -0.311 -0.909 -0.366 -0.326 -0.118 -0.023 0.282 0.114 -0.014 0.031 0.063 0.255 -0.110 -1936 24 24 24 24 24 24 24 24 24 24 24 24 -1937 -0.253 0.105 -0.689 -0.314 -0.029 0.085 0.002 0.215 0.245 0.322 0.112 -0.301 -0.042 -1937 24 24 24 24 24 24 24 24 24 24 24 24 -1938 0.157 0.162 0.369 0.331 0.045 -0.034 0.076 0.113 0.274 0.448 0.336 -0.357 0.160 -1938 24 24 24 24 24 24 24 24 24 24 24 24 -1939 0.128 0.058 -0.464 -0.175 0.052 0.109 -0.022 0.069 -0.139 -0.282 -0.008 0.697 0.002 -1939 25 25 25 25 25 25 25 25 25 25 25 25 -1940 -0.796 -0.111 -0.151 0.006 -0.160 0.152 0.119 0.024 0.158 -0.029 -0.038 0.279 -0.046 -1940 25 25 25 25 25 25 25 25 25 25 25 25 -1941 -0.212 0.085 -0.377 -0.036 -0.049 0.107 0.174 0.021 -0.190 0.192 -0.204 -0.164 -0.054 -1941 26 26 26 26 26 26 26 26 26 26 26 26 -1942 0.076 -0.428 -0.294 -0.111 -0.110 0.082 -0.099 -0.008 0.104 0.103 0.071 0.033 -0.049 -1942 27 27 27 27 27 27 26 27 27 27 27 27 -1943 -0.375 0.286 -0.247 0.176 0.056 -0.243 -0.003 -0.148 0.040 0.386 -0.010 0.283 0.017 -1943 27 27 27 27 27 27 27 27 27 27 27 27 -1944 0.719 0.299 0.035 -0.146 -0.042 -0.025 0.049 0.139 0.397 0.325 -0.055 -0.358 0.111 -1944 27 27 27 27 27 27 27 27 27 27 27 27 -1945 -0.256 -0.426 -0.186 0.189 -0.275 -0.118 -0.187 0.398 0.136 0.125 -0.130 -0.674 -0.117 -1945 27 27 27 27 27 27 27 27 27 27 27 27 -1946 0.121 0.134 -0.291 0.184 -0.098 -0.246 0.007 -0.026 0.053 -0.162 0.006 -0.553 -0.073 -1946 28 28 27 27 27 27 27 27 27 28 28 27 -1947 -0.293 -0.332 -0.003 0.147 -0.039 0.050 -0.056 0.076 0.060 0.413 0.072 -0.236 -0.011 -1947 28 28 28 28 28 28 28 28 28 28 28 28 -1948 0.416 -0.090 -0.399 -0.008 0.156 0.124 -0.104 -0.051 0.019 0.139 0.047 -0.223 0.002 -1948 28 28 29 29 29 29 28 28 28 29 29 29 -1949 0.366 -0.338 -0.392 -0.098 -0.110 -0.240 -0.173 -0.045 -0.084 0.125 -0.011 -0.349 -0.112 -1949 29 29 29 29 29 29 29 29 29 29 29 29 -1950 -0.626 -0.366 -0.244 -0.294 -0.151 -0.093 -0.230 -0.284 -0.145 -0.184 -0.670 -0.308 -0.300 -1950 30 30 30 31 31 30 30 31 31 31 31 30 -1951 -0.517 -0.795 -0.427 -0.078 -0.031 -0.102 -0.046 0.108 0.187 0.133 -0.033 0.393 -0.101 -1951 34 35 35 35 35 34 35 35 35 35 35 35 -1952 0.278 0.063 -0.408 -0.036 -0.054 -0.020 0.104 0.078 0.110 -0.147 -0.540 -0.208 -0.065 -1952 35 35 35 35 35 35 35 35 35 35 35 35 -1953 0.153 0.170 0.201 0.201 0.038 0.162 -0.018 0.130 0.099 0.181 -0.148 0.241 0.117 -1953 36 35 36 36 36 36 36 36 36 36 36 36 -1954 -0.494 -0.227 -0.361 -0.204 -0.297 -0.101 -0.178 -0.062 -0.004 0.071 0.231 -0.186 -0.151 -1954 36 36 36 36 36 36 36 36 36 36 36 36 -1955 0.393 -0.130 -0.773 -0.346 -0.196 -0.107 -0.200 0.033 -0.035 0.006 -0.373 -0.480 -0.184 -1955 36 37 37 37 37 37 37 37 36 37 36 36 -1956 -0.322 -0.686 -0.560 -0.518 -0.413 -0.249 -0.283 -0.371 -0.401 -0.305 -0.466 -0.359 -0.411 -1956 37 37 37 37 37 37 37 37 37 37 37 37 -1957 -0.265 -0.217 -0.449 -0.167 -0.114 0.055 -0.175 0.025 0.049 -0.066 0.130 0.345 -0.071 -1957 37 37 37 37 37 37 37 37 37 37 37 37 -1958 0.480 0.227 -0.088 0.008 0.040 -0.083 0.009 -0.040 -0.127 0.004 0.075 -0.016 0.041 -1958 37 37 37 37 37 37 37 37 37 37 37 37 -1959 0.120 0.057 0.193 0.060 -0.120 0.117 0.047 0.023 0.046 -0.164 -0.176 -0.017 0.016 -1959 38 37 38 37 38 37 38 38 38 38 37 37 -1960 0.029 0.333 -0.693 -0.292 -0.322 0.052 -0.092 -0.021 0.037 0.011 -0.213 0.299 -0.073 -1960 37 38 37 37 37 38 38 37 38 38 38 37 -1961 0.090 0.278 0.106 0.092 0.077 0.083 -0.034 -0.009 0.046 -0.050 0.006 -0.178 0.042 -1961 38 38 38 38 38 39 38 38 38 38 38 38 -1962 0.186 0.263 -0.018 0.046 -0.083 -0.145 -0.076 -0.071 -0.021 0.111 0.071 0.042 0.025 -1962 39 39 39 38 39 39 39 39 39 39 39 39 -1963 -0.086 0.388 -0.309 -0.228 -0.170 -0.135 0.059 0.042 0.132 0.384 0.266 -0.133 0.018 -1963 39 39 39 39 39 39 39 39 39 39 39 39 -1964 -0.036 -0.288 -0.515 -0.369 -0.178 -0.235 -0.194 -0.275 -0.279 -0.394 -0.282 -0.426 -0.289 -1964 39 39 39 39 39 39 39 39 39 39 39 39 -1965 -0.066 -0.443 -0.282 -0.455 -0.225 -0.106 -0.283 -0.218 -0.202 -0.062 -0.257 -0.151 -0.229 -1965 39 39 39 39 39 39 39 39 39 39 39 39 -1966 -0.230 -0.169 -0.134 -0.297 -0.193 0.032 -0.015 -0.028 -0.035 -0.180 -0.163 -0.391 -0.150 -1966 39 39 39 39 39 39 39 39 39 39 39 39 -1967 -0.276 -0.434 -0.097 -0.079 0.099 -0.242 -0.121 -0.091 -0.132 0.333 0.014 -0.124 -0.096 -1967 39 39 39 39 39 39 39 39 39 39 39 39 -1968 -0.359 -0.249 0.351 -0.203 -0.316 -0.294 -0.217 -0.203 -0.206 -0.134 -0.257 -0.441 -0.211 -1968 39 39 39 39 39 39 39 39 39 39 39 39 -1969 -0.648 -0.734 -0.331 -0.028 -0.035 -0.108 0.028 -0.033 -0.092 -0.070 0.171 0.226 -0.138 -1969 39 39 39 39 39 39 39 39 39 39 39 39 -1970 -0.068 0.279 -0.266 0.001 -0.113 0.076 -0.014 -0.078 0.003 -0.182 -0.118 -0.318 -0.067 -1970 39 39 39 39 39 39 39 39 39 39 39 39 -1971 -0.063 -0.421 -0.473 -0.365 -0.289 -0.301 -0.167 -0.108 -0.001 -0.064 0.006 -0.079 -0.194 -1971 38 38 38 38 38 38 38 38 39 38 38 38 -1972 -0.779 -0.630 -0.280 -0.168 -0.185 -0.039 -0.082 -0.083 -0.232 -0.108 -0.280 0.057 -0.234 -1972 38 38 38 38 38 38 38 38 38 38 38 38 -1973 0.150 0.505 0.362 0.245 0.140 0.237 0.069 0.021 -0.049 -0.011 -0.138 -0.062 0.122 -1973 38 38 38 38 38 38 38 38 38 38 38 38 -1974 -0.460 -0.605 -0.194 -0.155 -0.234 -0.188 -0.079 -0.120 -0.246 -0.347 -0.319 -0.326 -0.273 -1974 38 38 38 38 38 38 38 38 38 38 38 38 -1975 0.133 -0.034 0.028 0.030 0.015 -0.005 -0.018 -0.139 0.049 -0.173 -0.271 -0.261 -0.054 -1975 38 38 38 38 38 38 38 38 38 38 38 38 -1976 -0.060 -0.339 -0.718 -0.168 -0.307 -0.274 -0.297 -0.277 -0.222 -0.691 -0.447 -0.296 -0.341 -1976 38 38 38 38 38 38 38 38 38 38 38 38 -1977 -0.441 0.101 0.259 0.261 0.105 0.201 0.063 -0.009 0.018 -0.056 0.381 0.023 0.075 -1977 38 38 38 38 38 38 38 38 38 38 38 38 -1978 0.045 -0.046 0.136 0.003 -0.076 -0.165 -0.075 -0.258 -0.069 -0.060 0.149 -0.172 -0.049 -1978 38 38 38 38 38 38 38 38 38 38 38 38 -1979 -0.157 -0.301 0.019 -0.261 -0.076 0.097 0.028 0.054 0.088 0.187 0.124 0.630 0.036 -1979 38 38 38 38 38 38 38 38 38 38 38 38 -1980 0.101 0.252 -0.113 0.168 0.243 0.139 0.063 0.082 0.104 0.122 0.266 0.139 0.131 -1980 38 38 38 38 38 38 38 38 38 38 38 38 -1981 0.756 0.611 0.506 0.415 0.118 0.150 0.099 0.252 0.113 0.034 0.158 0.508 0.310 -1981 37 37 37 37 37 37 37 37 37 37 37 37 -1982 -0.200 0.102 -0.242 0.072 0.032 -0.114 0.041 0.027 0.106 -0.116 -0.093 0.461 0.006 -1982 36 36 36 36 36 36 36 37 37 37 37 37 -1983 0.744 0.490 0.368 0.138 0.142 0.107 0.246 0.333 0.379 0.216 0.560 0.139 0.322 -1983 36 37 37 37 37 37 37 37 37 37 37 37 -1984 0.179 -0.076 0.020 -0.133 0.218 0.035 0.005 0.044 -0.125 0.084 -0.203 -0.505 -0.038 -1984 37 36 37 37 37 37 37 37 36 36 37 37 -1985 -0.003 -0.309 -0.087 0.017 0.073 -0.051 -0.115 0.120 -0.020 0.048 -0.044 0.035 -0.028 -1985 37 37 37 37 37 37 37 37 37 37 37 37 -1986 0.448 0.192 0.237 0.293 0.203 0.161 -0.030 0.032 -0.053 0.011 -0.024 0.036 0.126 -1986 37 37 37 37 37 36 37 36 37 37 37 37 -1987 0.170 0.662 -0.118 0.183 0.171 0.212 0.356 0.184 0.298 0.227 0.098 0.516 0.247 -1987 37 37 37 37 37 37 37 37 37 37 37 37 -1988 0.621 0.250 0.402 0.324 0.360 0.418 0.342 0.326 0.381 0.373 0.122 0.507 0.369 -1988 37 37 37 37 37 36 36 37 36 37 37 36 -1989 0.180 0.455 0.424 0.249 0.191 0.167 0.222 0.220 0.217 0.304 0.076 0.314 0.252 -1989 36 36 36 36 36 36 36 36 36 36 36 36 -1990 0.393 0.536 1.160 0.603 0.400 0.418 0.283 0.351 0.196 0.477 0.571 0.413 0.483 -1990 36 36 36 36 36 36 36 36 36 36 36 36 -1991 0.510 0.523 0.247 0.581 0.328 0.566 0.438 0.364 0.381 0.362 0.260 0.117 0.390 -1991 35 34 34 34 33 34 33 34 34 34 33 33 -1992 0.723 0.568 0.431 0.080 0.173 0.063 -0.224 -0.058 -0.204 -0.192 -0.218 0.192 0.111 -1992 34 33 33 33 34 34 34 34 34 33 34 34 -1993 0.495 0.458 0.366 0.214 0.261 0.205 0.178 0.129 -0.109 0.111 -0.352 0.251 0.184 -1993 34 34 34 34 33 33 34 34 34 33 33 34 -1994 0.351 -0.176 0.338 0.424 0.387 0.440 0.280 0.278 0.350 0.508 0.500 0.420 0.342 -1994 34 34 34 34 34 34 34 34 34 35 35 34 -1995 0.879 1.223 0.511 0.513 0.362 0.522 0.423 0.605 0.477 0.678 0.572 0.221 0.582 -1995 34 34 35 35 34 35 34 35 34 35 34 34 -1996 0.202 0.455 0.138 0.008 0.307 0.284 0.303 0.242 0.090 0.184 0.245 0.411 0.239 -1996 34 34 34 34 34 34 34 34 34 34 34 34 -1997 0.428 0.657 0.554 0.422 0.359 0.528 0.405 0.501 0.576 0.656 0.535 0.626 0.521 -1997 34 34 34 34 34 33 34 34 33 34 34 33 -1998 0.722 1.318 0.768 0.927 0.810 0.858 0.954 0.894 0.687 0.682 0.275 0.854 0.812 -1998 34 34 34 34 34 34 34 35 34 34 34 34 -1999 0.767 1.163 0.218 0.492 0.404 0.476 0.514 0.445 0.573 0.555 0.390 0.819 0.568 -1999 34 34 34 34 34 34 34 34 34 34 34 34 -2000 0.347 0.984 0.778 0.946 0.516 0.448 0.405 0.560 0.415 0.292 0.073 0.204 0.497 -2000 34 34 34 34 34 34 34 34 34 34 33 33 -2001 0.725 0.520 0.789 0.673 0.706 0.576 0.652 0.804 0.556 0.612 1.019 0.542 0.681 -2001 31 31 32 32 31 31 31 32 32 32 31 31 -2002 1.198 1.358 1.144 0.692 0.600 0.664 0.784 0.624 0.697 0.498 0.734 0.278 0.773 -2002 32 32 33 32 32 32 32 32 33 33 33 33 -2003 0.988 0.573 0.578 0.613 0.845 0.713 0.640 0.799 0.753 0.915 0.668 1.113 0.766 -2003 33 32 32 32 32 33 32 32 32 33 32 33 -2004 0.805 1.060 0.911 0.706 0.371 0.543 0.422 0.481 0.573 0.775 1.065 0.475 0.682 -2004 33 33 33 32 32 32 32 32 32 32 32 32 -2005 0.957 0.562 0.853 1.009 0.762 0.867 0.829 0.716 0.902 1.077 1.186 0.782 0.875 -2005 32 32 32 32 32 32 32 32 32 31 32 32 -2006 0.618 0.907 0.767 0.580 0.523 0.869 0.801 0.736 0.772 0.950 0.874 1.303 0.808 -2006 31 32 32 32 32 32 33 32 32 32 33 31 -2007 1.670 0.886 1.009 1.109 0.790 0.594 0.633 0.792 0.832 0.987 0.784 0.804 0.907 -2007 32 32 32 32 32 33 33 32 33 32 32 32 -2008 0.418 0.531 1.194 0.572 0.495 0.624 0.741 0.589 0.576 1.001 1.057 0.671 0.706 -2008 32 32 32 33 33 32 32 33 33 33 33 32 -2009 0.778 0.790 0.639 0.826 0.598 0.665 0.596 0.867 0.965 0.718 0.853 0.484 0.731 -2009 32 33 32 31 32 32 33 32 32 32 32 32 -2010 0.767 0.804 1.087 1.040 0.917 0.982 0.965 0.859 0.730 0.874 1.170 0.454 0.887 -2010 33 32 33 32 32 32 33 33 33 33 33 33 -2011 0.482 0.398 0.596 0.836 0.536 0.777 0.778 0.775 0.851 0.933 0.621 0.832 0.701 -2011 31 32 31 31 31 31 31 31 31 31 31 31 -2012 0.589 0.352 0.679 1.019 0.953 0.934 0.799 0.884 0.896 0.834 0.838 0.290 0.756 -2012 31 30 29 29 29 29 29 29 30 29 29 29 -2013 0.919 0.938 0.663 0.588 0.768 0.818 0.680 0.706 0.771 0.810 0.969 0.786 0.784 -2013 28 28 28 28 29 29 28 28 28 28 29 28 -2014 0.930 0.421 0.921 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.757 -2014 27 28 28 0 0 0 0 0 0 0 0 0 diff --git a/examples/temperature/temperature.jl b/examples/temperature/temperature.jl deleted file mode 100644 index e041e46..0000000 --- a/examples/temperature/temperature.jl +++ /dev/null @@ -1,45 +0,0 @@ -""" -Reads a temperature data file as taken from -http://www.cru.uea.ac.uk/cru/data/temperature/ -""" -import PyPlot.plt -using LinearLeastSquares - -data = readdlm("CRUTEM4v-gl.dat", ' '); - -# We only want to extract the odd lines (1-indexed) -data = vcat([data[i, :] for i in 1 : 2 : size(data, 1)]...); -data = float(data); -years = data[:, 1]; -temperatures = vec(data[:, 2 : end - 1]'); - -num_samples = size(temperatures, 1) - -yearly_tread = Variable(num_samples); -a = Variable(); -b = Variable(); - - -equality_constraints = EqConstraint[]; -for i in 13 : num_samples - push!(equality_constraints, yearly_tread[i] == yearly_tread[i - 12]) -end - -t = [0 : num_samples - 1]; -t2 = t .^ 2; - - -smoothing = 0.01; -objective = sum_squares(yearly_tread + t * a + t2 * b - temperatures); -objective += smoothing * sum_squares(yearly_tread[1 : num_samples - 1] - yearly_tread[2 : num_samples]); - -optval = minimize!(objective, equality_constraints); - -plt.figure(0) -plt.plot(temperatures) -plt.plot(yearly_tread.value + t*a.value + t2*b.value, color="r") -plt.show() - -# plt.figure(1) -# plt.plot(yearly_tread.value + t*b.value + t2*a.value - temperatures) -# plt.show() diff --git a/examples/time_series/ar_fit.png b/examples/time_series/ar_fit.png new file mode 100644 index 0000000..5fc81fc Binary files /dev/null and b/examples/time_series/ar_fit.png differ diff --git a/examples/time_series/melbourne.png b/examples/time_series/melbourne.png new file mode 100644 index 0000000..c9ad5dd Binary files /dev/null and b/examples/time_series/melbourne.png differ diff --git a/examples/time_series/residuals.png b/examples/time_series/residuals.png new file mode 100644 index 0000000..4e98d7a Binary files /dev/null and b/examples/time_series/residuals.png differ diff --git a/examples/time_series/time_series.jl b/examples/time_series/time_series.jl index de50957..c3d49a8 100644 --- a/examples/time_series/time_series.jl +++ b/examples/time_series/time_series.jl @@ -1,34 +1,39 @@ using LinearLeastSquares -import PyPlot.plt -temps = readdlm("melbourne_temps.txt", ','); -plt.figure(0) -plt.plot(temps, color="b") -plt.title("Melbourne Daily Temperature") -n = size(temps, 1); -plt.xlim([0, n]) +using Gadfly -yearly = Variable(n) +temps = readdlm("melbourne_temps.txt", ',') +n = size(temps)[1] +p = plot( + x=1:1500, y=temps[1:1500], Geom.line, + Theme(panel_fill=color("white")) +) +# draw(PNG("melbourne.png", 16cm, 12cm), p) +yearly = Variable(n) eq_constraints = EqConstraint[] for i in 365 + 1 : n - eq_constraints += yearly[i] == yearly[i - 365]; + eq_constraints += yearly[i] == yearly[i - 365] end -smoothing = 100; -smooth_objective = sum_squares(yearly[1 : n - 1] - yearly[2 : n]); -optval = minimize!(sum_squares(temps - yearly) + smoothing * smooth_objective, eq_constraints); -residuals = temps - evaluate(yearly); +smoothing = 100 +smooth_objective = sum_squares(yearly[1 : n - 1] - yearly[2 : n]) +optval = minimize!(sum_squares(temps - yearly) + smoothing * smooth_objective, eq_constraints) +residuals = temps - evaluate(yearly) # Plot smooth fit -plt.figure(1) -plt.plot(temps) -plt.plot(evaluate(yearly), color="r") -plt.title("Smooth Fit of Data") -plt.xlim([0, n]) +p = plot( + layer(x=1:1500, y=evaluate(yearly)[1:1500], Geom.line, Theme(default_color=color("red"), line_width=2px)), + layer(x=1:1500, y=temps[1:1500], Geom.line), + Theme(panel_fill=color("white")) +) +# draw(PNG("yearly_fit.png", 16cm, 12cm), p) # Plot residuals for a few days -plt.figure(2) -plt.plot(residuals[1:100], color="g") +p = plot( + x=1:100, y=residuals[1:100], Geom.line, + Theme(default_color=color("green"), panel_fill=color("white")) +) +# draw(PNG("residuals.png", 16cm, 12cm), p) # Generate the residuals matrix ar_len = 5 @@ -41,19 +46,25 @@ end ar_coef = Variable(ar_len) optval2 = minimize!(sum_squares(residuals_mat * ar_coef - residuals[ar_len + 1 : end])) -# plot autoregressive fit of daily fluctuations for first few days -plt.figure(3) -plt.plot(residuals[ar_len + 1 : ar_len + 100], color="g") -plt.plot(residuals_mat[1:100, :] * evaluate(ar_coef), color="r") -plt.title("Autoregressive Fit of Residuals") +# plot autoregressive fit of daily fluctuations for a few days +ar_range = 1:145 +day_range = ar_range + ar_len +p = plot( + layer(x=day_range, y=residuals[day_range], Geom.line, Theme(default_color=color("green"))), + layer(x=day_range, y=residuals_mat[ar_range, :] * evaluate(ar_coef), Geom.line, Theme(default_color=color("red"))), + Theme(panel_fill=color("white")) +) +# draw(PNG("ar_fit.png", 16cm, 12cm), p) + -# plot final fit of data -plt.figure(4) -plt.plot(temps) total_estimate = evaluate(yearly) total_estimate[ar_len + 1 : end] += residuals_mat * evaluate(ar_coef) -plt.plot(total_estimate, color="r", alpha=0.5) -plt.title("Total Fit of Data") -plt.xlim([0, n]) -plt.show() +# plot final fit of data +p = plot( + layer(x=1:1500, y=total_estimate[1:1500], Geom.line, Theme(default_color=color("red"))), + layer(x=1:1500, y=temps[1:1500], Geom.line), + Theme(panel_fill=color("white")) +) +# draw(PNG("total_fit.png", 16cm, 12cm), p) + diff --git a/examples/time_series/total_fit.png b/examples/time_series/total_fit.png new file mode 100644 index 0000000..030ed31 Binary files /dev/null and b/examples/time_series/total_fit.png differ diff --git a/examples/time_series/yearly_fit.png b/examples/time_series/yearly_fit.png new file mode 100644 index 0000000..a763dfe Binary files /dev/null and b/examples/time_series/yearly_fit.png differ diff --git a/examples/tomography/tomography.jl b/examples/tomography/tomography.jl index f29a408..6b5e9c3 100644 --- a/examples/tomography/tomography.jl +++ b/examples/tomography/tomography.jl @@ -1,9 +1,10 @@ using LinearLeastSquares +using Gadfly -line_mat_x = readdlm("tux_sparse_x.txt"); -line_mat_y = readdlm("tux_sparse_y.txt"); -line_mat_val = readdlm("tux_sparse_val.txt"); -line_vals = readdlm("tux_sparse_lines.txt"); +line_mat_x = readdlm("tux_sparse_x.txt") +line_mat_y = readdlm("tux_sparse_y.txt") +line_mat_val = readdlm("tux_sparse_val.txt") +line_vals = readdlm("tux_sparse_lines.txt") # Form the sparse matrix from the data # Image is 50 x 50 @@ -11,21 +12,33 @@ img_size = 50 # The number of pixels in the image num_pixels = img_size * img_size -line_mat = spzeros(3300, num_pixels); +line_mat = spzeros(3300, num_pixels) num_vals = length(line_mat_val) for i in 1:num_vals - x = int(line_mat_x[i]); - y = int(line_mat_y[i]); - line_mat[x + 1, y + 1] = line_mat_val[i]; + x = int(line_mat_x[i]) + y = int(line_mat_y[i]) + line_mat[x + 1, y + 1] = line_mat_val[i] end x = Variable(num_pixels) -objective = sum_squares(line_mat * x - line_vals); -optval = minimize!(objective); +objective = sum_squares(line_mat * x - line_vals) +optval = minimize!(objective) -import PyPlot.plt -import PyPlot.cm -using PyPlot -plt.imshow(reshape(evaluate(x), img_size,img_size), cmap = get_cmaps()[29]) +rows = zeros(img_size*img_size) +cols = zeros(img_size*img_size) +for i = 1:img_size + for j = 1:img_size + rows[(i-1)*img_size + j] = i + cols[(i-1)*img_size + j] = img_size + 1 - j + end +end + +p = plot( + x=rows, y=cols, color=reshape(evaluate(x), img_size, img_size), Geom.rectbin, + Scale.ContinuousColorScale(Scale.lab_gradient(color("black"), color("white"))) +) +draw(PNG("tomography.png", 16cm, 14cm), p) + +# plt.imshow(reshape(evaluate(x), img_size, img_size), cmap = get_cmaps()[29]) diff --git a/examples/tomography/tomography.png b/examples/tomography/tomography.png new file mode 100644 index 0000000..6aa45f9 Binary files /dev/null and b/examples/tomography/tomography.png differ diff --git a/src/atoms/add_subtract.jl b/src/atoms/add_subtract.jl index 9dbb79b..ad022c6 100644 --- a/src/atoms/add_subtract.jl +++ b/src/atoms/add_subtract.jl @@ -12,7 +12,7 @@ function -(x::AffineExpr) vars_to_coeffs_map[v] = -c end this = AffineExpr(:-, (x,), vars_to_coeffs_map, -x.constant, x.size) - this.evaluate = ()-> -this.value; + this.evaluate = ()->(-x.evaluate()) return this end @@ -78,8 +78,10 @@ end function +(x::SumSquaresExpr, y::SumSquaresExpr) affines = copy(x.affines) + multipliers = copy(x.multipliers) append!(affines, y.affines) - this = SumSquaresExpr(:+, affines) + append!(multipliers, y.multipliers) + this = SumSquaresExpr(:+, affines, multipliers) return this end diff --git a/src/atoms/diag.jl b/src/atoms/diag.jl index 3677011..9cf21dd 100644 --- a/src/atoms/diag.jl +++ b/src/atoms/diag.jl @@ -14,10 +14,16 @@ function diag(x::AffineExpr, num::Int64) end len = min(x.size[1] - start_row + 1, x.size[2] - start_col + 1) start_ind = (start_col - 1) * x.size[1] + start_row - indexer = Constant(spzeros(len, x.size[1] * x.size[2])) - for i = 1 : len - indexer.value[i, start_ind + (i - 1) * (x.size[1] + 1)] = 1 - end + + indexer = Constant( + sparse( + 1 : len, + start_ind : x.size[1] + 1 : start_ind + (len - 1) * (x.size[1] + 1), + 1.0, + len, + x.size[1] * x.size[2] + ) + ) vars_to_coeffs_map = Dict{Uint64, Constant}() for (v, c) in x.vars_to_coeffs_map @@ -35,10 +41,16 @@ function diagm(x::AffineExpr) if x.size[2] != 1 error("Can only call diagm on column vectors") end - indexer = Constant(spzeros(x.size[1] * x.size[1], x.size[1])) - for i = 1 : x.size[1] - indexer.value[1 + (i - 1) * (x.size[1] + 1), i] = 1 - end + + indexer = Constant( + sparse( + 1 : x.size[1] + 1 : 1 + (x.size[1] - 1) * (x.size[1] + 1), + 1 : x.size[1], + 1.0, + x.size[1] * x.size[1], + x.size[1] + ) + ) vars_to_coeffs_map = Dict{Uint64, Constant}() for (v, c) in x.vars_to_coeffs_map diff --git a/src/atoms/index.jl b/src/atoms/index.jl index b5e0eb9..57ed627 100644 --- a/src/atoms/index.jl +++ b/src/atoms/index.jl @@ -6,19 +6,17 @@ getindex{T <: Real}(x::Constant, inds::AbstractArray{T, 1}) = Constant(getindex( getindex{T <: Real}(x::Constant, rows::AbstractArray{T, 1}, cols::AbstractArray{T, 1}) = Constant(getindex(x.value, rows, cols)) function getindex{T <: Real}(x::AffineExpr, inds::AbstractArray{T, 1}) - indexer = Constant(spzeros(length(inds), x.size[1] * x.size[2])) - k = 1 - for i in inds - indexer.value[k, i] = 1 - k += 1 - end + # number of rows/cols in the coefficient for x in our canonical form + num_rows_coeff = length(inds) + num_cols_coeff = x.size[1] * x.size[2] + indexer = Constant(sparse(1:length(inds), inds, 1.0, num_rows_coeff, num_cols_coeff)) vars_to_coeffs_map = Dict{Uint64, Constant}() for (v, c) in x.vars_to_coeffs_map vars_to_coeffs_map[v] = indexer * c end constant = indexer * x.constant - this = AffineExpr(:index, (x,), vars_to_coeffs_map, constant, (length(inds), 1)) + this = AffineExpr(:index, (x,), vars_to_coeffs_map, constant, (num_rows_coeff, 1)) this.evaluate = ()->x.evaluate()[inds] return this end @@ -28,18 +26,17 @@ function getindex{T <: Real}(x::AffineExpr, rows::AbstractArray{T, 1}, cols::Abs num_rows_coeff = length(rows) * length(cols) num_cols_coeff = x.size[1] * x.size[2] - indexer = Constant(spzeros(num_rows_coeff, num_cols_coeff)) - # Create the indexing matrix such that indexer * vec(x) = vec(x[rows, cols]) + J = Array(Int64, num_rows_coeff) k = 1 num_rows = x.size[1] for c in cols for r in rows - idx = num_rows * (convert(Int64, c) - 1) + convert(Int64, r) - indexer.value[k, idx] = 1 + J[k] = num_rows * (convert(Int64, c) - 1) + convert(Int64, r) k += 1 end end + indexer = Constant(sparse(1:num_rows_coeff, J, 1.0, num_rows_coeff, num_cols_coeff)) vars_to_coeffs_map = Dict{Uint64, Constant}() for (v, c) in x.vars_to_coeffs_map diff --git a/src/atoms/multiply_divide.jl b/src/atoms/multiply_divide.jl index c49d2a3..8c196c9 100644 --- a/src/atoms/multiply_divide.jl +++ b/src/atoms/multiply_divide.jl @@ -92,38 +92,40 @@ end # Only support scalar division function /(x::AffineExpr, y::Constant) - if y.size == (1, 1) - return x ./ y + if y.value == 0 + error("Cannot divide an affine expression by 0") else - error("Cannot divide by an expression whose size is not 1 by 1") + return x ./ y end end -./(x::AffineExpr, y::Numeric) = ./(x, convert(Constant, y)) -/(x::AffineExpr, y::Numeric) = /(x, convert(Constant, y)) +/(x::AffineExpr, y::Number) = /(x, convert(Constant, y)) ## Sum of squares expressions -function *(x::Constant, y::SumSquaresExpr) - if x.size != (1, 1) || x.value < 0 +function *(x::Number, y::SumSquaresExpr) + try + x = convert(Float64, x) + catch + error("Sum Squares expressions can only be multiplied by real numbers") + end + if x < 0 error("Sum Squares expressions can only be multiplied by nonegative scalars") end - x_sqrt = Constant(sqrt(x.value)) - affines = [x_sqrt * affine for affine in y.affines] - this = SumSquaresExpr(:*, affines) - return this + return SumSquaresExpr(:*, y.affines, x * y.multipliers) end -*(x::SumSquaresExpr, y::Constant) = *(y, x) -*(x::SumSquaresExpr, y::Numeric) = *(x, convert(Constant, y)) -*(x::Numeric, y::SumSquaresExpr) = *(convert(Constant, x), y) +*(x::SumSquaresExpr, y::Number) = *(y, x) -function /(x::SumSquaresExpr, y::Constant) - if y.size != (1, 1) || all(y.value .<= 0) +function /(x::SumSquaresExpr, y::Number) + try + y = convert(Float64, y) + catch + error("Sum Squares expressions can only be divided by real numbers") + end + if y <= 0 error("Sum Squares expressions can only be divided by positive scalars") end - return Constant(1 ./ y.value) * x + return SumSquaresExpr(:/, x.affines, x.multipliers / y) end - -/(x::SumSquaresExpr, y::Numeric) = /(x, convert(Constant, y)) diff --git a/src/atoms/stack.jl b/src/atoms/stack.jl index ea7d4a2..d545590 100644 --- a/src/atoms/stack.jl +++ b/src/atoms/stack.jl @@ -79,5 +79,11 @@ function hcat(args::AffineExpr...) end function hvcat(rows::(Int64...), args::AffineExpr...) - error("hvcat not yet implemented") + row_exprs = AffineExpr[] + index = 0 + for row_size in rows + push!(row_exprs, hcat(args[index + 1 : index + row_size]...)) + index += row_size + end + return vcat(row_exprs...) end diff --git a/src/atoms/transpose.jl b/src/atoms/transpose.jl index 7b026f8..7f7963b 100644 --- a/src/atoms/transpose.jl +++ b/src/atoms/transpose.jl @@ -6,19 +6,22 @@ ctranspose(x::Constant) = transpose(x) function transpose(x::AffineExpr) vec_sz = x.size[1] * x.size[2] - perm_matrix = spzeros(vec_sz, vec_sz) num_rows = x.size[1] num_cols = x.size[2] + I = Array(Int64, vec_sz) + J = Array(Int64, vec_sz) + + k = 1 for r = 1:num_rows for c = 1:num_cols - i = (c - 1) * num_rows + r - j = (r - 1) * num_cols + c - perm_matrix[i, j] = 1.0 + J[k] = (c - 1) * num_rows + r + I[k] = (r - 1) * num_cols + c + k += 1 end end - perm_constant = Constant(perm_matrix) + perm_constant = Constant(sparse(I, J, 1.0)) vars_to_coeffs_map = Dict{Uint64, Constant}() for (v, c) in x.vars_to_coeffs_map @@ -26,7 +29,7 @@ function transpose(x::AffineExpr) end constant = perm_constant * x.constant - this = AffineExpr(:transpose, (x,), vars_to_coeffs_map, constant, (x.size[2], x.size[1])) + this = AffineExpr(:transpose, (x,), vars_to_coeffs_map, constant, (num_cols, num_rows)) this.evaluate = ()->x.evaluate()' return this end diff --git a/src/solve.jl b/src/solve.jl index 81b1517..e0b8c5a 100644 --- a/src/solve.jl +++ b/src/solve.jl @@ -1,4 +1,4 @@ -export minimize!, satisfy! +export minimize!, solve! # Eventually, this can take in arguments to specify solving method. # For now, backslash is the only supported method. @@ -6,9 +6,8 @@ function solve!(problem_type::Symbol, objective::SumSquaresExpr, constraints::Ar p = Problem(problem_type, objective, constraints) # return backslash_solve!(p) backslash_solve!(p) - if p.status == "infeasible" - println("Problem was infeasible") - return nothing + if p.status == "KKT singular" + error("KKT system is singular. No unique solution") elseif p.status == "solved" return p.optval else @@ -19,8 +18,9 @@ end minimize!(objective::SumSquaresExpr, constraints::Array{EqConstraint}) = solve!(:minimize, objective, constraints) minimize!(objective::SumSquaresExpr, constraint::EqConstraint) = minimize!(objective, [constraint]) +minimize!(objective::SumSquaresExpr, constraints::EqConstraint...) = minimize!(objective, [constraints...]) minimize!(objective::SumSquaresExpr) = minimize!(objective, EqConstraint[]) -solve!(constraints::Array{EqConstraint}) = solve!(:satisfy, SumSquaresExpr(:empty, AffineExpr[]), constraints) -solve!(constraint::EqConstraint) = satisfy!([constraint]) -solve!(constraints::EqConstraint...) = satisfy!([constraints...]) +solve!(constraints::Array{EqConstraint}) = solve!(:solve, SumSquaresExpr(:empty, AffineExpr[], Float64[]), constraints) +solve!(constraint::EqConstraint) = solve!([constraint]) +solve!(constraints::EqConstraint...) = solve!([constraints...]) diff --git a/src/solvers/backslash_solver.jl b/src/solvers/backslash_solver.jl index 999482f..96c1d98 100644 --- a/src/solvers/backslash_solver.jl +++ b/src/solvers/backslash_solver.jl @@ -1,4 +1,4 @@ -export backslash_solve! +export backslash_solve!, build_kkt_system function reset_value_and_add_vars!(x::AffineOrConstant, unique_vars_map::Dict{Uint64, AffineExpr}) if x.head != :constant @@ -78,11 +78,14 @@ function build_kkt_system(A, b, C, d) constant[1 : size(A, 2)] = -2 * A' * b if size(C, 1) != 0 coefficient[1 : size(C, 2), size(A, 2) + 1 : end] = C' + coefficient[size(A, 2) + 1 : end, 1 : size(C, 2)] = C + constant[size(A, 2) + 1 : end] = -d end - end - if size(C, 1) != 0 - coefficient[size(A, 2) + 1 : end, 1 : size(C, 2)] = C - constant[size(A, 2) + 1 : end] = -d + elseif size(C, 1) != 0 + coefficient[:, 1 : size(C, 2)] = C + constant = -d + else + error("KKT system should not be empty") end return coefficient, constant @@ -104,19 +107,24 @@ function backslash_solve!(p::Problem) push!(canon_forms, constraint.canon_form) end C, d = coalesce_affine_exprs(vars_to_ranges_map, num_vars, canon_forms) - # return A, b, C, d - coefficient, constant = build_kkt_system(A, b, C, d) solution = nothing - try - solution = coefficient\full(constant) - catch - println("KKT system is singular") + # return A, b, C, d + if size(A, 1) > 0 || size(C, 1) > 0 + coefficient, constant = build_kkt_system(A, b, C, d) + try + # sparse solver doesnt error out some 2x2 singular KKT systems + solution = coefficient\full(constant) + end end if solution == nothing - p.status = "infeasible" + p.status = "KKT singular" else p.status = "solved" - p.optval = norm(A*solution[1:num_vars] + b)^2 + if size(A, 1) > 0 + p.optval = sum((A*solution[1:num_vars] + b).^2) + else + p.optval = 0.0 + end populate_vars!(unique_vars_map, vars_to_ranges_map, solution) end return p.optval diff --git a/src/types/expressions.jl b/src/types/expressions.jl index 8285ff8..953b59c 100644 --- a/src/types/expressions.jl +++ b/src/types/expressions.jl @@ -12,17 +12,18 @@ Numeric = Union(Number, AbstractArray) ValueOrNothing = Union(Value, Nothing) abstract AbstractExpr +abstract AffineOrConstant <: AbstractExpr function evaluate(x::AbstractExpr) if x.size == (1, 1) return (x.evaluate())[1] + elseif x.size[2] == 1 + return vec(full(x.evaluate())) else return full(x.evaluate()) end end -abstract AffineOrConstant <: AbstractExpr - type Constant <: AffineOrConstant head::Symbol value::ValueOrNothing @@ -72,6 +73,9 @@ function AffineConstant(value::Value) end function Variable(size::(Int64, Int64)) + if (size[1] < 1 || size[2] < 1) + error("invalid variable size") + end vec_sz = size[1]*size[2] this = AffineExpr(:variable, (), Dict{Uint64, Constant}(), Constant(spzeros(vec_sz, 1)), size) this.vars_to_coeffs_map[this.uid] = Constant(speye(vec_sz)) @@ -95,18 +99,19 @@ type SumSquaresExpr <: AbstractExpr head::Symbol value::ValueOrNothing affines::Array{AffineExpr} + multipliers::Array{Float64} uid::Uint64 evaluate::Function - function SumSquaresExpr(head::Symbol, affines::Array{AffineExpr}) - this = new(head, nothing, affines) + function SumSquaresExpr(head::Symbol, affines::Array{AffineExpr}, multipliers::Array{Float64}) + this = new(head, nothing, affines, multipliers) this.uid = object_id(this) this.evaluate = ()->begin - sum = 0 - for affine in affines - sum += norm(vec(affine.evaluate()))^2 + value = 0 + for i in 1:length(this.affines) + value += multipliers[i] * sum((this.affines[i].evaluate()).^2) end - return sum + return value end return this end @@ -115,7 +120,11 @@ end function sum_squares(affine::AffineExpr) affines = AffineExpr[] push!(affines, affine) - this = SumSquaresExpr(:sum_squares, affines) + this = SumSquaresExpr(:sum_squares, affines, ones(length(affines))) +end + +function evaluate(x::SumSquaresExpr) + return (x.evaluate())[1] end endof(x::AffineOrConstant) = x.size[1]*x.size[2] diff --git a/test/runtests.jl b/test/runtests.jl new file mode 100644 index 0000000..ab231ce --- /dev/null +++ b/test/runtests.jl @@ -0,0 +1 @@ +include("unit_tests/unit_tests.jl") diff --git a/test/unit_tests/add_subtract.jl b/test/unit_tests/add_subtract.jl new file mode 100644 index 0000000..92a7626 --- /dev/null +++ b/test/unit_tests/add_subtract.jl @@ -0,0 +1,19 @@ +a = create_var(1, 1) +b = eye(3) +c = a + b +test_expr(c, (a,)) + +p = b + c +test_expr(p, (a,)) + +x = create_var(3, 4) +d = a - x +test_expr(d, (a, x)) + +f = x + a +test_expr(f, (a, x)) + +g = x - a +test_expr(g, (a, x)) + +info("All add/subtract tests passed") diff --git a/test/unit_tests/compound.jl b/test/unit_tests/compound.jl new file mode 100644 index 0000000..db97015 --- /dev/null +++ b/test/unit_tests/compound.jl @@ -0,0 +1,22 @@ +# tests for sum, mean, and var, which implemented by calling other atoms + +x = create_var(5, 3) +s = sum(x) +m = mean(x) +test_expr(s, (x,)) +test_expr(m, (x,)) + +t = sum(x, 1) +n = mean(x, 1) +test_expr(t, (x,)) +test_expr(n, (x,)) + +u = sum(x, 2) +p = mean(x, 2) +test_expr(u, (x,)) +test_expr(p, (x,)) + +v = var(x) +test_expr(v) + +info("All sum/mean/var tests passed") diff --git a/test/unit_tests/constraints.jl b/test/unit_tests/constraints.jl new file mode 100644 index 0000000..4ff464c --- /dev/null +++ b/test/unit_tests/constraints.jl @@ -0,0 +1,9 @@ +x = Variable(4, 3) +y = Variable() + +x + y == 3 +3 == y +rand(4, 3) == x - 4*y +sum(y, 1) == rand(1, 3) + +info("All equality constraints tests passed") diff --git a/test/unit_tests/diag.jl b/test/unit_tests/diag.jl new file mode 100644 index 0000000..4a52e8c --- /dev/null +++ b/test/unit_tests/diag.jl @@ -0,0 +1,20 @@ +# diagm test +x = create_var(4, 1) +d = diagm(x) +test_expr(d, (x,)) + +# main diagonal +x = create_var(3, 4) +d = diag(x, 0) +test_expr(d, (x,)) + +# off diagonals +x = create_var(3, 4) +d = diag(x, 2) +test_expr(d, (x,)) + +x = create_var(3, 4) +d = diag(x, -2) +test_expr(d, (x,)) + +info("All diag tests passed") diff --git a/test/unit_tests/helpers.jl b/test/unit_tests/helpers.jl new file mode 100644 index 0000000..8f9d707 --- /dev/null +++ b/test/unit_tests/helpers.jl @@ -0,0 +1,46 @@ +export eval_internals, set_value!, test_expr, create_var + +function eval_internals(x::AffineExpr, vars::Tuple) + (rows, cols) = x.size + value = zeros(rows, cols) + for var in vars + coeff = x.vars_to_coeffs_map[object_id(var)] + value += reshape(coeff.value * vec(var.value), rows, cols) + end + value += reshape(x.constant.value, rows, cols) + return value +end + +function eval_internals(x::SumSquaresExpr) + value = 0 + for i in 1:length(x.affines) + value += x.multipliers[i]*vecnorm(evaluate(x.affines[i]))^2 + end + return value +end + +function set_value!(x::AffineExpr, value::Number) + x.value = convert(Float64, value) +end + +function set_value!(x::AffineExpr, value::AbstractArray) + x.value = convert(SparseMatrixCSC{Float64, Int64}, sparse(value)) +end + +function test_expr(x::SumSquaresExpr) + val1 = evaluate(x) + val2 = eval_internals(x) + @assert abs(val1 .- val2) <= 0.0001 +end + +function test_expr(x::AffineExpr, vars::Tuple) + val1 = eval_internals(x, vars) + val2 = evaluate(x) + @assert all(abs(val1 .- val2) .<= 0.0001) +end + +function create_var(m::Int64, n::Int64) + x = Variable(m, n) + set_value!(x, rand(m, n)) + return x +end diff --git a/test/unit_tests/index.jl b/test/unit_tests/index.jl new file mode 100644 index 0000000..6cc4bcc --- /dev/null +++ b/test/unit_tests/index.jl @@ -0,0 +1,8 @@ +# indexing tests +x = create_var(4, 6) +y = x[1:3, 5:end] +z = x[:] +test_expr(y, (x,)) +test_expr(z, (x,)) + +info("All index tests passed") diff --git a/test/unit_tests/multiply_divide.jl b/test/unit_tests/multiply_divide.jl new file mode 100644 index 0000000..3832aed --- /dev/null +++ b/test/unit_tests/multiply_divide.jl @@ -0,0 +1,19 @@ +a = create_var(1, 1) +b = eye(3) +c = a * b +test_expr(c, (a,)) + +d = create_var(5, 3) +f = a * 3.1 + 5.3 * d +test_expr(f, (a, d)) + +g = rand(4, 5) * d +test_expr(g, (d,)) + +h = d * rand(3, 4) +test_expr(h, (d,)) + +j = d / 4 +test_expr(j, (d,)) + +info("All multiply/divide tests passed") diff --git a/test/unit_tests/problems.jl b/test/unit_tests/problems.jl new file mode 100644 index 0000000..57bfce1 --- /dev/null +++ b/test/unit_tests/problems.jl @@ -0,0 +1,47 @@ +TOLERANCE = 0.0001 + +x = Variable() +solve!(x == 3) +@assert evaluate(x) == 3 + +y = Variable(3, 1) +A = rand(3, 3) +solve!(x == 3, A*y == 4) +@assert evaluate(x) == 3 +@assert all(abs(A*evaluate(y) - 4) .<= TOLERANCE) + +solve!([x == 5, y == 7]) +@assert evaluate(x) == 5 +@assert all(evaluate(y) .== 7) + +x = Variable(3) +y = sum_squares(x) +optval = minimize!(y) +@assert abs(optval - 0) < TOLERANCE +@assert all(abs(evaluate(x)) .<= TOLERANCE) + +x = Variable(3) +y = sum_squares(x) +optval = minimize!(y, x == 3) +@assert abs(optval - 27) < TOLERANCE +@assert all(abs(evaluate(x) .- 3) .<= TOLERANCE) + +x = Variable(2) +A = [1 0; 1 0] +y = sum_squares(x) +@test_throws ErrorException optval = minimize!(y, A * x == [1; 2]) + +A = randn(5, 5) +x_real = randn(5, 5) +x = Variable(5, 5) +y = A * x +b = A * x_real +constraints = EqConstraint[] +for i = 1:5 + constraints += y[:, i] == b[:, i] +end +optval = minimize!(sum_squares(x), constraints) +@assert abs(sum(x_real.^2) - optval) <= TOLERANCE +@assert all(abs(evaluate(x) .- x_real) .<= TOLERANCE) + +info("All solve/minimize tests passed") diff --git a/test/unit_tests/repmat.jl b/test/unit_tests/repmat.jl new file mode 100644 index 0000000..a8d8df1 --- /dev/null +++ b/test/unit_tests/repmat.jl @@ -0,0 +1,19 @@ +x = create_var(3, 2) +r = repmat(x, 4, 5) +test_expr(r, (x,)) + +x = create_var(1, 1) +r = repmat(x, 1, 1) +test_expr(r, (x,)) + +x = create_var(3, 3) +y = x + [1 2 3; 4 5 6; 7 8 9] +r = repmat(y, 1, 3) +test_expr(r, (x,)) + +x = create_var(3, 3) +y = x + [1 2 3; 4 5 6; 7 8 9] +r = repmat(y, 4, 1) +test_expr(r, (x,)) + +info("All repmat tests passed") diff --git a/test/unit_tests/reshape.jl b/test/unit_tests/reshape.jl new file mode 100644 index 0000000..695503f --- /dev/null +++ b/test/unit_tests/reshape.jl @@ -0,0 +1,10 @@ +# reshape test +x = create_var(6, 4) +y = reshape(x, 8, 3) +z = reshape(x, 4, 6) +v = vec(x) +test_expr(y, (x,)) +test_expr(z, (x,)) +test_expr(v, (x,)) + +info("All reshape tests passed") diff --git a/test/unit_tests/stack.jl b/test/unit_tests/stack.jl new file mode 100644 index 0000000..9fb2857 --- /dev/null +++ b/test/unit_tests/stack.jl @@ -0,0 +1,21 @@ +# hcat test +x = create_var(3, 1) +y = create_var(3, 2) +z = create_var(3, 3) +h = [x y z] +test_expr(h, (x, y, z)) + +# vcat test +x = create_var(1, 3) +y = create_var(2, 3) +z = create_var(3, 3) +v = [x; y; z] +test_expr(v, (x, y, z)) +v = [x, z, y] +test_expr(v, (x, y, z)) + +# hvcat test +hv = [x x; y y; z z] +test_expr(hv, (x, y, z)) + +info("All stack tests passed") diff --git a/test/unit_tests/sum_squares.jl b/test/unit_tests/sum_squares.jl new file mode 100644 index 0000000..4b87364 --- /dev/null +++ b/test/unit_tests/sum_squares.jl @@ -0,0 +1,15 @@ +x = create_var(1, 1) +y = create_var(3, 1) +s = sum_squares(x) + sum_squares(y) +test_expr(s) + +s = 4 * s +test_expr(s) + +s = s * 2 +test_expr(s) + +s = s / 4 +test_expr(s) + +info("All sum squares tests passed") diff --git a/test/unit_tests/transpose.jl b/test/unit_tests/transpose.jl new file mode 100644 index 0000000..de56572 --- /dev/null +++ b/test/unit_tests/transpose.jl @@ -0,0 +1,9 @@ +x = create_var(3, 4) +y = x' +test_expr(y, (x,)) + +x = create_var(10, 1) +y = x'' +test_expr(y, (x,)) + +info("All transpose tests passed") diff --git a/tests/unit_tests/unit_tests.jl b/test/unit_tests/unit_tests.jl similarity index 56% rename from tests/unit_tests/unit_tests.jl rename to test/unit_tests/unit_tests.jl index 0cc2715..4bae10b 100644 --- a/tests/unit_tests/unit_tests.jl +++ b/test/unit_tests/unit_tests.jl @@ -1,5 +1,6 @@ # Run all unit tests using LinearLeastSquares +using Base.Test include("helpers.jl") include("add_subtract.jl") @@ -8,3 +9,9 @@ include("stack.jl") include("diag.jl") include("repmat.jl") include("index.jl") +include("reshape.jl") +include("transpose.jl") +include("sum_squares.jl") +include("compound.jl") +include("constraints.jl") +include("problems.jl") diff --git a/tests/test.jl b/tests/test.jl deleted file mode 100644 index 5ed8b11..0000000 --- a/tests/test.jl +++ /dev/null @@ -1,46 +0,0 @@ -using LinearLeastSquares - -TOLERANCE = 1e-4; - -x = Variable(3); -y = sum_squares(x); -optval = minimize!(y); -@assert abs(optval - 0) < TOLERANCE - -x = Variable(3); -y = sum_squares(x); -optval = minimize!(y, x == 3); -@assert abs(optval - 27) < TOLERANCE - -x = Variable(2); -A = [1 0; 1 0]; -y = sum_squares(x); -optval = minimize!(y, A * x == [1; 2]); -@assert optval == nothing - -A = randn(5, 5); -x_real = randn(5, 5); -x = Variable(5, 5); -y = A * x; -b = A * x_real; -constraints = EqConstraint[]; -for i = 1:5 - constraints += y[:, i] == b[:, i]; -end -optval = minimize!(sum_squares(x), constraints); - -n = 200; -# Specify the true value of the variable -true_vect = [-1; 1]; -# Create data and labels -X = randn(n, 2) * 2; -b = randn(n, 1) * 0.1; -y = X * true_vect + b; -a = Variable(2); -optval = minimize!(sum_squares(X * a - y)); - - -# x = Variable(5); -# p = minimize!(sum_squares(sum(x)), x == 3); -# solve!(p); -# @assert abs(p.optval - 225) < TOLERANCE diff --git a/tests/unit_tests/add_subtract.jl b/tests/unit_tests/add_subtract.jl deleted file mode 100644 index a6cbaca..0000000 --- a/tests/unit_tests/add_subtract.jl +++ /dev/null @@ -1,9 +0,0 @@ -a = Variable(); -b = eye(3); -c = a + b; -set_value!(a, 1); -val1 = eval_internals(c, (a,)); -val2 = c.evaluate(); -@assert all(val1 .== val2) - -info("All add/subtract tests passed") diff --git a/tests/unit_tests/diag.jl b/tests/unit_tests/diag.jl deleted file mode 100644 index 920bd09..0000000 --- a/tests/unit_tests/diag.jl +++ /dev/null @@ -1,33 +0,0 @@ -# diagm test -x = Variable(4); -d = diagm(x); -set_value!(x, [1; 2; 3; 4]); -val1 = eval_internals(d, (x,)); -val2 = d.evaluate(); -@assert all(val1 .== val2); - -# diag test 1 -x = Variable(3, 4); -d = diag(x, 0); -set_value!(x, [1 2 3 4; 5 6 7 8; 9 10 11 12]); -val1 = eval_internals(d, (x,)) -val2 = d.evaluate(); -@assert all(val1 .== val2); - -# diag test 2 -x = Variable(3, 4); -d = diag(x, 2); -set_value!(x, [1 2 3 4; 5 6 7 8; 9 10 11 12]); -val1 = eval_internals(d, (x,)) -val2 = d.evaluate(); -@assert all(val1 .== val2); - -# diag test 3 -x = Variable(3, 4); -d = diag(x, -2); -set_value!(x, [1 2 3 4; 5 6 7 8; 9 10 11 12]); -val1 = eval_internals(d, (x,)) -val2 = d.evaluate(); -@assert all(val1 .== val2); - -info("All diag tests passed.") diff --git a/tests/unit_tests/helpers.jl b/tests/unit_tests/helpers.jl deleted file mode 100644 index c9a7564..0000000 --- a/tests/unit_tests/helpers.jl +++ /dev/null @@ -1,20 +0,0 @@ -export eval_internals, set_value! - -function eval_internals(x::AffineExpr, vars::Tuple) - (rows, cols) = x.size - value = zeros(rows, cols) - for var in vars - coeff = x.vars_to_coeffs_map[object_id(var)] - value += reshape(coeff.value * vec(var.value), rows, cols) - end - value += reshape(x.constant.value, rows, cols) - return value -end - -function set_value!(x::AffineExpr, value::Number) - x.value = convert(Float64, value) -end - -function set_value!(x::AffineExpr, value::AbstractArray) - x.value = convert(SparseMatrixCSC{Float64, Int64}, sparse(value)) -end diff --git a/tests/unit_tests/index.jl b/tests/unit_tests/index.jl deleted file mode 100644 index f0f8564..0000000 --- a/tests/unit_tests/index.jl +++ /dev/null @@ -1,13 +0,0 @@ -# indexing tests -x = Variable(4, 6); -set_value!(x, [1 2 3 4 5 6; 7 8 9 10 11 12; 13 14 15 16 17 18; 19 20 21 22 23 24]); -y = x[1:3, 5:end]; -z = x[:]; -val1 = eval_internals(y, (x,)); -val2 = eval_internals(z, (x,)); -val3 = y.evaluate(); -val4 = z.evaluate(); -@assert all(val1 .== val3) -@assert all(val2 .== val4) - -info("All index tests passed.") diff --git a/tests/unit_tests/multiply_divide.jl b/tests/unit_tests/multiply_divide.jl deleted file mode 100644 index df5f5e9..0000000 --- a/tests/unit_tests/multiply_divide.jl +++ /dev/null @@ -1,9 +0,0 @@ -a = Variable(); -b = eye(3); -c = a * b; -set_value!(a, 1); -val1 = eval_internals(c, (a,)); -val2 = c.evaluate(); -@assert all(val1 .== val2) - -info("All multiply/divide tests passed") diff --git a/tests/unit_tests/repmat.jl b/tests/unit_tests/repmat.jl deleted file mode 100644 index 1363ac1..0000000 --- a/tests/unit_tests/repmat.jl +++ /dev/null @@ -1,32 +0,0 @@ -# repmat test 1 -x = Variable(3, 2); -r = repmat(x, 4, 5); -set_value!(x, [5 8; 6 9; 7 0]); -val1 = eval_internals(r, (x,)); -val2 = r.evaluate(); -@assert all(val1 .== val2) - -x = Variable(); -r = repmat(x, 1, 1); -set_value!(x, 3); -val1 = eval_internals(r, (x,)); -val2 = r.evaluate(); -@assert all(val1 .== val2) - -x = Variable(3, 3); -y = x + [1 2 3; 4 5 6; 7 8 9]; -r = repmat(y, 1, 3); -set_value!(x, eye(3)); -val1 = eval_internals(r, (x,)); -val2 = r.evaluate(); -@assert all(val1 .== val2) - -x = Variable(3, 3); -y = x + [1 2 3; 4 5 6; 7 8 9]; -r = repmat(y, 4, 1); -set_value!(x, eye(3)); -val1 = eval_internals(r, (x,)); -val2 = r.evaluate(); -@assert all(val1 .== val2) - -info("All repmat tests passed.") diff --git a/tests/unit_tests/stack.jl b/tests/unit_tests/stack.jl deleted file mode 100644 index 96922db..0000000 --- a/tests/unit_tests/stack.jl +++ /dev/null @@ -1,29 +0,0 @@ -# hcat test -x = Variable(3); -y = Variable(3, 2); -z = Variable(3, 3); -h = [x y z]; -set_value!(x, [2; 3; 4]); -set_value!(y, [5 8; 6 9; 7 0]); -set_value!(z, eye(3)); -val1 = eval_internals(h, (x, y, z)); -val2 = h.evaluate(); -@assert all(val1 .== val2); - -# vcat test -x = Variable(1, 3); -y = Variable(2, 3); -z = Variable(3, 3); -v = [x; y; z]; -set_value!(x, [2 3 4]); -set_value!(y, [5 6 7; 8 9 0]); -set_value!(z, eye(3)); -val1 = eval_internals(v, (x, y, z)); -val2 = v.evaluate(); -@assert all(val1 .== val2); - -# hvcat test -# TODO: add when hvcat is implemented - - -info("All stack tests passed.")