text:
/*
 * Copyright 2017 Antonio Pastor anpastor{at}it.uc3m.es
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *  http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 *
 */

#include <cstdlib>
#include <cstdio>
#include <csignal>
#include <unistd.h>
#include <iostream>     // std::cout
#include <vector>
#include <tuple>        // std::tuple, std::get, std::tie, std::ignore
#include <map>
#include <unordered_map>
#include <gflags/gflags.h>
#include <boost/thread.hpp>
#include <pqxx/pqxx>

#include <fstream>

#include "zmqpp/zmqpp.hpp"

#include "../data-functions.h"

using std::cout;
using std::endl;
using std::string;

#define T_MSG   1 // microseconds
#define ONEEXPNINE 1000000000 //1e9

// TO-DO: define arguments for test running time and throughtput
DEFINE_string(day, "161201", "Day of the database to use for the test queries (in format YYMMDD)");
DEFINE_string(dbIP, "127.0.0.1", "IP address of the database");
DEFINE_int32(dbPORT, 5432, "Port of the database");
DEFINE_string(dbPWD, "", "password of the database");
DEFINE_string(dbUSER, "user", "database user");
DEFINE_string(dbNAME, "nameles", "database name");
//DEFINE_int32(test_duration, 60, "Test duration in seconds");
//DEFINE_int32(MPS, 30000, "Queries per second");
DEFINE_int32(port, 58510, "\"Send to\" port");
DEFINE_int32(nWorkers, 4, "Number of workers");
DEFINE_int32(max_msgs, 1e4, "Maximum number of messages per thread per database serialization");

typedef std::tuple<uint32_t,uint32_t,string> msg_tuple;
typedef std::vector<msg_tuple> msgs_list;

long time_diff(const timespec t0, const timespec t1);

void send_msgs(const string& sendToSocket, const string& dbConnect, const int n_msgs);

boost::thread senderThread, receiverThread;

void SIGINT_handler(int s){
					printf("Caught signal %d\n",s);
					receiverThread.interrupt();
					senderThread.interrupt();
}

int main (int argc, char *argv[]){
	google::ParseCommandLineFlags(&argc, &argv, true);

	int n_msgs = 8*FLAGS_max_msgs; //FLAGS_test_duration * FLAGS_MPS;
	string sendToSocket = "tcp://127.0.0.1:" + std::to_string(FLAGS_port);
	string receiveFromSocket = "tcp://*:" + std::to_string(FLAGS_port);
	string dbConnect = "dbname=" + FLAGS_dbNAME + " user="+ FLAGS_dbUSER +
										" host="+ FLAGS_dbIP + " port=" + std::to_string(FLAGS_dbPORT) + " password=" + FLAGS_dbPWD;

	senderThread = boost::thread(send_msgs, sendToSocket, dbConnect, n_msgs);

	MessagesQueue requests_queue;
	receiverThread = boost::thread(receive_msgs, &requests_queue, receiveFromSocket);

	boost::thread_group workers;
	for( int x=0; x<FLAGS_nWorkers; ++x ) {
		workers.add_thread(new boost::thread(queue_consumer, &requests_queue,
			 																		dbConnect, FLAGS_max_msgs));
	}
	struct sigaction sigIntHandler;

	sigIntHandler.sa_handler = SIGINT_handler;
	sigemptyset(&sigIntHandler.sa_mask);
	sigIntHandler.sa_flags = 0;

	sigaction(SIGINT, &sigIntHandler, NULL);
	senderThread.join();
	cout << "main all sent" << endl;
	receiverThread.join();
	cout << "main receiver down" << endl;
	workers.interrupt_all();
	workers.join_all();

  auto t = std::time(nullptr);
  auto tm = *std::localtime(&t);
  std::ostringstream date_oss;
  date_oss << std::put_time(&tm, "%C%m%d");

	pqxx::connection c(dbConnect);
	pqxx::read_transaction txn(c);
	pqxx::result r = txn.exec("SELECT referrer, host(ip), sum(cnt) AS cnt FROM tuples.temp_" + date_oss.str()
						+ " GROUP BY host(ip),referrer");

	std::ofstream file_db;
	file_db.open("db_counts.csv");
	for (auto row: r) {
		file_db << row[0].as<string>() << "," << row[1].as<string>() << "," << row[2].as<string>() << endl;
	}
	file_db.close();
	exit(0);
}

