text:
#include <rosbag/bag.h>
#include <rosbag/view.h>
#include <std_msgs/Int32.h>
#include <std_msgs/String.h>
#include <fstream>
#include <ros_control_boilerplate/MatchSpecificData.h>
#include <frc_msgs/PDPData.h>
#include <boost/lexical_cast.hpp>
#include <boost/foreach.hpp>
#define foreach BOOST_FOREACH

int main(int argc, char **argv)
{
	std::string bag_name = argv[1];

	std::string c_string = argv[2];
	int current_needed = boost::lexical_cast<int>(c_string);
	std::cout << current_needed;

	rosbag::Bag bag;
	bag.open(bag_name, rosbag::bagmode::Read);

	std::ofstream temp_file;
	temp_file.open ("temp_file.txt");

	std::ofstream times;
	times.open ("times.txt");

	std::vector<std::string> topics;
	topics.push_back(std::string("/frcrobot/pdp_states"));

	rosbag::View view(bag, rosbag::TopicQuery(topics));

	frc_msgs::PDPData::ConstPtr s = view[0].instantiate<frc_msgs::PDPData>();
	if (s != NULL)
		int start_time = floor(s->header.stamp.toSec());

	foreach(rosbag::MessageInstance const m, view)
	{
		frc_msgs::PDPData::ConstPtr s = m.instantiate<frc_msgs::PDPData>();
		if (s != NULL){
			temp_file << "current at " << current_needed << ": " <<  s->current[current_needed] << std::endl;
			if  (s->current[current_needed] > 0)
					times << "current at channel " << current_needed << " is significant at header " << floor(s->header.stamp.toSec()) - start_time << std::endl;
		}

	}
	bag.close();

	return 0;
}


text:
#ifndef STAN_MCMC_HMC_INTEGRATORS_IMPL_LEAPFROG_HPP
#define STAN_MCMC_HMC_INTEGRATORS_IMPL_LEAPFROG_HPP

#include <Eigen/Dense>
#include <stan/mcmc/hmc/integrators/base_leapfrog.hpp>

namespace stan {
  namespace mcmc {

    template <typename Hamiltonian>
    class impl_leapfrog: public base_leapfrog<Hamiltonian> {
    public:
      impl_leapfrog(): base_leapfrog<Hamiltonian>(),
                       max_num_fixed_point_(10),
                       fixed_point_threshold_(1e-8) {}

      void begin_update_p(typename Hamiltonian::PointType& z,
                          Hamiltonian& hamiltonian,
                          double epsilon,
                          callbacks::logger& logger) {
        hat_phi(z, hamiltonian, epsilon, logger);
        hat_tau(z, hamiltonian, epsilon, this->max_num_fixed_point_,
                logger);
      }

      void update_q(typename Hamiltonian::PointType& z,
                    Hamiltonian& hamiltonian,
                    double epsilon,
                    callbacks::logger& logger) {
        // hat{T} = dT/dp * d/dq
        Eigen::VectorXd q_init = z.q + 0.5 * epsilon * hamiltonian.dtau_dp(z);
        Eigen::VectorXd delta_q(z.q.size());

        for (int n = 0; n < this->max_num_fixed_point_; ++n) {
          delta_q = z.q;
          z.q.noalias() = q_init + 0.5 * epsilon * hamiltonian.dtau_dp(z);
          hamiltonian.update_metric(z, logger);

          delta_q -= z.q;
          if (delta_q.cwiseAbs().maxCoeff() < this->fixed_point_threshold_)
            break;
        }
        hamiltonian.update_gradients(z, logger);
      }

      void end_update_p(typename Hamiltonian::PointType& z,
                        Hamiltonian& hamiltonian,
                        double epsilon,
                        callbacks::logger& logger) {
        hat_tau(z, hamiltonian, epsilon, 1, logger);
        hat_phi(z, hamiltonian, epsilon, logger);
      }

      // hat{phi} = dphi/dq * d/dp
      void hat_phi(typename Hamiltonian::PointType& z,
                   Hamiltonian& hamiltonian,
                   double epsilon,
                   callbacks::logger& logger) {
        z.p -= epsilon * hamiltonian.dphi_dq(z, logger);
      }

      // hat{tau} = dtau/dq * d/dp
      void hat_tau(typename Hamiltonian::PointType& z,
                   Hamiltonian& hamiltonian,
                   double epsilon,
                   int num_fixed_point,
                   callbacks::logger& logger) {
        Eigen::VectorXd p_init = z.p;
        Eigen::VectorXd delta_p(z.p.size());

        for (int n = 0; n < num_fixed_point; ++n) {
          delta_p = z.p;
          z.p.noalias() = p_init - epsilon * hamiltonian.dtau_dq(z, logger);
          delta_p -= z.p;
          if (delta_p.cwiseAbs().maxCoeff() < this->fixed_point_threshold_)
            break;
        }
      }

      int max_num_fixed_point() {
        return this->max_num_fixed_point_;
      }

      void set_max_num_fixed_point(int n) {
        if (n > 0) this->max_num_fixed_point_ = n;
      }

      double fixed_point_threshold() {
        return this->fixed_point_threshold_;
      }

      void set_fixed_point_threshold(double t) {
        if (t > 0) this->fixed_point_threshold_ = t;
      }

    private:
      int max_num_fixed_point_;
      double fixed_point_threshold_;
    };

  }  // mcmc
}  // stan

#endif

ext: hpp

max_stars_repo_name: drezap/stan

max_issues_repo_name: drezap/stan

max_forks_repo_name: drezap/stan


text:
/*!
 * @file     mppi_new.cpp
 * @author   Giuseppe Rizzi
 * @date     01.07.2020
 * @version  1.0
 * @brief    description
 */

#include <algorithm>

#include <Eigen/Core>
#include <chrono>
#include <iomanip>
#include <iostream>
#include <map>
#include <random>
#include <vector>

#include "mppi/core/config.h"
#include "mppi/core/cost.h"
#include "mppi/core/dynamics.h"
#include "mppi/core/rollout.h"
#include "mppi/core/solver.h"

#include "mppi/utils/logging.h"
#include "mppi/utils/savgol_filter.h"

#ifdef SIGNAL_LOGGER
#include <signal_logger/signal_logger.hpp>
#endif

