C++ I1函数代码示例

void Postprocessing::initialise()
{
	// Read in the input maps
	fn_I1 = fn_in + "_half1_class001_unfil.mrc";
	fn_I2 = fn_in + "_half2_class001_unfil.mrc";

	if (verb > 0)
	{
		std::cout <<"== Reading input half-reconstructions: " <<std::endl;
		std::cout.width(35); std::cout << std::left <<"  + half1-map: "; std::cout << fn_I1 << std::endl;
		std::cout.width(35); std::cout << std::left <<"  + half2-map: "; std::cout << fn_I2 << std::endl;
	}

	I1.read(fn_I1);
	I2.read(fn_I2);
	I1().setXmippOrigin();
	I2().setXmippOrigin();

	if (!I1().sameShape(I2()))
	{
		std::cerr << " Size of half1 map: "; I1().printShape(std::cerr); std::cerr << std::endl;
		std::cerr << " Size of half2 map: "; I2().printShape(std::cerr); std::cerr << std::endl;
		REPORT_ERROR("Postprocessing::initialise ERROR: The two half reconstructions are not of the same size!");
	}

	if (do_auto_mask && fn_mask != "")
		REPORT_ERROR("Postprocessing::initialise ERROR: provide either --auto_mask OR --mask, but not both!");

	if (do_auto_bfac && ABS(adhoc_bfac) > 0.)
		REPORT_ERROR("Postprocessing::initialise ERROR: provide either --auto_bfac OR --adhoc_bfac, but not both!");
}

static void reset_locks(void)
{
	local_irq_disable();
	I1(A); I1(B); I1(C); I1(D);
	I1(X1); I1(X2); I1(Y1); I1(Y2); I1(Z1); I1(Z2);
	lockdep_reset();
	I2(A); I2(B); I2(C); I2(D);
	init_shared_classes();
	local_irq_enable();
}

void S2weaveCell::FindCellIndex(const P2& lptc)  
{
	// set the index for this cell 
	iu = FindCellParal(ps2w->ufibs, lptc.u); 
	iv = FindCellParal(ps2w->vfibs, lptc.v); 

	ASSERT(I1(ps2w->ufibs[iu - 1].wp, ps2w->ufibs[iu].wp).Contains(lptc.u)); 
	ASSERT(I1(ps2w->vfibs[iv - 1].wp, ps2w->vfibs[iv].wp).Contains(lptc.v)); 

	// rebuild the rest of the cell info
	ConstructCellBounds(); 
	CreateBoundList();
}

void S2weaveCellLinearCutTraverse::AdvanceToLamPos(double llamcp) 
{
	ASSERT(!bOnContour); 
	ASSERT(I1(lamcb, lamcf).Contains(lamcp)); 
	ASSERT(llamcp >= lamcb); 
	TOL_ZERO((ptcst + vbearing * lamcp - ptcp).Len()); 
	ASSERT((ibb == -1) || (lamcpbb > llamcp)); 

	lamcp = llamcp; 
	double clam = I1(lamcb, lamcf).InvAlong(lamcp); 
	ptcp = Along(clam, ptcb, ptcf); // or use an AlongAcc.  
	TOL_ZERO((ptcst + vbearing * lamcp - ptcp).Len()); 
}

void nemo_main ()
{
  setparams();	/* get cmdline stuff and compute x,y,u,v,etot,lz */

  optr = NULL;				/* make an orbit */
  allocate_orbit (&optr, 3, 1);
  Masso(optr)  = 1.0;                     /* and set Mass */
  Key(optr) = 0;
  Torb(optr,0) = tnow;
  Xorb(optr,0) = x;			/* .. positions */
  Yorb(optr,0) = y;
  Zorb(optr,0) = z;
  Uorb(optr,0) = u;			/* .. velocities */
  Vorb(optr,0) = v;
  Worb(optr,0) = w;
  I1(optr) = etot;			/*  energy (zero if not used) */
  I2(optr) = lz;                          /* angular momentum */
  
  dprintf(0,"pos: %f %f %f  \nvel: %f %f %f  \netot: %f\nlz=%f\n",
	  x,y,z,u,v,w,etot,lz);

  outstr = stropen (outfile,"w");		/* write to file */
  put_history(outstr);
  PotName(optr) = p.name;
  PotPars(optr) = p.pars;
  PotFile(optr) = p.file;
  write_orbit(outstr,optr);
  strclose(outstr);
}