void send_msgs(const string& sendToSocket, const string& dbConnect, const int n_msgs){
	zmqpp::context_t context;
	context.set(zmqpp::context_option::io_threads, 2);
	zmqpp::socket_t  sender(context, zmqpp::socket_type::push);
	sender.set(zmqpp::socket_option::linger, 1000); // 1 second
	try {
	sender.connect(sendToSocket);
	cout << "sender thread socket: " << sendToSocket << endl;
	} catch (zmqpp::zmq_internal_exception &e){
		cout << "Exception: " << e.what() << endl;
		sender.close();
		context.terminate();
		return;
	}

	strings_pair_counter_map bidTuplesCounter;

	pqxx::connection c(dbConnect);
	pqxx::read_transaction txn(c);
	pqxx::result r = txn.exec("SELECT referrer, host(ip) FROM tuples.ip_ref_" + FLAGS_day + " LIMIT " + std::to_string(n_msgs) );
	// cout << "Test queries loaded, starting latency test" << endl;
	uint32_t reqID = 1;
	struct timespec t0, t1;
	struct timespec ts0, ts1;
	timespec tspec_sleep = {0, 0};
	int msg_time = 1e2; // In nanoseconds
	string referrer, ip;
	clock_gettime(CLOCK_MONOTONIC_RAW, &t0);
	for (auto row: r) {
		clock_gettime(CLOCK_MONOTONIC_RAW, &ts0);
		zmqpp::message_t bid_req;
		referrer = row[0].as<string>();
		ip = row[1].as<string>();
		bid_req << reqID << referrer << ip;
		sender.send(bid_req, true);
		reqID++;
		bidTuplesCounter[strings_pair(referrer,ip)]++;
		clock_gettime(CLOCK_MONOTONIC_RAW, &ts1);
		int t_sleep = msg_time-time_diff(ts0,ts1);
		if (t_sleep>0){
			tspec_sleep.tv_nsec = t_sleep;
			clock_nanosleep(CLOCK_MONOTONIC, 0, &tspec_sleep, NULL);
		}
	}
	sender.close();
	context.terminate();

	std::vector<std::pair<strings_pair, unsigned int>> vector( bidTuplesCounter.begin(), bidTuplesCounter.end() );

	// Sort the vector according to the count in descending order.
	std::sort( vector.begin(), vector.end(),
						[]( const auto & lhs, const auto & rhs )
						{ return lhs.second > rhs.second; } );

	std::ofstream file_input;
	file_input.open("input_counts.csv");
	for ( const auto & it : vector){
		file_input << it.first.first << "," << it.first.second << "," << it.second << endl;
	}
	file_input.close();

	clock_gettime(CLOCK_MONOTONIC_RAW, &t1);
	cout << "Test duration: " << time_diff(t0, t1)/1e9 << endl;
	cout << "# of messages sent: " << std::dec << n_msgs << endl;
}

long time_diff(const timespec t0, const timespec t1){
	long sec = t1.tv_sec - t0.tv_sec;
	long nsec = t1.tv_nsec - t0.tv_nsec;
	if (nsec<0){
		nsec = ONEEXPNINE - nsec;
		sec--;
	}
	return sec*ONEEXPNINE + nsec;
}

max_stars_repo_name: Nameles-Org/Nameless-realtime

max_issues_repo_name: Nameles-Org/Nameless-realtime

max_forks_repo_name: Nameles-Org/Nameless-realtime


text:
#include <algorithm>
#include <fstream>

#include <boost/archive/text_oarchive.hpp>
#include <boost/archive/text_iarchive.hpp>

#include <Eigen/Dense>
#include "gmpxx.h"

#include "Expression.h"
#include "MonomialExpression.h"
#include "Scalar.h"
#include "ScalarSum.h"
#include "TensorMonomial.h"

template <typename Archive>
void Expression::serialize(Archive & ar, unsigned int const version) {
  if (version > 0) {
  } else {
  }
  ar & dimension;
  ar & summands;
}

template void Expression::serialize<boost::archive::text_oarchive>(boost::archive::text_oarchive&, unsigned int const);
template void Expression::serialize<boost::archive::text_iarchive>(boost::archive::text_iarchive&, unsigned int const);

typedef Eigen::Matrix<mpq_class, Eigen::Dynamic, Eigen::Dynamic> MatrixXq;

Expression::Expression() : dimension(4), summands(std::make_unique<std::vector<Summand>>()) { }

Expression::Expression(Expression const & other) : dimension(other.dimension), summands(std::make_unique<Sum>()) {
  std::for_each(other.summands->cbegin(), other.summands->cend(), [this](auto const & it) {
    summands->push_back(std::pair<std::unique_ptr<MonomialExpression>, std::unique_ptr<ScalarSum>>());
    auto monexp_new = std::make_unique<MonomialExpression>(*(it.first));
    auto scalar_new = std::make_unique<ScalarSum>(*(it.second));
    std::swap((*(--summands->end())).first, monexp_new);
    std::swap((*(--summands->end())).second, scalar_new);
  });
}

Expression & Expression::operator=(Expression const & other) {
  Expression tmp ( other );
  dimension = other.dimension;
  std::swap(summands, tmp.summands);
  return *this;
}

void Expression::set_dimension (unsigned int const dimension_set) {
  dimension = dimension_set;
}

unsigned int Expression::get_dimension (void) const {
  return dimension;
}

void Expression::AddSummand (MonomialExpression const & monomial_expression, ScalarSum const & scalar_sum) {
  summands->push_back(std::make_pair(std::make_unique<MonomialExpression>(monomial_expression), std::make_unique<ScalarSum>(scalar_sum)));
}

void Expression::AddSummand (MonomialExpression const & monomial_expression, Scalar const & scalar) {
  summands->push_back(std::make_pair(std::make_unique<MonomialExpression>(monomial_expression), std::make_unique<ScalarSum>(scalar)));
}

void Expression::AddSummand (MonomialExpression const & monomial_expression) {
  AddSummand(monomial_expression, Scalar());
}

void Expression::AddSummand (MonomialExpression const & monomial_expression, Rational const & rational) {
  AddSummand(monomial_expression, Scalar(rational));
}

void Expression::AddSummand (MonomialExpression const & monomial_expression, size_t variable) {
  AddSummand(monomial_expression, Scalar(variable));
}

void Expression::AddSummand (MonomialExpression const & monomial_expression, Rational const & rational, size_t variable) {
  AddSummand(monomial_expression, Scalar(rational, variable));
}

std::set<size_t> Expression::GetVariableSet () const {
  std::set<size_t> ret;
  std::for_each(summands->cbegin(), summands->cend(),
    [&ret] (auto const & summand) {
      ret.merge(summand.second->CoefficientSet());
    });

  return ret;
}

