//---------------------------------------------------------
DVec& NDG2D::PoissonIPDGbc2D
(DVec& ubc, //[in]
 DVec& qbc  //[in]
)
//---------------------------------------------------------
{
  // function [OP] = PoissonIPDGbc2D()
  // Purpose: Set up the discrete Poisson matrix directly
  //          using LDG. The operator is set up in the weak form

  // build DG derivative matrices
  int max_OP = (K*Np*Np*(1+Nfaces));

  // initialize parameters
  DVec faceR("faceR"), faceS("faceS");
  DMat V1D("V1D"), Dx("Dx"),Dy("Dy"), Dn1("Dn1"), mmE_Fm1("mmE(:,Fm1)");
  IVec Fm("Fm"), Fm1("Fm1"), fidM("fidM");
  double lnx=0.0,lny=0.0,lsJ=0.0,hinv=0.0,gtau=0.0;
  int i=0,k1=0,f1=0,id=0;
  IVec i1_Nfp = Range(1,Nfp);
  double N1N1 = double((N+1)*(N+1));

  // build local face matrices
  DMat massEdge[4]; // = zeros(Np,Np,Nfaces);
  for (i=1; i<=Nfaces; ++i) {
    massEdge[i].resize(Np,Np);
  }

  // face mass matrix 1
  Fm = Fmask(All,1); faceR = r(Fm); 
  V1D = Vandermonde1D(N, faceR);
  massEdge[1](Fm,Fm) = inv(V1D*trans(V1D));

  // face mass matrix 2
  Fm = Fmask(All,2); faceR = r(Fm); 
  V1D = Vandermonde1D(N, faceR);
  massEdge[2](Fm,Fm) = inv(V1D*trans(V1D));

  // face mass matrix 3
  Fm = Fmask(All,3); faceS = s(Fm); 
  V1D = Vandermonde1D(N, faceS); 
  massEdge[3](Fm,Fm) = inv(V1D*trans(V1D));
 
  // build DG right hand side
  DVec* pBC = new DVec(Np*K, "bc", OBJ_temp); 
  DVec& bc = (*pBC);  // reference, for syntax
  ////////////////////////////////////////////////////////////////

  umMSG(1, "\n ==> {OP} assembly [bc]: ");
  for (k1=1; k1<=K; ++k1)
  {
    if (! (k1%100)) { umMSG(1, "%d, ",k1); }

    // rows1 = outer(Range((k1-1)*Np+1,k1*Np), Ones(NGauss));

    // Build element-to-element parts of operator
    for (f1=1; f1<=Nfaces; ++f1)
    {
      if (BCType(k1,f1))
      {
        ////////////////////////added by Kevin ///////////////////////////////
        Fm1 = Fmask(All,f1); 
        fidM  = (k1-1)*Nfp*Nfaces + (f1-1)*Nfp + i1_Nfp;
        id = 1+(f1-1)*Nfp + (k1-1)*Nfp*Nfaces;

        lnx = nx(id); lny = ny(id); 
        lsJ = sJ(id); hinv = Fscale(id);

        Dx = rx(1,k1)*Dr + sx(1,k1)*Ds;  
        Dy = ry(1,k1)*Dr + sy(1,k1)*Ds;
        Dn1 = lnx*Dx + lny*Dy;

      //mmE = lsJ*massEdge(:,:,f1);
      //bc(All,k1) += (gtau*mmE(All,Fm1) - Dn1'*mmE(All,Fm1))*ubc(fidM);

        mmE_Fm1 = massEdge[f1](All,Fm1);  mmE_Fm1 *= lsJ;

        gtau = 10*N1N1*hinv; // set penalty scaling
        //bc(All,k1) += (gtau*mmE_Fm1 - trans(Dn1)*mmE_Fm1) * ubc(fidM);

        switch(BCType(k1,f1)){
	  case BC_Dirichlet: 
            bc(Np*(k1-1)+Range(1,Np)) += (gtau*mmE_Fm1 - trans(Dn1)*mmE_Fm1)*ubc(fidM);  
            break;
          case BC_Neuman:
            bc(Np*(k1-1)+Range(1,Np)) += mmE_Fm1*qbc(fidM);
            break;
	default:
	  std::cout<<"warning: boundary condition is incorrect"<<std::endl;
	}
      }
    }
  }
  return bc;
}

 inline oclMat oclMat::colRange(int startcol, int endcol) const
 {
     return oclMat(*this, Range::all(), Range(startcol, endcol));
 }

