C++ valarray::size方法代码示例

void save( Archive & ar, const STD::valarray<U> &t, const unsigned int /*file_version*/ )
{
    const collection_size_type count(t.size());
    ar << BOOST_SERIALIZATION_NVP(count);
    if (t.size()){
        // explict template arguments to pass intel C++ compiler
        ar << serialization::make_array<const U, collection_size_type>(
            static_cast<const U *>( boost::addressof(t[0]) ),
            count
        );
    }
}

template <class T> std::valarray<T> mmpbsa_utils::zip(const std::vector<T>& left,
            const std::valarray<T>& right)
{
    if(left.size() != right.size())
        throw mmpbsa::MMPBSAException("When using zip, the two arrays must have "
                "the same size.",mmpbsa::DATA_FORMAT_ERROR);

    size_t oldsize = left.size();
    std::valarray<T> returnMe(2*oldsize);
    size_t returnMeIndex = 0;
    for(size_t i = 0;i<oldsize;i++)
    {
        returnMe[returnMeIndex++] = left[i];
        returnMe[returnMeIndex++] = right[i];
    }
    return returnMe;
}

      template<class T> std::set<T> unique(const std::valarray<T> & rhs)
      {
         using namespace std;
         int i;
         
         valarray<T> newArray;

         set<T> newSet;
                  
         for (i=0; i<rhs.size(); i++)
         {
               //cout << rhs[i] << endl;
            newSet.insert(rhs[i]);
         }

         return newSet;
      }

    // method to retrieve the parafoveal magnocellular pathway response (no energy motion in fovea)
    const bool RetinaFilter::getMagnoParaFoveaResponse(std::valarray<float> &magnoParafovealResponse)
    {
        if (!_useMagnoOutput)
            return false;
        if (magnoParafovealResponse.size() != _MagnoRetinaFilter.getNBpixels())
            return false;

        register const float *magnoXOutputPTR= get_data(_MagnoRetinaFilter.getOutput());
        register float *parafovealMagnoResponsePTR=&magnoParafovealResponse[0];
        register float *hybridParvoMagnoCoefTablePTR=&_retinaParvoMagnoMapCoefTable[0]+1;

        for (unsigned int i=0 ; i<_photoreceptorsPrefilter.getNBpixels() ; ++i, hybridParvoMagnoCoefTablePTR+=2)
        {
            *(parafovealMagnoResponsePTR++)=*(magnoXOutputPTR++)**(hybridParvoMagnoCoefTablePTR);
        }

        return true;
    }

double ColaTopologyAddon::applyForcesAndConstraints(
        cola::ConstrainedFDLayout *layout, const vpsc::Dim dim,
        std::valarray<double>& g, vpsc::Variables& vs,
        vpsc::Constraints& cs, std::valarray<double> &coords, 
        cola::DesiredPositionsInDim& des, double oldStress)
{
    FILE_LOG(cola::logDEBUG1) << "applying topology preserving layout...";
    vpsc::Rectangle::setXBorder(0);
    vpsc::Rectangle::setYBorder(0);
    if(dim==vpsc::HORIZONTAL) {
        vpsc::Rectangle::setXBorder(0);
    }
    topology::setNodeVariables(topologyNodes,vs);
    topology::TopologyConstraints t(dim, topologyNodes, topologyRoutes,
            layout->clusterHierarchy, vs, cs);
    bool interrupted;
    int loopBreaker=100;
    cola::SparseMap HMap(layout->n);
    layout->computeForces(dim,HMap,g);
    std::valarray<double> oldCoords=coords;
    t.computeForces(g,HMap);
    cola::SparseMatrix H(HMap);
    layout->applyDescentVector(g,oldCoords,coords,oldStress,
            layout->computeStepSize(H,g,g));
    cola::setVariableDesiredPositions(vs,cs,des,coords);
    do {
        interrupted=t.solve();
        unsigned vptr=0;
        for(topology::Nodes::iterator i=topologyNodes.begin();
                i!=topologyNodes.end();++i,++vptr) {
            topology::Node* v=*i;
            coords[v->id]=v->rect->getCentreD(dim);
        }
        for(;vptr<coords.size();vptr++) {
            double d = vs[vptr]->finalPosition;
            coords[vptr]=d;
            layout->boundingBoxes[vptr]->moveCentreD(dim,d);
        }
        loopBreaker--;
    } while(interrupted&&loopBreaker>0);
    vpsc::Rectangle::setXBorder(0);
    vpsc::Rectangle::setYBorder(0);
    return layout->computeStress();
}