std::map<size_t, size_t> Expression::GetCoefficientMap () const {
  std::set<size_t> variables = GetVariableSet ();

  std::map<size_t, size_t> coefficient_map;
  std::for_each(variables.cbegin(), variables.cend(), [&coefficient_map,n=0] (auto const & a) mutable { coefficient_map[a] = n++; });

  return coefficient_map;
}

std::vector<std::vector<Rational>> Expression::GetPrefactorMatrix () const {
  std::map<size_t, size_t> coefficient_map = GetCoefficientMap ();

  std::vector<std::vector<Rational>> prefactor_matrix;
  std::for_each(summands->cbegin(), summands->cend(),
    [&prefactor_matrix, &coefficient_map] (auto const & summand) {
      prefactor_matrix.push_back(summand.second->CoefficientVector(coefficient_map));
    });

  std::sort(prefactor_matrix.begin(), prefactor_matrix.end());
  prefactor_matrix.erase(std::unique(prefactor_matrix.begin(), prefactor_matrix.end()), prefactor_matrix.end());

  return prefactor_matrix;
}

void Expression::SortMonomials() {
  auto summands_new = std::make_unique<Sum>();
  std::transform(summands->cbegin(), summands->cend(), std::back_inserter(*summands_new), [](auto const & a) { auto ret = std::make_pair(std::make_unique<MonomialExpression>(*a.first), std::make_unique<ScalarSum>(*a.second)); ret.first->Sort(); return ret; });
  std::swap(summands, summands_new);
}

void Expression::ApplyMonomialSymmetries() {
  for (auto & summand : *summands) {
    bool even = summand.first->ApplySymmetries();
    if (!even) {
      summand.second->Negate();
    }
  }
}

void Expression::ApplyMonomialSymmetriesToContractions() {
  bool repeat = true;

  while (repeat) {
    repeat = false;
    std::for_each(summands->begin(), summands->end(),
      [&repeat](auto & a) {
        Status ret = a.first->ApplySymmetriesToContractions();
        if (ret == SYM_EVEN) {
          repeat = true;
        } else if (ret == SYM_ODD) {
          a.second->Negate();
          repeat = true;
        } else if (ret != NO_ACTION) {
          assert(false);
        }
      });
  }
}

void Expression::CollectPrefactors() {
  if (summands->size() < 2) {
    return;
  }

  auto summands_new = std::make_unique<std::vector<std::pair<std::unique_ptr<MonomialExpression>, std::unique_ptr<ScalarSum>>>>();

  auto last = summands->begin();
  for (auto it = summands->begin() + 1; it < summands->end(); ++it) {
    if (*(last->first) == *(it->first)) {
      last->second->MergeWithOther(*(it->second));
      if (std::distance(it, summands->end()) == 1) {
        summands_new->push_back(std::pair<std::unique_ptr<MonomialExpression>, std::unique_ptr<ScalarSum>>());
        std::swap(*(--(summands_new->end())), *last);
      }
    } else {
        summands_new->push_back(std::pair<std::unique_ptr<MonomialExpression>, std::unique_ptr<ScalarSum>>());
        std::swap(*(--(summands_new->end())), *last);
      if (std::distance(it, summands->end()) == 1) {
        summands_new->push_back(std::pair<std::unique_ptr<MonomialExpression>, std::unique_ptr<ScalarSum>>());
        std::swap(*(--(summands_new->end())), *it);
      } else {
        last = it;
      }
    }
  }

  std::swap(summands, summands_new);

  CanonicalisePrefactors();
}

void Expression::CanonicalisePrefactors() {
  std::for_each(summands->begin(), summands->end(), [](auto & a) { a.second->Sort(); a.second->Collect(); });

  EliminateZeros();
}

void Expression::EliminateGamma() {
  bool repeat = true;

  while (repeat) {
    auto summands_new = std::make_unique<Sum>();
    repeat = false;
    std::transform(summands->cbegin(), summands->cend(), std::back_inserter(*summands_new),
      [&repeat](auto const & a) {
        auto summand_new = std::make_pair(std::make_unique<MonomialExpression>(*a.first), std::make_unique<ScalarSum>(*a.second));;
        Status ret = summand_new.first->EliminateGamma();
        if (ret == DELTA_OK ) {
          repeat = true;
        } else if (ret == DELTA_TO_TRACE) {
          summand_new.second->MultiplyCoefficient(Rational(3, 1));
          repeat = true;
        } else if (ret == ERROR) {
          assert(false);
        }
        return summand_new;
      });
    std::swap(summands, summands_new);
  }

}

void Expression::EliminateDelta() {
  bool repeat = true;

  while (repeat) {
    repeat = false;
    std::for_each(summands->begin(), summands->end(),
      [&repeat,this](auto & a) {
        Status ret = a.first->EliminateDelta();
        if (ret == DELTA_OK ) {
          repeat = true;
        } else if (ret == DELTA_TO_TRACE) {
          a.second->MultiplyCoefficient(Rational(dimension, 1));
          repeat = true;
        } else if (ret == ERROR) {
          assert(false);
        }
      });
  }

}

void Expression::EliminateEpsilon() {
  bool repeat = true;

  while (repeat) {
    repeat = false;
    std::for_each(summands->begin(), summands->end(),
      [&repeat,this](auto & a) {
        Status ret = a.first->EliminateEpsilon(dimension);
        if (ret == EPSILON_TO_ZERO) {
          a.second->MultiplyCoefficient(Rational(0, 1));
          repeat = true;
        } else if (ret == ERROR) {
          assert(false);
        }
      });
    EliminateZeros();
  }
}