void OpponentColorDescriptorExtractor::computeImpl( const Mat& bgrImage, std::vector<KeyPoint>& keypoints, Mat& descriptors ) const
{
    std::vector<Mat> opponentChannels;
    convertBGRImageToOpponentColorSpace( bgrImage, opponentChannels );

    const int N = 3; // channels count
    std::vector<KeyPoint> channelKeypoints[N];
    Mat channelDescriptors[N];
    std::vector<int> idxs[N];

    // Compute descriptors three times, once for each Opponent channel to concatenate into a single color descriptor
    int maxKeypointsCount = 0;
    for( int ci = 0; ci < N; ci++ )
    {
        channelKeypoints[ci].insert( channelKeypoints[ci].begin(), keypoints.begin(), keypoints.end() );
        // Use class_id member to get indices into initial keypoints vector
        for( size_t ki = 0; ki < channelKeypoints[ci].size(); ki++ )
            channelKeypoints[ci][ki].class_id = (int)ki;

        descriptorExtractor->compute( opponentChannels[ci], channelKeypoints[ci], channelDescriptors[ci] );
        idxs[ci].resize( channelKeypoints[ci].size() );
        for( size_t ki = 0; ki < channelKeypoints[ci].size(); ki++ )
        {
            idxs[ci][ki] = (int)ki;
        }
        std::sort( idxs[ci].begin(), idxs[ci].end(), KP_LessThan(channelKeypoints[ci]) );
        maxKeypointsCount = std::max( maxKeypointsCount, (int)channelKeypoints[ci].size());
    }

    std::vector<KeyPoint> outKeypoints;
    outKeypoints.reserve( keypoints.size() );

    int dSize = descriptorExtractor->descriptorSize();
    Mat mergedDescriptors( maxKeypointsCount, 3*dSize, descriptorExtractor->descriptorType() );
    int mergedCount = 0;
    // cp - current channel position
    size_t cp[] = {0, 0, 0};
    while( cp[0] < channelKeypoints[0].size() &&
           cp[1] < channelKeypoints[1].size() &&
           cp[2] < channelKeypoints[2].size() )
    {
        const int maxInitIdx = std::max( 0, std::max( channelKeypoints[0][idxs[0][cp[0]]].class_id,
                                                      std::max( channelKeypoints[1][idxs[1][cp[1]]].class_id,
                                                                channelKeypoints[2][idxs[2][cp[2]]].class_id ) ) );

        while( channelKeypoints[0][idxs[0][cp[0]]].class_id < maxInitIdx && cp[0] < channelKeypoints[0].size() ) { cp[0]++; }
        while( channelKeypoints[1][idxs[1][cp[1]]].class_id < maxInitIdx && cp[1] < channelKeypoints[1].size() ) { cp[1]++; }
        while( channelKeypoints[2][idxs[2][cp[2]]].class_id < maxInitIdx && cp[2] < channelKeypoints[2].size() ) { cp[2]++; }
        if( cp[0] >= channelKeypoints[0].size() || cp[1] >= channelKeypoints[1].size() || cp[2] >= channelKeypoints[2].size() )
            break;

        if( channelKeypoints[0][idxs[0][cp[0]]].class_id == maxInitIdx &&
            channelKeypoints[1][idxs[1][cp[1]]].class_id == maxInitIdx &&
            channelKeypoints[2][idxs[2][cp[2]]].class_id == maxInitIdx )
        {
            outKeypoints.push_back( keypoints[maxInitIdx] );
            // merge descriptors
            for( int ci = 0; ci < N; ci++ )
            {
                Mat dst = mergedDescriptors(Range(mergedCount, mergedCount+1), Range(ci*dSize, (ci+1)*dSize));
                channelDescriptors[ci].row( idxs[ci][cp[ci]] ).copyTo( dst );
                cp[ci]++;
            }
            mergedCount++;
        }
    }
    mergedDescriptors.rowRange(0, mergedCount).copyTo( descriptors );
    std::swap( outKeypoints, keypoints );
}

//\begin{>>PPP_2d_Jacobi_1d_Part.tex}{\subsubsection{runBenchmark}}
void PPP_2d_Jacobi_1d_Part::
runBenchmark(  )
//================================================================
// /Description:  Executes the benchmarking code the specified
//  number of times and stores the times in the array.
//
// /Return Values: None.
//
// /Errors:
//   None.
//
// /Author:  BJM
// /Date:  24 August 2000
//\end{PPP_2d_Jacobi_1d_Part.tex}
//================================================================
{
  double theTime,theStartTime,theEndTime;
  
  int theNumProcs = 1;
  int thisProcNum = 0;

// MPI_Comm_size(MPI_COMM_WORLD, &theNumProcs);
// MPI_Comm_rank(MPI_COMM_WORLD, &thisProcNum);
 
  Partitioning_Type Partitioning(Range(0,theNumProcs-1));
  Partitioning.SpecifyInternalGhostBoundaryWidths(1,1);
  Partitioning.SpecifyDecompositionAxes(1);

  int theArraySideLength = sqrt( mUnknownsPerProc * theNumProcs );

  printf ("theArraySideLength = %d \n",theArraySideLength);
  
  doubleArray U_old(theArraySideLength,theArraySideLength,Partitioning);
  U_old = 0.0;

  Index I (1,theArraySideLength-2);
  Index J (1,theArraySideLength-2);
  int i;

// printf ("Warming up ... \n");
  for(i=0;i< mNumberOfWarmupIterations;i++)
  {
    U_old(I,J) =
      (U_old(I+1,J+1) + U_old(I+1,J) + U_old(I+1,J-1) + U_old(I,J+1) +
       U_old(I,J-1) + U_old(I-1,J+1) + U_old(I-1,J) + U_old(I-1,J-1)) / 8.0;
  }
  
  // Now time the problem
// printf ("Running timing loop ... \n");
  for (i=0; i < mNumberOfTimingIterations; i++)
  {
//  printf ("Running timing loop iteration %d \n",i);
    theStartTime = Communication_Manager::Wall_Clock_Time();
//    theStartTime = clock();
    U_old(I,J) =
      (U_old(I+1,J+1) + U_old(I+1,J) + U_old(I+1,J-1) + U_old(I,J+1) +
       U_old(I,J-1) + U_old(I-1,J+1) + U_old(I-1,J) + U_old(I-1,J-1)) / 8.0;
    theEndTime = Communication_Manager::Wall_Clock_Time();
//    theEndTime = clock();
    theTime = theEndTime - theStartTime;
    
    mTimes[0][i] = theTime;
    printf("time= %f\n",theTime);
    
  }//end of loop over timed iterations
}