TargetInfo::TargetInfo(const std::valarray<double>& t) {
    _weights.resize(0);
    _targets.resize(t.size());
    _targets = t;
    
    _tmean = _targets.sum()/_targets.size();
    
    _tcov_part.resize(_targets.size());
    _tcov_part = _targets;
    _tcov_part -= _tmean;
	
    std::valarray<double> tmp = _tcov_part;
    tmp = _tcov_part;
    tmp *= tmp;
	
    _tvar = tmp.sum() / (tmp.size()-1);
    _tstd = sqrt(_tvar);
    _tmed = 0;
}

void stats(std:: valarray<float> map, float &mean, float &sigma, float &kurt, float &skew) {
    float sum = map.sum();
    int n = map.size();
    mean = map.sum()/float(n);
    std:: valarray <float> maps(n);
    valarray<float> maps2(n),maps3(n),maps4(n);
    for(int i=0; i<n; i++) {
        maps2[i] = gsl_pow_2(map[i] - mean);
        maps3[i] = gsl_pow_3(map[i] - mean);
        maps4[i] = gsl_pow_4(map[i] - mean);
    }
    sum = maps2.sum();
    sigma = sqrt(sum/(float(n)-1.));
    sum = maps3.sum();
    double mu3 = sum/(float(n)-1.);
    sum = maps4.sum();
    double mu4 = sum/(float(n)-1.);
    kurt = mu4/gsl_pow_4(sigma) -3;
    skew = mu3/gsl_pow_3(sigma);
}

static void
test_gslice (std::size_t  start,
             const char  *sizes,
             const char  *strides)
{
    const std::valarray<std::size_t> asizes (make_array (sizes));
    const std::valarray<std::size_t> astrides (make_array (strides));

    const std::gslice gsl (start, asizes, astrides);

    const std::valarray<std::size_t> indices = get_index_array (gsl);

    std::size_t maxinx = 0;

    for (std::size_t i = 0; i != indices.size (); ++i)
        if (maxinx < indices [i])
            maxinx = indices [i];

    std::valarray<std::size_t> va (maxinx + 1);
    for (std::size_t i = 0; i != va.size (); ++i)
        va [i] = i;

    for (int i = 0; i != 3; ++i) {
        // repeat each test three to verify that operator[](gslice)
        // doesn't change the observable state of its argument and
        // that the same result is obtained for a copy of gslice

        const std::valarray<std::size_t> array_slice =
            i < 2 ? va [gsl] : va [std::gslice (gsl)];

        bool equal = array_slice.size () == indices.size ();

        rw_assert (equal, 0, __LINE__,
                   "size() == %zu, got %zu\n",
                   indices.size (), array_slice.size ());

        if (equal) {
            for (std::size_t j = 0; j != array_slice.size (); ++j) {

                equal = array_slice [j] == va [indices [j]];

                rw_assert (equal, 0, __LINE__,
                           "mismatch at %u, index %u: expected %u, got %u\n",
                           j, indices [j], va [indices [j]],
                           array_slice [j]);
            }
        }
    }
}

// calculate the members, now in the context of a mask
void TargetInfo::set_training_mask(const std::valarray<bool>& tmask) {
    
    TargetInfo tmp;
    
    if (has_weights() ) {
	tmp = TargetInfo( _targets[tmask], _weights[tmask]);
    } else {
	tmp = TargetInfo( _targets[tmask] );
    }
    
    _tcov_part.resize(tmp._tcov_part.size());
    _tcov_part = tmp._tcov_part;

    _tmean = tmp._tmean;
    _tvar  = tmp._tvar;
    _tstd  = tmp._tstd;
    _tmed =  tmp._tmed;

    _training_mask.resize(tmask.size());
    _training_mask = tmask;
}