void Expression::EliminateEpsilonI() {
  bool repeat = true;

  while (repeat) {
    repeat = false;
    std::for_each(summands->begin(), summands->end(),
      [&repeat](auto & a) {
        Status ret = a.first->EliminateEpsilonI();
        if (ret == EPSILONI) {
          repeat = true;
        } else if (ret != NO_ACTION) {
          assert(false);
        }
      });
    EliminateZeros();
  }
}

static bool Sort(std::vector<size_t> vec) {
  if (vec.empty()) {
    return true;
  }
  unsigned int swap_counter = 0;
  bool swapped;
  do {
    swapped = false;
    for (size_t counter = 0; counter < vec.size() - 1; ++counter) {
      if (vec[counter] > vec[counter + 1]) {
        std::swap(vec[counter], vec[counter + 1]);
        swapped = true;
        swap_counter += 1;
      }
    }
  } while (swapped);
  return (swap_counter % 2 == 0) ? true : false;
}

void Expression::EliminateEpsilonEpsilonI() {
  bool repeat = true;

  while (repeat) {
    // new vector of summands
    auto summands_new = std::make_unique<Sum>();
    // don't repeat unless there has been some action
    repeat = false;
    std::for_each(summands->begin(), summands->end(),
      [&repeat,&summands_new,this](auto const & a) {
        // eliminate epsilon and epsilonI in summand
        // return indices of epsilon and epsilonI
        auto indices_pair = a.first->EliminateEpsilonEpsilonI(dimension);

        // if there was no action (returned indices are of size zero),
        // push back old summand to new vector of summands
        //
        // else, for every permutation, push back the (dimension) deltas and
        // then what remains from the original summand (if there is any tensor left)
        if ( indices_pair.first.size() == 0 && indices_pair.second.size() == 0 && !a.first->IsZero()) {
          summands_new->push_back(std::make_pair(std::make_unique<MonomialExpression>(*a.first), std::make_unique<ScalarSum>(*a.second)));
        } else if ( indices_pair.first.size() == dimension && indices_pair.second.size() == dimension ) {
          Tensor delta (2, (dimension == 3 ? "gamma" : "delta"));
          delta.SetSymmetric();

          TensorMonomial dim_deltas;

          std::vector<size_t> index_numbers (dimension);
          std::generate(index_numbers.begin(), index_numbers.end(),
            [&delta,&dim_deltas,n=0] () mutable {
              dim_deltas.AddFactorRight(delta);
              return n++;
            });

          do {
            auto index_number_copy = index_numbers;
            auto const index_number_copy_const = index_numbers;
            bool sign = Sort(index_number_copy);

            ScalarSum coeff = *a.second;
            if (dimension == 4) {
              if (sign) {
                coeff.Negate(); // Attention! Additional minus sign because epsilon epsilonI = - delta delta delta delta
              }
            } else if (dimension == 3) {
              if (!sign) {
                coeff.Negate();
              }
            } else {
              assert(false);
            }

            std::vector<Index> dim_deltas_indices (2*dimension);

            std::for_each(index_number_copy_const.cbegin(), index_number_copy_const.cend(),
              [&dim_deltas_indices,&indices_pair,&index_numbers] (auto const & a) {
                dim_deltas_indices[2*a] = indices_pair.first.at(a);
                dim_deltas_indices[2*a+1] = indices_pair.second.at(index_numbers.at(a));
              });

            MonomialExpression dim_deltas_me (dim_deltas, Indices(dim_deltas_indices));
            MonomialExpression new_summand = dim_deltas_me.MultiplyOther(*a.first);

            summands_new->push_back(std::make_pair(std::make_unique<MonomialExpression>(new_summand), std::make_unique<ScalarSum>(coeff)));
          } while (std::next_permutation(index_numbers.begin(), index_numbers.end()));

//          repeat = true;

        } else {
          assert(false);
        }
      });
    std::swap(summands, summands_new);
  }
}

void Expression::EliminateTracefree() {
  std::for_each(summands->begin(), summands->end(),
    [](auto & a) {
      if (a.first->EliminateTracefree() == TRACE_TO_ZERO) {
        a.second->MultiplyCoefficient(Rational(0, 1));
      }
    });

  EliminateZeros();
}

void Expression::EliminateGammaGamma() {
  bool repeat = true;

  while (repeat) {
    repeat = false;
    std::for_each(summands->begin(), summands->end(),
      [&repeat](auto & a) {
        Status ret = a.first->EliminateGammaGamma();
        if (ret == GAMMA_GAMMA_TO_TRACE) {
          a.second->MultiplyCoefficient(Rational(3, 1));
          repeat = true;
        } else if (ret == GAMMA_GAMMA_TO_DELTA) {
          repeat = true;
        } else if (ret == ERROR) {
          assert(false);
        }
      });
  }
}

void Expression::EliminateEtaEta() {
  bool repeat = true;

  while (repeat) {
    repeat = false;
    std::for_each(summands->begin(), summands->end(),
      [&repeat](auto & a) {
        Status ret = a.first->EliminateEtaEta();
        if (ret == ETA_ETA_TO_TRACE) {
          a.second->MultiplyCoefficient(Rational(4, 1));
          repeat = true;
        } else if (ret == ETA_ETA_TO_DELTA) {
          repeat = true;
        } else if (ret == ERROR) {
          assert(false);
        }
      });
  }
}