void Solver::DiagramI1(){

    double time0, time1;
    time0 = omp_get_wtime();

    #pragma omp parallel
    {
        int id = omp_get_thread_num();
        int threads = omp_get_num_threads();

        for (int n=id; n<Npphh; n+=threads) blockspphh[n]->SetUpMatrices_I1(t0);
    }
    time1 = omp_get_wtime(); //cout << "Inside I1. SetUpMatrices needs: " << time1-time0 << " seconds" << endl;

    time0 = omp_get_wtime();
    for (int n=0; n<Npphh; n++){

        mat I1 = 0.5* blockspphh[n]->T * blockspphh[n]->I1 / blockspphh[n]->epsilon;

        for (double I=0; I<blockspphh[n]->Holes.n_rows; I++){
            for (double A=0; A<blockspphh[n]->Particles.n_rows; A++){

                int i = blockspphh[n]->Holes(I,0); int j = blockspphh[n]->Holes(I,1);
                int a = blockspphh[n]->Particles(A,0); int b = blockspphh[n]->Particles(A,1);

                t( Index(a,b,i,j) ) += I1(A,I);
                //StoreT(Index(a,b,i,j), I1(A,I));
            }
        }
    }
    time1 = omp_get_wtime(); //cout << "Inside I1. Calculate matrices needs: " << time1-time0 << " seconds" << endl;
}

Mat Assignment2::computeH(const vector<Point2f> &query, const vector<Point2f> &train)
{
	Mat I1(3,4,CV_64FC1,0.0);
	Mat I2(3,4,CV_64FC1,0.0);
	Mat H(3,3,CV_64FC1,0.0);
	
	for(unsigned int i=0; i<(query.size()); i++)
	{
		// Assign co-ordinates to the matrix
		I1.at<double>(0,i)=query[i].x;
		I1.at<double>(1,i)=query[i].y;
		I1.at<double>(2,i)=1;
		
		I2.at<double>(0,i)=train[i].x;
		I2.at<double>(1,i)=train[i].y;
		I2.at<double>(2,i)=1;
	}
	
	// solve linear equations
	solve(I1.t(), I2.t(), H, DECOMP_SVD);

	H = H.t();

	return H;
}

static void reset_locks(void)
{
	local_irq_disable();
	lockdep_free_key_range(&ww_lockdep.acquire_key, 1);
	lockdep_free_key_range(&ww_lockdep.mutex_key, 1);

	I1(A); I1(B); I1(C); I1(D);
	I1(X1); I1(X2); I1(Y1); I1(Y2); I1(Z1); I1(Z2);
	I_WW(t); I_WW(t2); I_WW(o.base); I_WW(o2.base); I_WW(o3.base);
	lockdep_reset();
	I2(A); I2(B); I2(C); I2(D);
	init_shared_classes();

	ww_mutex_init(&o, &ww_lockdep); ww_mutex_init(&o2, &ww_lockdep); ww_mutex_init(&o3, &ww_lockdep);
	memset(&t, 0, sizeof(t)); memset(&t2, 0, sizeof(t2));
	memset(&ww_lockdep.acquire_key, 0, sizeof(ww_lockdep.acquire_key));
	memset(&ww_lockdep.mutex_key, 0, sizeof(ww_lockdep.mutex_key));
	local_irq_enable();
}

double K1(double z){//1st Modified Bessel Function of the second kind
  z=fabs(z/2);
  if(z<=1)
    return log(z)*I1(2*z)+1./(2*z)*
      (1.+0.15443144*z*z-0.67278579*z*z*z*z-0.18156897*pow(z,6)-
       0.01919402*pow(z,8)-0.00110404*pow(z,10)-0.00004686*pow(z,12));
  else
    return exp(-2*z)/sqrt(2*z)*
      (1.25331414+0.23498619/z-0.0365562/z/z+0.01504268/z/z/z-
       0.00780353/z/z/z/z+0.00325614/pow(z,5)-0.00038245/pow(z,6));
}

