C++ Compute函数代码示例

/*
 * Move a selected region up or down, reusing already computed values.
 */
static void MoveUD(int offset) {
  int width = (int) canvas_width;

  if (offset > 0) {
    char* p = pixmap;
    char* q = pixmap + offset * width;
    size_t size = (width - offset) * width;

    /* Move the invariant picture down by offset pixels. */
    memmove(q, p, size);
    /* Recompute the newly revealed top portion */
    Compute(0, width, x_left,
            0, offset, new_y_top,
            (int) canvas_width, pixmap);
  } else {
    char* p = pixmap - offset * width;
    char* q = pixmap;
    size_t size = (width + offset) * width;

    /* Move the invariant picture up by offset pixels. */
    memmove(q, p, size);
    /* Recompute the newly revealed bottom portion */
    Compute(0, width, x_left,
            width + offset, width, y_bottom,
            (int) canvas_width, pixmap);
  }
}

void NewWindowPlugin::domenu(const QString &menu_name, V3DPluginCallback2 &v3d, QWidget *parent)
{
    if (menu_name == tr("Invert Color (in current window)"))
    {
    	Invert(v3d, parent);
    }
    else
    if (menu_name == tr("Threshold..."))
    {
    	Threshold(v3d, parent);
    }
    else
    if (menu_name == tr("Add 2 Images..."))
    {
    	Compute('+', v3d, parent);
    }
    else
    if (menu_name == tr("Differ 2 Images..."))
    {
    	Compute('-', v3d, parent);
    }
    else
    if (menu_name == tr("Modulate 2 Images..."))
    {
    	Compute('*', v3d, parent);
    }
    else
    {
    	QMessageBox::information(parent, title, "V3DPluginInterface Demo version 1.2"
    			"\ndeveloped by Zongcai Ruan 2009-2010. (Janelia Farm Research Campus, HHMI)");
    }
}

  //  
  //         f(1) - f(0)
  // f'(0) = -----------      f(1) = f(0+h)
  //              h
  //
  Eigen::Vector3d OBFunction::NumericalDerivative(unsigned int index)
  {
    double e_orig, e_plus_delta, delta, dx, dy, dz;
    delta = 1.0e-5;

    const Eigen::Vector3d va = m_positions.at(index);
    Compute(OBFunction::Value);
    e_orig = GetValue();
    
    // X direction
    m_positions[index].x() += delta;
    Compute(OBFunction::Value);
    e_plus_delta = GetValue();
    dx = (e_plus_delta - e_orig) / delta;
    
    // Y direction
    m_positions[index].x() = va.x();
    m_positions[index].y() += delta;
    Compute(OBFunction::Value);
    e_plus_delta = GetValue();
    dy = (e_plus_delta - e_orig) / delta;
    
    // Z direction
    m_positions[index].y() = va.y();
    m_positions[index].z() += delta;
    Compute(OBFunction::Value);
    e_plus_delta = GetValue();
    dz = (e_plus_delta - e_orig) / delta;

    // reset coordinates to original
    m_positions[index].z() = va.z();

    return Eigen::Vector3d(-dx, -dy, -dz);
  }

/*
 * Move a selected region left or right, reusing already computed values.
 */
static void MoveLR(int offset) {
  int i;
  int j;
  int width = canvas_width;

  if (offset > 0) {
    /* Move the invariant picture left by offset pixels. */
    for (i = 0; i < width; ++i) {
      for (j = 0; j < width - offset; ++j) {
        pixmap[i * width + j] = pixmap[i * width + j + offset];
      }
    }
    /* Recompute the newly revealed right portion */
    Compute(width - offset, width, x_right,
            0, width, y_top,
            (int) canvas_width, pixmap);
  } else {
    /* Move the invariant picture right by offset pixels. */
    for (i = 0; i < width; ++i) {
      for (j = width - 1; j >= -offset; --j) {
        pixmap[i * width + j] = pixmap[i * width + j + offset];
      }
    }
    /* Recompute the newly revealed left portion */
    Compute(0, -offset, new_x_left,
            0, width, y_top,
            (int) canvas_width, pixmap);
  }
}

