text:
#include<stdio.h>
#include"mex.h"

#include "EdgeSE3ProjectDirectWithDirstortG2oLM.cpp"

#include <Eigen/Core>
#include <Eigen/Geometry>

// 李群李代数 库
#include "sophus/so3.hpp"
#include "sophus/se3.hpp"

#include <g2o/core/base_unary_edge.h>
#include <g2o/core/block_solver.h>
#include <g2o/core/optimization_algorithm_levenberg.h>
#include <g2o/solvers/dense/linear_solver_dense.h>
#include <g2o/core/robust_kernel.h>
#include <g2o/core/robust_kernel_impl.h>
#include <g2o/types/sba/types_six_dof_expmap.h>

using namespace g2o;

// 一次测量的值，包括一个世界坐标系下三维点,以及投影的对应的图像
struct Measurement
{
    Measurement ( Eigen::Vector3d p, float * im) : pos_world (p), image(im){}
    Eigen::Vector3d pos_world;
    float* image;
};

bool calibrationEstimationDirectG2OLM ( const vector<Measurement>& measurements, double *D, double *K, double *imageSize, Eigen::Isometry3d& Tcw, bool verbose, int max_iter)
{
    // 初始化g2o
    typedef g2o::BlockSolver<g2o::BlockSolverTraits<6,1>> DirectBlock;  // 求解的向量是6＊1的
    DirectBlock::LinearSolverType* linearSolver = new g2o::LinearSolverDense< DirectBlock::PoseMatrixType > ();
    DirectBlock* solver_ptr = new DirectBlock ( linearSolver );
    // g2o::OptimizationAlgorithmGaussNewton* solver = new g2o::OptimizationAlgorithmGaussNewton( solver_ptr ); // G-N
    g2o::OptimizationAlgorithmLevenberg* solver = new g2o::OptimizationAlgorithmLevenberg ( solver_ptr ); // L-M
    g2o::SparseOptimizer optimizer;
    optimizer.setAlgorithm ( solver );
    optimizer.setVerbose( verbose );

    // 添加顶点
    g2o::VertexSE3Expmap* pose = new g2o::VertexSE3Expmap();
    pose->setEstimate ( g2o::SE3Quat ( Tcw.rotation(), Tcw.translation() ) );
    pose->setId ( 0 );
    optimizer.addVertex ( pose );

    // 添加边
    int id=1;
    for ( Measurement m: measurements )
    {
        EdgeSE3ProjectDirectWithDirstortG2oLM* edge = new EdgeSE3ProjectDirectWithDirstortG2oLM (
            m.pos_world,
            D, K, imageSize, m.image
        );
        edge->setVertex (0, pose );
        edge->setMeasurement (0);
        edge->setInformation ( Eigen::Matrix<double,1,1>::Identity() );
        edge->setId ( id++ );
        optimizer.addEdge ( edge );
    }
    // cout<<"edges in graph: "<<optimizer.edges().size() <<endl;
    mexPrintf("edges in graph: %d\n", optimizer.edges().size());
    optimizer.initializeOptimization();
    optimizer.optimize ( max_iter );
    Tcw = pose->estimate();
}

void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]){
    // nlhs represent the number of parameters of the output
    // plhs is a array of the mxarray pointers, each pointing to the output
    // nrhs represents the number of parameters of the input
    // prhs is a array of the mxarray pointers, each pointing to the input

    // prhs[0], 6x1 matrix
    // prhs[1], Mx1 cell, each cell with NX3 points
    // prhs[2], Mx1 cell, each cell with PXQ single matrix
    // prhs[3], 1x2, or 1x3, or 1x4, or 1x5 matrix
    // prhs[4], 3x3 matrix
    // prhs[5], 1x2 matrix
    // prhs[6], 1x1 matrix
    // prhs[7], 1x1 matrix

    if(nrhs < 8){
        mexErrMsgIdAndTxt( "EdgeSE3ProjectDirectWithDirstortJacobian:invalidNumInputs", "at least 8 input arguments required");
        return;
    }

    // get the init xi
    const size_t *dimArrayOfInitXi = mxGetDimensions(prhs[0]);
    size_t sizeRowsInitXi = *(dimArrayOfInitXi + 0);
    size_t sizeColsInitXi = *(dimArrayOfInitXi + 1);
    if(sizeRowsInitXi != 6 || sizeColsInitXi != 1){
        mexErrMsgIdAndTxt( "EdgeSE3ProjectDirectWithDirstortJacobian:invalidInputs", "the 1st param should be 6x1");
        return;
    }
    double *ptrInitXi = (double *)(mxGetPr(prhs[0]));
    Eigen::Matrix<double, 6, 1> initXi;
    for(int i = 0; i < 6; i++){
        initXi(i, 0) = *(ptrInitXi + i);
    }
    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    // get the chessboard lidar points
    const size_t *dimArrayOfChessboardLidarPoints = mxGetDimensions(prhs[1]);
    size_t sizeRowsCellChessboardLidarPoints = *(dimArrayOfChessboardLidarPoints + 0);
    size_t sizeColsCellChessboardLidarPoints = *(dimArrayOfChessboardLidarPoints + 1);
    if(sizeColsCellChessboardLidarPoints != 1){
        mexErrMsgIdAndTxt( "EdgeSE3ProjectDirectWithDirstortJacobian:invalidInputs", "the 2st param should be Mx1 cell");
        return;
    }
    std::vector<std::vector<Eigen::Vector3d>>chessboardLidarPoints;
    size_t total_measurements = 0;
    for(int i = 0; i < sizeRowsCellChessboardLidarPoints; i++){
        // get the dimensions of each cell
        mxArray *cur_mArray_chessboardLidarPoints = mxGetCell(prhs[1], i);
        const size_t *dimArrayOfCurrentChessboardLidarPoints = mxGetDimensions(cur_mArray_chessboardLidarPoints);
        size_t sizeRowsCurrentChessboardLidarPoints = *(dimArrayOfCurrentChessboardLidarPoints + 0);
        size_t sizeColsCurrentChessboardLidarPoints = *(dimArrayOfCurrentChessboardLidarPoints + 1);
        if(sizeColsCurrentChessboardLidarPoints != 3){
            mexErrMsgIdAndTxt( "EdgeSE3ProjectDirectWithDirstortJacobian:invalidInputs", "the 2st param should be cell with Px3");
            return;
        }
        float *cur_ptr_chessboardLidarPoints = (float *)(mxGetPr(cur_mArray_chessboardLidarPoints));
        std::vector<Eigen::Vector3d>currentCellChessboardLidarPoints;
        total_measurements = total_measurements + sizeRowsCurrentChessboardLidarPoints;
        for(int j = 0; j < sizeRowsCurrentChessboardLidarPoints; j++){
            Eigen::Vector3d current_chessboardLidarPoint;
            current_chessboardLidarPoint[0] = cur_ptr_chessboardLidarPoints[0 * sizeRowsCurrentChessboardLidarPoints + j];
            current_chessboardLidarPoint[1] = cur_ptr_chessboardLidarPoints[1 * sizeRowsCurrentChessboardLidarPoints + j];
            current_chessboardLidarPoint[2] = cur_ptr_chessboardLidarPoints[2 * sizeRowsCurrentChessboardLidarPoints + j];
            currentCellChessboardLidarPoints.push_back(current_chessboardLidarPoint);
        }
        chessboardLidarPoints.push_back(currentCellChessboardLidarPoints);
    }

    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    // get the chessoard_Dt_Mask
    const size_t *dimArrayOfChessboardMaskDT = mxGetDimensions(prhs[2]);
    size_t sizeRowsCellChessboardMaskDT = *(dimArrayOfChessboardMaskDT + 0);
    size_t sizeColsCellChessboardMaskDT = *(dimArrayOfChessboardMaskDT + 1);
    if(sizeColsCellChessboardMaskDT != 1){
        mexErrMsgIdAndTxt( "EdgeSE3ProjectDirectWithDirstortJacobian:invalidInputs", "the 3st param should be Mx1 cell");
        return;
    }
    std::vector<float *>chessboardMaskDTPtr;
    for(int i = 0; i < sizeRowsCellChessboardMaskDT; i++){
        // get the dimensions of each cell
        mxArray *cur_mArray_chessboardMaskDT = mxGetCell(prhs[2], i);

        float *cur_ptr_chessboardMaskDT = (float *)(mxGetPr(cur_mArray_chessboardMaskDT));
        chessboardMaskDTPtr.push_back(cur_ptr_chessboardMaskDT);
    }

    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    // get the distortion parameters
    const size_t *dimArrayOfD = mxGetDimensions(prhs[3]);
    size_t sizeRowsD = *(dimArrayOfD + 0);
    size_t sizeColsD = *(dimArrayOfD + 1);
    if(sizeRowsD != 1 || sizeColsD != 5){
        mexErrMsgIdAndTxt( "EdgeSE3ProjectDirectWithDirstortJacobian:invalidInputs", "the 4st param should be 1x5 array");
        return;
    }
    double *ptrD = (double *)(mxGetPr(prhs[3]));

    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    // get the intrinsic parameters K
    const size_t *dimArrayOfK = mxGetDimensions(prhs[4]);
    size_t sizeRowsK = *(dimArrayOfK + 0);
    size_t sizeColsK = *(dimArrayOfK + 1);
    if(sizeRowsK != 3 || sizeColsK != 3){
        mexErrMsgIdAndTxt( "EdgeSE3ProjectDirectWithDirstortJacobian:invalidInputs", "the 5st param should be 3X3 array");
        return;
    }
    double *ptrK = (double *)(mxGetPr(prhs[4]));

    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    // get the imageSize
    const size_t *dimArrayOfImageSize = mxGetDimensions(prhs[5]);
    size_t sizeRowsImageSize = *(dimArrayOfImageSize + 0);
    size_t sizeColsImageSize = *(dimArrayOfImageSize + 1);
    if(sizeRowsImageSize != 1 || sizeColsImageSize != 2){
        mexErrMsgIdAndTxt( "EdgeSE3ProjectDirectWithDirstortJacobian:invalidInputs", "the 6st param should be 1x2 array");
        return;
    }
    double *ptrImageSize = (double *)(mxGetPr(prhs[5]));

    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    // get the verbose params
    const size_t *dimArrayOfVerbose = mxGetDimensions(prhs[6]);
    size_t sizeRowsVerbose= *(dimArrayOfVerbose + 0);
    size_t sizeColsVerbose = *(dimArrayOfVerbose + 1);
    if(sizeRowsVerbose != 1 || sizeColsVerbose != 1){
        mexErrMsgIdAndTxt( "EdgeSE3ProjectDirectWithDirstortJacobian:invalidInputs", "the 7st param should be 1x1 array");
        return;
    }
    double *ptrVerbose = (double *)(mxGetPr(prhs[6]));

    //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    // get the maxIter params
    const size_t *dimArrayOfMaxIter = mxGetDimensions(prhs[6]);
    size_t sizeRowsMaxIter = *(dimArrayOfMaxIter + 0);
    size_t sizeColsMaxIter = *(dimArrayOfMaxIter + 1);
    if(sizeRowsMaxIter != 1 || sizeColsMaxIter != 1){
        mexErrMsgIdAndTxt( "EdgeSE3ProjectDirectWithDirstortJacobian:invalidInputs", "the 8st param should be 1x1 array");
        return;
    }
    double *ptrMaxIter = (double *)(mxGetPr(prhs[7]));

    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
    std::vector<Measurement> measurements;
    std::vector<bool> isInners;
    Eigen::Isometry3d Tlc = Eigen::Isometry3d::Identity();

    // get the init tcl from the init Xi
    Sophus::SE3<double> SE3 = Sophus::SE3<double>::exp(initXi);
    Eigen::Matrix<double, 4, 4> SE3Matrix = SE3.matrix();
    Tlc.matrix() << SE3Matrix;

    for(int i = 0; i < chessboardMaskDTPtr.size(); i++){
        float *current_chessboardMaskDTPtr = chessboardMaskDTPtr[i];
        std::vector<Eigen::Vector3d> current_chessboardLidarPoints = chessboardLidarPoints[i];
        for(int j = 0; j < current_chessboardLidarPoints.size(); j++){
            measurements.push_back(Measurement(current_chessboardLidarPoints[j], current_chessboardMaskDTPtr));
        }
    }
    isInners.resize(measurements.size());

    calibrationEstimationDirectG2OLM ( measurements, ptrD, ptrK, ptrImageSize, Tlc, bool(ptrVerbose[0]), int(ptrMaxIter[0]));

    Sophus::SE3<double> SE3_Rt(Tlc.rotation(), Tlc.translation());

    // cout<<"SE3 = "<<endl<<SE3_Rt.matrix()<<endl;

    Eigen::Matrix<double,6,1> se3_out = SE3_Rt.log();

    // the output xi will be 1x6
    size_t dimArrayOfXi[2] = { 1, 6 };
    plhs[0] = mxCreateNumericArray(2, dimArrayOfXi, mxDOUBLE_CLASS, mxREAL);
    double *out_xi = (double *)mxGetData(plhs[0]);

    for(int i = 0; i < 6; i++){
        out_xi[i] = se3_out(i,0);
    }
}