namespace mppi {

Solver::Solver(dynamics_ptr dynamics, cost_ptr cost, policy_ptr policy,
               const Config& config)
    : dynamics_(std::move(dynamics)),
      cost_(std::move(cost)),
      policy_(std::move(policy)),
      config_(config) {
  init_data();
  init_threading();

  if (config_.display_update_freq) {
    start_time_ = std::chrono::high_resolution_clock::now();
  }

#ifdef SIGNAL_LOGGER
  signal_logger::add(min_cost_, "solver/rollouts/min_cost");
  signal_logger::add(max_cost_, "solver/rollouts/max_cost");
  signal_logger::add(rollouts_cost_, "solver/rollouts/costs");
  signal_logger::add(delay_steps_, "solver/delay_steps");
  signal_logger::add(rate_, "solver/rate");
  signal_logger::logger->updateLogger();
#endif
}

void Solver::init_data() {
  // TODO(giuseppe) this should be automatically computed when config is parsed
  steps_ = static_cast<int>(std::ceil(config_.horizon / config_.step_size));
  nx_ = dynamics_->get_state_dimension();
  nu_ = dynamics_->get_input_dimension();

  opt_roll_ = Rollout(steps_, nu_, nx_);
  opt_roll_cache_ = Rollout(steps_, nu_, nx_);
  nominal_.setZero(steps_, nu_);

  weights_.resize(config_.rollouts, 1.0 / config_.rollouts);
  rollouts_.resize(config_.rollouts, Rollout(steps_, nu_, nx_));
  cached_rollouts_ = std::ceil(config_.caching_factor * config_.rollouts);

  delay_steps_ = 0;
  step_count_ = 0;
  observation_set_ = false;
  reference_set_ = false;
}

void Solver::init_filter() {}

void Solver::init_threading() {
  if (config_.threads > 1) {
    std::cout << "Using multithreading. Number of threads: " << config_.threads
              << std::endl;
    pool_ = std::make_unique<ThreadPool>(config_.threads);
    futures_.resize(config_.threads);

    for (size_t i = 0; i < config_.threads; i++) {
      dynamics_v_.push_back(dynamics_->create());
      cost_v_.push_back(cost_->create());
    }
  }
}

void Solver::update_policy() {
  if (!observation_set_) {
    log_warning_throttle(1.0,
                         "Observation has never been set. Dropping update");
  } else if (!reference_set_) {
    log_warning_throttle(1.0, "Reference has never been set. Dropping update");
  } else {
    auto start = std::chrono::steady_clock::now();
    update_delay();
    copy_observation();

    for (size_t i = 0; i < config_.substeps; i++) {
      prepare_rollouts();
      update_reference();
      sample_trajectories();
      optimize();
      filter_input();

      stage_cost_ =
          cost_->get_stage_cost(x0_internal_, opt_roll_.uu[0], t0_internal_);
    }
    swap_policies();
    auto end = std::chrono::steady_clock::now();
    rate_ = std::chrono::duration_cast<std::chrono::nanoseconds>(end - start).count() / 1e9;
  }
}

void Solver::time_it() {
  auto stop = std::chrono::high_resolution_clock::now();
  auto duration =
      std::chrono::duration_cast<std::chrono::microseconds>(stop - start_time_)
          .count();
  auto timing_string =
      "Time since optimization start: " + std::to_string(duration) + "[μs], " +
      std::to_string(duration / 1000000.0) + "[s], " +
      std::to_string(1 / (duration / 1000000.0)) + "[Hz]";
  log_info(timing_string);
  start_time_ = std::chrono::high_resolution_clock::now();
}

void Solver::set_observation(const observation_t& x, const double t) {
  {
    std::unique_lock<std::shared_mutex> lock(state_mutex_);
    x0_ = x;
    reset_time_ = t;
  }

  // initialization of rollouts data
  if (first_step_) {
    copy_observation();
    initialize_rollouts();
    first_step_ = false;
  }

  observation_set_ = true;
}

void Solver::update_delay() {
  delay_steps_ = std::ceil((reset_time_ - t0_internal_) / config_.step_size);
  policy_->update_delay(delay_steps_);
}

void Solver::copy_observation() {
  std::shared_lock<std::shared_mutex> lock(state_mutex_);
  x0_internal_ = x0_;
  t0_internal_ = reset_time_;
}

void Solver::initialize_rollouts() {
  std::shared_lock<std::shared_mutex> lock_state(state_mutex_);
  opt_roll_.clear();
  std::fill(opt_roll_.uu.begin(), opt_roll_.uu.end(),
            dynamics_->get_zero_input(x0_internal_));

  std::shared_lock<std::shared_mutex> lock(rollout_cache_mutex_);
  std::fill(opt_roll_cache_.xx.begin(), opt_roll_cache_.xx.end(), x0_internal_);
  std::fill(opt_roll_cache_.uu.begin(), opt_roll_cache_.uu.end(),
            dynamics_->get_zero_input(x0_internal_));
  for (int i = 0; i < steps_; i++) {
    opt_roll_cache_.tt[i] = t0_internal_ + config_.step_size * i;
  }
}

void Solver::prepare_rollouts() {
  // cleanup
  rollouts_cost_.clear();
  for (auto& roll : rollouts_) {
    roll.clear_cost();
    roll.valid = true;
  }
}

void Solver::set_reference_trajectory(mppi::reference_trajectory_t& ref) {
  if (ref.rr.size() != ref.tt.size()) {
    std::stringstream error;
    error << "The reference trajectory state and time dimensions do not match: "
          << ref.rr.size() << " != " << ref.tt.size();
    throw std::runtime_error(error.str());
  }
  std::unique_lock<std::shared_mutex> lock(reference_mutex_);
  rr_tt_ref_ = ref;
  reference_set_ = true;
}

void Solver::update_reference() {
  std::shared_lock<std::shared_mutex> lock(reference_mutex_);
  cost_->set_reference_trajectory(rr_tt_ref_);

  if (config_.threads > 1) {
    for (auto& cost : cost_v_) cost->set_reference_trajectory(rr_tt_ref_);
  }

}

void Solver::sample_noise(input_t& noise) { sampler_->get_sample(noise); }

void Solver::sample_trajectories_batch(dynamics_ptr& dynamics, cost_ptr& cost,
                                       const size_t start_idx,
                                       const size_t end_idx) {
  observation_t x;
  for (size_t k = start_idx; k < end_idx; k++) {
    dynamics->reset(x0_internal_, t0_internal_);
    x = x0_internal_;
    double ts;
    for (int t = 0; t < steps_; t++) {
      ts = t0_internal_ + t * config_.step_size;
      rollouts_[k].tt[t] = ts;
      rollouts_[k].uu[t] = policy_->sample(ts, k);

      // compute input-state stage cost
      double cost_temp;
      cost_temp = std::pow(config_.discount_factor, t) *
                  cost->get_stage_cost(x, rollouts_[k].uu[t], ts);

      if (std::isnan(cost_temp)) {
        std::stringstream ss;
        ss << "Something went wrong ... dynamics diverged?\n" << std::endl;
        ss << "Rollout#" << k << std::endl;
        ss << "Time  where diverged: " << t << std::endl;
        ss << "State where diverged: " << x.transpose() << std::endl;
        ss << "Previous state: " << rollouts_[k].xx[t - 1].transpose()
           << std::endl;
        ss << "Input where diverged: " << rollouts_[k].uu[t].transpose();
        ss << "Previous input: " << rollouts_[k].uu[t - 1].transpose()
           << std::endl;
        rollouts_[k].valid = false;
        rollouts_[k].total_cost = std::numeric_limits<double>::infinity();
        log_warning(ss.str());
        break;
      }

      // store data
      rollouts_[k].xx[t] = x;
      rollouts_[k].cc(t) = cost_temp;
      rollouts_[k].total_cost += cost_temp;

      // integrate dynamics    auto start = std::chrono::steady_clock::now();
      x = dynamics->step(rollouts_[k].uu[t], config_.step_size);
    }
  }
}

void Solver::sample_trajectories() {
  policy_->shift(t0_internal_);
  policy_->update_samples(weights_, cached_rollouts_);

  if (config_.threads == 1) {
    sample_trajectories_batch(dynamics_, cost_, 0, config_.rollouts);
  } else {
    for (size_t i = 0; i < config_.threads; i++) {
      futures_[i] = pool_->enqueue(
          std::bind(&Solver::sample_trajectories_batch, this,
                    std::placeholders::_1, std::placeholders::_2,
                    std::placeholders::_3, std::placeholders::_4),
          dynamics_v_[i], cost_v_[i],
          (size_t)i * config_.rollouts / config_.threads,
          (size_t)(i + 1) * config_.rollouts / config_.threads);
    }

    for (size_t i = 0; i < config_.threads; i++) futures_[i].get();
  }
}

void Solver::compute_weights() {
  // keep all non diverged rollouts
  double min_cost = std::numeric_limits<double>::max();
  double max_cost = -min_cost;

  for (size_t k = 0; k < config_.rollouts; k++) {
    if (rollouts_[k].valid) {
      const double& cost = rollouts_[k].total_cost;
      min_cost = (cost < min_cost) ? cost : min_cost;
      max_cost = (cost > max_cost) ? cost : max_cost;
    }
  }

  min_cost_ = min_cost;   //*std::min_element(rollouts_cost_.begin(), rollouts_cost_.end());
  max_cost_ = max_cost;   //*std::max_element(rollouts_cost_.begin(), rollouts_cost_.end());

  double sum = 0.0;
  for (size_t k = 0; k < config_.rollouts; k++) {
    double modified_cost = config_.h * (rollouts_[k].total_cost - min_cost_) /
                           (max_cost_ - min_cost_);

    weights_[k] = rollouts_[k].valid ? std::exp(-modified_cost) : 0.0;
    sum += weights_[k];
  }
  std::transform(weights_.begin(), weights_.end(), weights_.begin(),
                 [&sum](double v) -> double { return v / sum; });
}

void Solver::optimize() {
  // get new rollouts weights
  compute_weights();

  // update policy according to new weights
  policy_->update(weights_, config_.alpha);

  // retrieve the nominal policy for each time step
  for (int t = 0; t < steps_; t++) {
    opt_roll_.tt[t] = t0_internal_ + t * config_.step_size;
    opt_roll_.uu[t] = policy_->nominal(t0_internal_ + t * config_.step_size);
  }

}

void Solver::filter_input() {
  // only if filter is available reset to the initial opt time
  if (filter_) {
    filter_->reset(x0_internal_, t0_internal_);
  }

  // reset the dynamics such that we rollout the dynamics again
  // with the input we filter step after step (sequentially)
  dynamics_->reset(x0_internal_, t0_internal_);
  opt_roll_.xx[0] = x0_internal_;

  // sequential filtering otherwise just nominal rollout
  for (int t = 0; t < steps_ - 1; t++) {
    if (filter_) {
      filter_->apply(opt_roll_.xx[t], opt_roll_.uu[t], opt_roll_.tt[t]);
      nominal_.row(t) = opt_roll_.uu[t].transpose();
    }
    opt_roll_.xx[t + 1] = dynamics_->step(opt_roll_.uu[t], config_.step_size);
  }

  // filter the last input in the sequence (if filter is available)
  if (filter_) {
    filter_->apply(opt_roll_.xx.back(), opt_roll_.uu.back(),
                   opt_roll_.tt.back());
    nominal_.bottomRows(1) = opt_roll_.uu.back().transpose();
    policy_->set_nominal(nominal_);
  }
}

void Solver::get_input(const observation_t& x, input_t& u, const double t) {
  static double coeff;
  static size_t idx;
  {
    std::shared_lock<std::shared_mutex> lock(rollout_cache_mutex_);
    if (t < opt_roll_cache_.tt.front()) {
      std::stringstream warning;
      warning << "Queried time " << t << " smaller than first available time "
              << opt_roll_cache_.tt.front();
      log_warning_throttle(1.0, warning.str());
      u = dynamics_->get_zero_input(x);
    }

    auto lower = std::lower_bound(opt_roll_cache_.tt.begin(),
                                  opt_roll_cache_.tt.end(), t);
    if (lower == opt_roll_cache_.tt.end()) {
      std::stringstream warning;
      warning << "Queried time " << t << " larger than last available time "
              << opt_roll_cache_.tt.back();
      log_warning_throttle(1.0, warning.str());
      u = opt_roll_cache_.uu.back();
      return;
    }

    idx = std::distance(opt_roll_cache_.tt.begin(), lower);
    // first index (time)
    if (idx == 0) {
      u = opt_roll_cache_.uu.front();
    }
    // last index (time larget then last step)
    else if (idx > opt_roll_cache_.steps_) {
      u = opt_roll_cache_.uu.back();
    }
    // interpolate
    else {
      coeff = (t - *(lower - 1)) / (*lower - *(lower - 1));
      u = (1 - coeff) * opt_roll_cache_.uu[idx - 1] +
          coeff * opt_roll_cache_.uu[idx];
    }
  }
}

void Solver::get_diagonal_variance(Eigen::VectorXd& var) const {
  var = sampler_->sigma().diagonal();
}

void Solver::get_input_state(const observation_t& x, observation_t& x_nom,
                             input_t& u_nom, const double t) {
  {
    std::shared_lock<std::shared_mutex> lock(rollout_cache_mutex_);
    if (t < opt_roll_cache_.tt.front()) {
      std::stringstream warning;
      warning << "Queried time " << t << " smaller than first available time "
              << opt_roll_cache_.tt.front();
      log_warning_throttle(1.0, warning.str());
      x_nom = opt_roll_cache_.xx.front();
      u_nom = dynamics_->get_zero_input(x);
      return;
    }

    auto lower = std::lower_bound(opt_roll_cache_.tt.begin(),
                                  opt_roll_cache_.tt.end(), t);
    if (lower == opt_roll_cache_.tt.end()) {
      std::stringstream warning;
      warning << "Queried time " << t << " larger than last available time "
              << opt_roll_cache_.tt.back();
      log_warning_throttle(1.0, warning.str());
      x_nom = opt_roll_cache_.xx.back();
      u_nom = dynamics_->get_zero_input(x);
      return;
    }

    size_t idx = std::distance(opt_roll_cache_.tt.begin(), lower);

    // first
    if (idx == 0) {
      x_nom = opt_roll_cache_.xx.front();
      u_nom = opt_roll_cache_.uu.front();
    }
    // last
    else if (idx > opt_roll_cache_.steps_) {
      x_nom = opt_roll_cache_.xx.back();
      u_nom = opt_roll_cache_.uu.back();
    }
    // interpolate
    else {
      double coeff = (t - *(lower - 1)) / (*lower - *(lower - 1));
      u_nom = (1 - coeff) * opt_roll_cache_.uu[idx - 1] +
              coeff * opt_roll_cache_.uu[idx];
      x_nom = opt_roll_cache_.xx[idx];  // TODO offer a way to also do
                                        // interpolation of the state
    }
  }
}

bool Solver::get_optimal_rollout(observation_array_t& xx, input_array_t& uu) {
  std::shared_lock<std::shared_mutex> lock(rollout_cache_mutex_);
  auto lower = std::lower_bound(opt_roll_cache_.tt.begin(),
                                opt_roll_cache_.tt.end(), reset_time_);
  if (lower == opt_roll_cache_.tt.end()) return false;
  size_t offset = std::distance(opt_roll_cache_.tt.begin(), lower);

  // fill with portion of vector starting from current time
  xx = observation_array_t(opt_roll_cache_.xx.begin() + offset,
                           opt_roll_cache_.xx.end());
  uu = input_array_t(opt_roll_cache_.uu.begin() + offset,
                     opt_roll_cache_.uu.end());
  return true;
}

void Solver::get_optimal_rollout(Rollout& r) {
  std::shared_lock<std::shared_mutex> lock(rollout_cache_mutex_);

  // fill with portion of vector starting from current reset time
  r.tt = time_array_t(opt_roll_cache_.tt.begin(),
                      opt_roll_cache_.tt.end());

  r.xx = observation_array_t(opt_roll_cache_.xx.begin(),
                             opt_roll_cache_.xx.end());

  r.uu = input_array_t(opt_roll_cache_.uu.begin(),
                       opt_roll_cache_.uu.end());
}

void Solver::swap_policies() {
  std::unique_lock<std::shared_mutex> lock(rollout_cache_mutex_);
  opt_roll_cache_ = opt_roll_;
}

template <typename T>
void Solver::shift_back(std::vector<T>& v_out, const T& fill,
                        const int offset) {
  std::rotate(v_out.begin(), v_out.begin() + offset, v_out.end());
  std::fill(v_out.end() - offset, v_out.end(), fill);
}

// explicit instantiation
template void Solver::shift_back<Eigen::VectorXd>(
    std::vector<Eigen::VectorXd>& v_out, const Eigen::VectorXd& fill,
    const int offset);

void Solver::print_cost_histogram() const {
  constexpr int nbins = 10;
  double delta = (max_cost_ - min_cost_) / nbins;

  std::cout << "Rollouts cost histogram " << std::endl;
  std::cout << "Max: " << max_cost_ << ", Min: " << min_cost_
            << ", delta: " << delta << std::endl;
  std::map<int, int> hist{};
  for (const auto& cost : rollouts_cost_) {
    ++hist[std::round((cost - min_cost_) / delta)];
  }
  for (auto p : hist) {
    std::cout << std::setw(2) << p.first << ' ' << std::string(p.second, '*')
              << '\n';
  }
  std::cout << std::endl << std::endl;
}

}  // namespace mppi