void Expression::EliminateEtaPartial() {
  bool repeat = true;

  while (repeat) {
    repeat = false;
    std::for_each(summands->begin(), summands->end(),
      [&repeat](auto & a) {
        Status ret = a.first->EliminateEtaPartial();
        if (ret == ETA_PARTIAL_CONTRACTION) {
          repeat = true;
        } else if (ret != NO_ACTION) {
          assert(false);
        }
      });
  }
}

void Expression::EliminateEtaRankOne() {
  bool repeat = true;

  while (repeat) {
    repeat = false;
    std::for_each(summands->begin(), summands->end(),
      [&repeat](auto & a) {
        Status ret = a.first->EliminateEtaRankOne();
        if (ret == ETA_CONTRACTION) {
          repeat = true;
        } else if (ret != NO_ACTION) {
          assert(false);
        }
      });
  }
}

void Expression::EliminateZeros() {
  auto summands_new = std::make_unique<std::vector<std::pair<std::unique_ptr<MonomialExpression>, std::unique_ptr<ScalarSum>>>>();

  for (auto & summand : *summands) {
    if (summand.second->IsZero()) {
      continue;
    }
    summands_new->push_back(std::pair<std::unique_ptr<MonomialExpression>, std::unique_ptr<ScalarSum>>());
    std::swap(*(--(summands_new->end())), summand);
  }

  std::swap(summands, summands_new);
}

void Expression::SubstituteIndices(Indices const & indices_old, Indices const & indices_new) {
  std::for_each(summands->begin(), summands->end(),
    [&indices_old,&indices_new] (auto & a) {
      a.first->ExchangeIndices(indices_old,indices_new);
    });
}

void Expression::SubstituteFreeIndices(Indices const & indices_old, Indices const & indices_new) {
  std::for_each(summands->begin(), summands->end(),
    [&indices_old,&indices_new] (auto & a) {
      a.first->ExchangeFreeIndices(indices_old,indices_new);
    });
}

void Expression::SubstituteVariable(size_t const variable_old, ScalarSum const & scalar_sum_new) {
  std::for_each(summands->begin(), summands->end(),
    [variable_old,scalar_sum_new](auto & a) {
      a.second->SubstituteVariable(variable_old, scalar_sum_new);
    });
}

void Expression::SubstituteVariables(std::map<size_t, ScalarSum> substitution_map) {
  std::for_each(summands->begin(), summands->end(),
    [substitution_map](auto & a) {
      a.second->SubstituteVariables(substitution_map);
    });
}

void Expression::SubstituteVariables(std::map<size_t, size_t> substitution_map) {
  std::for_each(summands->begin(), summands->end(),
    [substitution_map](auto & a) {
      a.second->SubstituteVariables(substitution_map);
    });
}

void Expression::EliminateVariable(size_t const variable) {
  std::for_each(summands->begin(), summands->end(), [variable](auto & a) { a.second->EliminateVariable(variable); });
}

void Expression::SortSummands() {
  std::sort(summands->begin(), summands->end());
}

bool CompareSummandsReversed(Summand const & a, Summand const & b);

void Expression::SortSummandsByLastFactors() {
  std::sort(summands->begin(), summands->end(), CompareSummandsReversed);
}

void Expression::ExchangeSymmetrise(Indices const & indices1, Indices const & indices2, bool anti) {
  auto summands_new = std::make_unique<Sum>();

  for (auto & summand : *summands) {
    summand.second->DivideByTwo();
    Summand summand_copy;
    summand_copy.first = std::make_unique<MonomialExpression>(*summand.first);
    summand_copy.second = std::make_unique<ScalarSum>(*summand.second);
    summands_new->push_back(std::make_pair(std::make_unique<MonomialExpression>(), std::make_unique<ScalarSum>()));
    std::swap(*(--(summands_new->end())), summand);
    if (anti) {
      summand_copy.second->Negate();
    }
    summand_copy.first->ExchangeFreeIndices(indices1, indices2);
    summands_new->push_back(std::make_pair(std::make_unique<MonomialExpression>(), std::make_unique<ScalarSum>()));
    std::swap(*(--(summands_new->end())), summand_copy);
  }

  std::swap(summands, summands_new);
}

void Expression::SortSummandsByPrefactors() {
  std::sort(summands->begin(), summands->end(), [](auto & a, auto & b) {
    if (*a.second < *b.second) {
      return true;
    } else if (*b.second < *a.second) {
      return false;
    } else if (*a.first < *b.first) {
      return true;
    } else {
      return false;
    }
  });
}

void Expression::ThreePlusOne(std::vector<Index> indices_to_zero) {
  auto summands_new = std::make_unique<Sum>();

  std::for_each(summands->begin(), summands->end(),
    [&indices_to_zero,&summands_new] (auto & a) {
      auto status = a.first->ThreePlusOne(indices_to_zero);
      if (status == ZERO_ZERO) {
        return;
      } else if (status == ZERO_POSITIVE) {
        summands_new->emplace_back(std::make_unique<MonomialExpression>(*a.first), std::make_unique<ScalarSum>(*a.second));
      } else if (status == ZERO_NEGATIVE) {
        summands_new->emplace_back(std::make_unique<MonomialExpression>(*a.first), std::make_unique<ScalarSum>(*a.second));
        summands_new->back().second->Negate();
      } else {
        assert(false);
      }
    });

  dimension = 3;
  std::swap(summands_new, summands);
}

bool Expression::ContainsMonomial (MonomialExpression const & monexpr) const {
  bool ret = false;
  for (auto & summand : *summands) {
    if (*(summand.first) == monexpr) {
      ret = true;
      break;
    }
  }

  return ret;
}