text:
// HEADER
#include <path_follower/controller/robotcontroller_kinematic_HBZ.h>

// THIRD PARTY
#include <nav_msgs/Path.h>
#include <visualization_msgs/Marker.h>
#include <visualization_msgs/MarkerArray.h>

// PROJECT
#include <path_follower/parameters/path_follower_parameters.h>
#include <path_follower/utils/cubic_spline_interpolation.h>
#include <cslibs_navigation_utilities/MathHelper.h>
#include <path_follower/utils/pose_tracker.h>
#include <path_follower/utils/visualizer.h>
#include <path_follower/utils/obstacle_cloud.h>
#include <path_follower/collision_avoidance/collision_avoider.h>
#include <path_follower/factory/controller_factory.h>

// SYSTEM
#include <cmath>
#include <deque>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <boost/algorithm/clamp.hpp>
#include <pcl_ros/point_cloud.h>
#include <path_follower/factory/controller_factory.h>

REGISTER_ROBOT_CONTROLLER(RobotController_Kinematic_HBZ, kinematic_hbz, default_collision_avoider);

using namespace Eigen;

RobotController_Kinematic_HBZ::RobotController_Kinematic_HBZ():
    RobotController(),
    cmd_(this),
    vn_(0),
    Ts_(0.02),
    ind_(0),
    proj_ind_(0),
    xe_(0),
    ye_(0),
    theta_e_(0),
    x_meas_(0),
    y_meas_(0),
    theta_meas_(0),
    angular_vel_(0),
    Vl_(0),
    Vr_(0),
    delta_(0)
{

    wheel_velocities_ = nh_.subscribe<std_msgs::Float64MultiArray>("wheel_velocities", 10,
                                                                   &RobotController_Kinematic_HBZ::WheelVelocities, this);

}

void RobotController_Kinematic_HBZ::stopMotion()
{

    cmd_.speed = 0;
    cmd_.direction_angle = 0;
    cmd_.rotation = 0;

    MoveCommand mcmd = cmd_;
    publishMoveCommand(mcmd);
}

void RobotController_Kinematic_HBZ::initialize()
{
    RobotController::initialize();

    //reset the index of the current point on the path
    ind_ = 0;

    theta_e_ = 0.0;

    // desired velocity
    vn_ = std::min(PathFollowerParameters::getInstance()->max_velocity(), velocity_);
    ROS_DEBUG_STREAM("velocity_: " << velocity_ << ", vn: " << vn_);

}

void RobotController_Kinematic_HBZ::WheelVelocities(const std_msgs::Float64MultiArray::ConstPtr& array)
{
    double frw = array->data[0];
    double flw = array->data[1];
    double brw = array->data[2];
    double blw = array->data[3];

    Vl_ = (flw + blw)/2.0;
    Vr_ = (frw + brw)/2.0;
}

void RobotController_Kinematic_HBZ::start()
{

}

void RobotController_Kinematic_HBZ::reset()
{
    RobotController::reset();
}

void RobotController_Kinematic_HBZ::setPath(Path::Ptr path)
{
    RobotController::setPath(path);
}

double RobotController_Kinematic_HBZ::computeSpeed()
{
    double v = vn_;

    double V1 = 1.0/2.0*(std::pow(xe_,2) + std::pow(ye_,2) + 1.0/opt_.lambda()*fabs(sin(theta_e_-delta_)));

    if(angular_vel_ > 0){

        if(V1 >= opt_.epsilon()){
            v = (-opt_.alpha_r()*opt_.y_ICR_l()*vn_)/(opt_.y_ICR_r() - opt_.y_ICR_l());
        }
        else if(V1 < opt_.epsilon()){
            v = (opt_.alpha_r()*vn_)/(1 + std::fabs(opt_.y_ICR_r()*path_interpol.curvature(ind_)));
        }
    }
    else if(angular_vel_ <= 0){

        if(V1 >= opt_.epsilon()){
            v = (opt_.alpha_l()*opt_.y_ICR_r()*vn_)/(opt_.y_ICR_r() - opt_.y_ICR_l());
        }
        else if(V1 < opt_.epsilon()){
            v = (opt_.alpha_l()*vn_)/(1 + std::fabs(opt_.y_ICR_l()*path_interpol.curvature(ind_)));
        }
    }

    return v;
}

RobotController::MoveCommandStatus RobotController_Kinematic_HBZ::computeMoveCommand(MoveCommand *cmd)
{
    // omni drive can rotate.
    *cmd = MoveCommand(true);

    if(path_interpol.n() < 2) {
        ROS_ERROR("[Line] path is too short (N = %d)", (int) path_interpol.n());

        stopMotion();
        return MoveCommandStatus::REACHED_GOAL;
    }

    /// get the pose as pose(0) = x, pose(1) = y, pose(2) = theta
    Eigen::Vector3d current_pose = pose_tracker_->getRobotPose();

    const geometry_msgs::Twist v_meas_twist = pose_tracker_->getVelocity();
    //get the robot's current angular velocity
    angular_vel_ = pose_tracker_->getVelocity().angular.z;

    double v_meas = getDirSign() * sqrt(v_meas_twist.linear.x * v_meas_twist.linear.x
            + v_meas_twist.linear.y * v_meas_twist.linear.y);

    x_meas_ = current_pose[0];
    y_meas_ = current_pose[1];
    theta_meas_ = current_pose[2];
    ///***///

    RobotController::findOrthogonalProjection();

    if(RobotController::isGoalReached(cmd)){
       return RobotController::MoveCommandStatus::REACHED_GOAL;
    }

    ///compute the control for the current point on the path

    //robot direction angle in path coordinates
    theta_e_ = MathHelper::AngleDelta(path_interpol.theta_p(ind_), theta_meas_);

    //robot position vector module
    double r = hypot(x_meas_ - path_interpol.p(ind_), y_meas_ - path_interpol.q(ind_));

    //robot position vector angle in world coordinates
    double theta_r = atan2(y_meas_ - path_interpol.q(ind_), x_meas_ - path_interpol.p(ind_));

    //robot position vector angle in path coordinates
    double delta_theta = MathHelper::AngleDelta(path_interpol.theta_p(ind_), theta_r);

    //current robot position in path coordinates
    xe_ = r * cos(delta_theta);
    ye_ = r * sin(delta_theta);

    ///***///

    ///Check the driving direction, and set the complementary angle in path coordinates, if driving backwards

    if (getDirSign() < 0.0) {
        theta_e_ = MathHelper::NormalizeAngle(M_PI + theta_e_);
        ROS_WARN_THROTTLE(1, "Driving backwards...");
    }

    ///***///

    ///Lyapunov-curvature speed control - skid steering

       double v = computeSpeed();

    ///***///

    ///compute the next point on the path

    double s_old = path_interpol.s_new();

    //calculate the speed of the "virtual vehicle"
    double s_prim_tmp = v * cos(theta_e_) + opt_.x_ICR()*angular_vel_*sin(theta_e_) + opt_.k1() * xe_;
    path_interpol.set_s_prim(s_prim_tmp > 0 ? s_prim_tmp : 0);

    //approximate the first derivative and calculate the next point
    double s_temp = Ts_*path_interpol.s_prim() + s_old;
    path_interpol.set_s_new(s_temp > 0 ? s_temp : 0);

    ///***///

    double exp_factor = RobotController::exponentialSpeedControl();
    v = v * exp_factor;

    //cmd_.speed = getDirSign()*std::max((double)PathFollowerParameters::getInstance()->min_velocity(), fabs(v));
    //allow the speed to be reduced to zero
    cmd_.speed = getDirSign() * v;

    ///***///

    ///Compute the delta_ and its derivative

    double delta_old = delta_;

    delta_ = MathHelper::AngleClamp(-opt_.theta_a()*tanh(ye_));

    double delta_prim = (delta_ - delta_old)/Ts_;
    ///***///

    ///Direction control

    cmd_.direction_angle = 0;

    double theta_diff = MathHelper::AngleDelta(delta_, theta_e_);

    double trig_ratio = (fabs(sin(theta_diff)))/(sin(theta_diff)*cos(theta_diff));
    double omega = delta_prim
            - opt_.lambda()*trig_ratio*ye_*v*sin(theta_e_)
            + opt_.lambda()*trig_ratio*ye_*opt_.x_ICR()*angular_vel_*cos(theta_e_)
            - opt_.k2()*trig_ratio*(theta_diff)*(theta_diff);

    omega += path_interpol.curvature(ind_)*path_interpol.s_prim();
    cmd_.rotation = boost::algorithm::clamp(omega, -opt_.max_angular_velocity(), opt_.max_angular_velocity());

    ///***///

    ///Plot the moving reference frame together with position vector and error components

    if (visualizer_->MarrayhasSubscriber()) {
        visualizer_->drawFrenetSerretFrame(getFixedFrame(), 0, current_pose, xe_, ye_, path_interpol.p(ind_),
                                           path_interpol.q(ind_), path_interpol.theta_p(ind_));
    }

    ///***///

    ///Compute the index of the new point

    double s_diff = std::numeric_limits<double>::max();
    uint old_ind = ind_;
    //if the robot is not moving, the Frenet-Serret frame should also not move
    if(fabs(v_meas) < 1e-3){
        s_diff = 0.0;
    }

    for (unsigned int i = old_ind; i < path_interpol.n(); i++){

        double s_diff_curr = std::abs(path_interpol.s_new() - path_interpol.s(i));

        if(s_diff_curr < s_diff){

            s_diff = s_diff_curr;
            ind_ = i;

        }

    }

    if(old_ind != ind_) {
        path_->fireNextWaypointCallback();
    }

    ///***///

    if (visualizer_->hasSubscriber()) {
        visualizer_->drawSteeringArrow(pose_tracker_->getFixedFrameId(), 1, pose_tracker_->getRobotPoseMsg(), cmd_.direction_angle, 0.2, 1.0, 0.2);
    }

    ///***///

    *cmd = cmd_;

    return MoveCommandStatus::OKAY;
}