ext: cpp


text:
/** This is free and unencumbered software released into the public domain.

The authors of ISIS do not claim copyright on the contents of this file.
For more details about the LICENSE terms and the AUTHORS, you will
find files of those names at the top level of this repository. **/

/* SPDX-License-Identifier: CC0-1.0 */

#include "SumFile.h"

#include <cmath>
#include <numeric>
#include <string>
#include <vector>

#include <QDebug>
#include <QFile>
#include <QTextStream>
#include <QtGlobal>

// boost library
#include <boost/foreach.hpp>

#include "Camera.h"
#include "IException.h"
#include "IString.h"
#include "Kernels.h"
#include "NaifStatus.h"
#include "Progress.h"

using namespace std;

namespace Isis {

  /**
   *  Default constructor creates empty SumFile object.
   */
  SumFile::SumFile()  { }

  /**
   * @brief Constructor reads and parses contents of SUMFILE
   *
   * @param sumFile Name of SUMFILE to read/parse
   */
  SumFile::SumFile(const QString &sumFile) {
    parseSumFile(sumFile);
  }

//  /**
//   *  Tests vailidty of the SumFile point.
//   *
//   *  @return @b bool
//   */
//  bool SumFile::isValid() const {
//    return (true);
//  }

  /**
   *  Returns the name of the SUMFILE (contents of first line).
   *
   *  @return @b QString Name of sumfile identifier
   */
  QString SumFile::name() const {
    return (m_id);

  }

  /**
   *  Returns the time found in the SUMFILE in UTC form
   *
   *  @return @b QString Time found in SUMFILE
   */
  QString SumFile::utc() const {
    return (m_obsTime.UTC());

  }

  /**
   *  Returns the SUMFILE time, in ET.
   *
   *  @return @b double Time converted to ET
   */
  double  SumFile::et() const {
    return (m_obsTime.Et());
  }

/**
 * @brief Returns the Time object found in SUMFILE
 *
 * @return const iTime& Time in SUMFILE
 */
  const iTime &SumFile::time() const {
    return ( m_obsTime );
  }

/**
 * @brief Update SPICE data in ISIS cube
 *
 * This method will update the SPICE blobs with the contents of the SUMSPICE.
 * The contents of the InstrumentPosition will be replaced with SUMSPICE SZ
 * vector and InstrumentPointing is replaced by the CX, CY, CZ matrix.
 *
 * @param cube    An intialized ISIS Cube object
 * @param camera  An option pointer to the camera. If not provided, it is
 *                retrieved from the Cube object
 *
 * @return bool   True if succesful, false if the operation fails
 */
  bool SumFile::updateSpice(Cube &cube, Camera *camera) const {
    bool good = updatePointing(cube, camera);
    if ( good ) {
     good = updatePosition(cube, camera);
    }
    return (good);
  }

  /**
   * @brief Convert the SUMFILE contents and apply pointing to ISIS camera
   *
   * It is up to the caller to ensure appropriate kernels are loaded. See
   * loadKernels();
   *
   * @param cube    Cube to update pointing for
   * @param camera  Camera object associated with cube for update
   *
   * @return bool  True if successful, false otherwise
   */
  bool SumFile::updatePointing(Cube &cube, Camera *camera) const {

    Camera *mycam = ( camera != NULL ) ? camera : cube.camera();

    // get new pointing quaternion from sum file

    // first, we get the rotattion between j2000 and target (i.e. body fixed frame)
    SpiceRotation *body = mycam->bodyRotation();
    Quaternion j2000ToTarget(body->Matrix());

    // next, get the constant rotation for the camera from the cube's table
    SpiceRotation *oldRotation = mycam->instrumentRotation();
    Quaternion oldConstantRotation(oldRotation->ConstantRotation()); // old TC rotation

    // Get the new rotation from the sum file.
    // this is the instrument frame (relative to the target, i.e. body-fixed)
    Quaternion newRotation = getPointing();

    // new target to instrument rotation (TC) is found by
    //     TC = inverse(rotation from sum file) * (old TC rotation)
    //
    // So, the new instrument frame relative to J2000 is
    //     CJ = inverse(TC) * TJ
    Quaternion newTimeBasedRotation = oldConstantRotation.Conjugate() * newRotation * j2000ToTarget; // new CJ

    Table table = oldRotation->Cache("InstrumentPointing");
    if ( table.Records() > 1 ) {
      QString message = "Expected/support only one InstrumentPointing record "
                     "(i.e., Framing camera) but got " +
                     toString(table.Records()) + " for file " + cube.fileName();
      throw IException(IException::User, message, _FILEINFO_);
    }

    TableRecord &rec0 = table[0];
    if ( rec0.Fields() < 5 ) {
      QString message = "Expected/support quaternion Table containing 5 or more fields "
                     "(i.e., Framing camera) but got " +
                     toString(rec0.Fields()) + " for file " + cube.fileName();
      throw IException(IException::User, message, _FILEINFO_);
    }

    for (int i = 0; i < 4; i++) {
      rec0[i] = newTimeBasedRotation[i];
    }

    // Put a comment in the record
    QString message = "Updated pointing with SUMFILE " + name() + " on " + iTime::CurrentLocalTime();
    PvlKeyword sumFileKeyword("SUMFILE", name());
    sumFileKeyword.addComment(message);
    table.Label().addKeyword(sumFileKeyword, PvlContainer::Replace);
    table.Update(rec0, 0);

    cube.write(table);

    // Delete polygon if it has one
    if ( cube.label()->hasObject("Polygon") ) {
      cube.label()->deleteObject("Polygon");
    }

    return (true);
  }

  /**
   * @brief Update spacecraft position with SUMFILE data
   *
   * @param cube    Cube to update position for
   * @param camera  Camera object associated with cube for update
   *
   * @return bool   True if successful, false otherwise
   */
  bool SumFile::updatePosition(Cube &cube, Camera *camera) const {

    Camera *mycam = ( !camera ) ? cube.camera() : camera;

    SpiceRotation *body = mycam->bodyRotation();
    Quaternion j2000ToTarget(body->Matrix());

    // Get body-fixed s/c vector - points from body to s/c
    vector<double> spacecraftPos = getPosition();

    // Have vector point from s/c to body
    vminus_c(&spacecraftPos[0], &spacecraftPos[0]);
    // Transform position vector from body-fixed to J2000
    spacecraftPos = j2000ToTarget.Conjugate().Qxv(spacecraftPos);

    // Retrieve s/c cache and update the position
    SpicePosition *oldPosition = mycam->instrumentPosition();
    Table table = oldPosition->Cache("InstrumentPosition");
    if ( table.Records() > 1 ) {
      QString message = "Expected/support only one InstrumentPosition record "
                     "(i.e., Framing camera) but got " +
                     toString(table.Records()) + " for file " + cube.fileName();
      throw IException(IException::User, message, _FILEINFO_);
    }

    TableRecord &rec0 = table[0];
    if ( rec0.Fields() < 4 ) {
      QString message = "Expected/support vector Table containing 4 or more fields "
                     "(i.e., Framing camera) but got " +
                     toString(rec0.Fields()) + " for file " + cube.fileName();
      throw IException(IException::User, message, _FILEINFO_);
    }

    // Now replace the quaternions in table
    for (unsigned int i = 0; i < spacecraftPos.size(); i++) {
      rec0[i] = spacecraftPos[i];
    }

    // Put a comment in the record
    QString message = "Updated position with SUMFILE " + name() + " on " + iTime::CurrentLocalTime();
    PvlKeyword sumFileKeyword("SUMFILE", name());
    sumFileKeyword.addComment(message);
    table.Label().addKeyword(sumFileKeyword, PvlContainer::Replace);
    table.Update(rec0, 0);

    cube.write(table);

    // Delete polygon if it has one
    if ( cube.label()->hasObject("Polygon") ) {
      cube.label()->deleteObject("Polygon");
    }

    return (true);

   }

  /**
   * @brief Return the point matrix, in body-fixed format.
   *
   * This method will return a quaternion, in body fixed coordinates,
   * corresponding to the 3x3 instrument pointing matrix found in the SUMFILE.
   * Note that there is no need to load kernels (or have transformations
   * determined) to provide this data.
   *
   * @return Quaternion The SUMFILE's pointing matrix
   */
  Quaternion SumFile::getPointing() const {
    vector<SpiceDouble> cmatrixBodyFixed(9);

    // copy the transposed matrix to a 9 unit vector to be converted to quaternion
    for (int i = 0; i < 3; i++) {
      for (int j = 0; j < 3; j++) {
        cmatrixBodyFixed[3*i+j] = m_pointingMatrix[i][j];// ji = instrument to target
      }
    }

    Quaternion quaternionBodyFixed(cmatrixBodyFixed);
    return quaternionBodyFixed;
  }

  /**
   * @brief Get spacecraft postion in body-fixed coordinates
   *
   * This method will return the vector found in the SUMFILE that represents
   * the coordinates of the spacecraft, in body fixed position.
   *
   * @return std::vector<double> The SUMFILE's spacecraft position
   */
  std::vector<double> SumFile::getPosition()  const{
    std::vector<double> spacecraftPosition;
    spacecraftPosition.reserve(3);
    spacecraftPosition.push_back(m_spacecraftPosition[0]);
    spacecraftPosition.push_back(m_spacecraftPosition[1]);
    spacecraftPosition.push_back(m_spacecraftPosition[2]);
    return (spacecraftPosition);
  }

  /**
   * @brief Get Sun postion in body-fixed coordinates
   *
   * This method will return the vector found in the SUMFILE that represents
   * the coordinates of the Sun, in body fixed position.
   *
   * @return std::vector<double> The SUMFILE's Sun position
   */
  std::vector<double> SumFile::getSunPosition()  const{
    std::vector<double> sunPosition;
    sunPosition.reserve(3);
    sunPosition.push_back(m_sunPosition[0]);
    sunPosition.push_back(m_sunPosition[1]);
    sunPosition.push_back(m_sunPosition[2]);
    return (sunPosition);
  }

  /**
   * Adds information about the SUMFILE to the output stream: name, start time
   * (UTC), start time (ET), number of lines, number of samples, DN minimum, DN
   * maximum.
   *
   * @param outs The output Stream to be filled.
   *
   * @return @b std::ostream& The filled output stream.
   */
  std::ostream &SumFile::brief(std::ostream &outs) const {
    outs << name() << ", " << utc() << ", "  << std::fixed << et() << ", "
         << m_numLines << ", " << m_numSamples
         << ", " << m_dnMin << ", " << m_dnMax;
    return (outs);
  }

