eigen_samples.cpp

#include <stdio>
#include <vector>
#include <Eigen>

int initializeABunch(MatrixXd &A)
{
    int n = 2;
    A = MatrixXd::Zero(n, n);
    A = MatrixXd::Random(n, n);
    A = MatrixXd::Identity(n, n);
    A = MatrixXd::Constant(n, n, 0.6);
    A = MatrixXd::One(n, n);
	MatrixXd B(n+1, n+1);
	B << A, 0, 1, 1;
	MatrixXd F(A); //constr F from A;
	F.resizeLike(A); //resize to same size as A
}

int accessABunch(MatrixXd &A, VectorXd &v)
{
	MatrixXd B = A.block(0,0,2,1); //get two rows, 1 wide from the top left corner
	VectorXd v2 = v.head(2); //get two first elemnts, also use: segement(0,2) and tail(4)
}

int calculateABunch(VectorXd x1, VectorXd x2)
{
	VectorXd a = VectorXd::Random(10);
	VectorXd b = VectorXd::Random(10);
	VectorXd c;

	c = a.dot(b); //scalar
	int sum = a.sum(); //calculate sum of elements of array
	int norm = a.norm(); //calculate norm of vector
	
	MatrixXd A = MatrixXd::Random(10,10);
	MatrixXd B = MatrixXd::Random(10,10);
	MatrixXd C;
	
	C = A*B;
	C = A.cwiseProduct(B); //cell wise product, like dot product for two 1:n matrixes
	C = A.transpose().inverse().determinant(); //have fun figuring out what to use this for
	
	//cross breeding
	C = a.asDiagonal(); //create matrix with diagonal entries from a;
	C = a.matrix();
	a = C.array();
	
	a = b.normalize();
	b.normalized() //does the normalize in situ (so now: a == b)
}

int normalizations()
{
    VectorXd n = VectorXd::Random(10);
    assert(n.normalized() == n / n.norm()) //or use normalize() for in-situ
}

int solveEigen()
{
    MatrixXd A = MatrixXd::Random(10, 10);
    EigenSolver<MatrixXd> eig = EigenSolver<MatrixXd>(A); //O(n^3)
    const VectorXd &V = eig.eigenvalues();
    const MatrixXd &M = eig.eigenvectors();
}

int solveLinearSystem()
{
    MatrixXd A = MatrixXd::Random(10, 10);
    VectorXd b = VectorXd::Random(10);
    PartialPivLU<MatrixXd>> save = A.partialPivLu();
    VectorXd x = save.solve(b);
    x = A.llt().solve(b); //solve with llt
}

int sparseMatrix()
{
    VectorXd b = VectorXd::Random(10);
    MatrixXd C = MatrixXd::Random(10, 10);
    SparseMatrix<double> A(10, 10);
    A.reserve(2);

    std::vector<Triplet<double>> vec;
    Triplet<double> trip(0, 0, 10);
    vec.push_back(trip);

    //insert triplets
    A.setFromTriplets(vec.begin(), vec.end());
    //or:
    for (auto it = vec.beginn(); it != vec.end(); it++)
        A.insert(it->row(), it->col(), it->value());
    //or:
    A = C.sparseView();
    A.makeCompressed();

    SparseLU<SparseMatrix<double>> solver;
    solver.analyzePattern(A);
	
	//use either factorize() or compute()
	solver.factorize(A);
	if (solver.inf() != Success)
		throw 1; //lu failed, A is not invertible
    solver.compute(A);
    VexctorXd x = solver.solve(b);
}

int ellpackMult(const std::vector<Triplet<double>> triplets, const int maxCols, const int rows)
{
    //ellpack format prepare
    std::vector<double> vals(rows * maxCols);
    std::fill(vals.begin(), vals.end(), 0);
    std::vector<double> columnPtr(rows * maxCols);
    std::fill(vals.begin(), vals.end(), -1);

	//create ellpack from triples
    for (auto it = tr : triplets)
    {
        for (int l = tr->row() * maxCols; l < (tr->row() + 1) * maxCols; l++)
        {
            if (columnPtr[l] == -1)
            {
                columnPtr[l] = tr->col();
                vals[l] = tr->value();
                goto found;
            }
        }
        throw 1; //not enought space in array
    found:
    }

    //multiply ellpack with vector x = A*v
    VectorXd v = VectorXd::Random(10);
    VectorXd x(10);
    for (int i = 0; i < v.size(); i++)
    {
        for (int j = i * maxCols; j < i * (maxCols + 1); i++)
        {
            if (columnPtr(j) != -1)
                x(i) += v(columnPtr(j)) * vals(j);
        }
    }

    //multiply ellpack with vector x = A.transpose()*v
    VectorXd v = VectorXd::Random(10);
    VectorXd x(10);
    for (int i = 0; i < v.size(); i++)
    {
        for (int j = i * maxCols; j < i * (maxCols + 1); i++)
        {
            if (columnPtr(j) != -1)
                x(columnPtr(j)) += v(i) * vals(j);
        }
    }    
}