void RobotController_Kinematic_HBZ::publishMoveCommand(const MoveCommand &cmd) const
{
    geometry_msgs::Twist msg;
    msg.linear.x  = cmd.getVelocity();
    msg.linear.y  = 0;
    msg.angular.z = cmd.getRotationalVelocity();

    cmd_pub_.publish(msg);
}

ext: cpp

max_stars_repo_name: sunarditay/gerona

max_issues_repo_name: sunarditay/gerona

max_forks_repo_name: sunarditay/gerona


text:
/*
 * $Id: vector.cc,v 1.7 2005/10/06 23:28:03 julianc Exp $
 *
 * Copyright (C) 1997 Todd Veldhuizen <tveldhui@oonumerics.org>
 * All rights reserved.  Please see <blitz/blitz.h> for terms and
 * conditions of use.
 *
 */

#ifndef BZ_VECTOR_CC
#define BZ_VECTOR_CC

#include <blitz/vector.h>
#include <blitz/update.h>

BZ_NAMESPACE(blitz)

template<typename P_numtype>
Vector<P_numtype> Vector<P_numtype>::copy() const
{
    Vector<P_numtype> newCopy(length_);

    if (stride_ == 1)
    {
        memcpy(newCopy.data(), data(), length_ * sizeof(P_numtype));
    }
    else {
        for (int i=0; i < length_; ++i)
        {
            // Can assume that newCopy has unit stride, and hence use
            // the operator(), which assumes unit stride.  Since this
            // vector doesn't have unit stride, use [].
            newCopy(i) = (*this)[i];
        }
    }

    return newCopy;
}

template<typename P_numtype>
void Vector<P_numtype>::makeUnique()
{
    if ((stride_ != 1) || (this->numReferences() > 1))
    {
        Vector<P_numtype> tmp = copy();
        reference(tmp);
    }
}

template<typename P_numtype>
void Vector<P_numtype>::reference(Vector<P_numtype>& x)
{
    MemoryBlockReference<P_numtype>::changeBlock(x);
    length_ = x.length_;
    stride_ = x.stride_;
}

template<typename P_numtype>
void Vector<P_numtype>::resize(int length)
{
    if (length != length_)
    {
        MemoryBlockReference<P_numtype>::newBlock(length);
        length_ = length;
        stride_ = 1;
    }
}

template<typename P_numtype>
void Vector<P_numtype>::resizeAndPreserve(int newLength)
{
    Vector<P_numtype> newVector(newLength);

    if (newLength > length_)
        newVector(Range(0,length_-1)) = (*this);
    else
        newVector(Range(0,newLength-1)) = (*this);

    reference(newVector);
}

/*****************************************************************************
 * Assignment operators with vector expression operand
 */

template<typename P_numtype> template<typename P_expr, typename P_updater>
inline
void Vector<P_numtype>::_bz_assign(P_expr expr, P_updater)
{
    BZPRECONDITION(expr.length(length_) == length_);

    // If all vectors in the expression, plus the vector to which the
    // result is being assigned have unit stride, then avoid stride
    // calculations.
    if ((stride_ == 1) && (expr._bz_hasFastAccess()))
    {
#ifndef BZ_PARTIAL_LOOP_UNROLL
        for (int i=0; i < length_; ++i)
            P_updater::update(data_[i], expr._bz_fastAccess(i));
#else
        // Unwind the inner loop, four elements at a time.
        int leftover = length_ & 0x03;

        int i=0;
        for (; i < leftover; ++i)
            P_updater::update(data_[i], expr._bz_fastAccess(i));

        for (; i < length_; i += 4)
        {
            // Common subexpression elimination: avoid doing i+1, i+2, i+3
            // multiple times (compilers *won't* do this cse automatically)
            int t1 = i+1;
            int t2 = i+2;
            int t3 = i+3;

            _bz_typename P_expr::T_numtype tmp1, tmp2, tmp3, tmp4;

            // Reading all the results first avoids aliasing
            // ambiguities, and provides more flexibility in
            // register allocation & instruction scheduling
            tmp1 = expr._bz_fastAccess(i);
            tmp2 = expr._bz_fastAccess(BZ_NO_PROPAGATE(t1));
            tmp3 = expr._bz_fastAccess(BZ_NO_PROPAGATE(t2));
            tmp4 = expr._bz_fastAccess(BZ_NO_PROPAGATE(t3));

            P_updater::update(data_[i], BZ_NO_PROPAGATE(tmp1));
            P_updater::update(data_[t1], tmp2);
            P_updater::update(data_[t2], tmp3);
            P_updater::update(data_[t3], tmp4);
        }
#endif
    }
    else {
        // Not all unit strides -- have to access all the vector elements
        // as data_[i*stride_], which is slower.
        for (int i=0; i < length_; ++i)
            P_updater::update((*this)[i], expr[i]);
    }
}

template<typename P_numtype> template<typename P_expr>
inline Vector<P_numtype>& Vector<P_numtype>::operator=(_bz_VecExpr<P_expr> expr)
{
    BZPRECONDITION(expr.length(length_) == length_);

    // If all vectors in the expression, plus the vector to which the
    // result is being assigned have unit stride, then avoid stride
    // calculations.
    if ((stride_ == 1) && (expr._bz_hasFastAccess()))
    {
#ifndef BZ_PARTIAL_LOOP_UNROLL
        for (int i=0; i < length_; ++i)
            data_[i] = (P_numtype)expr._bz_fastAccess(i);
#else
        // Unwind the inner loop, four elements at a time.
        int leftover = length_ & 3;

        int i=0;
        for (; i < leftover; ++i)
            data_[i] = (P_numtype)expr._bz_fastAccess(i);

        for (; i < length_; i += 4)
        {
            // Common subexpression elimination: avoid doing i+1, i+2, i+3
            // multiple times (compilers *won't* do this cse automatically)
            int t1 = i+1;
            int t2 = i+2;
            int t3 = i+3;

            _bz_typename P_expr::T_numtype tmp1, tmp2, tmp3, tmp4;

            // Reading all the results first avoids aliasing
            // ambiguities, and provides more flexibility in
            // register allocation & instruction scheduling
            tmp1 = expr._bz_fastAccess(i);
            tmp2 = expr._bz_fastAccess(BZ_NO_PROPAGATE(t1));
            tmp3 = expr._bz_fastAccess(BZ_NO_PROPAGATE(t2));
            tmp4 = expr._bz_fastAccess(BZ_NO_PROPAGATE(t3));

            data_[i] = (P_numtype)BZ_NO_PROPAGATE(tmp1);
            data_[t1] = (P_numtype)tmp2;
            data_[t2] = (P_numtype)tmp3;
            data_[t3] = (P_numtype)tmp4;
        }
#endif
    }
    else {
        // Not all unit strides -- have to access all the vector elements
        // as data_[i*stride_], which is slower.
        for (int i=0; i < length_; ++i)
            (*this)[i] = (P_numtype)expr[i];
    }
    return *this;
}

#ifdef BZ_PARTIAL_LOOP_UNROLL