  /**
   * @brief Opens and parses the contents of a SUMFILE.
   *
   * THis method will open a Gaskell SPC SUMFILE and parse the contents. The
   * SUMFILE contents are internalized in the SumFile object.
   *
   * @param sumFile Name of a SPC SUMFILE
   */
  void SumFile::parseSumFile(const QString &sumFile) {
    FileName sfile(sumFile);
    QFile sumF(sfile.expanded());
    if ( !sumF.open(QIODevice::ReadOnly | QIODevice::Text) ) {
      QString message = "Failed to open FROM file \"" + sfile.original() + "\"";
      throw IException(IException::User, message, _FILEINFO_);
    }

    // Open the text stream
    QTextStream sumin(&sumF);

    // Get the image id
    QStringList values = getSumLine(sumin.readLine(), 1);
    m_id  = values.first();

    // Get the time. We expect spaces in the time specification
    values = getSumLine(sumin.readLine(), 4);
    m_obsTime = iTime(values.join("-"));

    // Get image lines/samples and lower/upper dn boundaries
    values = getSumLine(sumin.readLine(), 7, "THRSH");
    m_numSamples = cvtDouble(values[0]);
    m_numLines = cvtDouble(values[1]);
    m_dnMin  = cvtDouble(values[2]);
    m_dnMax  = cvtDouble(values[3]);

    // Get pixel size and boresight line/sample
    values = getSumLine(sumin.readLine(), 5, "CTR");
    m_pxlSize = cvtDouble(values[0]);
    m_centerSample   = cvtDouble(values[1]);
    m_centerLine   = cvtDouble(values[2]);

    // Get spacecraft position
    values = getSumLine(sumin.readLine(), 4, "SCOBJ");
    for (int i = 0; i < 3; i++) {
      m_spacecraftPosition[i] = cvtDouble(values[i]);
    }

    // Get X component of pointing matrix
    values = getSumLine(sumin.readLine(), 4, "CX");
    for (int i = 0; i < 3; i++) m_pointingMatrix[0][i] = cvtDouble(values[i]);

    // Get Y component of pointing matrix
    values = getSumLine(sumin.readLine(), 4,"CY");
    for (int i = 0; i < 3; i++) m_pointingMatrix[1][i] = cvtDouble(values[i]);

    // Get Z component of pointing matrix
    values = getSumLine(sumin.readLine(), 4, "CZ");
    for (int i = 0; i < 3; i++) m_pointingMatrix[2][i] = cvtDouble(values[i]);

    // Get sun position
    values = getSumLine(sumin.readLine(), 4, "SZ");
    for (int i = 0; i < 3; i++) m_sunPosition[i] = cvtDouble(values[i]);

    // Get K-matrix
    values = getSumLine(sumin.readLine(), 7, "K-MATRIX");
    for (int i = 0; i < 6; i++) m_kmatrix[i] = cvtDouble(values[i]);

    // Get distortion parameters or SIGPXL values.
    values = getSumLine(sumin.readLine(), 0);

    if ( values.last() == "DISTORTION" ) {
      for (int i = 0; i < 4; i++) m_distortion[i] = cvtDouble(values[i]);
    }

    // Get formal spacecraft position uncertainty
    values = getSumLine(sumin.readLine(), 4, "SIGMA_VSO");
    for (int i = 0; i < 3; i++) m_sigmaSCPos[i] = cvtDouble(values[i]);

    // Get formal pointing uncertainty
    values = getSumLine(sumin.readLine(), 4, "SIGMA_PTG");
    for (int i = 0; i < 3; i++) m_sigmaSCPos[i] = cvtDouble(values[i]);

    // Should be at the landmarks line in the file now.  This test ensures it
    values = getSumLine(sumin.readLine(), 1, "LANDMARKS");

    // Read all landmarks (ControlMeasures) in this section until we encounter
    //limb fits...

    // Read to end of file
    while ( !sumin.atEnd() ) {
      values = getSumLine(sumin.readLine());
    }

   return;
  }

  /**
   * @brief Get referred input line from SUMFILE line
   *
   * This method will tokenize a simple string with optional error checking for
   * known line content. The string is expected to be line of a Gaskell SUMFILE.
   *
   * Each data string has control characters removed, commas converted to spaces
   * and multiple spaces compressed to a single space. The string is then
   * tokenized with each value separated by a space.
   *
   * Callers can use the nexpected parameter to ensure the expected fields are
   * read from the input stream. If 0, it will accept any number of tokens
   * (including 0). If the expected number of tokens are not met, an exception is
   * thrown to indicate unexpected input.
   *
   * The tag parameter is provided to check for a specific value in the last
   * column of a parsed SUMFILE line. The tag string is compared to the last
   * field string parsed from the processed input data string using a case
   * insensitve comparison. If the values are not equivalent, an exception is
   * thrown. If a value is not provided, this check is not applied.
   *
   * @param data       Line from a Gaskell SUMFILE
   * @param nexpected  Optional number of columns/fields after tokenization. a
   *                   value of 0 will not apply this check.
   * @param tag        If specified, is the expect string in last column of a
   *                   parsed SUMFILE line. String comparison is case
   *                   insensitive.
   *
   * @return QStringList List of tokenized string fields
   */
  QStringList SumFile::getSumLine(const QString &data, const int &nexpected,
                                  const QString &tag) const {

    // Copy string and parse results
    QString line = data.simplified();
    line.replace(",", " ");
    QStringList values = line.split(" ", QString::SkipEmptyParts);

    if ( nexpected > 0 ) {
      if ( values.size() != nexpected ) {
        QString message = "Expected " + toString(nexpected) + " but got "
                       + toString(values.size());
        throw IException(IException::User, message, _FILEINFO_);
      }
    }

    // Check for tags
    if ( !tag.isEmpty() ) {
      if ( values.size() > 0 ) {
        if ( values.last().toLower() != tag.toLower() ) {
          QString message = "Expected line tag given (" + tag +
                         ") does not match contents (" + values.last() + ")";
          throw IException(IException::User, message, _FILEINFO_);
        }
      }
      else {
        QString message = "Line tag given (" + tag + ") but line has not values";
        throw IException(IException::User, message, _FILEINFO_);
      }
    }

    return (values);
  }

  /**
   * @brief Convert text to double value
   *
   * This method will convert special forms of double text value formats to binary
   * form. It is particulary useful when handle text of the form 2.440D+3". This
   * format is not directly supported by conversion facilities.  The "D" is
   * converted to "e" in these cases.
   *
   * @param value Text value to convert
   *
   * @return double Converted text value
   */
  double SumFile::cvtDouble(const QString &value) const {
    QString dblstr(value);
    return (dblstr.replace("D","e").toDouble());
  }

  /**
   * @brief Load a list of SUMFILEs from a file
   *
   * This static method loads a list of SUMFILEs contained in a file. Each line in
   * the file is assumsed to be a SUMFILE and will be initialized using a new
   * SumFile object for each file contained in the list.
   *
   * @param sumFiles     Name of file containing SUMFILEs
   *
   * @return SumFileList List of SumFile objects created from SUMFILE list
   */
  SumFileList loadSumFiles(const FileList &sumFiles) {
    Progress progress;
    progress.SetText("Loading Sum File List");
    progress.SetMaximumSteps(sumFiles.size());
    SumFileList sum_list;

    for (int f = 0; f < sumFiles.size(); f++) {
      progress.CheckStatus();
      SharedSumFile sumF(new SumFile(sumFiles[f].original()));
      sum_list.append(sumF);
    }

    return (sum_list);
  }

}
// namespace Isis


text: #include <boost/geometry/core/access.hpp>

max_stars_repo_name: miathedev/BoostForArduino

max_issues_repo_name: miathedev/BoostForArduino

max_forks_repo_name: miathedev/BoostForArduino


text:
/******************************************************************************

  Implementation of mean.

  Copyright (c) 2013
  Alexander Rukletsov <rukletsov@gmail.com>
  All rights reserved.

  Redistribution and use in source and binary forms, with or without
  modification, are permitted provided that the following conditions
  are met:
  1.  Redistributions of source code must retain the above copyright
      notice, this list of conditions and the following disclaimer.
  2.  Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.

  THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS "AS IS" AND
  ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
  FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
  OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
  OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
  SUCH DAMAGE.

*******************************************************************************/

#ifndef MEAN_HPP_BB12C8E0_E48D_11E2_83DA_9FFD99890B3C
#define MEAN_HPP_BB12C8E0_E48D_11E2_83DA_9FFD99890B3C

#include <vector>
#include <algorithm>
#include <functional>
#include <stdexcept>

// Suppress annoying MSVC's C4244 conversion warnings popping out from boost::accumulators
// library (boost::numeric::functional namespace) mainly due to division operation
// in mean computation. Suppress MSVC's C4512 warning for boost auxiliary classes as well.
#ifdef _MSC_VER
#   pragma warning(push)
#   pragma warning(disable:4244)
#   pragma warning(disable:4512)
#endif // _MSC_VER

#include <boost/accumulators/accumulators.hpp>
#include <boost/accumulators/statistics/stats.hpp>
#include <boost/accumulators/statistics/mean.hpp>

#ifdef _MSC_VER
#   pragma warning(pop)
#endif // _MSC_VER

namespace bo {
namespace math {

// Returns the mean of the given samples.
template <typename SampleType>
SampleType mean(const std::vector<SampleType>& data)
{
    // The mean of an empty set is meaningless.
    if (data.size() == 0)
        throw std::logic_error("The mean of an empty set is meaningless.");

    // Initialize boost accumulator.
    namespace accs = boost::accumulators;
    typedef accs::accumulator_set<SampleType, accs::stats<accs::tag::mean> > Acc;
    Acc acc;

    // Fill accumulator with data.
    acc = std::for_each(data.begin(), data.end(), acc);

    // Request and return mean.
    return accs::mean(acc);
}

} // namespace math
} // namespace bo

#endif // MEAN_HPP_BB12C8E0_E48D_11E2_83DA_9FFD99890B3C

max_stars_repo_name: rukletsov/bo

max_issues_repo_name: rukletsov/bo

max_forks_repo_name: rukletsov/bo


text: #include <boost/preprocessor/repetition/enum_params_with_defaults.hpp>


text:
/*!@file
 * @copyright This code is licensed under the 3-clause BSD license.
 *   Copyright ETH Zurich, Laboratory of Physical Chemistry, Reiher Group.
 *   See LICENSE.txt for details.
 */

#include <boost/test/unit_test.hpp>

#include "Molassembler/Modeling/ShapeInference.h"
#include "Molassembler/Shapes/Data.h"

#include "boost/optional.hpp"

#include <iostream>

using namespace Scine;
using namespace Molassembler;
using namespace ShapeInference;

using TestCaseType = std::tuple<
  std::string, // a Name for the current compound
  Utils::ElementType, // The central atom type
  std::vector<BindingSite>, // a list of ligand types
  int // charge centered on the central atom
>;

// Runs a specific test case type
void testVsepr(
  const Shapes::Shape& expectedShapeName,
  const std::vector<TestCaseType>& testCases
) {
  std::string complexName;
  int charge;
  std::vector<BindingSite> ligands;
  Utils::ElementType centerAtomType;

  for(const auto& testCase : testCases) {

    std::tie(complexName, centerAtomType, ligands, charge) = testCase;

    boost::optional<Shapes::Shape> shape;

    try {
      shape = vsepr(
        centerAtomType,
        ligands,
        charge
      );
    } catch(std::exception& e) {
      std::cout << "Caught an exception when trying "
        << complexName << "." << std::endl;
      throw e;
    }

    BOOST_CHECK_MESSAGE(
      shape,
      "VSEPR could not find a shape for this case"
    );

    if(shape) {
      BOOST_CHECK_MESSAGE(
        shape.value() == expectedShapeName,
        complexName << " has been determined as "
          << Shapes::name(shape.value())
          << ", expected: "
          << Shapes::name(expectedShapeName)
      );
    }
  }
}

// Helper function to compose ligand situations
auto makeLigand(
  const unsigned L,
  const unsigned X,
  const Utils::ElementType& type,
  const BondType& bty
) {
  return BindingSite {
    L,
    X,
    {type},
    bty
  };
}

// Helper function to compose ligand situations
std::vector<BindingSite> repeat(
  const BindingSite& ligand,
  const unsigned N
) {
  return std::vector<BindingSite> (N, ligand);
}

// Helper function to compose ligand situations
std::vector<BindingSite> merge(
  const std::vector<BindingSite>& a,
  const std::vector<BindingSite>& b
) {
  std::vector<BindingSite> ret = a;
  std::copy(
    b.begin(),
    b.end(),
    std::back_inserter(ret)
  );

  return ret;
}