void Expression::RedefineScalars() {
  if (summands->size() == 0) {
    return;
  }

  std::map<size_t, size_t> coefficient_map = GetCoefficientMap ();

  std::vector<std::vector<Rational>> prefactor_matrix = GetPrefactorMatrix ();

  MatrixXq mq(prefactor_matrix.size(), prefactor_matrix.front().size());

  for (size_t row_counter = 0; row_counter < prefactor_matrix.size(); ++row_counter) {
    for (size_t column_counter = 0; column_counter < prefactor_matrix.front().size(); ++column_counter) {
      Fraction frac = prefactor_matrix[row_counter][column_counter].get_fraction();
      mq(row_counter, column_counter) = mpq_class(frac.first, frac.second);
    }
  }

  Eigen::FullPivLU<MatrixXq> lu_decompq(mq);

  MatrixXq kq = lu_decompq.kernel();

  std::set<size_t> coeff_removed;

  for (int column_counter = 0; column_counter < kq.cols(); ++column_counter) {
    for (int row_counter = kq.rows() - 1; row_counter >= 0; --row_counter) {
      if (kq(row_counter, column_counter) == 0 ) {
        continue;
      } else if (std::find(coeff_removed.begin(), coeff_removed.end(), row_counter) != coeff_removed.end()) {
        continue;
      } else {
        coeff_removed.insert(row_counter);
        break;
      }
    }
  }

  std::map<size_t, size_t> coefficient_rmap;
  std::for_each(coefficient_map.begin(), coefficient_map.end(), [&coefficient_rmap](auto a) { coefficient_rmap.insert(std::make_pair(a.second, a.first)); });

  Expression ret (*this);
  std::for_each(coeff_removed.cbegin(), coeff_removed.cend(), [&ret, &coefficient_rmap, this] (auto a) { EliminateVariable(coefficient_rmap.at(a)); });

  CanonicalisePrefactors();

  std::map<size_t, size_t> substitution_map;

  std::for_each(coefficient_map.cbegin(), coefficient_map.cend(),
    [&substitution_map,&coeff_removed,n=1] (auto const & a) mutable {
      if (coeff_removed.find(a.second) == coeff_removed.end()) {
        substitution_map[a.first] = n++;
      }
    });

  SubstituteVariables(substitution_map);

/*
  auto summands_new = std::make_unique<Sum>();
  summands_new->push_back(std::make_pair(std::make_unique<MonomialExpression>(), std::make_unique<ScalarSum>(Scalar(1))));
  std::swap(summands->front().first, summands_new->back().first);

  size_t variable_counter = 0;

  auto it_last = summands->begin();

  for (auto it = summands->begin() + 1; it != summands->end(); ++it) {
    Rational coefficient(1,1);
    if (*(it->second) == *((it - 1)->second)) {
      coefficient = coefficient.DivideOther((it_last)->second->Ratio(*(it->second)));
    } else {
      ++variable_counter;
      it_last = it;
    }
    summands_new->push_back(std::make_pair(std::make_unique<MonomialExpression>(), std::make_unique<ScalarSum>(Scalar(coefficient, variable_counter + 1))));
    std::swap(it->first, summands_new->back().first);
  }

  std::swap(summands, summands_new);
*/
}

void Expression::RenameDummies() {
  std::for_each(summands->begin(), summands->end(),
    [](auto & a) {
      a.first->RenameDummies();
    });
}

void Expression::MultiplyCoefficient (Rational const & q) {
  std::for_each(summands->begin(), summands->end(),
    [q] (auto & a) {
      a.second->MultiplyCoefficient(q);
    });
}

ScalarSum Expression::EvaluateIndices(Indices const & indices, std::vector<size_t> const & numbers) const {
  std::map<Index, size_t> evaluation_map;

  for (size_t counter = 0; counter < indices.size(); ++counter) {
    evaluation_map[indices.at(counter)] = numbers.at(counter);
  }
  ScalarSum ret;
  std::for_each(summands->begin(), summands->end(),
    [&ret, &evaluation_map](auto & a) {
      ScalarSum new_summand(*(a.second));
      new_summand.MultiplyCoefficient(a.first->EvaluateIndices(evaluation_map));
      ret.MergeWithOther(new_summand);
    });

  ret.Sort();
  ret.Collect();

  return ret;
}

void Expression::AddOther(Expression const & other) {
  auto summands_new = std::make_unique<Sum>();

  summands_new->reserve(this->summands->size() + other.summands->size());
  std::for_each(summands->cbegin(), summands->cend(),
    [&summands_new] (auto const & a) {
      summands_new->emplace_back(std::make_pair(std::make_unique<MonomialExpression>(*a.first), std::make_unique<ScalarSum>(*a.second)));
    });
  std::for_each(other.summands->cbegin(), other.summands->cend(),
    [&summands_new] (auto const & a) {
      summands_new->emplace_back(std::make_pair(std::make_unique<MonomialExpression>(*a.first), std::make_unique<ScalarSum>(*a.second)));
    });

  assert(summands_new->size() == this->summands->size() + other.summands->size());

  std::swap(summands, summands_new);
}

void Expression::MultiplyOther(Expression const & other) {
  auto summands_new = std::make_unique<Sum>();

  if (!(IsZero() || other.IsZero())) {

  std::for_each(summands->cbegin(), summands->cend(),
    [&summands_new, &other] (auto const & a) {
      std::for_each(other.summands->cbegin(), other.summands->cend(),
        [&summands_new, &a] (auto const & b) {
          summands_new->push_back(std::make_pair(std::make_unique<MonomialExpression>(a.first->MultiplyOther(*(b.first))), std::make_unique<ScalarSum>(a.second->MultiplyOther(*(b.second)))));
        });
      });

  }

  assert (summands_new->size() == this->summands->size() * other.summands->size());

  std::swap(summands, summands_new);
}