void HiggsPlot::PlotBF(int iDecay, double tBmH_max, double BF_max){
  if(iDecay<0 || iDecay>2) {cout<<"iDecay must be 0, 1 or 2"<<endl; return;}
  styles style; style.setPadsStyle(-1); style.setDefaultStyle();
  int nBins = 1000;
  TString hName, epsName = "public_html/Higgs_BF_TEMP_BaBar.eps", label, Llabel, Rlabel;
  TString yTitle[] = {"BF(B#rightarrow#tau#nu) (10^{-5})", "R(D)", "R(D*)"};
  TString TagDecay[] = {"BF", "R(D)", "R(D*)"};
  TCanvas can;
  TH1F *hBF[2];
  for(int his=0; his<1; his++) {
    hName = "hBF"; hName += his;
    hBF[his] = new TH1F(hName,"",nBins,0,tBmH_max);
  }
  double tBmH, BF[2];
  for(int bin=1; bin<=nBins; bin++){
    tBmH = hBF[0]->GetBinCenter(bin);
    Compute(tBmH,BF,iDecay);
    hBF[0]->SetBinContent(bin, BF[0]);
    hBF[0]->SetBinError(bin, BF[1]);
  }
  hBF[0]->SetFillColor(2); hBF[0]->SetLineColor(2);
  hName += "1";
  hBF[1] = static_cast<TH1F*>(hBF[0]->Clone(hName));
  for(int bin=1; bin<=nBins; bin++)hBF[1]->SetBinError(bin,0);
  hBF[1]->SetFillColor(0);hBF[1]->SetLineColor(1); hBF[1]->SetLineWidth(2);

  TBox box; box.SetLineColor(4);box.SetFillColor(4);
  TLine line; line.SetLineStyle(1); line.SetLineColor(4); line.SetLineWidth(3);
  if(BF_max>0) hBF[0]->SetMaximum(BF_max);
  hBF[0]->Draw("e3");
  box.SetFillStyle(3002); 
  box.DrawBox(0,Measurement[iDecay][0]-Measurement[iDecay][1],
	      tBmH_max,Measurement[iDecay][0]+Measurement[iDecay][1]);
  line.DrawLine(0,Measurement[iDecay][0],tBmH_max,Measurement[iDecay][0]);
  hBF[0]->Draw("e3 same");
  hBF[1]->Draw("c same");

  Compute(0,BF,iDecay);
  label = "#splitline{"; label += TagDecay[iDecay]; label += "_{SM} = ";
  if(iDecay==0){
    label+="(";label+=RoundNumber(BF[0],1); label+=" #pm ";
    label+=RoundNumber(BF[1],1); label+=")#times10^{-5}}{BF_{exp} = (";
    label+=RoundNumber(Measurement[iDecay][0],1); label+=" #pm ";
    label+=RoundNumber(Measurement[iDecay][1],1); label+=")#times10^{-5}}";
    Llabel = ""; Rlabel = label;
  } else {
    label+=RoundNumber(BF[0],3); label+=" #pm ";
    label+=RoundNumber(BF[1],3); label+="}{"; label += TagDecay[iDecay]; label += "_{exp} = ";
    label+=RoundNumber(Measurement[iDecay][0],3); label+=" #pm ";
    label+=RoundNumber(Measurement[iDecay][1],3); label+="}";
    Rlabel = ""; Llabel = label;
  }
  style.setTitles(hBF[0],"tan#beta/m_{H^{+}} (GeV^{-1})",yTitle[iDecay],Llabel,Rlabel);
  epsName.ReplaceAll("TEMP",DecayName[iDecay]);
  can.SaveAs(epsName);
  for(int his=0; his<2; his++) hBF[his]->Delete();
}