BOOST_AUTO_TEST_CASE(VseprTests, *boost::unit_test::label("Molassembler")) {
  using Element = Utils::ElementType;

  testVsepr( // AX2E0
    Shapes::Shape::Line,
    {
      TestCaseType {
        "CO2",
        Element::C,
        repeat(
          makeLigand(0, 2, Element::O, BondType::Double),
          2
        ),
        0
      },
      TestCaseType {
        "BeF2",
        Element::Be,
        repeat(
          makeLigand(0, 1, Element::F, BondType::Single),
          2
        ),
        0
      },
      TestCaseType {
        "C≡C-H",
        Element::C,
        {
          makeLigand(0, 3, Element::C, BondType::Triple),
          makeLigand(0, 1, Element::H, BondType::Single)
        },
        0
      }
    }
  );

  testVsepr( // AX3E0
    Shapes::Shape::EquilateralTriangle,
    {
      TestCaseType {
        "BCl3",
        Element::B,
        repeat(
          makeLigand(0, 1, Element::Cl, BondType::Single),
          3
        ),
        0
      },
      TestCaseType {
        "H2C=O",
        Element::C,
        {
          makeLigand(0, 1, Element::H, BondType::Single),
          makeLigand(0, 1, Element::H, BondType::Single),
          makeLigand(0, 2, Element::O, BondType::Double),
        },
        0
      },
      TestCaseType {
        "SO3",
        Element::S,
        repeat(
          makeLigand(0, 2, Element::O, BondType::Double),
          3
        ),
        0
      }
    }
  );

  testVsepr( // AX2E1 (V-shape / Bent)
    Shapes::Shape::Bent,
    {
      TestCaseType {
        "SO2",
        Element::S,
        repeat(
          makeLigand(0, 2, Element::O, BondType::Double),
          2
        ),
        0
      },
      TestCaseType {
        "CH2··",
        Element::C,
        repeat(
          makeLigand(0, 1, Element::H, BondType::Single),
          2
        ),
        0
      },
      TestCaseType {
        "NH2^+",
        Element::N,
        repeat(
          makeLigand(0, 1, Element::H, BondType::Single),
          2
        ),
        +1
      },
      TestCaseType {
        /* Lewis structure is
         *
         *  O = N⁺ – O⁻  < – >  ⁻O – N⁺ = O
         *
         */
        "NO2·",
        Element::N,
        {
          makeLigand(0, 1, Element::O, BondType::Single),
          makeLigand(0, 2, Element::O, BondType::Double),
        },
        +1
      }
    }
  );

  testVsepr( // AX4E0
    Shapes::Shape::Tetrahedron,
    {
      TestCaseType {
        "CH4",
        Element::C,
        repeat(
          makeLigand(0, 1, Element::H, BondType::Single),
          4
        ),
        0
      },
      TestCaseType {
        "XeO4",
        Element::Xe,
        repeat(
          makeLigand(0, 2, Element::O, BondType::Double),
          4
        ),
        0
      },
      TestCaseType {
        "ClO4-",
        Element::Cl,
        {
          makeLigand(0, 2, Element::O, BondType::Double),
          makeLigand(0, 2, Element::O, BondType::Double),
          makeLigand(0, 2, Element::O, BondType::Double),
          makeLigand(0, 1, Element::O, BondType::Single)
        },
        0 // charge is not on central chlorine atom!
      },
      TestCaseType {
        "N≡SF3",
        Element::S,
        {
          makeLigand(0, 1, Element::F, BondType::Single),
          makeLigand(0, 1, Element::F, BondType::Single),
          makeLigand(0, 1, Element::F, BondType::Single),
          makeLigand(0, 3, Element::N, BondType::Triple)
        },
        0
      }
    }
  );

  testVsepr( // AX3E1
    Shapes::Shape::VacantTetrahedron,
    {
      TestCaseType {
        "NH3",
        Element::N,
        repeat(
          makeLigand(0, 1, Element::H, BondType::Single),
          3
        ),
        0
      },
      TestCaseType {
        "XeO3",
        Element::Xe,
        repeat(
          makeLigand(0, 2, Element::O, BondType::Double),
          3
        ),
        0
      }
    }
  );

  testVsepr( // AX2E2
    Shapes::Shape::Bent,
    {
      TestCaseType {
        "OH2",
        Element::O,
        repeat(
          makeLigand(0, 1, Element::H, BondType::Single),
          2
        ),
        0
      }
    }
  );

  testVsepr(
    Shapes::Shape::TrigonalBipyramid,
    {
      TestCaseType {
        "PCl5",
        Element::P,
        repeat(
          makeLigand(0, 1, Element::Cl, BondType::Single),
          5
        ),
        0
      }
    }
  );

  testVsepr(
    Shapes::Shape::Seesaw,
    {
      TestCaseType {
        "SF4",
        Element::S,
        repeat(
          makeLigand(0, 1, Element::F, BondType::Single),
          4
        ),
        0
      }
    }
  );

  testVsepr(
    Shapes::Shape::T,
    {
      TestCaseType {
        "ClF3",
        Element::Cl,
        repeat(
          makeLigand(0, 1, Element::F, BondType::Single),
          3
        ),
        0
      }
    }
  );

  testVsepr( // AX2E3
    Shapes::Shape::Line,
    {
      TestCaseType {
        "I3-",
        Element::I,
        {
          makeLigand(0, 1, Element::I, BondType::Single),
          makeLigand(1, 0, Element::I, BondType::Single)
        },
        -1
      },
      TestCaseType {
        "XeF2",
        Element::Xe,
        repeat(
          makeLigand(0, 1, Element::F, BondType::Single),
          2
        ),
        0
      }
    }
  );

  testVsepr(
    Shapes::Shape::Octahedron,
    {
      TestCaseType {
        "SF6",
        Element::S,
        repeat(
          makeLigand(0, 1, Element::F, BondType::Single),
          6
        ),
        0
      }
    }
  );

  testVsepr(
    Shapes::Shape::SquarePyramid,
    {
      TestCaseType {
        "BrF5",
        Element::Br,
        repeat(
          makeLigand(0, 1, Element::F, BondType::Single),
          5
        ),
        0
      }
    }
  );

  testVsepr(
    Shapes::Shape::Square,
    {
      TestCaseType {
        "XeF4",
        Element::Xe,
        repeat(
          makeLigand(0, 1, Element::F, BondType::Single),
          4
        ),
        0
      }
    }
  );

  testVsepr(
    Shapes::Shape::PentagonalBipyramid,
    {
      TestCaseType {
        "IF7",
        Element::I,
        repeat(
          makeLigand(0, 1, Element::F, BondType::Single),
          7
        ),
        0
      }
    }
  );

  testVsepr(
    Shapes::Shape::PentagonalPyramid,
    {
      TestCaseType {
        "XeOF5^-",
        Element::Xe,
        merge(
          repeat(
            makeLigand(0, 1, Element::F, BondType::Single),
            5
          ),
          repeat(
            makeLigand(0, 2, Element::O, BondType::Double),
            1
          )
        ),
        -1
      }
    }
  );

  testVsepr(
    Shapes::Shape::Pentagon,
    {
      TestCaseType {
        "XeF5-",
        Element::Xe,
        repeat(
          makeLigand(0, 1, Element::F, BondType::Single),
          5
        ),
        -1
      }
    }
  );

  testVsepr(
    Shapes::Shape::SquareAntiprism,
    {
      TestCaseType {
        "XeF8^2-",
        Element::Xe,
        repeat(
          makeLigand(0, 1, Element::F, BondType::Single),
          8
        ),
        0 // charge is not located on Xe!
      }, TestCaseType {
        "IF8^-",
        Element::I,
        repeat(
          makeLigand(0, 1, Element::F, BondType::Single),
          8
        ),
        -1
      }
    }
  );

  // dative bonding
  testVsepr(
    Shapes::Shape::Tetrahedron,
    {
      TestCaseType {
        "H3N -> BF3, on N",
        Element::N,
        {
          makeLigand(0, 1, Element::H, BondType::Single),
          makeLigand(0, 1, Element::H, BondType::Single),
          makeLigand(0, 1, Element::H, BondType::Single),
          makeLigand(0, 0, Element::B, BondType::Single)
        },
        +1
      },
      TestCaseType {
        "H3N -> BF3, on B",
        Element::B,
        {
          makeLigand(0, 1, Element::H, BondType::Single),
          makeLigand(0, 1, Element::H, BondType::Single),
          makeLigand(0, 1, Element::H, BondType::Single),
          makeLigand(1, 0, Element::N, BondType::Single)
        },
        -1
      }
    }
  );
}


text:
#include "connection.hpp"
#include <serverpp/tcp_server.hpp>
#include <boost/range/adaptor/filtered.hpp>
#include <boost/range/adaptor/transformed.hpp>
#include <boost/range/algorithm/for_each.hpp>
#include <boost/range/algorithm_ext/erase.hpp>
#include <algorithm>
#include <functional>
#include <iostream>
#include <memory>
#include <utility>
#include <vector>

namespace {

serverpp::byte_storage rot13_encode(serverpp::bytes data)
{
    auto const &encode_byte = [](serverpp::byte datum)
    {
        if (datum >= 'A' && datum <= 'Z')
        {
            return serverpp::byte('A' + (((datum - 'A') + 13) % 26));
        }
        else if (datum >= 'a' && datum <= 'z')
        {
            return serverpp::byte('a' + (((datum - 'a') + 13) % 26));
        }
        else if (datum == '\0')
        {
            return serverpp::byte('\n');
        }
        else
        {
            return datum;
        }

    };

    auto const encoded_data = data | boost::adaptors::transformed(encode_byte);
    return serverpp::byte_storage{encoded_data.begin(), encoded_data.end()};
}

}

namespace rot13 {

class server
{
public:
    server()
      : work_(boost::asio::make_work_guard(io_context_)),
        tcp_server_(
            io_context_,
            0,
            [this](auto &&new_socket)
            {
                this->new_connection(std::forward<decltype(new_socket)>(new_socket));
            })
    {
        std::cout << "TCP Server started up on port " << tcp_server_.port() << "\n";
    }

    void run()
    {
        io_context_.run();
    }

private:
    void new_connection(serverpp::tcp_socket &&new_socket)
    {
        std::cout << "Accepted new socket\n";

        connections_.emplace_back(new connection(std::move(new_socket)));
        auto &connection = connections_.back();

        connection->async_get_terminal_type(
            [](std::string const &ttype)
            {
                std::cout << "Terminal type = " << ttype << "\n";
            });

        connection->on_window_size_changed(
            [](std::uint16_t width, std::uint16_t height)
            {
                std::cout << "Window size is now " << width << "x" << height << "\n";
            });

        schedule_read(*connection);
    }

    void schedule_read(connection &cnx)
    {
        cnx.async_read(
            [this, &cnx](serverpp::bytes data)
            {
                read_handler(cnx, data);
            },
            [this, &cnx]()
            {
                if (cnx.is_alive())
                {
                    schedule_read(cnx);
                }
                else
                {
                    connection_death_handler(cnx);
                }
            });
    }

    void read_handler(connection &cnx, serverpp::bytes data)
    {
        cnx.write(rot13_encode(data));
    }

    void connection_death_handler(connection &dead_connection)
    {
        std::cout << "Connection died\n";

        const auto is_dead_connection = [&dead_connection](auto const &connection)
        {
            return connection.get() == &dead_connection;
        };

        boost::for_each(
            connections_ | boost::adaptors::filtered(is_dead_connection),
            [](auto &connection)
            {
                connection.reset();
            });

        boost::remove_erase_if(
            connections_,
            [](auto const &connection)
            {
                return !connection;
            });
    }

    boost::asio::io_context io_context_;
    boost::asio::executor_work_guard<
        boost::asio::io_context::executor_type> work_;

    serverpp::tcp_server tcp_server_;

    std::vector<std::unique_ptr<connection>> connections_;
};

}

int main()
{
    rot13::server server;
    server.run();
}

max_stars_repo_name: KazDragon/telnetpp

max_issues_repo_name: KazDragon/telnetpp