int test_string_cast_vector()
{
	int Error = 0;
	
	glm::vec2 A1(1, 2);
	std::string A2 = glm::to_string(A1);
	Error += A2 != std::string("fvec2(1.000000, 2.000000)") ? 1 : 0;
	
	glm::vec3 B1(1, 2, 3);
	std::string B2 = glm::to_string(B1);
	Error += B2 != std::string("fvec3(1.000000, 2.000000, 3.000000)") ? 1 : 0;

	glm::vec4 C1(1, 2, 3, 4);
	std::string C2 = glm::to_string(C1);
	Error += C2 != std::string("fvec4(1.000000, 2.000000, 3.000000, 4.000000)") ? 1 : 0;
	
	glm::ivec2 D1(1, 2);
	std::string D2 = glm::to_string(D1);
	Error += D2 != std::string("ivec2(1, 2)") ? 1 : 0;
	
	glm::ivec3 E1(1, 2, 3);
	std::string E2 = glm::to_string(E1);
	Error += E2 != std::string("ivec3(1, 2, 3)") ? 1 : 0;
	
	glm::ivec4 F1(1, 2, 3, 4);
	std::string F2 = glm::to_string(F1);
	Error += F2 != std::string("ivec4(1, 2, 3, 4)") ? 1 : 0;
	
	glm::hvec2 G1(1, 2);
	std::string G2 = glm::to_string(G1);
	Error += G2 != std::string("hvec2(1.0000, 2.0000)") ? 1 : 0;
	
	glm::hvec3 H1(1, 2, 3);
	std::string H2 = glm::to_string(H1);
	Error += H2 != std::string("hvec3(1.0000, 2.0000, 3.0000)") ? 1 : 0;
	
	glm::hvec4 I1(1, 2, 3, 4);
	std::string I2 = glm::to_string(I1);
	Error += I2 != std::string("hvec4(1.0000, 2.0000, 3.0000, 4.0000)") ? 1 : 0;
	
	glm::dvec2 J1(1, 2);
	std::string J2 = glm::to_string(J1);
	Error += J2 != std::string("dvec2(1.000000, 2.000000)") ? 1 : 0;
	
	glm::dvec3 K1(1, 2, 3);
	std::string K2 = glm::to_string(K1);
	Error += K2 != std::string("dvec3(1.000000, 2.000000, 3.000000)") ? 1 : 0;
	
	glm::dvec4 L1(1, 2, 3, 4);
	std::string L2 = glm::to_string(L1);
	Error += L2 != std::string("dvec4(1.000000, 2.000000, 3.000000, 4.000000)") ? 1 : 0;
	
	return Error;
}

double IN(int n,double z){
  double i0=I0(z),i1=I1(z),in=0;
  int i;
  if(n==0)return i0;
  if(n==1)return i1;
  for(i=1;i<n;i++){
    in=i0-(2./z)*i*i1;
    i0=i1;
    i1=in;
  }
  return in;
}

/* decode binex mesaage 0x01-02: decoded glonass ephmemeris ------------------*/
static int decode_bnx_01_02(raw_t *raw, unsigned char *buff, int len)
{
    geph_t geph={0};
    unsigned char *p=buff;
    double tod,tof,tau_gps;
    int prn,day,leap;
    
    trace(4,"binex 0x01-02: len=%d\n",len);
    
    if (len>=119) {
        prn        =U1(p)+1;   p+=1;
        day        =U2(p);     p+=2;
        tod        =U4(p);     p+=4;
        geph.taun  =-R8(p);    p+=8;
        geph.gamn  =R8(p);     p+=8;
        tof        =U4(p);     p+=4;
        geph.pos[0]=R8(p)*1E3; p+=8;
        geph.vel[0]=R8(p)*1E3; p+=8;
        geph.acc[0]=R8(p)*1E3; p+=8;
        geph.pos[1]=R8(p)*1E3; p+=8;
        geph.vel[1]=R8(p)*1E3; p+=8;
        geph.acc[1]=R8(p)*1E3; p+=8;
        geph.pos[2]=R8(p)*1E3; p+=8;
        geph.vel[2]=R8(p)*1E3; p+=8;
        geph.acc[2]=R8(p)*1E3; p+=8;
        geph.svh   =U1(p)&0x1; p+=1;
        geph.frq   =I1(p);     p+=1;
        geph.age   =U1(p);     p+=1;
        leap       =U1(p);     p+=1;
        tau_gps    =R8(p);     p+=8;
        geph.dtaun =R8(p);
    }
    else {
        trace(2,"binex 0x01-02: length error len=%d\n",len);
        return -1;
    }
    if (!(geph.sat=satno(SYS_GLO,prn))) {
        trace(2,"binex 0x01-02: satellite error prn=%d\n",prn);
        return -1;
    }
    if (raw->time.time==0) return 0;
    geph.toe=utc2gpst(adjday(raw->time,tod-10800.0));
    geph.tof=utc2gpst(adjday(raw->time,tof-10800.0));
    geph.iode=(int)(fmod(tod+10800.0,86400.0)/900.0+0.5);
    
    if (!strstr(raw->opt,"-EPHALL")) {
        if (fabs(timediff(geph.toe,raw->nav.geph[prn-MINPRNGLO].toe))<1.0&&
            geph.svh==raw->nav.geph[prn-MINPRNGLO].svh) return 0; /* unchanged */
    }
    raw->nav.geph[prn-1]=geph;
    raw->ephsat=geph.sat;
    return 2;
}