#define BZ_VECTOR_ASSIGN(op)                                            \
template<typename P_numtype> template<typename P_expr>                        \
inline Vector<P_numtype>& Vector<P_numtype>::                           \
    operator op (_bz_VecExpr<P_expr> expr)                              \
{                                                                       \
    BZPRECONDITION(expr.length(length_) == length_);                    \
    if ((stride_ == 1) && (expr._bz_hasFastAccess()))                   \
    {                                                                   \
        int leftover = length_ & 0x03;                                  \
                                                                        \
        int i=0;                                                        \
        for (; i < leftover; ++i)                                       \
            data_[i] op expr._bz_fastAccess(i);                         \
                                                                        \
        for (; i < length_; i += 4)                                     \
        {                                                               \
            int t1 = i+1;                                               \
            int t2 = i+2;                                               \
            int t3 = i+3;                                               \
                                                                        \
            _bz_typename P_expr::T_numtype tmp1, tmp2, tmp3, tmp4;      \
                                                                        \
            tmp1 = expr._bz_fastAccess(i);                              \
            tmp2 = expr._bz_fastAccess(t1);                             \
            tmp3 = expr._bz_fastAccess(t2);                             \
            tmp4 = expr._bz_fastAccess(t3);                             \
                                                                        \
            data_[i] op tmp1;                                           \
            data_[t1] op tmp2;                                          \
            data_[t2] op tmp3;                                          \
            data_[t3] op tmp4;                                          \
        }                                                               \
    }                                                                   \
    else {                                                              \
        for (int i=0; i < length_; ++i)                                 \
            (*this)[i] op expr[i];                                      \
    }                                                                   \
    return *this;                                                       \
}

#else   // not BZ_PARTIAL_LOOP_UNROLL

#ifdef BZ_ALTERNATE_FORWARD_BACKWARD_TRAVERSALS

/*
 * NEEDS_WORK: need to incorporate BZ_NO_PROPAGATE here.  This
 * will require doing away with the macro BZ_VECTOR_ASSIGN, and
 * adopting an approach similar to that used in arrayeval.cc.
 */

#define BZ_VECTOR_ASSIGN(op)                                            \
template<typename P_numtype> template<typename P_expr>                        \
inline Vector<P_numtype>& Vector<P_numtype>::                           \
    operator op (_bz_VecExpr<P_expr> expr)                              \
{                                                                       \
    BZPRECONDITION(expr.length(length_) == length_);                    \
    static int traversalOrder = 0;                                      \
    if ((stride_ == 1) && (expr._bz_hasFastAccess()))                   \
    {                                                                   \
        if (traversalOrder & 0x01)                                      \
            for (int i=length_-1; i >= 0; --i)                          \
                data_[i] op expr._bz_fastAccess(i);                     \
        else                                                            \
            for (int i=0; i < length_; ++i)                             \
                data_[i] op expr._bz_fastAccess(i);                     \
    }                                                                   \
    else {                                                              \
        for (int i=0; i < length_; ++i)                                 \
            (*this)[i] op expr[i];                                      \
    }                                                                   \
    traversalOrder ^= 0x01;                                             \
    return *this;                                                       \
}

#else   // not BZ_ALTERNATE_FORWARD_BACKWARD_TRAVERSALS

#define BZ_VECTOR_ASSIGN(op)                                            \
template<typename P_numtype> template<typename P_expr>                        \
inline Vector<P_numtype>& Vector<P_numtype>::                           \
    operator op (_bz_VecExpr<P_expr> expr)                              \
{                                                                       \
    BZPRECONDITION(expr.length(length_) == length_);                    \
    if ((stride_ == 1) && (expr._bz_hasFastAccess()))                   \
    {                                                                   \
        for (int i=0; i < length_; ++i)                                 \
            data_[i] op expr._bz_fastAccess(i);                         \
    }                                                                   \
    else {                                                              \
        for (int i=0; i < length_; ++i)                                 \
            (*this)[i] op expr[i];                                      \
    }                                                                   \
    return *this;                                                       \
}
#endif // BZ_ALTERNATE_FORWARD_BACKWARD_TRAVERSALS
#endif // BZ_PARTIAL_LOOP_UNROLL

BZ_VECTOR_ASSIGN(+=)
BZ_VECTOR_ASSIGN(-=)
BZ_VECTOR_ASSIGN(*=)
BZ_VECTOR_ASSIGN(/=)
BZ_VECTOR_ASSIGN(%=)
BZ_VECTOR_ASSIGN(^=)
BZ_VECTOR_ASSIGN(&=)
BZ_VECTOR_ASSIGN(|=)
BZ_VECTOR_ASSIGN(>>=)
BZ_VECTOR_ASSIGN(<<=)

#if NOT_DEFINED

template<typename P_numtype> template<typename P_expr>
inline Vector<P_numtype>& Vector<P_numtype>::operator+=(_bz_VecExpr<P_expr> expr)
{
    _bz_assign(expr, _bz_plus_update<P_numtype,
        _bz_typename P_expr::T_numtype>());
    return *this;
}

template<typename P_numtype> template<typename P_expr>
inline Vector<P_numtype>& Vector<P_numtype>::operator-=(_bz_VecExpr<P_expr> expr)
{
    _bz_assign(expr, _bz_minus_update<P_numtype,
        _bz_typename P_expr::T_numtype>());
    return *this;
}

template<typename P_numtype> template<typename P_expr>
inline Vector<P_numtype>& Vector<P_numtype>::operator*=(_bz_VecExpr<P_expr> expr)
{
    _bz_assign(expr, _bz_multiply_update<P_numtype,
        _bz_typename P_expr::T_numtype>());
    return *this;
}

template<typename P_numtype> template<typename P_expr>
inline Vector<P_numtype>& Vector<P_numtype>::operator/=(_bz_VecExpr<P_expr> expr)
{
    _bz_assign(expr, _bz_divide_update<P_numtype,
        _bz_typename P_expr::T_numtype>());
    return *this;
}

template<typename P_numtype> template<typename P_expr>
inline Vector<P_numtype>& Vector<P_numtype>::operator%=(_bz_VecExpr<P_expr> expr)
{
    _bz_assign(expr, _bz_mod_update<P_numtype,
        _bz_typename P_expr::T_numtype>());
    return *this;
}

template<typename P_numtype> template<typename P_expr>
inline Vector<P_numtype>& Vector<P_numtype>::operator^=(_bz_VecExpr<P_expr> expr)
{
    _bz_assign(expr, _bz_xor_update<P_numtype,
        _bz_typename P_expr::T_numtype>());
    return *this;
}

template<typename P_numtype> template<typename P_expr>
inline Vector<P_numtype>& Vector<P_numtype>::operator&=(_bz_VecExpr<P_expr> expr)
{
    _bz_assign(expr, _bz_bitand_update<P_numtype,
        _bz_typename P_expr::T_numtype>());
    return *this;
}

template<typename P_numtype> template<typename P_expr>
inline Vector<P_numtype>& Vector<P_numtype>::operator|=(_bz_VecExpr<P_expr> expr)
{
    _bz_assign(expr, _bz_bitor_update<P_numtype,
        _bz_typename P_expr::T_numtype>());
    return *this;
}

template<typename P_numtype> template<typename P_expr>
inline Vector<P_numtype>& Vector<P_numtype>::operator>>=(_bz_VecExpr<P_expr> expr)
{
    _bz_assign(expr, _bz_shiftr_update<P_numtype,
        _bz_typename P_expr::T_numtype>());
    return *this;
}

template<typename P_numtype> template<typename P_expr>
inline Vector<P_numtype>& Vector<P_numtype>::operator<<=(_bz_VecExpr<P_expr> expr)
{
    _bz_assign(expr, _bz_shiftl_update<P_numtype,
        _bz_typename P_expr::T_numtype>());
    return *this;
}
#endif   // NOT_DEFINED