max_forks_repo_name: KazDragon/telnetpp


text:
//==============================================================================
//         Copyright 2003 - 2012   LASMEA UMR 6602 CNRS/Univ. Clermont II
//         Copyright 2009 - 2012   LRI    UMR 8623 CNRS/Univ Paris Sud XI
//
//          Distributed under the Boost Software License, Version 1.0.
//                 See accompanying file LICENSE.txt or copy at
//                     http://www.boost.org/LICENSE_1_0.txt
//==============================================================================
#ifndef BOOST_SIMD_UNIT_ROLL_TEST_HPP_INCLUDED
#define BOOST_SIMD_UNIT_ROLL_TEST_HPP_INCLUDED

#include <boost/simd/swar/include/functions/shuffle.hpp>
#include <nt2/sdk/unit/tests.hpp>

//==============================================================================
template<class T, int IA, int IB>
struct roll_unary_test_2
{
  static void call( T a )
  {
    using boost::simd::shuffle;
    T ref, res;
    ref[0] = (IA == -1) ? 0 : a[IA];
    ref[1] = (IB == -1) ? 0 : a[IB];
    res = shuffle<IA,IB>(a);
    NT2_TEST_EQUAL(res,ref);
    roll_unary_test_2<T,IA+1,IB>::call(a);
  }
};

template<class T, int IB>
struct roll_unary_test_2<T,2,IB>
{
  static void call( T a ) { roll_unary_test_2<T,-1,IB+1>::call(a); }
};

template<class T>
struct roll_unary_test_2<T,-1,2> { static void call( T a ){} };

//==============================================================================
template<class T, int I0, int I1, int I2, int I3>
struct roll_unary_test_4
{
  static void call( T a )
  {
    using boost::simd::shuffle;
    T ref, res;
    ref[0] = (I0 == -1) ? 0 : a[I0];
    ref[1] = (I1 == -1) ? 0 : a[I1];
    ref[2] = (I2 == -1) ? 0 : a[I2];
    ref[3] = (I3 == -1) ? 0 : a[I3];
    res = shuffle<I0,I1,I2,I3>(a);
    NT2_TEST_EQUAL(res,ref);
    roll_unary_test_4<T,I0+1,I1,I2,I3>::call(a);
  }
};

template<class T, int I1, int I2, int I3>
struct roll_unary_test_4<T,4,I1,I2,I3>
{
  static void call( T a ) { roll_unary_test_4<T,-1,I1+1,I2,I3>::call(a); }
};

template<class T, int I2, int I3>
struct roll_unary_test_4<T,-1,4,I2,I3>
{
  static void call( T a ) { roll_unary_test_4<T,-1,-1,I2+1,I3>::call(a); }
};

template<class T, int I3>
struct roll_unary_test_4<T,-1,-1,4,I3>
{
  static void call( T a ) { roll_unary_test_4<T,-1,-1,-1,I3+1>::call(a); }
};

template<class T>
struct roll_unary_test_4<T,-1,-1,-1,4> { static void call( T a ){} };

//==============================================================================
template<class T, int IA, int IB>
struct roll_binary_test_2
{
  static void call( T a, T b)
  {
    using boost::simd::shuffle;
    T ref, res;
    ref[0] = (IA == -1) ? 0 : (IA<2 ? a[IA] : b[IA-2]);
    ref[1] = (IB == -1) ? 0 : (IB<2 ? a[IB] : b[IB-2]);
    res = shuffle<IA,IB>(a,b);
    NT2_TEST_EQUAL(res,ref);
    roll_binary_test_2<T,IA+1,IB>::call(a,b);
  }
};

template<class T, int IB>
struct roll_binary_test_2<T,4,IB>
{
  static void call( T a, T b)
  {
     roll_binary_test_2<T,-1,IB+1>::call(a,b);
  }
};

template<class T>
struct roll_binary_test_2<T,-1,4>
{
  static void call( T a, T b) {}
};

//==============================================================================
template<class T, int I0, int I1, int I2, int I3>
struct roll_binary_test_4
{
  static void call( T a, T b )
  {
    using boost::simd::shuffle;
    T ref, res;
    ref[0] = (I0 == -1) ? 0 : ((I0<4)?a[I0]:b[I0-4]);
    ref[1] = (I1 == -1) ? 0 : ((I1<4)?a[I1]:b[I1-4]);
    ref[2] = (I2 == -1) ? 0 : ((I2<4)?a[I2]:b[I2-4]);
    ref[3] = (I3 == -1) ? 0 : ((I3<4)?a[I3]:b[I3-4]);
    res = shuffle<I0,I1,I2,I3>(a,b);
    NT2_TEST_EQUAL(res,ref);
    roll_binary_test_4<T,I0+1,I1,I2,I3>::call(a,b);
  }
};

template<class T, int I1, int I2, int I3>
struct roll_binary_test_4<T,8,I1,I2,I3>
{
  static void call( T a, T b ) { roll_binary_test_4<T,-1,I1+1,I2,I3>::call(a,b); }
};

template<class T, int I2, int I3>
struct roll_binary_test_4<T,-1,8,I2,I3>
{
  static void call( T a, T b ) { roll_binary_test_4<T,-1,-1,I2+1,I3>::call(a,b); }
};

template<class T, int I3>
struct roll_binary_test_4<T,-1,-1,8,I3>
{
  static void call( T a, T b ) { roll_binary_test_4<T,-1,-1,-1,I3+1>::call(a,b); }
};

template<class T>
struct roll_binary_test_4<T,-1,-1,-1,8> { static void call( T a, T b ){} };

#endif


text:
// system includes -----------------------------------------------
#include <Eigen/Geometry>
#include <boost/math/constants/constants.hpp>
#include <boost/program_options.hpp>
#include <iostream>
// own includes --------------------------------------------------
#include <base/eigen2hdf.hpp>
#include <base/init.hpp>
#include <base/timer.hpp>
#include <fft/fft2.hpp>
#include <fft/fft2_r2c.hpp>
#include <ridgelet/ridgelet_cell_array.hpp>
#include <ridgelet/ridgelet_frame.hpp>
#include <ridgelet/rt.hpp>

#include <boundary_conditions/inflow_bc.hpp>
#include <operators/operators.hpp>
#include <solver/cg.hpp>
#include <solver/ridgelet_solver.hpp>

using namespace std;

typedef RT<> RT_t;
typedef RT_t::array_t array_t;
typedef RT_t::complex_array_t complex_array_t;
typedef RT_t::rt_coeff_t rt_coeff_t;
typedef RidgeletCellArray<rt_coeff_t> rca_t;

namespace po = boost::program_options;

const double Lx = 1.2;
const double Ly = 1.2;

class RidgeletInflowTest
{
 public:
  RidgeletInflowTest(po::variables_map vm);

  void init_sigma(array_t& Sigma);
  /// initialize fft(f(t=0, x))
  void init_fh0(complex_array_t& Fh);
  /// prescribed f(x) in absorption layer
  void init_fbc(array_t& fbc, double vx, double vy, double ql, double qt);
  /// apply boundary conditions in absorption layer to fft(f(x))
  void apply_bc(complex_array_t& Fh, const array_t& fbc);
  /**
   *
   * @param vx
   * @param vy
   * @param ql    inflow left
   * @param qt    inflow top
   * @param rho   density
   */
  void run(double vx, double vy, double ql, double qt, double rho, int ntsteps);

  RidgeletSolver<rt_coeff_t>& solver() { return rt_solver; }

  void log_cg(bool flag) { log_hist = flag; }

 private:
  RidgeletFrame rf;
  RT_t rt;
  FFTr2c<PlannerR2C> fft;
  RidgeletSolver<rt_coeff_t> rt_solver;
  // parameters
  po::variables_map vm;
  bool log_hist = false;
};

RidgeletInflowTest::RidgeletInflowTest(po::variables_map vm_)
    : vm(vm_)
{
  unsigned int Jx = vm["Jx"].as<unsigned int>();
  unsigned int Jy = vm["Jy"].as<unsigned int>();
  unsigned int rho_x = vm["rx"].as<unsigned int>();
  unsigned int rho_y = vm["ry"].as<unsigned int>();

  rf = RidgeletFrame(Jx, Jy, rho_x, rho_y);
  rt = RT_t(rf);

  const unsigned int Nx = rf.Nx();  // #cols
  const unsigned int Ny = rf.Ny();  // #rows
  init_fftw(fft, FFTW_MEASURE, rf);
  fft.get_plan().create_and_get_plan(Ny / 2, Nx / 2, PlannerR2C::INV);
  fft.get_plan().create_and_get_plan(Ny / 2, Nx / 2, PlannerR2C::FWD);
}

void
RidgeletInflowTest::init_sigma(array_t& Sigma)
{
  double sigma = vm["sigma"].as<double>();
  unsigned int Nx = rf.Nx();
  unsigned int Ny = rf.Ny();

  Eigen::ArrayXd xi = Eigen::ArrayXd::LinSpaced(Nx + 1, 0, Lx).segment(0, Nx);
  Eigen::ArrayXd yi = Eigen::ArrayXd::LinSpaced(Ny + 1, 0, Ly).segment(0, Ny);

  Sigma.setZero();

  auto sigmaf_const = [sigma](double x, double y) {
    double d = std::max(x, y);
    if (d > 1) {
      return sigma;
    } else {
      return 0.0;
    }
  };
  auto sigmaf_lin = [sigma](double x, double y) {
    double d = std::max(x, y);
    if (d > 1 && d <= 1.1)
      return (d - 1) / 0.1 * sigma;
    else
      return sigma;
  };
  auto sigmaf_smooth = [sigma](double x, double y) {
    double d = std::max(x, y);
    double z = d - 1;
    if (d > 1 && d <= 1.2)
      return 400 * std::pow(z, 2) - 4000 * std::pow(z, 3) + 10000 * std::pow(z, 4);
    else
      return 0.;
  };
  std::function<double(double, double)> sigmaf;
  if (vm["sigmaf"].as<std::string>() == "const")
    sigmaf = sigmaf_const;
  else if (vm["sigmaf"].as<std::string>() == "lin")
    sigmaf = sigmaf_lin;
  else if (vm["sigmaf"].as<std::string>() == "smooth")
    sigmaf = sigmaf_smooth;
  else {
    std::cerr << "Option --sigmaf=" << vm["sigmaf"].as<std::string>() << " not found." << std::endl;
  }

  // initialize sigma
  for (unsigned int i = 0; i < Ny; ++i) {
    for (unsigned int j = 0; j < Nx; ++j) {
      if (yi[j] < 1 and yi[i] < 1) continue;
      Sigma(i, j) = sigmaf(xi[j], yi[i]);
    }
  }
}

void
RidgeletInflowTest::init_fh0(complex_array_t& Fh)
{
  unsigned int Nx = rf.Nx();
  unsigned int Ny = rf.Ny();

  Eigen::ArrayXd xi = Eigen::ArrayXd::LinSpaced(Nx / 2 + 1, 0, Lx).segment(0, Nx / 2);
  Eigen::ArrayXd yi = Eigen::ArrayXd::LinSpaced(Ny / 2 + 1, 0, Ly).segment(0, Ny / 2);

  array_t X = xi.transpose().replicate(Ny / 2, 1);
  array_t Y = yi.replicate(1, Nx / 2);
  array_t F = array_t::Ones(Ny / 2, Nx / 2);
  F = (X < 1 && Y < 1).select(F, array_t::Zero(Ny / 2, Nx / 2));
  complex_array_t fh(Ny / 2, Nx / 2);
  fft.ft(fh, F, false);
  Fh.setZero();
  ftcut(Fh, Ny / 2, Nx / 2) = 4 * fh;
}