double HiggsPlot::ProbChi2(int iDecay, double tBmH, double rL){
  double RD[2], RDs[2], chi2;
  double valMeas[3][2] = {{Measurement[0][0], Measurement[1][0]},
			  {MeasuredRD[0]->Eval(tBmH), MeasuredRD[2]->Eval(tBmH)},
			  {MeasuredRD[1]->Eval(tBmH), MeasuredRD[3]->Eval(tBmH)}};
  int ndof=1;
  if(iDecay<3) {
    Compute(tBmH,RD,iDecay);
    chi2 = Chi2(RD,valMeas[iDecay]);
  } else {
    if(iDecay==3) {Compute(tBmH,RD,1); Compute(tBmH,RDs,2);}
    else {
      double mb = 4.2, rRrL = tBmH + rL, tBmH_Eff = sqrt(fabs(rRrL/(mTau*mb))); 
      if(rRrL > 0) RDCoef[0][0][1] = fabs(RDCoef[0][0][1]);
      Compute(tBmH_Eff,RD,1);
      valMeas[1][0] = MeasuredRD[0]->Eval(tBmH_Eff);
      valMeas[1][1] = MeasuredRD[2]->Eval(tBmH_Eff);
      rRrL = tBmH - rL; tBmH_Eff = sqrt(fabs(rRrL/(mTau*mb)));
      if(rRrL > 0) RDCoef[1][0][1] = fabs(RDCoef[1][0][1]);
      Compute(tBmH_Eff,RDs,2);
      RDCoef[0][0][1] = -fabs(RDCoef[0][0][1]); RDCoef[1][0][1] = -fabs(RDCoef[1][0][1]);
      valMeas[2][0] = MeasuredRD[1]->Eval(tBmH_Eff);
      valMeas[2][1] = MeasuredRD[3]->Eval(tBmH_Eff);
//       cout<<RoundNumber(RD[0],3)<<" +- "<<RoundNumber(RD[1],3)<<", "
// 	  <<RoundNumber(RDs[0],3)<<" +- "<<RoundNumber(RDs[1],3)<<" at tBmH = "<<tBmH_Eff<<endl;
    }
    //if(tBmH==0) {RD[0] = 0.316; RD[1] = 0.014;}  // MILC 2012
    //if(tBmH==0) {RD[0] = 0.302; RD[1] = 0.016;}  // Tanaka 2010
    //if(tBmH==0) {RD[0] = 0.310; RD[1] = 0.020;}  // Nierste 2008
    TMatrixT<double> CovM(2,2), DiffRD(1,2), MDiff(1,2); //CovM is the covariance matrix
    if(_varyRD){
      CovM(0,0) = pow(RD[1],2) +pow(valMeas[1][1],2);
      CovM(1,1) = pow(RDs[1],2)+pow(valMeas[2][1],2);
      CovM(1,0) = Measurement[3][0]*valMeas[1][1]*valMeas[2][1];
      DiffRD(0,0) = valMeas[1][0]-RD[0]; DiffRD(0,1) = valMeas[2][0]-RDs[0];
    } else {
      CovM(0,0) = pow(RD[1],2) +pow(Measurement[1][1],2);
      CovM(1,1) = pow(RDs[1],2)+pow(Measurement[2][1],2);
      CovM(1,0) = Measurement[3][0]*Measurement[1][1]*Measurement[2][1];
      DiffRD(0,0) = Measurement[1][0]-RD[0]; DiffRD(0,1) = Measurement[2][0]-RDs[0];
    }
    CovM(0,1) = CovM(1,0);
    CovM.Invert();
    MDiff.Mult(DiffRD,CovM);
    chi2 = MDiff(0,0)*DiffRD(0,0) + MDiff(0,1)*DiffRD(0,1);
    ndof++;
//     if(tBmH==0 && _varyRD) 
//       cout<<endl<<"SM chi2 is "<<RoundNumber(chi2,2)<<" with a p-value of "<<TMath::Prob(chi2,ndof)<<". This is "
// 	  <<RoundNumber(sqrt(2)*TMath::ErfInverse(1 - TMath::Prob(chi2,ndof)),2)<<" sigma away"<<endl<<endl;
  }

  return 1 - TMath::Prob(chi2,ndof);
}

  /**
   * Computes and returns the number of valid pixels in the boxcar at pixel index
   * (zero based). No error checks are done for out of array bounds conditions.
   * The routine works the fastest when sequentially accessing the counts
   * (e.g., index = 0,1,2,...).
   *
   * @param index Zero based sample position
   *
   * @return double
   */
  int QuickFilter::Count(const int index) {
    // Move the boxcar if necessary
    Compute(index);

    // Return the valid count
    return p_lastCount;
  }

void Machine::Start() {
	if ( true ) {
		i = 0;
		s = 0;
		Compute();
	}
}

bool
DebugReplayIGC::Next()
{
  last_basic = basic;
  last_calculated = calculated;

  const char *line;
  while ((line = reader->read()) != NULL) {
    if (line[0] == 'B') {
      IGCFix fix;
      if (IGCParseFix(line, extensions, fix) && fix.gps_valid) {
        CopyFromFix(fix);

        Compute();
        return true;
      }
    } else if (line[0] == 'H') {
      BrokenDate date;
      if (memcmp(line, "HFDTE", 5) == 0 &&
          IGCParseDateRecord(line, date)) {
        (BrokenDate &)basic.date_time_utc = date;
        basic.date_available = true;
      }
    } else if (line[0] == 'I') {
      IGCParseExtensions(line, extensions);
    }
  }

  return false;
}