bool Expression::IsZero() const { return summands->empty(); }

void Expression::SaveToFile(std::string const & filename) const {
  std::ofstream file;
  file.open(filename);
  boost::archive::text_oarchive oa {file};
  oa << *this;
}

void Expression::LoadFromFile(std::string const & filename) {
  std::ifstream file;
  file.open(filename);
  boost::archive::text_iarchive ia {file};
  ia >> *this;
}

std::string const Expression::GetLatexString(std::string base_name, bool upper) const {
  std::stringstream ss;

  bool first = true;

  if (summands->empty()) {
    return "0";
  }

  for (auto const & summand : *summands) {
    if (!first) {
      ss << std::endl;
    }
    ss << summand.second->ToString(base_name, !first) << " ";
    first = false;

    ss << summand.first->GetLatexString(upper);
  }

  return ss.str();
}

bool operator< (Summand const & lop, Summand const & rop) {
  if (*lop.first < *rop.first) {
    return true;
  } else if (*rop.first < *lop.first) {
    return false;
  } else if (*lop.second < *rop.second) {
    return true;
  } else {
    return false;
  }
}

bool CompareSummandsReversed (Summand const & lop, Summand const & rop) {
  if (MonomialExpression::CompareReversed(*lop.first, *rop.first)) {
    return true;
  } else if (MonomialExpression::CompareReversed(*rop.first, *lop.first)) {
    return false;
  } else if (*lop.second < *rop.second) {
    return true;
  } else {
    return false;
  }
}

ext: cpp

max_stars_repo_name: nilsalex/tensor-algebra

max_issues_repo_name: nilsalex/tensor-algebra

max_forks_repo_name: nilsalex/tensor-algebra


text:
/* =========================================================================
   Copyright (c) 2010-2016, Institute for Microelectronics,
                            Institute for Analysis and Scientific Computing,
                            TU Wien.
   Portions of this software are copyright by UChicago Argonne, LLC.

                            -----------------
                  ViennaCL - The Vienna Computing Library
                            -----------------

   Project Head:    Karl Rupp                   rupp@iue.tuwien.ac.at

   (A list of authors and contributors can be found in the PDF manual)

   License:         MIT (X11), see file LICENSE in the base directory
============================================================================= */

/** \example eigen-with-viennacl.cpp
*
*   This tutorial shows how data can be directly transferred from the Eigen Library to ViennaCL objects using the built-in convenience wrappers.
*
*   The first step is to include the necessary headers and activate the Eigen convenience functions in ViennaCL:
**/

// System headers
#include <iostream>

// Eigen headers
#include <Eigen/Core>
#include <Eigen/Sparse>

// IMPORTANT: Must be set prior to any ViennaCL includes if you want to use ViennaCL algorithms on Eigen objects
#define VIENNACL_WITH_EIGEN 1

// ViennaCL includes
#include "viennacl/vector.hpp"
#include "viennacl/matrix.hpp"
#include "viennacl/compressed_matrix.hpp"
#include "viennacl/linalg/prod.hpp"

// Helper functions for this tutorial:
#include "vector-io.hpp"

/**
*   The following is a set of auxiliary dispatchers for obtaining the right Eigen types for a given floating point type.
*   This is merely an implementation detail, so feel free to skip over it.
**/

//dense matrix:
template<typename T>
struct Eigen_dense_matrix
{
  typedef typename T::ERROR_NO_EIGEN_TYPE_AVAILABLE   error_type;
};

template<>
struct Eigen_dense_matrix<float>
{
  typedef Eigen::MatrixXf  type;
};

template<>
struct Eigen_dense_matrix<double>
{
  typedef Eigen::MatrixXd  type;
};

//sparse matrix
template<typename T>
struct Eigen_vector
{
  typedef typename T::ERROR_NO_EIGEN_TYPE_AVAILABLE   error_type;
};

template<>
struct Eigen_vector<float>
{
  typedef Eigen::VectorXf  type;
};

template<>
struct Eigen_vector<double>
{
  typedef Eigen::VectorXd  type;
};