Range<Scalar,Dim> Range<Scalar,Dim>::unitRange()
{
    return Range(Vector<Scalar,Dim>(0),Vector<Scalar,Dim>(1.0));
}

 static Range Invalid() {
     return Range(std::numeric_limits<size_t>::max(),
                  std::numeric_limits<size_t>::min());
 }

 //! return shifted Range
 Range operator + (const size_t& shift) const {
     return Range(begin + shift, end + shift);
 }

//.........这里部分代码省略.........
			}
		}*/
#pragma endregion
	}

#pragma region // --- OpenGLにデータ渡す ---
	
	//メモリ割り当て
	allpoints = (GLfloat *)malloc(sizeof(float)*rows*fileTotal*cols);
	for(fileCount=0;fileCount<fileTotal;fileCount++)
	{
		//座標値をallpointsに入れる
		for(int i=0;i<rows;i++){
			for(int j=0;j<cols;j++){
				allpoints[fileCount*rows*cols+i*cols+j] = shape_fixed[fileCount].at<float>(i,j);
			}
		}
	}
#pragma endregion

	
#pragma region // --- カメラRTの計算 ---
	Mat cameraRT[fileTotal];
	Mat cameraR[fileTotal];
	Mat cameraT[fileTotal];
	cameraRT[0] = Mat::eye(4,4,CV_32F);
	cameraR[0] = Mat::eye(3,3,CV_32F);
	cameraT[0] = Mat::zeros(1,3,CV_32F);
	for(i=1;i<fileTotal;i++){
		cameraRT[i] = Mat::eye(4,4,CV_32F);
		cameraR[i] = Mat::eye(3,3,CV_32F);
		cameraT[i] = Mat::zeros(1,3,CV_32F);
		
		Mat r = my_rt[i].operator()(Range(0,3),Range(0,3));
		cameraR[i] = cameraR[i-1]*r.t();
		Mat t = my_rt[i].operator()(Range(3,4),Range(0,3));
		cameraT[i] = t*cameraR[i-1].t() + cameraT[i-1];
		
		cameraRT[i].at<float>(0,0) = cameraR[i].at<float>(0,0);
		cameraRT[i].at<float>(0,1) = cameraR[i].at<float>(0,1);
		cameraRT[i].at<float>(0,2) = cameraR[i].at<float>(0,2);
		cameraRT[i].at<float>(1,0) = cameraR[i].at<float>(1,0);
		cameraRT[i].at<float>(1,1) = cameraR[i].at<float>(1,1);
		cameraRT[i].at<float>(1,2) = cameraR[i].at<float>(1,2);
		cameraRT[i].at<float>(2,0) = cameraR[i].at<float>(2,0);
		cameraRT[i].at<float>(2,1) = cameraR[i].at<float>(2,1);
		cameraRT[i].at<float>(2,2) = cameraR[i].at<float>(2,2);
		
		cameraRT[i].at<float>(3,0) = cameraT[i].at<float>(0,0);
		cameraRT[i].at<float>(3,1) = cameraT[i].at<float>(0,1);
		cameraRT[i].at<float>(3,2) = cameraT[i].at<float>(0,2);
	}
#pragma endregion