void ScaledGeometricalMoments<VoxelT,MomentT>::Init (
        const VoxelT* _voxels,
        int _xDim, int _yDim, int _zDim, 
        double _xCOG, double _yCOG, double _zCOG,
        double _scale,
        int _maxOrder)
{
    xDim_ = _xDim;
    yDim_ = _yDim;
    zDim_ = _zDim;

    maxOrder_ = _maxOrder;

    size_t totalSize = xDim_ * yDim_ * zDim_;
    voxels_.resize (totalSize);
    for (int i=0; i<totalSize; ++i)
    {
        voxels_[i] = _voxels[i];
    }

    moments_.resize (maxOrder_ + 1);
    for (int i=0; i<=maxOrder_; ++i)
    {
        moments_[i].resize (maxOrder_ - i + 1);
        for (int j=0; j<=maxOrder_ - i; ++j)
        {
            moments_[i][j].resize (maxOrder_ - i - j + 1);
        }
    }

    ComputeSamples (_xCOG, _yCOG, _zCOG, _scale);

    Compute ();
}   

void ControlReactivo::SetCommand(float vela, float velg,double dist)
{
    va=vela;
    vg=velg;
    disttray=dist;
    Compute();
}

static sci_tab_t *SelectN(
  unsigned int 	B,
  double 	Phi )
{
  unsigned int i ;
  struct sci_tab* best = NULL ;

  for ( i = 0 ; i < 4 ; i++ )
  {
    double err ;

    Compute( i, B, Phi, &test_array[i] );
    err = fabs( test_array[i].err );

    if ( best )
    {
      if ( err < fabs( best->err ) )
        best = &test_array[i] ;
    }
    else
      best = &test_array[i] ;
  }

  return best ;
}

void Pdb::Explorer()
{
	VectorMap<String, Value> prev = DataMap(explorer);
	explorer.Clear();
	try {
		String x = ~expexp;
		if(!IsNull(x)) {
			CParser p(x);
			Val v = Exp(p);
			Vis(explorer, "=", prev, Visualise(v));
			if(v.type >= 0 && v.ref == 0 && !v.rvalue)
				Explore(v, prev);
			if(v.ref > 0 && GetRVal(v).address)
				for(int i = 0; i < 20; i++)
					Vis(explorer, Format("[%d]", i), prev, Visualise(DeRef(Compute(v, RValue(i), '+'))));
		}
	}
	catch(CParser::Error e) {
		Visual v;
		v.Cat(e, LtRed);
		explorer.Add("", RawPickToValue(v));
	}
	exback.Enable(exprev.GetCount());
	exfw.Enable(exnext.GetCount());
}

//==============================================================================
Ifpack2_Amesos::Ifpack2_Amesos(const Ifpack2_Amesos& rhs) :
  Matrix_(Teuchos::rcp( &rhs.Matrix(), false )),
  Label_(rhs.Label()),
  IsInitialized_(false),
  IsComputed_(false),
  NumInitialize_(rhs.NumInitialize()),
  NumCompute_(rhs.NumCompute()),
  NumApplyInverse_(rhs.NumApplyInverse()),
  InitializeTime_(rhs.InitializeTime()),
  ComputeTime_(rhs.ComputeTime()),
  ApplyInverseTime_(rhs.ApplyInverseTime()),
  ComputeFlops_(rhs.ComputeFlops()),
  ApplyInverseFlops_(rhs.ApplyInverseFlops()),
  Condest_(rhs.Condest())
{

  Problem_ = Teuchos::rcp( new Tpetra_LinearProblem );

  // copy the RHS list in *this.List
  Teuchos::ParameterList RHSList(rhs.List());
  List_ = RHSList;

  // I do not have a copy constructor for Amesos,
  // so Initialize() and Compute() of this object 
  // are called if the rhs did so
  if (rhs.IsInitialized()) {
    IsInitialized_ = true;
    Initialize();
  }
  if (rhs.IsComputed()) {
    IsComputed_ = true;
    Compute();
  }

}

bool
DebugReplayNMEA::Next()
{
  last_basic = computed_basic;

  const char *line;
  while ((line = reader->ReadLine()) != NULL) {
    raw_basic.clock = clock.NextClock(raw_basic.time_available
                                      ? raw_basic.time
                                      : -1.);

    if (!device || !device->ParseNMEA(line, raw_basic))
      parser.ParseLine(line, raw_basic);

    if (raw_basic.location_available != last_basic.location_available) {
      Compute();
      return true;
    }
  }

  if (computed_basic.time_available)
    flying_computer.Finish(calculated.flight, computed_basic.time);

  return false;
}