int matrixPow(const MatrixXd &A, int k)
{
    MatrixXd X = MatrixXd::Indetity(A.rows(), A.cols());
    unsigned int p = 1;
    for (unsigned int j = 1; j <= std::ceil(std::log2(k)); j++)
    {
        if ((~p & j) == 0) //check if k has a 1 at pos j, if so, multiply
            X = X * A;
        A = A * A;
        p <<= 1;
    }
}

int cooMult(std::vector<Triplet<double>> &A, std::vector<Triplet<double>> & B, size_t n)
{
	MatrixXd C = MatrixXd::Zero(n,n);
	//"simple"
	// C = A * B
	for (auto it = tr : A)
	{
		for (auto it_2 = tr : B)
		{
			if (it->col() == it_2->row()) {
				C(it->col(), it_2->row()) += it->value() * it_2->value();
			}
		}
	}
	
	//efficient
	//first: sort triplets by col / row
	std::sort(A.begin(), A.end(), [](const Triplet<double> &a, const Triplet<double> &b) { return a->col() < b->col();});
	std::sort(B.begin(), B.end(), [](const Triplet<double> &a, const Triplet<double> &b) { return a->row() < b->row();});
	
	//second: find intersections
	std::vector<int> intersects;
	//go only once cause we're smart
	auto a_it = A.begin();
	auto b_it = B.begin();
	while (a_it != A.end() && b_it != B.end())
	{
		if (a_it->col() == b_it->row())
		{
			intersects.push_back(a_it->col());
			while (a_it->col() == b_it->row()) {
				a_it++;
			}
		} else if (a->col() < b_it->row()) {
			a_it++;
		} else {
			b_it++;
		}
	}
	
	//third: multiply! such smart wow many work
	auto A_it = A.begin();
	auto B_it = B.begin();
	std::vector<Triplet<double>> res_vec;
	for (auto i : intersects)
	{
		//search for iterators to start at
		A_it = std::find_if(A_it, A.end(), [i](const Triplet<double> a1) { return a1.col() == i; });
		B_it = std::find_if(B_it, B.end(), [i](const Triplet<double> b1) { return b1.row() == i; });
		
		for (; A_it != A.end(); A_it++)
		{
			if (A_it->col() != i)
				break;
			for (auto B_itt = B_it;B_itt != B.end(); B_itt++)
			{
				if (B_itt->row() != i)
					break;
				else
					res_vec.push_back(Triplet<double>(A_it->row(), B_itt->col(), A_it->value() * B_itt->value()));
			}
		}
	}
}

int maps()
{
	MatrixXd A;
	Map<VectorXd>(A.data(), A.size()> map = Map<VectorXd>(A.data(), A.size()); //A.size() is A.cols() * A.rows()
}

int grahmSchmidtFun(const MatrixXd &A)
{
    //create Q from A
    MatrixXd Q(A);
    //normalize q
    Q.col(0).normalize();
	//do tha loooop
    for (unsigned int i = 1; i < Q.cols(); i++)
    {
        Q.col(i) -= Q.leftCols(i) * (Q.leftCols(i).transpose() * A.col(j));
        double eps = std::numeric_limits<double>::denorm_min();
        if (Q.col(i).norm() <= eps * Q.col(i).norm())
            throw 1; //precicion is not good enough, or columns almost lineary dependet
        Q.col(i).normalize();
    }
}

int testOrthogonality(const MatrixXd &A)
{
    double val = (A * A.transpose() - MatrixXd::Identity(A.rows(), A.cols())).norm();
    double eps = std::numeric_limits<double>::denorm_min();
    if (val < eps)
        throw 1;
}

int normSolve(const MatrixXd &A, const Vectorxd &b)
{
	MatrixXd at = A.transpose() * A;
	MatrixXd right = A.transpose() * b;
	VectorXd solution = at.llt().solve(right);
}

int sparseStuff()
{
    SparseMatrix<double> matr;

    VectorXd v = matr.col(0).normalize();
}

int svd_solve()
{
	MatrixXd A = MatrixXd::Random(10,10);
	JacobiSVD<MatrixXd> svd(A, ComputeFullU | ComputeFullV); //ComputerThinV ... etc
	MatrixXd U = svd.matrixU();
	MatrixXd V = sdv.matrixV(),
	Matrix E = svd.singularValues().asDiagonal();
}

int qr_decomp()
{
	MatrixXd A = MatrixXd::Random(5,10);
	HouseholderQR<MatrixXd> hh(A);
	VectorXd f = vectorXd::Random(10);
	auto a = f.asDiagonal() * A; //exception
	
	a = d.asDiagonal();
	auto b = a * A; //no exception
	
	MatrixXd Q = (hh.householderQ() * Identity(5, 10)); //Identity part weglassen damit schneller calculiert wird 
	MatrixXd R = hh.matrixQR().block(0,0,10,10).template TriangularView<Upper>();
}

int main()
{
    std::cout << std::scientific << std::setprecision(3) << std::setw(15) << "hallo welt" << std::setw(2) << "n" << std::endl;
    Timer timer;

    for (unsigned int i = 0; i < 20; i++)
    {
        timer.start();

        timer.stop();
    }

    std::cout << timer.min() << std::endl;
}