// --- データ出力 ---
#if FILEOUTPUT

	///////////////////////////////
	// 全ての点群(shape_fixed)をまとめて書き出し
	// pcd
	//
	FILE *outfp;
	char outfilename[100];
	sprintf(outfilename,"%s/%s/result_xyz.pcd",outdir,dir);
	outfp = fopen(outfilename,"w");
	if(outfp == NULL){

//! given a global range [0,global_size) and p PEs to split the range, calculate
//! the [local_begin,local_end) index range assigned to the PE i.
static inline Range CalculateLocalRange(
    size_t global_size, size_t p, size_t i) {
    return Range(0, global_size).Partition(i, p);
}

void X_Comp::EvalProducts()
  {
  if (NJoins()>0)
    switch (SolveMode())
      {
      case PBMODE:
        DoBreak();
        break;
      case SSMODE:
        DoBreak();
        break;
      case DYNMODE:
        {
        int JoinId=0;
        double Press=Joins[JoinId].Pressure();

        int I[MaxIOList+1];
        FillIOIndexList(JoinId, I);

        flag dbg=0;//((Debug|dbgDerivs) & DBG_Derivs);

        if (SolveMode()==PBMODE)
          Press=GetPBInputPressure(JoinId);

        StkSpConduit SdLcl("SdLcl", chLINEID(), this);
        SpConduit &Sd = SdLcl();
        Sd.QZero();
        Sd.SetPress(Press);
        double  Qot=0.0;
        int NFeeds=0;
        for (int i, ic = 0; (i=I[ic])>=0; ic++)
          if (IO_In(i))
            {
            Sd.SelectModel(IOConduit(i), NFeeds==0);
            Sd.QAddF(*IOConduit(i), som_ALL, 1.0);
            NFeeds++;
            }
          else
            Qot += IOQmEst_Out(i);

        if (dbg)
          for (ic = 0; (i=I[ic])>=0; ic++)
            if (IO_In(i))
              {
              char t[128];
              sprintf(t, "%s.%s", FullObjTag(), IODesc_Self(i)->pName);
              dbgpln("i-Xfr %-12.12s : %14.6g, %14.6g | %14.6g",t,
                     IOConduit(i)->QMass(som_SL), IOConduit(i)->QMass(som_Vap), K_2_C(IOConduit(i)->Temp()));
              }

        double denom;
        double  P1 = IOP_Rmt (IOWithId_Self(ioid_In));
        double  P2 = IOP_Self(IOWithId_Self(ioid_In));

        // hss 9/3/98 - Use data base for specific heat ratio.
        P_SpecificHeatRatio = Max(1.1, Sd.CpCv());

        if( Sd.Temp() > 0.01 )
          Tin = Sd.Temp();

        // hss Calc Polytropic Efficiency
        /*double FlowMin  = pI->QVolume() * 60.0;
        double rpm      = Max(SpeedRatio, 10000.0);
        double EffTemp  = Efficiency.Zxy(FlowMin,rpm);
        if (EffTemp >= 0.6)
          P_PolytropicEff = EffTemp;
        else
          P_PolytropicEff = 0.6;*/

        denom = P_SpecificHeatRatio * P_PolytropicEff;
        if( fabs(denom) < 1.0e-30)
          denom = 1.0e-30;
        // hss Try a fix to prevent crash when P1 is negative
        if ((P1 > 0.0) && (P2 > P1))
          Tout = Tin*pow(P2/P1,(P_SpecificHeatRatio - 1.0)/denom);
        else
          Tout = Tin;      // end of kluge

        Sd.SetTemp(Tout);

        double Qin=Sd.QMass(som_ALL);
        // What Comes in must go out
        double Scl=Range(0.0, Qin/GTZ(Qot), 1000.0);
        for (ic = 0; (i=I[ic])>=0; ic++)
          if (IO_Out(i))
            IOConduit(i)->QSetM(Sd, som_ALL, IOQmEst_Out(i)*Scl, Press);//Joins[JoinId].Pressure());//PMax);
        }
        break;
      }
  };

void X_Comp::EvalDiscrete()
  {
  if (Control==CTSpeed)
    {
    SpConduit &C=*IOConduit(IOWithId_Self(ioid_In));
    //double Cflow = GEZ(IOFB(IOWithId_Self(ioid_In), 0)->GetQm());
    //double VFlow = GEZ(IOFB(IOWithId_Self(ioid_In), 0)->GetQv());
    double Qv=C.QVolume();
    double Qm=C.QMass();

    MaxSpeed=Max(100.0, MaxSpeed);
    MinSpeed=Range(0.0, MinSpeed, MaxSpeed);
    double RqdS=(P_Status ? Range(MinSpeed/MaxSpeed, SpdSpt,1.0) : 0);
    SpdFbk+=(RqdS-SpdFbk)*ICGetTimeInc()/Max(1.0, SpdTau);
    ActSpeed=SpdFbk*MaxSpeed;
    SpeedRatio=SpdFbk;

    double Hd=-166.466+1.581213*ActSpeed+(-17.5093-0.15871*ActSpeed)*Qv+(0.0003296+0.014441*ActSpeed)*Qv*Qv;
//    double SrgVolFlw=-0.96597+0.010344*ActSpeed;
    double SrgSpd=Range(166.0, ActSpeed, 275.0);
    QvSrg=7E-05*SrgSpd*SrgSpd - 0.0199*SrgSpd + 2.2831;
    QvIn=C.QVolume();

    //Qs=Range(0.0, Qm*SrgVolFlw/GTZ(Qm), 150.0);

    Pin  = Max(100.0,IOP_Rmt (IOWithId_Self(ioid_In)));
//    double PCalc=(7990.392-21.47525359*SpdFbk*266.67+26.75703018*Range(0.0,Cflow,100.0));
//    PCalc=Max(1000.0, PCalc);
//    double PRatio=7990/PCalc;

    //PressIn=GTZ(PressIn);
//    if ((Qv>1.0e-6) && (C.QMass()>1.0e-6))
//      { // Only Change if Direction is Forward
    double GammaIn=1.262;//Range(1.2, C.CpCv(), 1.5);
    double MoleWtFIn=C.MoleWt();//*Rho;
    double TempIn=C.Temp();
//      }


    double R=8.314/MoleWtFIn;
    double X=GammaIn/(GammaIn-1.0);
    // PHT_Isentropic:

    // Units of Hd Expected to be kNm/kg == kJ/kg;
    // See Perry P6-17 & WMC-KNS Acid Plant Curves.
    double dEfficiency=0.9;
    double P2=Pin*Pow(GEZ(Hd*dEfficiency/GTZ(X*R*TempIn)+1.0), X);
    double Bst=GEZ(P2-Pin);

    double TheBoostDmp=1.0-(1.0/Max(ICGetTimeInc(), SpdTau));
    TheBoost=TheBoostDmp*TheBoost+(1-TheBoostDmp)*GEZ(Bst);
    //TheBoost=TheBoostDmp*TheBoost+(1-TheBoostDmp)*GEZ(Pin*PRatio-Pin);

    }
  else
    {
    double Q     = IOQm_In(IOWithId_Self(ioid_In));
    double P1    = IOP_Rmt(IOWithId_Self(ioid_In));
    double P2    = IOP_Self(IOWithId_Self(ioid_In));
    double pwr   = Pwr_Curve(Q,P1,P2);
    double goodp = P2_Curve(Q,P1,P_MaxPower);
    double Qc    = Flw_Curve(P_MaxPower,P1,P_SetpointPressure);

    if( fabs(Q - Qprv)/Max(1.0e-60,Qprv) < 0.001 && ((P2 < P1) || (Q < 0.001)) )
      {
      X_SetpointPressure = 0.0;
      }
    else if( Q > 1.0e-50 && fabs(Q - Qprv)/Max(1.0e-60,Qprv) < 0.001 && P1 < P2 )
      {
      if( Q > Qc )
        {
        double newval = (P_SetpointPressure - goodp);
        double test   = (X_SetpointPressure*0.5) + (newval*0.5);
        if( P2-test < P1 )
          {// This will fail at unrealistically high flow
          X_SetpointPressure = (X_SetpointPressure*0.5) + 0.5*(P2-P1);
          }
        else
          {
          X_SetpointPressure = test;
          }
        //TRACE("EVAL DISCREET CONVERGED %f %f %f %f\n\n",Q,pwr,goodp,X_SetpointPressure);
        }
      else
        {
        X_SetpointPressure /= 2.0;
        }
      }
    Qprv = Q;
    }
  }

    // Test object detection network from Darknet framework.
    void testDarknetModel(const std::string& cfg, const std::string& weights,
                          const std::vector<std::vector<int> >& refClassIds,
                          const std::vector<std::vector<float> >& refConfidences,
                          const std::vector<std::vector<Rect2d> >& refBoxes,
                          double scoreDiff, double iouDiff, float confThreshold = 0.24, float nmsThreshold = 0.4)
    {
        checkBackend();

        Mat img1 = imread(_tf("dog416.png"));
        Mat img2 = imread(_tf("street.png"));
        std::vector<Mat> samples(2);
        samples[0] = img1; samples[1] = img2;

        // determine test type, whether batch or single img
        int batch_size = refClassIds.size();
        CV_Assert(batch_size == 1 || batch_size == 2);
        samples.resize(batch_size);

        Mat inp = blobFromImages(samples, 1.0/255, Size(416, 416), Scalar(), true, false);

        Net net = readNet(findDataFile("dnn/" + cfg, false),
                          findDataFile("dnn/" + weights, false));
        net.setPreferableBackend(backend);
        net.setPreferableTarget(target);
        net.setInput(inp);
        std::vector<Mat> outs;
        net.forward(outs, getOutputsNames(net));

        for (int b = 0; b < batch_size; ++b)
        {
            std::vector<int> classIds;
            std::vector<float> confidences;
            std::vector<Rect2d> boxes;
            for (int i = 0; i < outs.size(); ++i)
            {
                Mat out;
                if (batch_size > 1){
                    // get the sample slice from 3D matrix (batch, box, classes+5)
                    Range ranges[3] = {Range(b, b+1), Range::all(), Range::all()};
                    out = outs[i](ranges).reshape(1, outs[i].size[1]);
                }else{
                    out = outs[i];
                }
                for (int j = 0; j < out.rows; ++j)
                {
                    Mat scores = out.row(j).colRange(5, out.cols);
                    double confidence;
                    Point maxLoc;
                    minMaxLoc(scores, 0, &confidence, 0, &maxLoc);

                    if (confidence > confThreshold) {
                        float* detection = out.ptr<float>(j);
                        double centerX = detection[0];
                        double centerY = detection[1];
                        double width = detection[2];
                        double height = detection[3];
                        boxes.push_back(Rect2d(centerX - 0.5 * width, centerY - 0.5 * height,
                                            width, height));
                        confidences.push_back(confidence);
                        classIds.push_back(maxLoc.x);
                    }
                }
            }

            // here we need NMS of boxes
            std::vector<int> indices;
            NMSBoxes(boxes, confidences, confThreshold, nmsThreshold, indices);

            std::vector<int> nms_classIds;
            std::vector<float> nms_confidences;
            std::vector<Rect2d> nms_boxes;

            for (size_t i = 0; i < indices.size(); ++i)
            {
                int idx = indices[i];
                Rect2d box = boxes[idx];
                float conf = confidences[idx];
                int class_id = classIds[idx];
                nms_boxes.push_back(box);
                nms_confidences.push_back(conf);
                nms_classIds.push_back(class_id);
            }

            normAssertDetections(refClassIds[b], refConfidences[b], refBoxes[b], nms_classIds,
                             nms_confidences, nms_boxes, format("batch size %d, sample %d\n", batch_size, b).c_str(), confThreshold, scoreDiff, iouDiff);
        }
    }

void SubmatrixQueriesTest::multipleBenchmarksRowQueries(size_t n, bench_time_t *naiveTime, bench_time_t *queryTime, bench_time_t *cascadingTime, bench_time_t *simpleCascadingTime)
{
    bench_time_t clock1, clock2;
    for (size_t i = 0; i < n; i++) {        
        size_t row = rand() % (_testMatrix->rows());
        
        size_t c1,c2;
        c1 = rand() % (_testMatrix->cols());
        c2 = rand() % (_testMatrix->cols());
        
        Range c = Range(min(c1,c2),max(c1,c2));

        clock1 = now();
        _queryDS->columnTree()->cascadingMaxInRange(row, c);
        clock2 = now();
        
        *cascadingTime = add(*cascadingTime,diff(clock2, clock1));
    }
    for (size_t i = 0; i < n; i++) {

        size_t row = rand() % (_testMatrix->rows());
        
        size_t c1,c2;
        c1 = rand() % (_testMatrix->cols());
        c2 = rand() % (_testMatrix->cols());
        
        Range c = Range(min(c1,c2),max(c1,c2));
        
        clock1 = now();
        _queryDS->columnTree()->simpleCascadingMaxInRange(row, c);
        clock2 = now();
        
        *simpleCascadingTime = add(*simpleCascadingTime,diff(clock2, clock1));
    }
    for (size_t i = 0; i < n; i++) {
        size_t row = rand() % (_testMatrix->rows());
        
        size_t c1,c2;
        c1 = rand() % (_testMatrix->cols());
        c2 = rand() % (_testMatrix->cols());
        
        Range c = Range(min(c1,c2),max(c1,c2));
        
        clock1 = now();
        _queryDS->columnTree()->maxForRowInRange(row, c.min, c.max);
        clock2 = now();
        
        *queryTime = add(*queryTime,diff(clock2, clock1));
    }
    if(SubmatrixQueriesTest::benchmarkNaiveQueries){
        for (size_t i = 0; i < n; i++) {
            size_t row = rand() % (_testMatrix->rows());
            
            size_t c1,c2;
            c1 = rand() % (_testMatrix->cols());
            c2 = rand() % (_testMatrix->cols());
            
            Range c = Range(min(c1,c2),max(c1,c2));

            clock1 = now();
            SubmatrixQueriesTest::naiveMaximumInRow(_testMatrix, c, row);
            clock2 = now();
            
            *naiveTime = add(*naiveTime,diff(clock2, clock1));
        }
    }
}

 //! calculate a partition range [begin,end) by taking the current Range
 //! splitting it into p parts and taking the i-th one.
 Range Partition(size_t i, size_t parts) const {
     assert(i < parts);
     return Range(CalculateBeginOfPart(i, parts),
                  CalculateBeginOfPart(i + 1, parts));
 }

void my_init(vector<bool> keepcol){
  Partial.extend(num_var_tape_ * 1);
  arg_mark_.resize(play_.op_arg_rec_.size());
  for(size_t i=0;i<arg_mark_.size();i++)arg_mark_[i]=false;
  /* Run a reverse test-sweep to store pointers once */
  tape_point tp;
  play_.reverse_start(tp.op, tp.op_arg, tp.op_index, tp.var_index);
  tp_.resize(tp.op_index+1);
  var2op_.resize(tp.var_index+1);
  op_mark_.resize(tp.op_index+1);
  for(size_t i=0;i<op_mark_.size();i++)op_mark_[i]=0;
  user_region_mark_.resize(tp.op_index+1);
  for(size_t i=0;i<user_region_mark_.size();i++)user_region_mark_[i]=0;
  tp_[tp.op_index]=tp;
  /* 1. We need to be able to find out, for a given variable, what operator created 
     the variable. This is easiest done by looping through the _operators_ because for a 
     given op we have access to all the resulting variables it creates.
     2. We precompute the a vector of "tape_points" so that instead of calling 
     "reverse_next", we simply get the next tape entry by tp_[i-1].
  */
  while(tp.op != BeginOp ){ /* tp.op_index is decremented by one in each iteration ... */
    // printTP(tp); /* For debugging */
    play_.reverse_next(tp.op, tp.op_arg, tp.op_index, tp.var_index);
    /* Csum is special case - see remarks in player.hpp and reverse_sweep.hpp */
    if(tp.op == CSumOp)play_.reverse_csum(tp.op, tp.op_arg, tp.op_index, tp.var_index);
    for(size_t i=0;i<NumRes(tp.op);i++)var2op_[tp.var_index-i]=tp.op_index;
    tp_[tp.op_index]=tp;
    markArgs(tp);
  }
  /* Lookup table: is tape_point within a UserOp region? */
  bool user_within=false;
  user_region_.resize(tp_.size());
  for(size_t i=0;i<tp_.size();i++){
    if(tp_[i].op==UserOp){
      user_region_[i]=true;
      user_within=!user_within;	
    } else {
      user_region_[i]=user_within;
    }
  }

  /* Lookup table: is tape_point a constant (=only fixed effect dependent) ? */
  constant_tape_point_.resize(tp_.size());
  int indep_var_number=0;
  for(size_t i=0;i<tp_.size();i++){
    if(tp_[i].op==InvOp){ /* All independent variables are marked according to being
			     random or fixed effect */
      constant_tape_point_[i]=!keepcol[indep_var_number];
      indep_var_number++;
    } else { /* Mark operator as constant if _all_ arguments are constant */
      constant_tape_point_[i] = is_tape_point_constant(i);
    }

    //std::cout << constant_tape_point_[i] << " "; printTP(tp_[i]);

  }
  // std::cout << "Total:   " << constant_tape_point_.size() << "\n";
  // int sum=0; for(int i=0;i<constant_tape_point_.size();i++)sum+=constant_tape_point_[i];
  // std::cout << "Constant:" << sum << "\n";


  // Calculate pattern
  int m=Range();
  colpattern.resize(m);
  for(int i=0;i<m;i++)my_pattern(i);
  for(size_t i=0;i<op_mark_.size();i++)op_mark_[i]=0; /* remember to reset marks */
  for(size_t i=0;i<user_region_mark_.size();i++)user_region_mark_[i]=0; /* remember to reset marks */
}

Range getOverlapped(const Range& a, const Range& b){
    if(a.first > b.first) return getOverlapped(b, a);
    if(b.first > a.second) return NULL_RANGE;
    return Range(b.first, (a.second > b.second)? b.second : a.second);
}

  void TaskManager :: CreateJob (const function<void(TaskInfo&)> & afunc,
                                 int antasks)
  {
    if (num_threads == 1 || !task_manager || func)
      {
        if (startup_function) (*startup_function)();
        
        TaskInfo ti;
        ti.ntasks = antasks;
        ti.thread_nr = 0; ti.nthreads = 1;
        // ti.node_nr = 0; ti.nnodes = 1;
        for (ti.task_nr = 0; ti.task_nr < antasks; ti.task_nr++)
          afunc(ti);

        if (cleanup_function) (*cleanup_function)();        
        return;
      }

    
    trace->StartJob(jobnr, afunc.target_type());

    func = &afunc;

    ntasks.store (antasks); // , memory_order_relaxed);
    ex = nullptr;


    nodedata[0]->start_cnt.store (0, memory_order_relaxed);

    jobnr++;
    
    for (int j = 0; j < num_nodes; j++)
      nodedata[j]->participate |= 1;

    if (startup_function) (*startup_function)();
    
    int thd = 0;
    int thds = GetNumThreads();
    int mynode = num_nodes * thd/thds;

    IntRange mytasks = Range(int(ntasks)).Split (mynode, num_nodes);
    NodeData & mynode_data = *(nodedata[mynode]);

    TaskInfo ti;
    ti.nthreads = thds;
    ti.thread_nr = thd;
    // ti.nnodes = num_nodes;
    // ti.node_nr = mynode;

    try
      {
        while (1)
          {
            int mytask = mynode_data.start_cnt++;
            if (mytask >= mytasks.Size()) break;
            
            ti.task_nr = mytasks.First()+mytask;
            ti.ntasks = ntasks;

            {
              RegionTracer t(ti.thread_nr, jobnr, RegionTracer::ID_JOB, ti.task_nr);
              (*func)(ti); 
            }
          }

      }
    catch (Exception e)
      {
        {
          lock_guard<mutex> guard(copyex_mutex);
          delete ex;
          ex = new Exception (e);
          mynode_data.start_cnt = mytasks.Size();
        }
      }

    if (cleanup_function) (*cleanup_function)();
    
    for (int j = 0; j < num_nodes; j++)
      if (workers_on_node[j])
        {
          while (complete[j] != jobnr)
            _mm_pause();
        }

    func = nullptr;
    if (ex)
      throw Exception (*ex);

    trace->StopJob();
  }

static void computeShapeByReshapeMask(const MatShape &srcShape,
                                      const MatShape &maskShape,
                                      Range srcRange /*= Range::all()*/,
                                      MatShape& dstShape)
{
    int srcShapeSize = (int)srcShape.size();
    int maskShapeSize = (int)maskShape.size();

    if (srcRange == Range::all())
        srcRange = Range(0, srcShapeSize);
    else
    {
        int sz = srcRange.size();
        srcRange.start = clamp(srcRange.start, srcShapeSize);
        srcRange.end = srcRange.end == INT_MAX ? srcShapeSize : srcRange.start + sz;
    }

    bool explicitMask = !maskShape.empty();  // All mask values are positive.
    for (int i = 0, n = maskShape.size(); i < n && explicitMask; ++i)
    {
        explicitMask = maskShape[i] > 0;
    }
    // Working range of source shape is a range where area(src) == area(mask).
    if (explicitMask)
    {
        int maskTotal = total(maskShape);
        // Go from the end of mask until we collect required total.
        bool matched = false;
        for (int i = srcRange.end - 1; i >= srcRange.start; --i)
        {
            if (matched)
            {
                if (i == 0 || total(srcShape, i, srcRange.end) != maskTotal)
                {
                    srcRange.start = i + 1;
                    break;
                }
            }
            else
            {
                matched = total(srcShape, i, srcRange.end) == maskTotal;
            }
        }
        CV_Assert(total(srcShape, srcRange.start, srcRange.end) == maskTotal);
    }

    CV_Assert(0 <= srcRange.start && srcRange.start <= srcRange.end && srcRange.end <= srcShapeSize);
    int dstShapeSize = srcShapeSize - srcRange.size() + maskShapeSize;
    dstShape.resize(dstShapeSize);

    std::copy(srcShape.begin(), srcShape.begin() + srcRange.start, dstShape.begin());
    std::copy(srcShape.begin() + srcRange.end, srcShape.begin() + srcShapeSize, dstShape.begin() + srcRange.start + maskShapeSize);

    int inferDim = -1;
    for (int i = 0; i < maskShapeSize; i++)
    {
        if (maskShape[i] > 0)
        {
            dstShape[srcRange.start + i] = maskShape[i];
        }
        else if (maskShape[i] == 0)
        {
            if (srcRange.start + i >= srcShapeSize)
                CV_Error(Error::StsBadArg, format("Copy dim[%d] (which has zero size) is out of the source shape bounds", srcRange.start + i));
            dstShape[srcRange.start + i] = srcShape[srcRange.start + i];
        }
        else if (maskShape[i] == -1)
        {
            if (inferDim != -1)
                CV_Error(Error::StsAssert, "Duplicate of inferred dim (which is denoted by -1)");
            inferDim = srcRange.start + i;
            dstShape[inferDim] = 1;
        }
        else
            CV_Error(Error::StsBadArg, "maskShape[i] >= -1");
    }

    size_t srcTotal = total(srcShape);
    size_t dstTotal = total(dstShape);

    if (inferDim != -1)
    {
        if (srcTotal % dstTotal != 0)
            CV_Error(Error::StsBackTrace, "Can't infer a dim denoted by -1");

        dstShape[inferDim] = (int)(srcTotal / dstTotal);
    }
    else
    {
        CV_Assert(srcTotal == dstTotal);
    }
}

  void TaskManager :: Loop(int thd)
  {
    /*
    static Timer tADD("add entry counter");
    static Timer tCASready1("spin-CAS ready tick1");
    static Timer tCASready2("spin-CAS ready tick2");
    static Timer tCASyield("spin-CAS yield");
    static Timer tCAS1("spin-CAS wait");
    static Timer texit("exit zone");
    static Timer tdec("decrement");
    */
    thread_id = thd;

    int thds = GetNumThreads();

    int mynode = num_nodes * thd/thds;

    NodeData & mynode_data = *(nodedata[mynode]);



    TaskInfo ti;
    ti.nthreads = thds;
    ti.thread_nr = thd;
    // ti.nnodes = num_nodes;
    // ti.node_nr = mynode;

      
#ifdef USE_NUMA
    numa_run_on_node (mynode);
#endif
    active_workers++;
    workers_on_node[mynode]++;
    int jobdone = 0;


#ifdef USE_MKL
    auto mkl_max = mkl_get_max_threads();
    mkl_set_num_threads_local(1);
#endif

    
    while (!done)
      {
        if (complete[mynode] > jobdone)
          jobdone = complete[mynode];

        if (jobnr == jobdone)
          {
            // RegionTracer t(ti.thread_nr, tCASyield, ti.task_nr);            
            if(sleep)
              this_thread::sleep_for(chrono::microseconds(sleep_usecs));
            else
              {
#ifdef WIN32
                this_thread::yield();
#else  // WIN32
                sched_yield();
#endif // WIN32
              }
            continue;
          }

        {
          // RegionTracer t(ti.thread_nr, tADD, ti.task_nr);

          // non-atomic fast check ...
          if ( (mynode_data.participate & 1) == 0) continue;

          int oldval = mynode_data.participate += 2;
          if ( (oldval & 1) == 0)
            { // job not active, going out again
              mynode_data.participate -= 2;
              continue;
            }
        }

        if (startup_function) (*startup_function)();
        
        IntRange mytasks = Range(int(ntasks)).Split (mynode, num_nodes);
          
        try
          {
            
            while (1)
              {
                if (mynode_data.start_cnt >= mytasks.Size()) break;
		int mytask = mynode_data.start_cnt.fetch_add(1, memory_order_relaxed);
                if (mytask >= mytasks.Size()) break;
                
                ti.task_nr = mytasks.First()+mytask;
                ti.ntasks = ntasks;
                
                {
                  RegionTracer t(ti.thread_nr, jobnr, RegionTracer::ID_JOB, ti.task_nr);
                  (*func)(ti);
                }
              }

          }
//.........这里部分代码省略.........

// Preprocess() - build a segment tree for O(log n) queries
void regProp2::compSeeds(void)
{
	Datareg2& reg2 = (Datareg2&)data;
	int i, j;
	int xdim, ydim;
	float val[4];
	Range* _prop_x, *prop_x;
	Range prop_y;
	Range propagated;
	Range c_prop;
	Range responsibility, c_respons;
	Range delay;
	Range y_comp;
	float min_x, min_y, max_x, max_y;
	float min_in, max_in, min4, max4;
	int nseed;
	xdim = reg2.dim[0];
	ydim = reg2.dim[1];
	_prop_x = new Range[ydim];
	// proceed through the slices computing seeds
	nseed=0;
	// process the k'th slab
	for(i=0; i<xdim-1; i++)
		for(j=0; j<ydim-1; j++)
		{
			prop_x = &_prop_x[j];
			// load the voxel data
			reg2.getCellValues(i, j, val);
			min_x = MIN2(val[0], val[3]);
			max_x = MAX2(val[0], val[3]);
			min_y = MIN2(val[0], val[1]);
			max_y = MAX2(val[0], val[1]);
			// set the incoming values if on a border
			if(i==0)
			{
				prop_x->Set(min_x, max_x);
			}
			if(j==0)
			{
				prop_y.Set(min_y, max_y);
			}
			// merge incoming information
			y_comp = prop_y.Complement(min_y, max_y);
			pro