/**
*    The following function contains the main code for this tutorial.
*    It consists of the following steps:
*      - Creates Eigen matrices and vectors
*      - Initializes them with data
*      - Create ViennaCL objects
*      - Copy them over to the respective ViennaCL objects
*      - Compute matrix-vector products in both Eigen and ViennaCL and compare results.
*
**/
template<typename ScalarType>
void run_tutorial()
{
  /**
  * Get Eigen matrix and vector types for the provided ScalarType.
  * Involves a little bit of template-metaprogramming.
  **/
  typedef typename Eigen_dense_matrix<ScalarType>::type  EigenMatrix;
  typedef typename Eigen_vector<ScalarType>::type        EigenVector;

  /**
  * Create and fill dense matrices from the Eigen library:
  **/
  EigenMatrix eigen_densemat(6, 5);
  EigenMatrix eigen_densemat2(6, 5);
  eigen_densemat(0,0) = 2.0;   eigen_densemat(0,1) = -1.0;
  eigen_densemat(1,0) = -1.0;  eigen_densemat(1,1) =  2.0;  eigen_densemat(1,2) = -1.0;
  eigen_densemat(2,1) = -1.0;  eigen_densemat(2,2) = -1.0;  eigen_densemat(2,3) = -1.0;
  eigen_densemat(3,2) = -1.0;  eigen_densemat(3,3) =  2.0;  eigen_densemat(3,4) = -1.0;
                               eigen_densemat(5,4) = -1.0;  eigen_densemat(4,4) = -1.0;
  Eigen::Map<EigenMatrix> eigen_densemat_map(eigen_densemat.data(), 6, 5); // same as eigen_densemat, but emulating user-provided buffer

  /**
  * Create and fill sparse matrices from the Eigen library:
  **/
  Eigen::SparseMatrix<ScalarType, Eigen::RowMajor> eigen_sparsemat(6, 5);
  Eigen::SparseMatrix<ScalarType, Eigen::RowMajor> eigen_sparsemat2(6, 5);
  eigen_sparsemat.reserve(5*2);
  eigen_sparsemat.insert(0,0) = 2.0;   eigen_sparsemat.insert(0,1) = -1.0;
  eigen_sparsemat.insert(1,1) = 2.0;   eigen_sparsemat.insert(1,2) = -1.0;
  eigen_sparsemat.insert(2,2) = -1.0;  eigen_sparsemat.insert(2,3) = -1.0;
  eigen_sparsemat.insert(3,3) = 2.0;   eigen_sparsemat.insert(3,4) = -1.0;
  eigen_sparsemat.insert(5,4) = -1.0;
  //eigen_sparsemat.endFill();

  /**
  * Create and fill a few vectors from the Eigen library:
  **/
  EigenVector eigen_rhs(5);
  Eigen::Map<EigenVector> eigen_rhs_map(eigen_rhs.data(), 5);
  EigenVector eigen_result(6);
  EigenVector eigen_temp(6);

  eigen_rhs(0) = 10.0;
  eigen_rhs(1) = 11.0;
  eigen_rhs(2) = 12.0;
  eigen_rhs(3) = 13.0;
  eigen_rhs(4) = 14.0;

  /**
  * Create the corresponding ViennaCL objects:
  **/
  viennacl::vector<ScalarType> vcl_rhs(5);
  viennacl::vector<ScalarType> vcl_result(6);
  viennacl::matrix<ScalarType> vcl_densemat(6, 5);
  viennacl::compressed_matrix<ScalarType> vcl_sparsemat(6, 5);

  /**
  * Directly copy the Eigen objects to ViennaCL objects
  **/
  viennacl::copy(&(eigen_rhs[0]), &(eigen_rhs[0]) + 5, vcl_rhs.begin());  // Method 1: via iterator interface (cf. std::copy())
  viennacl::copy(eigen_rhs, vcl_rhs);                                     // Method 2: via built-in wrappers (convenience layer)
  viennacl::copy(eigen_rhs_map, vcl_rhs);                                 // Same as method 2, but for a mapped vector

  viennacl::copy(eigen_densemat, vcl_densemat);
  viennacl::copy(eigen_densemat_map, vcl_densemat); //same as above, using mapped matrix
  viennacl::copy(eigen_sparsemat, vcl_sparsemat);
  std::cout << "VCL sparsematrix dimensions: " << vcl_sparsemat.size1() << ", " << vcl_sparsemat.size2() << std::endl;

  // For completeness: Copy matrices from ViennaCL back to Eigen:
  viennacl::copy(vcl_densemat, eigen_densemat2);
  viennacl::copy(vcl_sparsemat, eigen_sparsemat2);

  /**
  * Run dense matrix-vector products and compare results:
  **/
  eigen_result = eigen_densemat * eigen_rhs;
  vcl_result = viennacl::linalg::prod(vcl_densemat, vcl_rhs);
  viennacl::copy(vcl_result, eigen_temp);
  std::cout << "Difference for dense matrix-vector product: " << (eigen_result - eigen_temp).norm() << std::endl;
  std::cout << "Difference for dense matrix-vector product (Eigen->ViennaCL->Eigen): "
            << (eigen_densemat2 * eigen_rhs - eigen_temp).norm() << std::endl;

  /**
  * Run sparse matrix-vector products and compare results:
  **/
  eigen_result = eigen_sparsemat * eigen_rhs;
  vcl_result = viennacl::linalg::prod(vcl_sparsemat, vcl_rhs);
  viennacl::copy(vcl_result, eigen_temp);
  std::cout << "Difference for sparse matrix-vector product: " << (eigen_result - eigen_temp).norm() << std::endl;
  std::cout << "Difference for sparse matrix-vector product (Eigen->ViennaCL->Eigen): "
            << (eigen_sparsemat2 * eigen_rhs - eigen_temp).norm() << std::endl;
}

/**
*   In the main() routine we only call the worker function defined above with both single and double precision arithmetic.
**/
int main(int, char *[])
{
  std::cout << "----------------------------------------------" << std::endl;
  std::cout << "## Single precision" << std::endl;
  std::cout << "----------------------------------------------" << std::endl;
  run_tutorial<float>();

#ifdef VIENNACL_HAVE_OPENCL
  if ( viennacl::ocl::current_device().double_support() )
#endif
  {
    std::cout << "----------------------------------------------" << std::endl;
    std::cout << "## Double precision" << std::endl;
    std::cout << "----------------------------------------------" << std::endl;
    run_tutorial<double>();
  }

  /**
  *   That's it. Print a success message and exit.
  **/
  std::cout << std::endl;
  std::cout << "!!!! TUTORIAL COMPLETED SUCCESSFULLY !!!!" << std::endl;
  std::cout << std::endl;

}

max_stars_repo_name: yuchengs/viennacl-dev

max_issues_repo_name: yuchengs/viennacl-dev

max_forks_repo_name: yuchengs/viennacl-dev