std::valarray<double> SourceTermsWindRMHD::dUdt(const std::valarray<double> &Uin)
{
  this->prepare_integration();

  std::valarray<double> Lsrc(Uin.size());

  double C = 2.0; // compactness of source region (larger value is more compact)
  double L = 0.125 / C; // cut-off radius, C=1 gives a source size of 1/8

  for (int i=0; i<stride[0]; i+=NQ) {

    int N[3];
    absolute_index_to_3d(i/NQ, N);

    double x[3] = {
      Mara->domain->x_at(N[0]),
      Mara->domain->y_at(N[1]),
      Mara->domain->z_at(N[2]) };

    double L0[8] = {0, 0, 0, 0, 0, 0, 0, 0};
    double R = sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]);
    double r = sqrt(x[0]*x[0] + x[1]*x[1]);

    double phi_hat[2] = { -x[1]/r, x[0]/r };

    double ramp = C * (1.0 - R/L > 0.0 ? 1.0 - R/L : 0.0);
    //double ramp = C * (R < L ? exp(L*L / (L*L - r*r)) : 0.0);
    L0[ddd] = ddot * ramp;
    L0[tau] = edot * ramp;
    L0[Sx] = sdot * phi_hat[0] * ramp;
    L0[Sy] = sdot * phi_hat[1] * ramp;

    for (int q=0; q<NQ; ++q) {
      Lsrc[i+q] += L0[q];
    }
  }

  return Lsrc;
}

pdf_redshift::pdf_redshift(std::valarray< double >& zSampling, std::valarray< double >& pdz)
{
  
  // Allocate arrays to store pdf samples
  nPoints = zSampling.size();
  x = (double*) malloc(sizeof(double) * nPoints);
  y = (double*) malloc(sizeof(double) * nPoints);
  
  // Copy array data
  for(int i =0; i < nPoints; i++){
      x[i] = zSampling[i];
      y[i] = pdz[i];
  }
  
  interpolator = gsl_interp_alloc(gsl_interp_cspline, nPoints);
  accelerator = gsl_interp_accel_alloc();
  
  gsl_interp_init(interpolator, x, y, nPoints);
  
  // Compute normalization factor
  normalizationFactor = 1.0 / gsl_interp_eval_integ(interpolator, x, y, x[0], x[nPoints-1], accelerator);
}