void
RidgeletInflowTest::init_fbc(array_t& fbc, double vx, double vy, double ql, double qt)
{
  unsigned int Nx = rf.Nx();
  unsigned int Ny = rf.Ny();

  fbc.resize(Ny / 2, Nx / 2);
  make_inflow_source(fbc, vx, vy, Lx, Ly, ql, qt);

  Eigen::ArrayXd xi = Eigen::ArrayXd::LinSpaced(Nx / 2 + 1, 0, Lx).segment(0, Nx / 2);
  Eigen::ArrayXd yi = Eigen::ArrayXd::LinSpaced(Ny / 2 + 1, 0, Ly).segment(0, Ny / 2);

  array_t X = xi.transpose().replicate(Ny / 2, 1);
  array_t Y = yi.replicate(1, Nx / 2);
  fbc = (X < 1 && Y < 1).select(array_t::Zero(Ny / 2, Nx / 2), fbc);
}

void
RidgeletInflowTest::apply_bc(complex_array_t& Fh, const array_t& fbc)
{
  unsigned int Nx = rf.Nx();
  unsigned int Ny = rf.Ny();

  complex_array_t fh = ftcut(Fh, Ny / 2, Nx / 2);
  array_t F(Ny / 2, Nx / 2);
  fft.ift(F, fh);
  F /= 4;
  F = (fbc > 0).select(fbc, F);
  fft.ft(fh, F, false);
  Fh.setZero();
  ftcut(Fh, Ny / 2, Nx / 2) = 4 * fh;
}