int decode_navgpstime(tUbxRawData *raw)
{
   unsigned char *p = raw->buff + UBX_MSGSTART_SHIFT;
   tExternalGNSSGPSTimePtr gpstime = Msg_Get_ExternalGnssGpsTime();

   /* do not reset this struct due to week and leapsecond */
   gpstime->iTow = U4(p);
   gpstime->fTow = I4(p+4);
   gpstime->gpsWeek = (short)(I2(p + 8));
   gpstime->leapSec = (short)(I1(p + 10));

   return 0; 
}

void bench_env_klip(Pt2di sz)
{
     Im2D_U_INT1 I1(sz.x,sz.y);
     Im2D_U_INT1 I2(sz.x,sz.y);
     INT vmax = 30;




     ELISE_COPY
     (
         I1.all_pts(),
         Min
         (
             vmax,
                1
             +  frandr()*5
             +  unif_noise_4(3) * 30
         ),
         I1.out() 
     );
     ELISE_COPY
     (
         I1.border(1),
         0,
         I1.out()
     );
     ELISE_COPY
     (
         I1.all_pts(),
         EnvKLipshcitz_32(I1.in(0),vmax),
         I2.out() 
     );


     U_INT1 ** i1 = I1.data();
     U_INT1 ** i2 = I2.data();

     for (INT x=0; x<sz.x ; x++)
         for (INT y=0; y<sz.y ; y++)
            if (i1[y][x])
            {
                 INT v  = std::min3 
                      (
                          std::min3(i2[y+1][x-1]+3,i2[y+1][x]+2,i2[y+1][x+1]+3),
                          std::min3(i2[y][x-1]+2,(INT)i1[y][x],i2[y][x+1]+2),
                          std::min3(i2[y-1][x-1]+3,i2[y-1][x]+2,i2[y-1][x+1]+3)
                      );
                 BENCH_ASSERT(v==i2[y][x]);
           }
}

void Postprocessing::makeGuinierPlot(MultidimArray<Complex > &FT, std::vector<fit_point2D> &guinier)
{

	MultidimArray<int> radial_count(XSIZE(FT));
	MultidimArray<DOUBLE> lnF(XSIZE(FT));
	fit_point2D      onepoint;

	FOR_ALL_ELEMENTS_IN_FFTW_TRANSFORM(FT)
	{
    	int r2 = kp * kp + ip * ip + jp * jp;
    	int ires = ROUND(sqrt((DOUBLE)r2));
		if (ires < XSIZE(radial_count))
		{

	        lnF(ires) += abs(DIRECT_A3D_ELEM(FT, k, i, j));
	        radial_count(ires)++;
		}
	}

	DOUBLE xsize = XSIZE(I1());
	guinier.clear();
	FOR_ALL_ELEMENTS_IN_ARRAY1D(radial_count)
	{

		DOUBLE res = (xsize * angpix)/(DOUBLE)i; // resolution in Angstrom
		if (res >= angpix * 2.) // Apply B-factor sharpening until Nyquist, then low-pass filter later on (with a soft edge)
        {
            onepoint.x = 1. / (res * res);
            if (DIRECT_A1D_ELEM(lnF, i) > 0.)
            {
                onepoint.y = log ( DIRECT_A1D_ELEM(lnF, i) / DIRECT_A1D_ELEM(radial_count, i) );
                if (res <= fit_minres && res >= fit_maxres)
                {
                    onepoint.w = 1.;
                }
                else
                {
                    onepoint.w = 0.;
                }
            }
            else
            {
                onepoint.y = -99.;
                onepoint.w = 0.;
            }
            //std::cerr << " onepoint.x= " << onepoint.x << " onepoint.y= " << onepoint.y << " onepoint.w= " << onepoint.w << std::endl;
            guinier.push_back(onepoint);
        }
	}

}