/*****************************************************************************
 * Assignment operators with scalar operand
 */

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::initialize(P_numtype x)
{
    typedef _bz_VecExprConstant<P_numtype> T_expr;
    (*this) = _bz_VecExpr<T_expr>(T_expr(x));
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator+=(P_numtype x)
{
    typedef _bz_VecExprConstant<P_numtype> T_expr;
    (*this) += _bz_VecExpr<T_expr>(T_expr(x));
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator-=(P_numtype x)
{
    typedef _bz_VecExprConstant<P_numtype> T_expr;
    (*this) -= _bz_VecExpr<T_expr>(T_expr(x));
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator*=(P_numtype x)
{
    typedef _bz_VecExprConstant<P_numtype> T_expr;
    (*this) *= _bz_VecExpr<T_expr>(T_expr(x));
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator/=(P_numtype x)
{
    typedef _bz_VecExprConstant<P_numtype> T_expr;
    (*this) /= _bz_VecExpr<T_expr>(T_expr(x));
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator%=(P_numtype x)
{
    typedef _bz_VecExprConstant<P_numtype> T_expr;
    (*this) %= _bz_VecExpr<T_expr>(T_expr(x));
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator^=(P_numtype x)
{
    typedef _bz_VecExprConstant<P_numtype> T_expr;
    (*this) ^= _bz_VecExpr<T_expr>(T_expr(x));
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator&=(P_numtype x)
{
    typedef _bz_VecExprConstant<P_numtype> T_expr;
    (*this) &= _bz_VecExpr<T_expr>(T_expr(x));
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator|=(P_numtype x)
{
    typedef _bz_VecExprConstant<P_numtype> T_expr;
    (*this) |= _bz_VecExpr<T_expr>(T_expr(x));
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator>>=(int x)
{
    typedef _bz_VecExprConstant<int> T_expr;
    (*this) >>= _bz_VecExpr<T_expr>(T_expr(x));
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator<<=(int x)
{
    typedef _bz_VecExprConstant<P_numtype> T_expr;
    (*this) <<= _bz_VecExpr<T_expr>(T_expr(x));
    return *this;
}

/*****************************************************************************
 * Assignment operators with vector operand
 */

// This version is for two vectors with the same template parameter.
// Does not appear to be supported by the current C++ standard; or
// is it?
#if 0
template<typename P_numtype> template<>
inline Vector<P_numtype>&
Vector<P_numtype>::operator=(const Vector<P_numtype>& x)
{
    // NEEDS_WORK: if unit strides, use memcpy or something similar.

    typedef VectorIterConst<P_numtype> T_expr;
    (*this) = _bz_VecExpr<T_expr>(T_expr(*this));
    return *this;
}
#endif

// This version is for two vectors with *different* template
// parameters.
template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>&
Vector<P_numtype>::operator=(const Vector<P_numtype2>& x)
{
    (*this) = _bz_VecExpr<VectorIterConst<P_numtype2> >(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>&
Vector<P_numtype>::operator+=(const Vector<P_numtype2>& x)
{
    (*this) += _bz_VecExpr<VectorIterConst<P_numtype2> >(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>&
Vector<P_numtype>::operator-=(const Vector<P_numtype2>& x)
{
    (*this) -= _bz_VecExpr<VectorIterConst<P_numtype2> >(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>&
Vector<P_numtype>::operator*=(const Vector<P_numtype2>& x)
{
    (*this) *= _bz_VecExpr<VectorIterConst<P_numtype2> >(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>&
Vector<P_numtype>::operator/=(const Vector<P_numtype2>& x)
{
    (*this) /= _bz_VecExpr<VectorIterConst<P_numtype2> >(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>&
Vector<P_numtype>::operator%=(const Vector<P_numtype2>& x)
{
    (*this) %= _bz_VecExpr<VectorIterConst<P_numtype2> >(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>&
Vector<P_numtype>::operator^=(const Vector<P_numtype2>& x)
{
    (*this) ^= _bz_VecExpr<VectorIterConst<P_numtype2> >(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>&
Vector<P_numtype>::operator&=(const Vector<P_numtype2>& x)
{
    (*this) &= _bz_VecExpr<VectorIterConst<P_numtype2> >(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>&
Vector<P_numtype>::operator|=(const Vector<P_numtype2>& x)
{
    (*this) |= _bz_VecExpr<VectorIterConst<P_numtype2> >(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>&
Vector<P_numtype>::operator<<=(const Vector<P_numtype2>& x)
{
    (*this) <<= _bz_VecExpr<VectorIterConst<P_numtype2> >(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>&
Vector<P_numtype>::operator>>=(const Vector<P_numtype2>& x)
{
    (*this) >>= _bz_VecExpr<VectorIterConst<P_numtype2> >(x.beginFast());
    return *this;
}

/*****************************************************************************
 * Assignment operators with Range operand
 */

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator=(Range r)
{
    (*this) = _bz_VecExpr<Range>(r);
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator+=(Range r)
{
    (*this) += _bz_VecExpr<Range>(r);
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator-=(Range r)
{
    (*this) -= _bz_VecExpr<Range>(r);
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator*=(Range r)
{
    (*this) *= _bz_VecExpr<Range>(r);
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator/=(Range r)
{
    (*this) /= _bz_VecExpr<Range>(r);
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator%=(Range r)
{
    (*this) %= _bz_VecExpr<Range>(r);
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator^=(Range r)
{
    (*this) ^= _bz_VecExpr<Range>(r);
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator&=(Range r)
{
    (*this) &= _bz_VecExpr<Range>(r);
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator|=(Range r)
{
    (*this) |= _bz_VecExpr<Range>(r);
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator>>=(Range r)
{
    (*this) >>= _bz_VecExpr<Range>(r);
    return *this;
}

template<typename P_numtype>
inline Vector<P_numtype>& Vector<P_numtype>::operator<<=(Range r)
{
    (*this) <<= _bz_VecExpr<Range>(r);
    return *this;
}

/*****************************************************************************
 * Assignment operators with VectorPick operand
 */

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>& Vector<P_numtype>::operator=(const
    VectorPick<P_numtype2>& x)
{
    typedef VectorPickIterConst<P_numtype2> T_expr;
    (*this) = _bz_VecExpr<T_expr>(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>& Vector<P_numtype>::operator+=(const
    VectorPick<P_numtype2>& x)
{
    typedef VectorPickIterConst<P_numtype2> T_expr;
    (*this) += _bz_VecExpr<T_expr>(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>& Vector<P_numtype>::operator-=(const
    VectorPick<P_numtype2>& x)
{
    typedef VectorPickIterConst<P_numtype2> T_expr;
    (*this) -= _bz_VecExpr<T_expr>(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>& Vector<P_numtype>::operator*=(const
    VectorPick<P_numtype2>& x)
{
    typedef VectorPickIterConst<P_numtype2> T_expr;
    (*this) *= _bz_VecExpr<T_expr>(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>& Vector<P_numtype>::operator/=(const
    VectorPick<P_numtype2>& x)
{
    typedef VectorPickIterConst<P_numtype2> T_expr;
    (*this) /= _bz_VecExpr<T_expr>(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>& Vector<P_numtype>::operator%=(const
    VectorPick<P_numtype2>& x)
{
    typedef VectorPickIterConst<P_numtype2> T_expr;
    (*this) %= _bz_VecExpr<T_expr>(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>& Vector<P_numtype>::operator^=(const
    VectorPick<P_numtype2>& x)
{
    typedef VectorPickIterConst<P_numtype2> T_expr;
    (*this) ^= _bz_VecExpr<T_expr>(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>& Vector<P_numtype>::operator&=(const
    VectorPick<P_numtype2>& x)
{
    typedef VectorPickIterConst<P_numtype2> T_expr;
    (*this) &= _bz_VecExpr<T_expr>(x.beginFast());
    return *this;
}

template<typename P_numtype> template<typename P_numtype2>
inline Vector<P_numtype>& Vector<P_numtype>::operator|=(const
    VectorPick<P_numtype2>& x)
{
    typedef VectorPickIterConst<P_numtype2> T_expr;
    (*this) |= _bz_VecExpr<T_expr>(x.beginFast());
    return *this;
}

/*****************************************************************************
 * Assignment operators with Random operand
 */

template<typename P_numtype> template<typename P_distribution>
Vector<P_numtype>& Vector<P_numtype>::operator=(Random<P_distribution>& rand)
{
    for (int i=0; i < length_; ++i)
        (*this)[i] = rand.random();
    return *this;
}

#ifdef BZ_PECULIAR_RANDOM_VECTOR_ASSIGN_BUG

template<typename P_numtype> template<typename P_distribution>
Vector<P_numtype>& Vector<P_numtype>::operator=(Random<P_distribution>& rand)
{
    (*this) = _bz_VecExpr<_bz_VecExprRandom<P_distribution> >
        (_bz_VecExprRandom<P_distribution>(rand));
    return *this;
}

template<typename P_numtype> template<typename P_distribution>
Vector<P_numtype>& Vector<P_numtype>::operator+=(Random<P_distribution>& rand)
{
    (*this) += _bz_VecExpr<_bz_VecExprRandom<P_distribution> >
        (_bz_VecExprRandom<P_distribution>(rand));
    return *this;
}

template<typename P_numtype> template<typename P_distribution>
Vector<P_numtype>& Vector<P_numtype>::operator-=(Random<P_distribution>& rand)
{
    (*this) -= _bz_VecExpr<_bz_VecExprRandom<P_distribution> >
        (_bz_VecExprRandom<P_distribution>(rand));
    return *this;
}

template<typename P_numtype> template<typename P_distribution>
Vector<P_numtype>& Vector<P_numtype>::operator*=(Random<P_distribution>& rand)
{
    (*this) *= _bz_VecExpr<_bz_VecExprRandom<P_distribution> >
        (_bz_VecExprRandom<P_distribution>(rand));
    return *this;
}

template<typename P_numtype> template<typename P_distribution>
Vector<P_numtype>& Vector<P_numtype>::operator/=(Random<P_distribution>& rand)
{
    (*this) /= _bz_VecExpr<_bz_VecExprRandom<P_distribution> >
        (_bz_VecExprRandom<P_distribution>(rand));
    return *this;
}

template<typename P_numtype> template<typename P_distribution>
Vector<P_numtype>& Vector<P_numtype>::operator%=(Random<P_distribution>& rand)
{
    (*this) %= _bz_VecExpr<_bz_VecExprRandom<P_distribution> >
        (_bz_VecExprRandom<P_distribution>(rand));
    return *this;
}

template<typename P_numtype> template<typename P_distribution>
Vector<P_numtype>& Vector<P_numtype>::operator^=(Random<P_distribution>& rand)
{
    (*this) ^= _bz_VecExpr<_bz_VecExprRandom<P_distribution> >
        (_bz_VecExprRandom<P_distribution>(rand));
    return *this;
}

template<typename P_numtype> template<typename P_distribution>
Vector<P_numtype>& Vector<P_numtype>::operator&=(Random<P_distribution>& rand)
{
    (*this) &= _bz_VecExpr<_bz_VecExprRandom<P_distribution> >
        (_bz_VecExprRandom<P_distribution>(rand));
    return *this;
}

template<typename P_numtype> template<typename P_distribution>
Vector<P_numtype>& Vector<P_numtype>::operator|=(Random<P_distribution>& rand)
{
    (*this) |= _bz_VecExpr<_bz_VecExprRandom<P_distribution> >
        (_bz_VecExprRandom<P_distribution>(rand));
    return *this;
}

#endif // BZ_PECULIAR_RANDOM_VECTOR_ASSIGN_BUG

BZ_NAMESPACE_END

#endif // BZ_VECTOR_CC

ext: cc

max_stars_repo_name: HuXiaoling/TopoNet

max_issues_repo_name: HuXiaoling/TopoNet

max_forks_repo_name: HuXiaoling/TopoNet


text:
// Boost.Geometry

// Copyright (c) 2018, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle

// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_GEOMETRY_SRS_PROJECTIONS_PAR_DATA_HPP
#define BOOST_GEOMETRY_SRS_PROJECTIONS_PAR_DATA_HPP

#include <string>
#include <vector>

#include <boost/geometry/core/assert.hpp>
#include <boost/geometry/core/config.hpp>

namespace boost { namespace geometry { namespace srs
{

#ifndef DOXYGEN_NO_DETAIL
namespace detail
{

struct nadgrids
    : std::vector<std::string>
{
    typedef std::vector<std::string> base_t;

    nadgrids()
    {}

    template <typename It>
    nadgrids(It first, It last)
        : base_t(first, last)
    {}

#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
    nadgrids(std::initializer_list<std::string> l)
        : base_t(l)
    {}
#endif

    nadgrids(std::string const& g0)
        : base_t(1)
    {
        base_t& d = *this;
        d[0] = g0;
    }
    nadgrids(std::string const& g0, std::string const& g1)
        : base_t(2)
    {
        base_t& d = *this;
        d[0] = g0; d[1] = g1;
    }
    nadgrids(std::string const& g0, std::string const& g1, std::string const& g2)
        : base_t(3)
    {
        base_t& d = *this;
        d[0] = g0; d[1] = g1; d[2] = g2;
    }
    nadgrids(std::string const& g0, std::string const& g1, std::string const& g2, std::string const& g3)
        : base_t(4)
    {
        base_t& d = *this;
        d[0] = g0; d[1] = g1; d[2] = g2; d[3] = g3;
    }
    nadgrids(std::string const& g0, std::string const& g1, std::string const& g2, std::string const& g3, std::string const& g4)
        : base_t(5)
    {
        base_t& d = *this;
        d[0] = g0; d[1] = g1; d[2] = g2; d[3] = g3; d[4] = g4;
    }
};

template <typename T = double>
struct towgs84
{
    static const std::size_t static_capacity = 7;

    typedef std::size_t size_type;
    typedef T value_type;
    typedef T* iterator;
    typedef const T* const_iterator;
    typedef T& reference;
    typedef const T& const_reference;

    towgs84()
        : m_size(0)
#ifdef BOOST_GEOMETRY_CXX11_ARRAY_UNIFIED_INITIALIZATION
        , m_data{0, 0, 0, 0, 0, 0, 0}
    {}
#else
    {
        std::fill(m_data, m_data + 7, T(0));
    }
#endif

    template <typename It>
    towgs84(It first, It last)
    {
        assign(first, last);
    }

#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
    towgs84(std::initializer_list<T> l)
    {
        assign(l.begin(), l.end());
    }
#endif

    towgs84(T const& v0, T const& v1, T const& v2)
        : m_size(3)
    {
        m_data[0] = v0;
        m_data[1] = v1;
        m_data[2] = v2;
    }

    towgs84(T const& v0, T const& v1, T const& v2, T const& v3, T const& v4, T const& v5, T const& v6)
        : m_size(7)
    {
        m_data[0] = v0;
        m_data[1] = v1;
        m_data[2] = v2;
        m_data[3] = v3;
        m_data[4] = v4;
        m_data[5] = v5;
        m_data[6] = v6;
    }

    void push_back(T const& v)
    {
        BOOST_GEOMETRY_ASSERT(m_size < static_capacity);
        m_data[m_size] = v;
        ++m_size;
    }

    template <typename It>
    void assign(It first, It last)
    {
        for (m_size = 0 ; first != last && m_size < 7 ; ++first, ++m_size)
            m_data[m_size] = *first;
    }

#ifndef BOOST_NO_CXX11_HDR_INITIALIZER_LIST
    void assign(std::initializer_list<T> l)
    {
        assign(l.begin(), l.end());
    }
#endif

    const_reference operator[](size_type i) const
    {
        BOOST_GEOMETRY_ASSERT(i < m_size);
        return m_data[i];
    }

    reference operator[](size_type i)
    {
        BOOST_GEOMETRY_ASSERT(i < m_size);
        return m_data[i];
    }

    size_type size() const
    {
        return m_size;
    }

    bool empty() const
    {
        return m_size == 0;
    }

    void clear()
    {
        m_size = 0;
    }

    iterator begin() { return m_data; }
    iterator end() { return m_data + m_size; }
    const_iterator begin() const { return m_data; }
    const_iterator end() const { return m_data + m_size; }

private:
    size_type m_size;
    T m_data[7];
};

} // namespace detail
#endif // DOXYGEN_NO_DETAIL

}}} // namespace boost::geometry::srs

#endif // BOOST_GEOMETRY_SRS_PROJECTIONS_SPAR_HPP

ext: hpp


text:
#ifndef _INTRUSIVE_WEAK_PTR_HPP_INCLUDED
#define _INTRUSIVE_WEAK_PTR_HPP_INCLUDED

//
// Based on:
//
//  intrusive_ptr.hpp
//
//  Copyright (c) 2001, 2002 Peter Dimov
//
// Distributed under the Boost Software License, Version 1.0. (See
// accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
//
//  See http://www.boost.org/libs/smart_ptr/intrusive_ptr.html for documentation.
//

#include <boost/config.hpp>

#include <boost/assert.hpp>
#include <boost/detail/workaround.hpp>
#include <boost/smart_ptr/detail/sp_convertible.hpp>
#include <boost/smart_ptr/detail/sp_nullptr_t.hpp>

#include <boost/config/no_tr1/functional.hpp>           // for std::less

#include <boost/intrusive_ptr.hpp>

#if !defined(BOOST_NO_IOSTREAM)
#if !defined(BOOST_NO_IOSFWD)
#include <iosfwd>               // for std::basic_ostream
#else
#include <ostream>
#endif
#endif

//
//  intrusive_weak_ptr
//
//  A smart pointer that uses intrusive reference counting.
//
//  Relies on unqualified calls to
//
//      void intrusive_weak_ptr_add_ref(T * p);
//      void intrusive_weak_ptr_release(T * p);
//      bool intrusive_weak_ptr_alive(T * p);
//
//          (p != 0)
//
//  The object is responsible for destroying itself.
//

template<class T> class intrusive_weak_ptr
{
private:

	typedef intrusive_weak_ptr this_type;

public:

	typedef T element_type;

	intrusive_weak_ptr() BOOST_NOEXCEPT : px( 0 )
	{
	}

	intrusive_weak_ptr( T * p, bool add_ref = true ): px( p )
	{
		if( px != 0 && add_ref ) intrusive_weak_ptr_add_ref( px );
	}

#if !defined(BOOST_NO_MEMBER_TEMPLATES) || defined(BOOST_MSVC6_MEMBER_TEMPLATES)

	template<class U>
#if !defined( BOOST_SP_NO_SP_CONVERTIBLE )

	intrusive_weak_ptr( intrusive_weak_ptr<U> const & rhs, typename boost::detail::sp_enable_if_convertible<U,T>::type = boost::detail::sp_empty() )

#else

	intrusive_weak_ptr( intrusive_weak_ptr<U> const & rhs )

#endif
	: px( rhs.get() )
	{
		if( px != 0 ) intrusive_weak_ptr_add_ref( px );
	}

#endif

	intrusive_weak_ptr(intrusive_weak_ptr const & rhs): px( rhs.px )
	{
		if( px != 0 ) intrusive_weak_ptr_add_ref( px );
	}

	template<class U>
	intrusive_weak_ptr( boost::intrusive_ptr<U> const & rhs ) // non-explicit
	: px( rhs.get() )
	{
		if( px != 0 ) intrusive_weak_ptr_add_ref( px );
	}

	~intrusive_weak_ptr()
	{
		if( px != 0 ) intrusive_weak_ptr_release( px );
	}

#if !defined(BOOST_NO_MEMBER_TEMPLATES) || defined(BOOST_MSVC6_MEMBER_TEMPLATES)

	template<class U> intrusive_weak_ptr & operator=(intrusive_weak_ptr<U> const & rhs)
	{
		this_type(rhs).swap(*this);
		return *this;
	}

#endif

// Move support

#if !defined( BOOST_NO_CXX11_RVALUE_REFERENCES )

	intrusive_weak_ptr(intrusive_weak_ptr && rhs) BOOST_NOEXCEPT : px( rhs.px )
	{
		rhs.px = 0;
	}

	intrusive_weak_ptr & operator=(intrusive_weak_ptr && rhs) BOOST_NOEXCEPT
	{
		this_type( static_cast< intrusive_weak_ptr && >( rhs ) ).swap(*this);
		return *this;
	}

#endif

	intrusive_weak_ptr & operator=(intrusive_weak_ptr const & rhs)
	{
		this_type(rhs).swap(*this);
		return *this;
	}

	intrusive_weak_ptr & operator=(T * rhs)
	{
		this_type(rhs).swap(*this);
		return *this;
	}

	void reset() BOOST_NOEXCEPT
	{
		this_type().swap( *this );
	}

	void reset( T * rhs )
	{
		this_type( rhs ).swap( *this );
	}

	T * get() const BOOST_NOEXCEPT
	{
		if (px && !intrusive_weak_ptr_alive(px)) {
			intrusive_weak_ptr_release( px );
			px = 0;
		}
		return px;
	}

	/*~ Forbidden! Use get() and check for NULL.
	T & operator*() const
	{
		T * p = get();
		BOOST_ASSERT( p != 0 );
		return *p;
	}

	T * operator->() const
	{
		T * p = get();
		BOOST_ASSERT( p != 0 );
		return p;
	}
	*/

#if !defined( BOOST_NO_CXX11_EXPLICIT_CONVERSION_OPERATORS ) && !defined( BOOST_NO_CXX11_NULLPTR )

	explicit operator bool () const BOOST_NOEXCEPT
	{
		return get();
	}

#elif defined( _MANAGED )

	static void unspecified_bool( this_type*** )
	{
	}

	typedef void (*unspecified_bool_type)( this_type*** );

	operator unspecified_bool_type() const BOOST_NOEXCEPT
	{
		return get() ? unspecified_bool : 0;
	}

#else

	typedef element_type * this_type::*unspecified_bool_type;

	operator unspecified_bool_type() const BOOST_NOEXCEPT
	{
		return get() ? &this_type::px : 0;
	}

#endif

	// operator! is redundant, but some compilers need it
	bool operator! () const BOOST_NOEXCEPT
	{
		return !get();
	}

	void swap(intrusive_weak_ptr & rhs) BOOST_NOEXCEPT
	{
		T * tmp = px;
		px = rhs.px;
		rhs.px = tmp;
	}

private:

	mutable T * px;
};

template<class T, class U> inline bool operator==(intrusive_weak_ptr<T> const & a, intrusive_weak_ptr<U> const & b)
{
	return a.get() == b.get();
}

template<class T, class U> inline bool operator!=(intrusive_weak_ptr<T> const & a, intrusive_weak_ptr<U> const & b)
{
	return a.get() != b.get();
}

template<class T, class U> inline bool operator==(intrusive_weak_ptr<T> const & a, U * b)
{
	return a.get() == b;
}

template<class T, class U> inline bool operator!=(intrusive_weak_ptr<T> const & a, U * b)
{
	return a.get() != b;
}

template<class T, class U> inline bool operator==(T * a, intrusive_weak_ptr<U> const & b)
{
	return a == b.get();
}

template<class T, class U> inline bool operator!=(T * a, intrusive_weak_ptr<U> const & b)
{
	return a != b.get();
}

#if __GNUC__ == 2 && __GNUC_MINOR__ <= 96

// Resolve the ambiguity between our op!= and the one in rel_ops

template<class T> inline bool operator!=(intrusive_weak_ptr<T> const & a, intrusive_weak_ptr<T> const & b)
{
	return a.get() != b.get();
}

#endif

#if !defined( BOOST_NO_CXX11_NULLPTR )

template<class T> inline bool operator==( intrusive_weak_ptr<T> const & p, boost::detail::sp_nullptr_t ) BOOST_NOEXCEPT
{
	return p.get() == 0;
}

template<class T> inline bool operator==( boost::detail::sp_nullptr_t, intrusive_weak_ptr<T> const & p ) BOOST_NOEXCEPT
{
	return p.get() == 0;
}

template<class T> inline bool operator!=( intrusive_weak_ptr<T> const & p, boost::detail::sp_nullptr_t ) BOOST_NOEXCEPT
{
	return p.get() != 0;
}

template<class T> inline bool operator!=( boost::detail::sp_nullptr_t, intrusive_weak_ptr<T> const & p ) BOOST_NOEXCEPT
{
	return p.get() != 0;
}

#endif

template<class T> inline bool operator<(intrusive_weak_ptr<T> const & a, intrusive_weak_ptr<T> const & b)
{
	return std::less<T *>()(a.get(), b.get());
}

template<class T> void swap(intrusive_weak_ptr<T> & lhs, intrusive_weak_ptr<T> & rhs)
{
	lhs.swap(rhs);
}

// mem_fn support

template<class T> T * get_pointer(intrusive_weak_ptr<T> const & p)
{
	return p.get();
}

template<class T, class U> intrusive_weak_ptr<T> static_pointer_cast(intrusive_weak_ptr<U> const & p)
{
	return static_cast<T *>(p.get());
}

template<class T, class U> intrusive_weak_ptr<T> const_pointer_cast(intrusive_weak_ptr<U> const & p)
{
	return const_cast<T *>(p.get());
}

template<class T, class U> intrusive_weak_ptr<T> dynamic_pointer_cast(intrusive_weak_ptr<U> const & p)
{
	return dynamic_cast<T *>(p.get());
}

// operator<<

#if !defined(BOOST_NO_IOSTREAM)

#if defined(BOOST_NO_TEMPLATED_IOSTREAMS) || ( defined(__GNUC__) &&  (__GNUC__ < 3) )

template<class Y> std::ostream & operator<< (std::ostream & os, intrusive_weak_ptr<Y> const & p)
{
	os << p.get();
	return os;
}

#else

// in STLport's no-iostreams mode no iostream symbols can be used
#ifndef _STLP_NO_IOSTREAMS

# if defined(BOOST_MSVC) && BOOST_WORKAROUND(BOOST_MSVC, < 1300 && __SGI_STL_PORT)
// MSVC6 has problems finding std::basic_ostream through the using declaration in namespace _STL
using std::basic_ostream;
template<class E, class T, class Y> basic_ostream<E, T> & operator<< (basic_ostream<E, T> & os, intrusive_weak_ptr<Y> const & p)
# else
template<class E, class T, class Y> std::basic_ostream<E, T> & operator<< (std::basic_ostream<E, T> & os, intrusive_weak_ptr<Y> const & p)
# endif
{
	os << p.get();
	return os;
}

#endif // _STLP_NO_IOSTREAMS

#endif // __GNUC__ < 3

#endif // !defined(BOOST_NO_IOSTREAM)

// hash_value

template< class T > struct hash;

template< class T > std::size_t hash_value( intrusive_weak_ptr<T> const & p )
{
	return boost::hash< T* >()( p.get() );
}

#endif  // #ifndef _INTRUSIVE_WEAK_PTR_HPP_INCLUDED

max_stars_repo_name: Vbif/geometric-diversity

max_issues_repo_name: Vbif/geometric-diversity

max_forks_repo_name: Vbif/geometric-diversity


text:
// Copyright (c) Microsoft Corporation. All rights reserved.
// Licensed under the MIT license.

#include "conv.h"
#include "common.h"
#include "gemm.h"
#include "utils/shape_utils.h"
#include <Eigen/Core>
#include <Eigen/Dense>
#include <vector>

// Wasm interop method
void conv_f32(void *data) {
  uint32_t *dataIndex = static_cast<uint32_t *>(data);
  uint32_t const argc = dataIndex[0];
  conv2D_f32_imp(
      PARAM_FLOAT_PTR(data, dataIndex[1]), PARAM_INT32_PTR(data, dataIndex[2]),
      PARAM_FLOAT_PTR(data, dataIndex[3]), PARAM_INT32_PTR(data, dataIndex[4]),
      PARAM_FLOAT_PTR(data, dataIndex[5]), PARAM_INT32_PTR(data, dataIndex[6]),
      PARAM_FLOAT_PTR(data, dataIndex[7]), PARAM_INT32_PTR(data, dataIndex[8]),
      PARAM_INT32(data, dataIndex[9]), PARAM_INT32_PTR(data, dataIndex[10]),
      PARAM_INT32_PTR(data, dataIndex[11]), PARAM_INT32(data, dataIndex[12]));
}

// TODO: There's lots of room for optimization here
// Core operator implementation
void conv2D_f32_imp(float *X, int *X_shape, float *W, int *W_shape, float *Y,
                    int *Y_shape, float *bias, int *dilations, int group,
                    int *pads, int *strides, int dimension) {
  const int N = X_shape[0]; // Batch size
  const int C = X_shape[1]; // Number of input channels

  const int M = W_shape[0]; // Output channels

  const int CG = C / group;

  std::vector<int32_t> inputSizes(X_shape,
                                  X_shape + dimension +
                                      2); // Input dimension sizes
  std::vector<int32_t> kernelSizes(W_shape,
                                   W_shape + dimension +
                                       2); // Kernel dimension sizes
  std::vector<int32_t> outputSizes(Y_shape,
                                   Y_shape + dimension +
                                       2); // Output dimension sizes

  const size_t kernelSize = ShapeUtils::size_from_dims(kernelSizes);
  const size_t inputSize = ShapeUtils::size_from_dims(inputSizes);
  const size_t outputSize = ShapeUtils::size_from_dims(outputSizes);

  const size_t outputDataSize = ShapeUtils::size_from_dims(outputSizes, 2);
  const size_t kernelDataSize = ShapeUtils::size_from_dims(kernelSizes, 2);

  std::vector<int32_t> outputStrides = ShapeUtils::compute_strides(outputSizes);
  std::vector<int32_t> inputStrides = ShapeUtils::compute_strides(inputSizes);
  std::vector<int32_t> kernelStrides = ShapeUtils::compute_strides(kernelSizes);

  std::vector<int32_t> outputIndices(dimension + 2, 0);
  size_t outputIndex = 0;

  std::vector<int32_t> inputIndices(dimension + 2, 0);
  size_t inputIndex = 0;

  std::vector<int32_t> kernelIndices(dimension + 2, 0);
  size_t kernelIndex = 0;

  // Iterate over all batches
  for (int32_t n = 0; n < N; n++) {
    outputIndices[0] = n;
    inputIndices[0] = n;
    // Iterate over all output channels
    for (int32_t m = 0; m < M; m++) {
      outputIndices[1] = m;
      kernelIndices[0] = m;
      if (bias != nullptr) {
        const float b = bias[m];

        size_t outIx = n * outputStrides[0] + m * outputStrides[1];
        const size_t end = outIx + outputDataSize;

        while (outIx < end) {
          Y[outIx] = b;
          outIx++;
        }
      }

      for (int32_t cg = 0; cg < CG; cg++) {
        const int32_t c = m * CG + cg;
        kernelIndices[1] = cg;
        inputIndices[1] = c;

        std::fill(outputIndices.begin() + 2, outputIndices.end(), 0);

        size_t outputBase = outputStrides[0] * n + outputStrides[1] * m;
        outputIndex = 0;
        while (outputIndex < outputDataSize) {

          std::fill(kernelIndices.begin() + 2, kernelIndices.end(), 0);
          kernelIndex = 0;
          size_t kernelBase = kernelStrides[0] * m + kernelStrides[1] * cg;
          while (kernelIndex < kernelDataSize) {
            bool skip = false;

            inputIndex = n * inputStrides[0] + c * inputStrides[1];
            for (size_t axis = 0; axis < dimension; axis++) {
              const int32_t ix = outputIndices[axis + 2] * strides[axis] -
                                 pads[axis] +
                                 kernelIndices[axis + 2] * dilations[axis];

              if (ix < 0 || ix >= inputSizes[axis + 2]) {
                skip = true;
                break;
              }

              inputIndices[axis + 2] = ix;
              inputIndex += inputStrides[axis + 2] * ix;
            }

            if (!skip) {
              Y[outputBase + outputIndex] +=
                  W[kernelBase + kernelIndex] * X[inputIndex];
            }

            ShapeUtils::increment_index(kernelIndices, kernelSizes,
                                        dimension + 2);
            kernelIndex++;
          }

          outputIndex++;
          ShapeUtils::increment_index(outputIndices, outputSizes,
                                      dimension + 2);
        }
      }
    }
  }
}


text:
// Boost.TypeErasure library
//
// Copyright 2012 Steven Watanabe
//
// Distributed under the Boost Software License Version 1.0. (See
// accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
//
// $Id$

#if !defined(BOOST_PP_IS_ITERATING)

#ifndef BOOST_TYPE_ERASURE_DETAIL_CHECK_CALL_HPP_INCLUDED
#define BOOST_TYPE_ERASURE_DETAIL_CHECK_CALL_HPP_INCLUDED

#include <boost/mpl/eval_if.hpp>
#include <boost/mpl/bool.hpp>
#include <boost/mpl/and.hpp>
#include <boost/type_traits/is_convertible.hpp>
#include <boost/type_traits/remove_reference.hpp>
#include <boost/type_traits/remove_cv.hpp>
#include <boost/type_traits/is_reference.hpp>
#include <boost/type_traits/is_const.hpp>
#include <boost/type_traits/function_traits.hpp>
#include <boost/preprocessor/iteration/iterate.hpp>
#include <boost/preprocessor/repetition/enum_params.hpp>
#include <boost/preprocessor/repetition/enum_trailing_params.hpp>
#include <boost/preprocessor/repetition/repeat.hpp>
#include <boost/preprocessor/cat.hpp>
#include <boost/type_erasure/placeholder_of.hpp>

namespace boost {
namespace type_erasure {
namespace detail {

template<class Sig, class Args>
struct check_call : ::boost::mpl::false_ {};

template<class T, class Enable = void>
struct qualified_placeholder
{
    typedef void type;
};

template<class T>
struct qualified_placeholder<T&, typename T::_boost_type_erasure_is_any>
{
    typedef typename ::boost::type_erasure::placeholder_of<T>::type placeholder;
    typedef typename ::boost::remove_reference<placeholder>::type unref;
    typedef typename ::boost::mpl::if_< ::boost::is_const<T>,
        const unref,
        unref
    >::type add_const;
    typedef typename ::boost::mpl::if_< ::boost::is_reference<placeholder>,
        placeholder,
        add_const&
    >::type type;
};

#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES

template<class T>
struct qualified_placeholder<T&&, typename T::_boost_type_erasure_is_any>
{
    typedef typename ::boost::type_erasure::placeholder_of<T>::type placeholder;
    typedef placeholder&& type;
};

#endif

template<class P, class A>
struct check_placeholder_arg_impl : ::boost::mpl::false_ {};

template<class P>
struct check_placeholder_arg_impl<P, P&> : ::boost::mpl::true_ {};

template<class P>
struct check_placeholder_arg_impl<P, const P&> : ::boost::mpl::true_ {};

#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES

template<class P>
struct check_placeholder_arg_impl<P, P&&> : ::boost::mpl::true_ {};

#endif

template<class P>
struct check_placeholder_arg_impl<P&, P&> : ::boost::mpl::true_ {};

template<class P>
struct check_placeholder_arg_impl<const P&, P&> : ::boost::mpl::true_ {};

template<class P>
struct check_placeholder_arg_impl<const P&, const P&> : ::boost::mpl::true_ {};

#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES

template<class P>
struct check_placeholder_arg_impl<const P&, P&&> : ::boost::mpl::true_ {};

template<class P>
struct check_placeholder_arg_impl<P&&, P&&> : ::boost::mpl::true_ {};

#endif

template<class P, class Arg>
struct check_placeholder_arg :
    check_placeholder_arg_impl<
        P,
        typename ::boost::type_erasure::detail::qualified_placeholder<Arg>::type
    >::type
{};

#if !defined(BOOST_NO_CXX11_RVALUE_REFERENCES) && \
    ((defined(__GNUC__) && !(__GNUC__ > 4 || (__GNUC__ == 4 && __GNUC_MINOR__ > 6))) || \
    defined(__MINGW32__) || defined(__MINGW64__))
#define BOOST_TYPE_ERASURE_BROKEN_RVALUE_IS_CONVERTIBLE
#endif

template<class FormalArg, class ActualArg>
struct check_arg
{
    typedef typename ::boost::mpl::eval_if<
        is_placeholder<
            typename ::boost::remove_cv<
                typename ::boost::remove_reference<FormalArg>::type
            >::type
        >,
        ::boost::type_erasure::detail::check_placeholder_arg<FormalArg, ActualArg>,
#ifdef BOOST_TYPE_ERASURE_BROKEN_RVALUE_IS_CONVERTIBLE
        ::boost::mpl::true_
#else
        ::boost::is_convertible<ActualArg, FormalArg>
#endif
    >::type type;
};

#define BOOST_PP_FILENAME_1 <boost/type_erasure/detail/check_call.hpp>
#define BOOST_PP_ITERATION_LIMITS (0, BOOST_TYPE_ERASURE_MAX_ARITY)
#include BOOST_PP_ITERATE()

}
}
}

#endif

#else

#define N BOOST_PP_ITERATION()

#define BOOST_TYPE_ERASURE_CHECK_ARG(z, n, data)                \
    typedef typename ::boost::type_erasure::detail::check_arg<  \
        BOOST_PP_CAT(T, n),                                     \
        BOOST_PP_CAT(U, n)                                      \
    >::type BOOST_PP_CAT(check, n);                             \
    typedef typename ::boost::mpl::and_<                        \
        BOOST_PP_CAT(type, n),                                  \
        BOOST_PP_CAT(check, n)                                  \
    >::type BOOST_PP_CAT(type, BOOST_PP_INC(n));

template<
    class R
    BOOST_PP_ENUM_TRAILING_PARAMS(N, class T)
    BOOST_PP_ENUM_TRAILING_PARAMS(N, class U)
>
struct check_call<R(BOOST_PP_ENUM_PARAMS(N, T)), void(BOOST_PP_ENUM_BINARY_PARAMS(N, U, u))> {
    typedef ::boost::mpl::true_ type0;
    BOOST_PP_REPEAT(N, BOOST_TYPE_ERASURE_CHECK_ARG, ~)
    typedef BOOST_PP_CAT(type, N) type;
};

#undef N

#endif

max_stars_repo_name: rudylee/expo

max_issues_repo_name: rudylee/expo

max_forks_repo_name: rudylee/expo


text:
//
//  Copyright (c) 2010 Athanasios Iliopoulos
//
//  Distributed under the Boost Software License, Version 1.0. (See
//  accompanying file LICENSE_1_0.txt or copy at
//  http://www.boost.org/LICENSE_1_0.txt)
//

#include <boost/numeric/ublas/assignment.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <boost/timer.hpp>

using namespace boost::numeric::ublas;

int main() {

    boost::timer timer;

    unsigned int iterations = 1000000000;
    double elapsed_exp, elapsed_assigner;

    std::cout << "Ublas vector<double> Benchmarks------------------------ " << "\n";

    {
    std::cout << "Size 2 vector: " << "\n";
    vector<double> a(2);

    timer.restart();
    for(unsigned int i=0; i!=iterations; i++) {
        a(0)=0; a(1)=1;
    }
    elapsed_exp = timer.elapsed();
    std::cout << "Explicit element assign time: " << elapsed_exp << " secs" << "\n";

    timer.restart();
    for(unsigned int i=0; i!=iterations; i++)
        a <<= 0, 1;
    elapsed_assigner = timer.elapsed();
    std::cout << "Assigner time: " << elapsed_assigner << " secs" << "\n";
    std::cout << "Difference: " << (elapsed_assigner/elapsed_exp-1)*100 << "%" << std::endl;
    }

    {
    std::cout << "Size 3 vector: " << "\n";
    vector<double> a(3);

    timer.restart();
    for(unsigned int i=0; i!=iterations; i++) {
        a(0)=0; a(1)=1; a(2)=2;
    }
    elapsed_exp = timer.elapsed();
    std::cout << "Explicit element assign time: " << elapsed_exp << " secs" << "\n";

    timer.restart();
    for(unsigned int i=0; i!=iterations; i++)
        a <<= 0, 1, 2;
    elapsed_assigner = timer.elapsed();
    std::cout << "Assigner time: " << elapsed_assigner << " secs" << "\n";
    std::cout << "Difference: " << (elapsed_assigner/elapsed_exp-1)*100 << "%" << std::endl;
    }

    iterations = 100000000;

    {
    std::cout << "Size 8 vector: " << "\n";
    vector<double> a(8);

    timer.restart();
    for(unsigned int i=0; i!=iterations; i++) {
        a(0)=0; a(1)=1; a(2)=2; a(3)=3; a(4)=4; a(5)=5; a(6)=6; a(7)=7;
    }
    elapsed_exp = timer.elapsed();
    std::cout << "Explicit element assign time: " << elapsed_exp << " secs" << "\n";

    timer.restart();
    for(unsigned int i=0; i!=iterations; i++)
        a <<= 0, 1, 2, 3, 4, 5, 6, 7;
    elapsed_assigner = timer.elapsed();
    std::cout << "Assigner time: " << elapsed_assigner << " secs" << "\n";
    std::cout << "Difference: " << (elapsed_assigner/elapsed_exp-1)*100 << "%" << std::endl;
    }

    std::cout << "Ublas matrix<double> Benchmarks------------------------ " << "\n";

    iterations = 200000000;
    {
    std::cout << "Size 3x3 matrix: " << "\n";
    matrix<double> a(3,3);

    timer.restart();
    for(unsigned int i=0; i!=iterations; i++) {
        a(0,0)=0; a(0,1)=1; a(0,2)=2;
        a(1,0)=3; a(1,1)=4; a(1,2)=5;
        a(2,0)=6; a(2,1)=7; a(2,2)=8;
    }
    elapsed_exp = timer.elapsed();
    std::cout << "Explicit element assign time: " << elapsed_exp << " secs" << "\n";

    timer.restart();
    for(unsigned int i=0; i!=iterations; i++)
        a <<= 0, 1, 2, 3, 4, 5, 6, 7, 8;
    elapsed_assigner = timer.elapsed();
    std::cout << "Assigner time: " << elapsed_assigner << " secs" << "\n";
    std::cout << "Difference: " << (elapsed_assigner/elapsed_exp-1)*100 << "%" << std::endl;
    }

    std::cout << "Size 2x2 matrix: " << "\n";
    iterations = 500000000;
    {
    matrix<double> a(2,2);

    timer.restart();
    for(unsigned int i=0; i!=iterations; i++) {
        a(0,0)=0; a(0,1)=1;
        a(1,0)=3; a(1,1)=4;
    }
    elapsed_exp = timer.elapsed();
    std::cout << "Explicit element assign time: " << elapsed_exp << " secs" << "\n";

    timer.restart();
    for(unsigned int i=0; i!=iterations; i++)
        a <<= 0, 1, 3, 4;
    elapsed_assigner = timer.elapsed();
    std::cout << "Assigner time: " << elapsed_assigner << " secs" << "\n";

    std::cout << "Difference: " << (elapsed_assigner/elapsed_exp-1)*100 << "%" << std::endl;

    timer.restart();
    for(unsigned int i=0; i!=iterations; i++)
        a <<= traverse_policy::by_row_no_wrap(), 0, 1, next_row(), 3, 4;
    elapsed_assigner = timer.elapsed();
    std::cout << "Assigner time no_wrap: " << elapsed_assigner << " secs" << "\n";
    std::cout << "Difference: " << (elapsed_assigner/elapsed_exp-1)*100 << "%" << std::endl;
    }

    return 0;
}