void TransientAreasSegmentationModuleImpl::_run(const std::valarray<float> &inputToSegment, const int channelIndex)
{
#ifdef SEGMENTATIONDEBUG
    std::cout<<"Input length vs internal buffers length = "<<inputToSegment.size()<<", "<<_localMotion.size()<<std::endl;
#endif
    // preliminary basic error check
    // FIXME validate basic tests
    //if (inputToSegment.size() != _localMotion.size())
    //    throw cv::Exception(-1, "Input matrix size does not match instance buffers setup !", "SegmentationModule::run", "SegmentationModule.cpp", 0);

    // first square the input in order to increase the signal to noise ratio
    // get motion local energy
    _squaringSpatiotemporalLPfilter(&inputToSegment[channelIndex*getNBpixels()], &_localMotion[0]);

    // second low pass filter: access to the neighborhood motion energy
    _spatiotemporalLPfilter(&_localMotion[0], &_neighborhoodMotion[0], 1);

    // third low pass filter: access to the background motion energy
    _spatiotemporalLPfilter(&_localMotion[0], &_contextMotionEnergy[0], 2);

    // compute the ON and OFF ways (positive and negative values of the difference of the two filterings)
    float*localMotionPTR=&_localMotion[0], *neighborhoodMotionPTR=&_neighborhoodMotion[0], *contextMotionPTR=&_contextMotionEnergy[0];

    // float meanEnergy=LPfilter2.sum()/(float)_LPfilter2.size();
    register bool *segmentationPicturePTR= &_segmentedAreas[0];
    for (unsigned int index=0; index<_filterOutput.getNBpixels() ; ++index, ++segmentationPicturePTR, ++localMotionPTR, ++neighborhoodMotionPTR, contextMotionPTR++)
    {
        float generalMotionContextDecision=*neighborhoodMotionPTR-*contextMotionPTR;

        if (generalMotionContextDecision>0) // local maximum should be detected in this case
        {
            /* apply segmentation on local motion superior to its neighborhood
             * => to segment objects moving faster than their neighborhood
             */
            if (generalMotionContextDecision>_segmentationParameters.thresholdON)// && meanEnergy*1.1<*neighborhoodMotionPTR)
            {
                *segmentationPicturePTR=((*localMotionPTR-*neighborhoodMotionPTR)>_segmentationParameters.thresholdON);
            }
            else
                *segmentationPicturePTR=false;
        }
#ifdef USE_LOCALMINIMUMS
        else  // local minimum should be detected in this case
        {
            /* apply segmentation for non moving objects
             * only if the wide area around motion energy is high
             * => to segment object moving slower than the neighborhood
             */
            if (-1.0*generalMotionContextDecision>_segmentationParameters.thresholdOFF && meanEnergy*0.9>*neighborhoodMotionPTR)
            {
                /* in order to segment non moving objects in a camera motion case
                 * we focus on local energy which is much lower than the wide neighborhood
                 */
                *segmentationPicturePTR+=(*neighborhoodMotionPTR-*localMotionPTR>_segmentationParameters.thresholdOFF)*127;
            }
        }
#else
        else
            *segmentationPicturePTR=false;
#endif
    }

void CreateChildProcess(const std::vector<mmpbsa::atom_t>& atoms,
						  const std::valarray<mmpbsa::Vector>& crds,
						  const std::map<std::string,mead_data_t>& radii,
						  int *error_flag)
// Create a child process that uses the previously created pipes for STDIN and STDOUT.
{ 
  const size_t max_chars = 4096;
  CHAR cmdstr[max_chars], appname[max_chars];
  PROCESS_INFORMATION piProcInfo; 
  STARTUPINFO siStartInfo;
  BOOL bSuccess = FALSE;

  if(mmpbsa_io::resolve_filename("./xyzr2sas.exe",appname,max_chars) != 0)
    ErrorExit(TEXT("Could not resolve the area application filename"));
  
  snprintf(cmdstr,max_chars,"\"%s\" %lu",appname,crds.size());			      
 
// Set up members of the PROCESS_INFORMATION structure. 
 
   ZeroMemory( &piProcInfo, sizeof(PROCESS_INFORMATION) );
 
// Set up members of the STARTUPINFO structure. 
// This structure specifies the STDIN and STDOUT handles for redirection.
 
   ZeroMemory( &siStartInfo, sizeof(STARTUPINFO) );
   siStartInfo.cb = sizeof(STARTUPINFO); 
   //siStartInfo.hStdError = g_hChildStd_OUT_Wr;
   siStartInfo.hStdOutput = g_hChildStd_OUT_Wr;
   siStartInfo.hStdInput = g_hChildStd_IN_Rd;
   siStartInfo.dwFlags |= STARTF_USESTDHANDLES;
 
// Create the child process. 
   UINT WINAPI orig_err_mask = SetErrorMode(SEM_NOGPFAULTERRORBOX);
   SetErrorMode(orig_err_mask | SEM_NOGPFAULTERRORBOX);

   bSuccess = CreateProcess(appname, 
			    cmdstr,     // command line
      NULL,          // process security attributes 
      NULL,          // primary thread security attributes 
      TRUE,          // handles are inherited 
			    CREATE_NO_WINDOW,			    //original: 0,             // creation flags 
      NULL,          // use parent's environment 
      NULL,          // use parent's current directory 
      &siStartInfo,  // STARTUPINFO pointer 
      &piProcInfo);  // receives PROCESS_INFORMATION 
   
   SetErrorMode(orig_err_mask);

   // If an error occurs, exit the application. 
   if ( ! bSuccess ) 
     {
       if(error_flag != NULL)
	 *error_flag  = 1;
       ErrorExit(TEXT("Could not start child process"));
     }
   else 
   {
      // Close handles to the child process and its primary thread.
      // Some applications might keep these handles to monitor the status
      // of the child process, for example. 
      CloseHandle(piProcInfo.hProcess);
      CloseHandle(piProcInfo.hThread);
   }
}

 valarray(const std::valarray<T> &valarray,
          const context &context = system::default_context())
     : m_buffer(context, valarray.size() * sizeof(T))
 {
     copy(&valarray[0], &valarray[valarray.size()], begin());
 }

 void get_gradient(std::valarray<double>& rv) const {rv.resize(vg.size()); rv=vg; }