void
RidgeletInflowTest::run(double vx, double vy, double ql, double qt, double rho, int ntsteps)
{
  rca_t x(rf);
  double dt = vm["dt"].as<double>();
  int cg_maxit = vm["cg_maxit"].as<int>();
  double cg_tol = vm["cg_tol"].as<double>();

  unsigned int Nx = rf.Nx();
  unsigned int Ny = rf.Ny();

  array_t Sigma(Ny, Nx);
  this->init_sigma(Sigma);

  array_t fbc(Ny / 2, Nx / 2);
  this->init_fbc(fbc, vx, vy, ql, qt);

  hid_t file, group_id;
  if (vm.count("save"))
    file = H5Fcreate("advection_inflow.h5", H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);

  array_t q(Ny, Nx);
  double r = sqrt(vx * vx + vy * vy);
  q.setZero();
  make_inflow_source(q, vx, vy, Lx, Ly, ql, qt);
  q *= (Sigma * r);

  if (vm.count("save")) {
    eigen2hdf::save(file, "q", q);
    eigen2hdf::save(file, "Sigma", Sigma);
    eigen2hdf::save(file, "fbc", fbc);
  }
  complex_array_t qh(Ny, Nx);
  fft.ft(qh, q, false);

  // make F decay smoothly in the absorption layer
  // F *= (array_t::Ones(Ny, Nx) -Sigma/sigma);
  complex_array_t Fh(Ny, Nx);
  this->init_fh0(Fh);

  auto& x_coeffs = x.coeffs();
  rt.rt(x_coeffs, Fh);

  // initialize operators
  // AhA = T'T
  AhAOpsigma AhA(Sigma, vx, vy, Lx, Ly, Nx, Ny, dt);
  // preconditioned operator
  RDTSCTimer timer;
  complex_array_t Bh(Ny, Nx);
  complex_array_t Rhsh(Ny, Nx);  // right hand side
  rca_t b(rf);
  rt_solver = RidgeletSolver<rt_coeff_t>(rf, vx, vy);
  PTransportOpBC<RT_t> A(rt, AhA, vx, vy);
  TransportOperatorBC T(Sigma, vx, vy, Lx, Ly, Nx, Ny, dt);  // required for rhs

  if (vm.count("save")) {
    complex_array_t solh(Ny / 2, Nx / 2);
    solh = ftcut(Fh, Ny / 2, Nx / 2);
    array_t sol(Ny / 2, Nx / 2);
    // scale fourier coefficients by grid_out factor
    // ftcut(solh, Ny, Nx) = Fh;
    fft.ift(sol, solh);
    sol /= 4;
    group_id = H5Gcreate(
        file, boost::lexical_cast<string>(0).c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
    eigen2hdf::save(group_id, "sol", sol);
    H5Gclose(group_id);
  }
  // do not print log to stdout
  Logger::GetInstance().detach_stdout();
  for (int i = 1; i <= ntsteps; ++i) {
    this->apply_bc(Fh, fbc);
    // T'*Fh
    Rhsh = Fh + dt * qh;  // rhs in fourier space
    T.apply(Bh, Rhsh, true /* hermitian transpose */);
    rt.rt(b.coeffs(), Bh);

    timer.start();
    Logger::GetInstance().push_prefix(boost::lexical_cast<std::string>(i) + " CGSOLVER");
    rt_solver.solve(x, A, b, cg_tol, cg_maxit, log_hist);
    Logger::GetInstance().pop_prefix();
    auto nc_solve = timer.stop();

    cout << "RidgeletSolver took: " << nc_solve / 1e9 << " Gcycles\n";
    // write arrays to hdf5
    // complex_array_t tmp(Ny, Nx);
    rt.irt(Fh, x.coeffs());

    // eigen output
    if (vm.count("save")) {
      complex_array_t solh(Ny / 2, Nx / 2);
      solh = ftcut(Fh, Ny / 2, Nx / 2);
      array_t sol(Ny / 2, Nx / 2);
      // scale fourier coefficients by grid_out factor
      // ftcut(solh, Ny, Nx) = Fh;
      fft.ift(sol, solh);
      sol /= 4;
      group_id = H5Gcreate(
          file, boost::lexical_cast<string>(i).c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
      eigen2hdf::save(group_id, "sol", sol);
      H5Gclose(group_id);
    }
  }
  if (vm.count("save")) H5Fclose(file);
}


text:
//=======================================================================
// Copyright 2001 Jeremy G. Siek, Andrew Lumsdaine, Lie-Quan Lee,
//
// This file is part of the Boost Graph Library
//
// You should have received a copy of the License Agreement for the
// Boost Graph Library along with the software; see the file LICENSE.
// If not, contact Office of Research, Indiana University,
// Bloomington, IN 47405.
//
// Permission to modify the code and to distribute the code is
// granted, provided the text of this NOTICE is retained, a notice if
// the code was modified is included with the above COPYRIGHT NOTICE
// and with the COPYRIGHT NOTICE in the LICENSE file, and that the
// LICENSE file is distributed with the modified code.
//
// LICENSOR MAKES NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED.
// By way of example, but not limitation, Licensor MAKES NO
// REPRESENTATIONS OR WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY
// PARTICULAR PURPOSE OR THAT THE USE OF THE LICENSED SOFTWARE COMPONENTS
// OR DOCUMENTATION WILL NOT INFRINGE ANY PATENTS, COPYRIGHTS, TRADEMARKS
// OR OTHER RIGHTS.
//=======================================================================
#include <boost/config.hpp>
#include <iostream>
#include <fstream>
#include <boost/graph/adjacency_list.hpp>
using namespace boost;

template < typename VertexDescriptor, typename VertexNameMap > void
print_vertex_name(VertexDescriptor v, VertexNameMap name_map)
{
  std::cout << get(name_map, v);
}

template < typename Graph, typename TransDelayMap, typename VertexNameMap > void
print_trans_delay(typename graph_traits < Graph >::edge_descriptor e,
                  const Graph & g, TransDelayMap delay_map,
                  VertexNameMap name_map)
{
  std::cout << "trans-delay(" << get(name_map, source(e, g)) << ","
    << get(name_map, target(e, g)) << ") = " << get(delay_map, e);
}

template < typename Graph, typename VertexNameMap > void
print_vertex_names(const Graph & g, VertexNameMap name_map)
{
  std::cout << "vertices(g) = { ";
  typedef typename graph_traits < Graph >::vertex_iterator iter_t;
  for (std::pair < iter_t, iter_t > p = vertices(g); p.first != p.second;
       ++p.first) {
    print_vertex_name(*p.first, name_map);
    std::cout << ' ';
  }
  std::cout << "}" << std::endl;
}

template < typename Graph, typename TransDelayMap, typename VertexNameMap > void
print_trans_delays(const Graph & g, TransDelayMap trans_delay_map,
                   VertexNameMap name_map)
{
  typename graph_traits < Graph >::edge_iterator first, last;
  for (tie(first, last) = edges(g); first != last; ++first) {
    print_trans_delay(*first, g, trans_delay_map, name_map);
    std::cout << std::endl;
  }
}

template < typename Graph, typename VertexNameMap, typename TransDelayMap > void
build_router_network(Graph & g, VertexNameMap name_map,
                     TransDelayMap delay_map)
{
  typename graph_traits < Graph >::vertex_descriptor a, b, c, d, e;
  a = add_vertex(g);
  name_map[a] = 'a';
  b = add_vertex(g);
  name_map[b] = 'b';
  c = add_vertex(g);
  name_map[c] = 'c';
  d = add_vertex(g);
  name_map[d] = 'd';
  e = add_vertex(g);
  name_map[e] = 'e';

  typename graph_traits < Graph >::edge_descriptor ed;
  bool inserted;

  tie(ed, inserted) = add_edge(a, b, g);
  delay_map[ed] = 1.2;
  tie(ed, inserted) = add_edge(a, d, g);
  delay_map[ed] = 4.5;
  tie(ed, inserted) = add_edge(b, d, g);
  delay_map[ed] = 1.8;
  tie(ed, inserted) = add_edge(c, a, g);
  delay_map[ed] = 2.6;
  tie(ed, inserted) = add_edge(c, e, g);
  delay_map[ed] = 5.2;
  tie(ed, inserted) = add_edge(d, c, g);
  delay_map[ed] = 0.4;
  tie(ed, inserted) = add_edge(d, e, g);
  delay_map[ed] = 3.3;

}

int
main()
{
  typedef adjacency_list < listS, listS, directedS,
    property < vertex_name_t, char >,
    property < edge_weight_t, double > > graph_t;
  graph_t g;

  property_map < graph_t, vertex_name_t >::type name_map =
    get(vertex_name, g);
  property_map < graph_t, edge_weight_t >::type delay_map =
    get(edge_weight, g);

  build_router_network(g, name_map, delay_map);
  print_vertex_names(g, name_map);
  print_trans_delays(g, delay_map, name_map);
}


text:
/*
 * Copyright 2013 Matthew Harvey
 *
 * 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.
 */

#ifndef GUARD_interval_type_hpp_5555643181109635
#define GUARD_interval_type_hpp_5555643181109635

#include <boost/date_time/gregorian/gregorian.hpp>
#include <wx/string.h>
#include <ostream>

namespace dcm
{

/**
 * Enumerated type representing different date intervals.
 *
 * Should be self-explanatory; except note that
 * \e month_ends represents an interval between the
 * last day of one month and the last day of the next,
 * whereas \e months simply represents an interval of
 * one month with no particular reference to the end
 * of the month.
 *
 * Note the numbering of the enumeration is significant.
 * The class DcmDatabaseConnection relies on the
 * numbering shown here.
 */
enum class IntervalType: unsigned char
{
    days = 1,
    weeks,
    months,
    month_ends
};

/**
 * @returns a short phrase describing the interval.
 * E.g. "days", "month", etc.. The phrase is plural if and
 * only if \e is_plural is true.
 */
wxString phrase(IntervalType x, bool is_plural);

/**
 * @returns \e true if and only if \e p_date can occur within
 * an indefinitely long sequence of dates separated by intervals of
 * \e p_interval_type. This is particularly relevant for \e
 * IntervalType::months and \e IntervalType::month_ends. E.g. a date
 * of 29 Feb. 2012 cannot be accommodated by IntervalType::months because
 * not all months have a 29th; and a date of 30 July 2014 cannot be
 * accommodated by IntervalType::month_ends because it is not the last
 * day of a month.
 */
bool is_valid_date_for_interval_type
(   boost::gregorian::date const& p_date,
    IntervalType p_interval_type
);

/**
 * Output IntervalType as "user unfriendly string"
 * suitable for developer-facing diagnostic messages.
 */
std::ostream&
operator<<(std::ostream& os, IntervalType);

}  // namespace dcm

#endif  // GUARD_interval_type_hpp_5555643181109635

max_stars_repo_name: skybaboon/dailycashmanager

max_issues_repo_name: skybaboon/dailycashmanager

max_forks_repo_name: skybaboon/dailycashmanager


text:
/*
 * Copyright (c) 2016 Spotify AB
 *
 * 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 <memory>
#include <string>

#include <boost/test/unit_test.hpp>

#include <spotify/json/codec/empty_as.hpp>
#include <spotify/json/codec/eq.hpp>
#include <spotify/json/codec/null.hpp>
#include <spotify/json/codec/number.hpp>
#include <spotify/json/codec/object.hpp>
#include <spotify/json/codec/omit.hpp>
#include <spotify/json/codec/smart_ptr.hpp>
#include <spotify/json/codec/string.hpp>
#include <spotify/json/decode.hpp>
#include <spotify/json/encode.hpp>

BOOST_AUTO_TEST_SUITE(spotify)
BOOST_AUTO_TEST_SUITE(json)
BOOST_AUTO_TEST_SUITE(codec)

namespace {

struct Val {
  std::string a;
  std::string b;
};

template <typename Codec>
typename Codec::object_type test_decode(const Codec &codec, const std::string &json) {
  decode_context c(json.c_str(), json.c_str() + json.size());
  auto obj = codec.decode(c);
  BOOST_CHECK_EQUAL(c.position, c.end);
  return obj;
}

template <typename Codec>
void test_decode_fail(const Codec &codec, const std::string &json) {
  decode_context c(json.c_str(), json.c_str() + json.size());
  BOOST_CHECK_THROW(codec.decode(c), decode_exception);
}

}  // namespace

/*
 * Constructing
 */

BOOST_AUTO_TEST_CASE(json_codec_empty_as_should_default_construct) {
  empty_as_t<null_t<std::string>, string_t> codec;
}

BOOST_AUTO_TEST_CASE(json_codec_empty_as_should_construct) {
  empty_as_t<null_t<std::string>, string_t> codec((null_t<std::string>()), string_t());
}

BOOST_AUTO_TEST_CASE(json_codec_empty_as_should_construct_with_helper) {
  empty_as(null<std::string>(), string());
}

BOOST_AUTO_TEST_CASE(json_codec_empty_as_null_should_construct) {
  empty_as_null(string());
}

BOOST_AUTO_TEST_CASE(json_codec_empty_as_omit_should_construct) {
  empty_as_omit(string());
}

/*
 * Decoding
 */

BOOST_AUTO_TEST_CASE(json_codec_empty_as_should_decode) {
  const auto codec = empty_as_omit(string());
  BOOST_CHECK_EQUAL(test_decode(codec, "\"\""), "");
  BOOST_CHECK_EQUAL(test_decode(codec, "\"hello\""), "hello");
}

BOOST_AUTO_TEST_CASE(json_codec_empty_as_should_decode_default) {
  const auto codec = empty_as_null(string());
  BOOST_CHECK_EQUAL(test_decode(codec, "null"), "");
}

BOOST_AUTO_TEST_CASE(json_codec_empty_as_should_fail_on_invalid_input) {
  const auto codec = empty_as_omit(string());
  test_decode_fail(codec, "null");
  test_decode_fail(codec, "e");
  test_decode_fail(codec, "[{},true]");
}

BOOST_AUTO_TEST_CASE(json_codec_empty_as_with_eq) {
  const auto codec = empty_as(eq(123), number<int>());
  BOOST_CHECK_EQUAL(encode(codec, 0), "123");
  BOOST_CHECK(test_decode(codec, "123") == 123);
}

/*
 * Encoding
 */

BOOST_AUTO_TEST_CASE(json_codec_empty_as_should_encode_default) {
  const auto codec = empty_as_null(string());
  BOOST_CHECK(codec.should_encode(""));
  BOOST_CHECK_EQUAL(encode(codec, ""), "null");
}

BOOST_AUTO_TEST_CASE(json_codec_empty_as_should_not_encode_default) {
  const auto codec = empty_as_omit(string());
  BOOST_CHECK(!codec.should_encode(""));
}

BOOST_AUTO_TEST_CASE(json_codec_empty_as_should_encode_non_default) {
  const auto codec = empty_as_omit(string());
  BOOST_CHECK(codec.should_encode("abc"));
  BOOST_CHECK_EQUAL(encode(codec, "abc"), "\"abc\"");
}

BOOST_AUTO_TEST_CASE(json_codec_empty_as_with_shared_ptr) {
  const auto codec = empty_as_null(default_codec<std::shared_ptr<std::string>>());
  BOOST_CHECK_EQUAL(encode(codec, std::shared_ptr<std::string>()), "null");
  BOOST_CHECK_EQUAL(encode(codec, std::make_shared<std::string>("abc")), "\"abc\"");
}

BOOST_AUTO_TEST_CASE(json_codec_empty_as_with_object_and_omit) {
  auto codec = object<Val>();
  codec.optional("a", &Val::a);
  codec.optional("b", &Val::b, empty_as_omit(string()));
  BOOST_CHECK_EQUAL(encode(codec, Val()), R"({"a":""})");  // no "b"
}

BOOST_AUTO_TEST_SUITE_END()  // codec
BOOST_AUTO_TEST_SUITE_END()  // json
BOOST_AUTO_TEST_SUITE_END()  // spotify

max_stars_repo_name: ivangalkin/spotify-json

max_issues_repo_name: ivangalkin/spotify-json

max_forks_repo_name: ivangalkin/spotify-json


text:
#define BOOST_TEST_DYN_LINK
#define BOOST_TEST_MODULE parse_literal
#include <boost/test/unit_test.hpp>

#include "../src/parse_literal.hpp"
#include "../src/parse_literal.hpp"

#include "../src/error/compile_exception.hpp"

#include "state_utils.hpp"

using std::string;

BOOST_AUTO_TEST_CASE(standard)
{
    state s = make_state("\"abc\"");

    lit_node* got = parse_literal(s);
    BOOST_CHECK(got);
    string expected = "abc";

    BOOST_CHECK(rangeify(expected) == rangeify(*got));
    BOOST_CHECK(remaining(s) == "");
}

BOOST_AUTO_TEST_CASE(digits)
{
    state s = make_state("912fas31");

    lit_node* got = parse_literal(s);
    BOOST_CHECK(got);
    string expected = "912";

    BOOST_CHECK(rangeify(*got) == rangeify(expected));
    BOOST_CHECK(remaining(s) == "fas31");
}

BOOST_AUTO_TEST_CASE(unterminated_quote)
{
    state s = make_state("\"abc");

    BOOST_CHECK_THROW(parse_literal(s), compile_exception);
}

BOOST_AUTO_TEST_CASE(quote_terminated_next_line)
{
    state s = make_state("\"abc\nasdf\"");

    BOOST_CHECK_THROW(parse_literal(s), compile_exception);
}

BOOST_AUTO_TEST_CASE(no_literal)
{
    state s = make_state("abc");

    lit_node* got = parse_literal(s);
    BOOST_CHECK(!got);
    BOOST_CHECK(remaining(s) == "abc");
}

max_stars_repo_name: mbmaier/asm-lisp

max_issues_repo_name: mbmaier/asm-lisp

max_forks_repo_name: mbmaier/asm-lisp


text:
#include "stats/rmsd.h"
#include <xmd/utils/convert.h>
#include <Eigen/Core>
#include <Eigen/SVD>

namespace xmd {
    void compute_rmsd::operator()() const {
        using matrix_t = Eigen::Matrix<real, Eigen::Dynamic, 3>;
        matrix_t P = matrix_t::Zero(num_particles, 3);
        matrix_t Q = matrix_t::Zero(num_particles, 3);

        for (int idx = 0; idx < num_particles; ++idx) {
            P.row(idx) = convert<real>(r[idx]);
            Q.row(idx) = convert<real>(ref_r[idx]);
        }

        Eigen::Matrix<real, 3, 3> H = P.transpose() * Q;
        auto svd = Eigen::JacobiSVD<decltype(H)>(H,
            Eigen::ComputeFullU & Eigen::ComputeFullV);

        auto U = svd.matrixU();
        auto V = svd.matrixV();
        auto d = (V * U.transpose()).determinant() > 0 ? 1 : -1;
        auto R = V * Eigen::DiagonalMatrix<real, 3>(1, 1, d) * U.transpose();

        auto D = (P - Q) * R;
        real rmsd_ = 0.0;
        for (int idx = 0; idx < num_particles; ++idx) {
            rmsd_ += D.row(idx).squaredNorm();
        }
        *rmsd = sqrt(rmsd_ / num_particles);
    }
}

max_stars_repo_name: vitreusx/xmd

max_issues_repo_name: vitreusx/xmd

max_forks_repo_name: vitreusx/xmd


text: #include <boost/mpl/aux_/preprocessed/msvc60/iter_fold_impl.hpp>


text:
/**
 * Copyright Soramitsu Co., Ltd. 2018 All Rights Reserved.
 * http://soramitsu.co.jp
 *
 * 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.
 */

#ifndef IROHA_VISITOR_HPP
#define IROHA_VISITOR_HPP

#include <boost/variant/apply_visitor.hpp>  // for boost::apply_visitor
#include <type_traits>                      // for std::decay
#include <utility>                          // for std::forward

namespace iroha {

  template <typename... Lambdas>
  struct lambda_visitor;

  template <typename Lambda, typename... Lambdas>
  struct lambda_visitor<Lambda, Lambdas...>
      : public Lambda, public lambda_visitor<Lambdas...> {
    using Lambda::operator();
    using lambda_visitor<Lambdas...>::operator();

    lambda_visitor(Lambda lambda, Lambdas... lambdas)
        : Lambda(lambda), lambda_visitor<Lambdas...>(lambdas...) {}
  };

  template <typename Lambda>
  struct lambda_visitor<Lambda> : public Lambda {
    using Lambda::operator();

    lambda_visitor(Lambda lambda) : Lambda(lambda) {}
  };

  /**
   * @brief Convenient in-place compile-time visitor creation, from a set of
   * lambdas
   *
   * @code
   * make_visitor([](int a){ return 1; },
   *              [](std::string b) { return 2; });
   * @nocode
   *
   * is essentially the same as
   *
   * @code
   * struct visitor : public boost::static_visitor<int> {
   *   int operator()(int a) { return 1; }
   *   int operator()(std::string b) { return 2; }
   * }
   * @nocode
   *
   * @param lambdas
   * @return visitor
   */
  template <class... Fs>
  constexpr auto make_visitor(Fs &&... fs) {
    using visitor_type = lambda_visitor<std::decay_t<Fs>...>;
    return visitor_type(std::forward<Fs>(fs)...);
  }

  /**
   * @brief Inplace visitor for boost::variant.
   * @code
   *   boost::variant<int, std::string> value = "1234";
   *   ...
   *   visit_in_place(value,
   *                  [](int v) { std::cout << "(int)" << v; },
   *                  [](std::string v) { std::cout << "(string)" << v;}
   *                  );
   * @nocode
   *
   * @param variant
   * @param lambdas
   * @param lambdas
   */
  template <typename TVariant, typename... TVisitors>
  constexpr decltype(auto) visit_in_place(TVariant &&variant, TVisitors &&... visitors) {
    return boost::apply_visitor(make_visitor(visitors...), std::forward<TVariant>(variant));
  }
}  // namespace iroha

#endif  // IROHA_VISITOR_HPP

max_stars_repo_name: truongnmt/iroha

max_issues_repo_name: truongnmt/iroha

max_forks_repo_name: truongnmt/iroha