Commit f28f5de7 authored by ralfh's avatar ralfh

New homework problems based on CodeExpert files

parent 96279a72
#ifndef MATRIXCLASS_HPP
#define MATRIXCLASS_HPP
// The purpose of this exercise is to introduce the Eigen::Matrix class.
#include <iostream>
// TO DO: Include the "Eigen/Dense" header file with #include <Eigen/Dense>
// START
#include <Eigen/Dense>
// END
/* SAM_LISTING_BEGIN_0 */
Eigen::Matrix2d smallTriangular( double a, double b, double c ){
/*
* This functions returns a 2 by 2 triangular matrix of doubles a, b, c.
*/
// We know in advance that we want to create a matrix of doubles with a fixed size of 2 by 2.
// Therefore, we pass the parameters <double, 2, 2> to the template class Eigen::Matrix.
Eigen::Matrix< double, 2, 2 > A;
// We have declared the variable A of the type Eigen::Matrix<double,2,2>,
// but we have not initialized its entries.
// We can do this using comma-initialization:
A << a, b, 0, c;
// Question: Is A now an upper triangular matrix, or a lower triangular matrix?
return A;
}
/* SAM_LISTING_END_0 */
Eigen::MatrixXd constantTriangular( int n, double val ){
/*
* This function returns an n by n upper triangular matrix with the constant value val in the upper triangular part.
*/
// Now we do not know the size of our matrix at compile time.
// Hence, we use the special value Eigen::Dynamic to set the size of A.
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > A;
// Eigen::Matrix has a method Zero(n_rows,n_cols) which returns the n_rows by n_cols matrix whose entries are all equal to zero.
A = Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic >::Zero( n, n );
// To get a triangular matrix, we must set the values above the diagonal to val.
// This we can do by a nested for-loop.
for(int row = 0; row < n; row++){
// TO DO: Write the inner for-loop.
// Hint: We can access and change the entries of A using parentheses, i.e. A(row,col) = val;
// Note that indexing starts at the value 0 (as usual), and not 1 (as in Matlab).
// START
for(int col=row; col<n; col++){
A( row, col ) = val;
}
// END
}
//std::cout << A << std::endl;
return A;
}
Eigen::VectorXd arithmetics(int m, int n){
/*
* This function does not do anything meaningful.
* It is only meant to show some simple Eigen arithmetics and matrix initializations.
*/
// Because the syntax Eigen::Matrix< type, n_rows, n_cols > is very cumbersome,
// Eigen provides convenience classes that work as shorthands.
// For example, Eigen::Matrix2d is shorthand for Eigen::Matrix< double, 2, 2 >.
// This declares dynamic size (signified by the letter 'X') matrices of doubles (signified by the letter 'd').
Eigen::MatrixXd A, B, C, I;
// We initialize the matrices arbitrarily.
A = constantTriangular( n, 5.0 ).transpose(); // an n by n lower triangular matrix
B = Eigen::MatrixXd::Random( m, n ); // a random m by n matrix with values between -1 and 1.
I = Eigen::MatrixXd::Identity( n, n ); // The n by n identity matrix.
// We can perform arithmetics on matrices: +, -, *
// Note that for + and -, the left hand matrix and the right hand matrix have to be the same size,
// and that matrix multiplication * is only defined if the number of columns
// in the left hand matrix is the same as the number of rows in the right hand matrix.
C = B*(A - 5.0*I);
// Eigen::VectorXd is shorthand for Eigen::Matrix< double, Eigen::Dynamic, 1 >.
// Hence, all the same arithmetics work for vectors.
Eigen::VectorXd u, v;
// TO DO: Initialize the entries of u as the integers from 1 to n, multiply u by the matrix C, and store the result in v.
// START
u = Eigen::VectorXd::LinSpaced( n, 1, n ); // Or use a for-loop.
v = C*u;
// END
//std::cout<<v<<std::endl;
return v;
}
double casting(){
/*
* This function does not do anything meaningful.
* It is only meant to introduce a few more typedefs and how to cast them.
*/
// RowVectorXi is a shorthand for Eigen::Matrix< int, 1, Eigen::Dynamic >
// Constant(2,1) creates a vector of size 2 and initializes the entries with the value 1.
Eigen::RowVectorXi u = Eigen::RowVectorXi::Constant(2,1);
// std::complex is a template class for complex numbers.
// Here we declare two complex numbers, with real and imaginary parts represented by doubles.
// z0 = 1 - i
// z1 = 5 + i
std::complex<double> z0(1,-1), z1(5,1);
// Declare and initialize a size 2 vector of std::complex<double>.
Eigen::Vector2cd v;
v(0) = z0;
v(1) = z1;
double x;
// TO DO: Calculate u*v and store the result in x.
// Hint: First use u.cast< NEW TYPE >() to cast the "int" vector u to a "std::complex<double>" vector.
// The result will be u*v = 1*(1-i) + 1*(5+i) = 6 + 0i, a real number. You can get the real part of a std::complex<double> by .real().
// START
std::complex<double> z = u.cast<std::complex<double>>()*v;
x = z.real();
// END
return x;
}
#endif
# Code Expert Generic Project Configuration.
#
# In case this file is omitted, the execution is done with the default values.
# You can either define a cmd or a script file to be executed.
# Default to script file setting a cmd overrides it.
cmd:
# compile: echo "COMPILE"
# run: echo "RUN"
# test: echo "TEST"
# submit: echo "SUBMIT"
# All scripts are set by default to "/scripts/${ACTION}.sh"
script:
compile: /scripts/compile.sh
run: /scripts/run.sh
test: /scripts/test.sh
submit: /scripts/submit.sh
\ No newline at end of file
descriptions:
default:
name: Default
default: true
# Optionally add mulitple versions here.
# Use the file name without the .md extension as key and add a `name:` to it.
\ No newline at end of file
# Tutorial: part 1/3
## The Eigen::Matrix template class
> The purpose of this exercise is to learn to declare and initialize [Eigen::Matrix objects](https://eigen.tuxfamily.org/dox/group__TutorialMatrixClass.html), and to get familiar with some handy typedefs and methods.
> Open "MatrixClass.hpp" and fill in the missing code in between the delimiters `// START` and `// END` according to the instructions preceded by `// TO DO:`.
> Read each function of the program carefully. The code will not compile until after the first task is finished, which is to include Eigen/Dense.
\ No newline at end of file
# Tests of MatrixClass.hpp
> There are four functions in MatrixClass.hpp and main.cpp has one test for each function.
***
> `smallTriangular(double a, double b, double c)`: This function does not need to be edited, so the first test is passed automatically as long as the code compiles. The function takes as arguments three doubles and creates a 2 by 2 upper triangular matrix from those numbers.
The test prints `smallTriangular(1,2,3)`, and the correct result is:
```
1 2
0 3
```
***
> `constantTriangular(int n, double val)`: This function creates an n by n triangular MatrixXd with the entries in upper triangular part equal to `val`.
The test prints `constantTriangular(3,20)`, and the correct result is:
```
20 20 20
0 20 20
0 0 20
```
***
> `arithmetics(int m, int n)`: This function performs some arbitrary arithmetics on a random m by n matrix (the random seed is fixed).
The test prints `arithmetics(2,5)`, and the correct result is:
```
-18.983
19.4972
```
***
> `casting()`: This function calculates the inner product of the integer vector `[2 1]` and the complex vector `[1-i, 5+i]` and returns the real part as a double.
The test prints `casting()`, and the correct result is:
```
6
```
#include "MatrixClass.hpp"
int main(){
std::cout<< "\nTest of smallTriangular():\n";
std::cout<< smallTriangular(1,2,3) << std::endl;
std::cout<< "\nTest of constantTriangular():\n";
std::cout<< constantTriangular(3,20) << std::endl;
std::cout<< "\nTest of arithmetics():\n";
std::srand(5); // So that random behaviour is predictable.
std::cout<< arithmetics(2,5) << std::endl;
std::cout<< "\nTest of casting():\n";
std::cout<< casting() << std::endl;
return 0;
}
#!/bin/bash
echo "Compiling ..."
# compile c code
# find . -iname '*.cpp' | sort | xargs g++ -fdiagnostics-color=always --pedantic -Wextra -Wall -std=c++14 -I/usr/include/python3.7m -llibpython -o bin/a.out || exit 1
# Compiling for 'matplotlibcpp.h'
mkdir -p bin
find . -iname '*.cpp' | sort | xargs \
g++ -fdiagnostics-color=always -std=c++11 \
-I/usr/local/include/python3.7m \
-I/usr/local/lib/python3.7/site-packages/numpy/core/include \
-I/usr/include/eigen3/ \
-lpython3.7m \
-lpthread -lutil -ldl \
-Xlinker -export-dynamic \
-o bin/a.out
echo "Compilation successful"
\ No newline at end of file
#!/bin/bash
export PYTHONIOENCODING="UTF-8"
# compile (call compile script)
bash "${WORKDIR}/scripts/compile.sh"
# run
bin/a.out
\ No newline at end of file
#!/bin/bash
echo "Not implemented"
\ No newline at end of file
#!/bin/bash
# run (call run script)
bash "${WORKDIR}/scripts/run.sh"
\ No newline at end of file
// "testname", "input", "expected output", "points", "execution time limit (ms)"
"smallTriangular", "(any input)", "1 2\n0 3", "1", "5000",
"constantTriangular", "(any input)", "20 20 20\n 0 20 20\n 0 0 20", "1", "5000",
"arithmetics", "(any input)", "-18.983\n19.4972", "1", "5000",
"casting", "(any input)", "6", "1", "5000",
\ No newline at end of file
## CodeExpert setup for homework problem "ArrowMatrix"
Compilation requires the following files:
* `Python,h`: must reside in the include path /usr/local/include/python3.7m
* `libpython3.7.m`: must be found in /usr/local/lib/python3.7
#ifndef ARROWMATRIX_FUNCTIONS_HPP
#define ARROWMATRIX_FUNCTIONS_HPP
#include <iostream>
#include <iomanip>
#include <Eigen/Dense>
#include <vector>
#include "timer.h"
#include "matplotlibcpp.h"
#include "plot.hpp"
using namespace Eigen;
namespace plt = matplotlibcpp;
/* @brief Build an "arrow matrix" and compute A*A*y
* Given vectors $a$ and $d$, returns A*A*x in $y$, where A is built from a, d
* @param[in] d An n-dimensional vector
* @param[in] a An n-dimensional vector
* @param[in] x An n-dimensional vector
* @param[out] y The vector y = A*A*x
*/
/* SAM_LISTING_BEGIN_0 */
void arrow_matrix_2_times_x(const VectorXd &d, const VectorXd &a,
const VectorXd &x, VectorXd &y) {
assert(d.size() == a.size() && a.size() == x.size() &&
"Vector size must be the same!");
int n = d.size();
// In this lines, we extract the blocks used to construct the matrix A.
VectorXd d_head = d.head(n-1);
VectorXd a_head = a.head(n-1);
MatrixXd d_diag = d_head.asDiagonal();
MatrixXd A(n,n);
// We build the matrix A using the "comma initialization": each expression separated
// by a comma is a "block" of the matrix we are building.
// d\_diag is the top left (n-1)x(n-1) block
// a\_head is the top right vertical vector
// a\_head.transpose() is the bottom left horizontal vector
// d(n-1) is a single element (a 1x1 matrix), on the bottom right corner
// This is how the matrix looks like:
// A = | D | a |
// |-----+--------|
// | a\^T | d(n-1) |
A << d_diag, a_head,
a_head.transpose(), d(n-1);
y = A*A*x;
}
/* SAM_LISTING_END_0 */
/* @brief Build an "arrow matrix"
* Given vectors $a$ and $b$, returns A*A*x in $y$, where A is build from a,d
* @param[in] d An n-dimensional vector
* @param[in] a An n-dimensional vector
* @param[in] x An n-dimensional vector
* @param[out] y The vector y = A*A*x
*/
/* SAM_LISTING_BEGIN_1 */
void efficient_arrow_matrix_2_times_x(const VectorXd &d,
const VectorXd &a,
const VectorXd &x,
VectorXd &y) {
assert(d.size() == a.size() && a.size() == x.size() &&
"Vector size must be the same!");
int n = d.size();
// Notice that we can compute (A*A)*x more efficiently using
// A*(A*x). This is, in fact, performing two matrix vector
// multiplications
// instead of a more expensive matrix-matrix multiplication.
// Therefore, as first step, we need a way to efficiently
// compute A*x
// This function computes A*x. you can use it
// by calling A\_times\_x(x).
// This is the syntax for lambda functions: notice the extra
// [variables] code. Each variable written within [] brackets
// will be captured (i.e. seen) inside the lambda function.
// Without \&, a copy of the variable will be performed.
// Notice that A = D + H + V, s.t. A*x = D*x + H*x + V*x
// D*x can be rewritten as d*x componentwise
// H*x is zero, except at the last component
// V*x is only affected by the last component of x
auto A_times_x = [&a, &d, n] (const VectorXd & x) {
// This takes care of the diagonal (D*x)
// Notice: we use d.array() to tell Eigen to treat
// a vector as an array. As a result: each operation
// is performed componentwise.
VectorXd Ax = ( d.array() * x.array() ).matrix();
// H*x only affects the last component of A*x
// This is a dot product between a and x with the last
// component removed
Ax(n-1) += a.head(n - 1).dot(x.head(n-1));
// V*x is equal to the vector
// (a(0)*x(n-1), ..., a(n-2)*x(n-1), 0)
Ax.head(n-1) += x(n-1) * a.head(n-1);
return Ax;
};
// <=> y = A*(A*x)
y = A_times_x(A_times_x(x));
}
/* SAM_LISTING_END_1 */
/* \brief Compute the runtime of arrow matrix multiplication.
* Repeat tests 10 times, and output the minimal runtime
* amongst all times. Test both the inefficient and the efficient
* versions.
*/
void runtime_arrow_matrix() {
/* SAM_LISTING_BEGIN_3 */
//Memory allocation for plot
std::vector<double> vec_size;
std::vector<double> elap_time, elap_time_eff;
//header for result print out
std::cout << std::setw(8) << "n"
<< std::setw(15) << "original"
<< std::setw(15) << "efficient"
<< std::endl;
for (unsigned n = 8; n <= 128; n *= 2) {
//save vector size (for plot)
vec_size.push_back(n);
// Number of repetitions
unsigned int repeats = 10;
Timer timer, timer_eff;
// Repeat test many times
for (unsigned int r = 0; r < repeats; ++r) {
// Create test input using random vectors
Eigen::VectorXd a = Eigen::VectorXd::Random(n),
d = Eigen::VectorXd::Random(n),
x = Eigen::VectorXd::Random(n),
y;
// Compute times
timer.start();
arrow_matrix_2_times_x(d, a, x, y);
timer.stop();
// Compute times for efficient implementation
timer_eff.start();
efficient_arrow_matrix_2_times_x(d, a, x, y);
timer_eff.stop();
}
// Print results (for grading): inefficient
std::cout << std::setw(8) << n
<< std::scientific << std::setprecision(3)
<< std::setw(15) << timer.min()
<< std::setw(15) << timer_eff.min()
<< std::endl;
//time needed
elap_time.push_back(timer.min());
elap_time_eff.push_back(timer_eff.min());
}
/* DO NOT CHANGE */
//create plot
plot(vec_size, elap_time, elap_time_eff, "./cx_out/text.png");
/* SAM_LISTING_END_3 */
}
// Additional
// Rename long variable name to duration_t (easy to change)
using duration_t = std::chrono::nanoseconds;
/* \brief Compute runtime of $F$, repeating test "repeats" times
* Will return minimal runtime.
* This function uses "crhono".
* \tparam Function type of F, must have an operator()
* \param[in] F Function for which you want to measure runtime.
* \param[in] repeats Number of repetitions.
*/
template <class Function>
duration_t timing(const Function & F, int repeats = 10) {
// Shortcut for time_point
using time_point_t = std::chrono::high_resolution_clock::time_point;
// Loop many times
duration_t min_elapsed;
for(int r = 0; r < repeats; r++) {
// Start clock (MATLAB: tic)
time_point_t start = std::chrono::high_resolution_clock::now();
// Run function
F();
// Stop clock (MATLAB: toc) and measure difference
duration_t elapsed = std::chrono::duration_cast<duration_t>(
std::chrono::high_resolution_clock::now() - start);
// Compute min between all runs
min_elapsed = r == 0 ? elapsed : std::min(elapsed, min_elapsed);
}
return min_elapsed;
}
/* \brief Compute timing using chrono
* Also demonstrate use of lambda functions
*/
void runtime_arrow_matrix_with_chrono() {
// Table header
std::cout << std::setw(8) << "n"
<< std::scientific << std::setprecision(3)
<< std::setw(15) << "original"
<< std::setw(15) << "efficient"
<< std::endl;
// Run from $2^5$ to $2^11$ with powers of two
for(unsigned int n = (1 << 5); n < (1 << 12); n = n << 1) {
// Create random vectors
VectorXd d = VectorXd::Random(n);
VectorXd a = VectorXd::Random(n);
VectorXd x = VectorXd::Random(n);
VectorXd y(n);
// Call "timing", using a lambda function for F
// Remember: we cannot pass arrow\_matrix\_2\_times\_x directly to timing
// the timing function expects a n object with operator()(void)
duration_t elapsed = timing([&a, &d, &x, &y] () {
arrow_matrix_2_times_x(d, a, x, y);
}, 10);
// Call "timing", using a lambda function for F
duration_t elapsed_efficient = timing([&a, &d, &x, &y] () {
efficient_arrow_matrix_2_times_x(d, a, x, y);
}, 10);
// Output timings (not part of the exercice)
std::cout << std::setw(8)<< n
<< std::scientific << std::setprecision(3)
<< std::setw(15) << elapsed.count() * 1e-9 // ns to s
<< std::setw(15) << elapsed_efficient.count() * 1e-9 // ns to s
<< std::endl;
}
}
#endif
\ No newline at end of file
# Code Expert Generic Project Configuration.
#
# In case this file is omitted, the execution is done with the default values.
# You can either define a cmd or a script file to be executed.
# Default to script file setting a cmd overrides it.
cmd:
# compile: echo "COMPILE"
# run: echo "RUN"
# test: echo "TEST"
# submit: echo "SUBMIT"
# All scripts are set by default to "/scripts/${ACTION}.sh"
script:
compile: /scripts/compile.sh
run: /scripts/run.sh
test: /scripts/test.sh
submit: /scripts/submit.sh
\ No newline at end of file
descriptions:
default:
name: Default
default: true
# Optionally add mulitple versions here.
# Use the file name without the .md extension as key and add a `name:` to it.
\ No newline at end of file
# Problem 2-1: Arrow matrix-vector mutliplication
> For the task description of this exercise, please refer to [NCSE19_Problems.pdf](
https://www.sam.math.ethz.ch/~grsam/NCSE19/HOMEWORK/NCSE19_Problems.pdf).
> Open "ArrowMatrix_functions.hpp" and fill in the missing code in between the delimiters `// START` and `// END` according to the instructions preceded by `// TO DO:`.
> To view plots, click the "Files" button on the toolbar at the bottom of the page.
\ No newline at end of file
# Tests of ArrowMatrix_functions.hpp
> There is no test for this exercise as zzzzz_test_runner.cpp does not work if there are plots created in the program.
\ No newline at end of file
#include <iostream>
#include <iomanip>
#include <Eigen/Dense>
#include <vector>
#include "timer.h"
#include "ArrowMatrix_functions.hpp"
using namespace Eigen;
int main(void) {
// Test vectors
VectorXd a(5);
a << 1., 2., 3., 4., 5.;
VectorXd d(5);
d <<1., 3., 4., 5., 6.;
VectorXd x(5);
x << -5., 4., 6., -8., 5.;
VectorXd yi;
// Run both functions
arrow_matrix_2_times_x(a,d,x,yi);
VectorXd ye(yi.size());
efficient_arrow_matrix_2_times_x(a,d,x,ye);
// Compute error
double err = (yi - ye).norm();
// Output error
std::cout << "--> Correctness test." << std::endl;
std::cout << "Error: " << err << std::endl;
// Print out runtime
std::cout << "--> Runtime test." << std::endl;
runtime_arrow_matrix();
std::cout << "Plot was created. See 'Files'." << std::endl;
}
#ifndef PLOT_HPP
#define PLOT_HPP
#include <cmath>
#include <Eigen/Dense>
#include <vector>
# include "matplotlibcpp.h"
namespace plt = matplotlibcpp;