C++ drand48函数代码示例

void gen_SSCA2_graph_MPI(graph_t* g)
{
    uint32_t local_TotVertices, TotVertices;
	uint32_t* clusterSizes;
	uint32_t* firstVsInCluster;
	uint32_t estTotClusters, local_totClusters;
	
	uint32_t *startVertex, *endVertex;
	uint32_t local_numEdges;
	weight_t* weights;
    weight_t MinWeight,MaxWeight;
    vertex_id_t u, v; 
    edge_id_t offset;
    uint32_t MaxCliqueSize;
    uint32_t MaxParallelEdges = 1;
    double ProbUnidirectional = 1.0;
    double ProbIntercliqueEdges = 0.1;
	uint32_t i_cluster, currCluster;
	uint32_t *startV, *endV, *d;
	uint32_t estNumEdges, edgeNum;

	uint32_t i, j, k, t, t1, t2, dsize;
	double p;
    int seed;
	uint32_t* permV;
    int size, rank, lgsize;
    

    g->directed = false;
    g->min_weight = 0;
    g->max_weight = 1;
    g->filename[0] = '\0';
    g->n = (vertex_id_t)1 << g->scale;
    if (strlen(g->filename) == 0) sprintf(g->filename, "ssca2-%d", g->scale);

    /*----------------------------------------------*/
	/*		initialize RNG 		*/
	/*----------------------------------------------*/
    seed = 2387;   
    srand48(seed);
   
    MinWeight = g->min_weight;
    MaxWeight = g->max_weight;
    TotVertices = g->n;
    size = g->nproc;
    g->local_n = g->n/size;
    local_TotVertices = g->local_n;
    rank = g->rank;
    for (lgsize = 0; lgsize < size; ++lgsize) {
        if ((1 << lgsize) == size) break;
    }
    
    MPI_Datatype MPI_VERTEX_ID_T;                                        
    MPI_Type_contiguous(sizeof(vertex_id_t), MPI_BYTE, &MPI_VERTEX_ID_T);
    MPI_Type_commit(&MPI_VERTEX_ID_T);                                   


	/*----------------------------------------------*/
	/*		generate local clusters		*/
	/*----------------------------------------------*/
    MaxCliqueSize = 49;

  	/* Estimate number of clusters required to make up 
	 * TotVertices and pad by 25% */
	estTotClusters = 1.25 * TotVertices / (size * MaxCliqueSize/2);
	
	/* Allocate a block of memory for this cluster-size array */
	clusterSizes = (uint32_t *) malloc(estTotClusters*sizeof(uint32_t));

	/* Generate random cluster sizes. */
    
	for(i = 0; i < estTotClusters; i++) {
		clusterSizes[i] = 1 + ( drand48() * MaxCliqueSize);
	}

	local_totClusters = 0;
  
	/* Allocate memory for storing the first vertex in each cluster */
	firstVsInCluster = (uint32_t *) malloc(estTotClusters*sizeof(uint32_t));

	/* Compute the first vertex in each cluster */
	firstVsInCluster[0] = 0;
	for (i=1; i<estTotClusters; i++) {
		firstVsInCluster[i] = firstVsInCluster[i-1] + clusterSizes[i-1];
		if (firstVsInCluster[i] > local_TotVertices-1)
			break;
	}

	local_totClusters = i;

	/* Fix the size of the last cluster */
	clusterSizes[local_totClusters-1] = local_TotVertices -\
		firstVsInCluster[local_totClusters-1];

	/*------------------------------------------------------*/
	/*		generate intra-cluster edges		*/
	/*------------------------------------------------------*/
 
	/* Roughly estimate the total number of edges */
	estNumEdges = (uint32_t) ((local_TotVertices * (double) MaxCliqueSize * (2-ProbUnidirectional)/2) +\
//.........这里部分代码省略.........

int main(int argc, char* argv[]) {
  if (argc != 4 ) {
    fprintf( stderr, "need csr-filename N reps!\n" );
    exit(-1);
  }

  char* l_csr_file;
  REALTYPE* l_a_sp; 
  unsigned int* l_rowptr;
  unsigned int* l_colidx;
  unsigned int l_rowcount, l_colcount, l_elements;
  
  REALTYPE* l_a_dense;
  REALTYPE* l_b;
  REALTYPE* l_c;
  REALTYPE* l_c_gold;
  REALTYPE* l_c_dense;
  REALTYPE l_max_error = 0.0;
  unsigned int l_m;
  unsigned int l_n;
  unsigned int l_k;

  unsigned int l_i;
  unsigned int l_j;
  unsigned int l_z;
  unsigned int l_elems;
  unsigned int l_reps;

  struct timeval l_start, l_end;
  double l_total;

  /* read sparse A */
  l_csr_file = argv[1];
  l_n = atoi(argv[2]);
  l_reps = atoi(argv[3]);
  if (my_csr_reader(  l_csr_file,
                 &l_rowptr,
                 &l_colidx,
                 &l_a_sp,
                 &l_rowcount, &l_colcount, &l_elements ) != 0 ) 
  {
    exit(-1);
  }
  l_m = l_rowcount;
  l_k = l_colcount;
  printf("CSR matrix data structure we just read:\n");
  printf("rows: %u, columns: %u, elements: %u\n", l_rowcount, l_colcount, l_elements);

  /* allocate dense matrices */
  l_a_dense = (REALTYPE*)_mm_malloc(l_k * l_m * sizeof(REALTYPE), 64);
  l_b = (REALTYPE*)_mm_malloc(l_k * l_n * sizeof(REALTYPE), 64);
  l_c = (REALTYPE*)_mm_malloc(l_m * l_n * sizeof(REALTYPE), 64);
  l_c_gold = (REALTYPE*)_mm_malloc(l_m * l_n * sizeof(REALTYPE), 64);
  l_c_dense = (REALTYPE*)_mm_malloc(l_m * l_n * sizeof(REALTYPE), 64);

  /* touch B */
  for ( l_i = 0; l_i < l_k*l_n; l_i++) {
    l_b[l_i] = (REALTYPE)drand48();
  }
  
  /* touch dense A */
  for ( l_i = 0; l_i < l_k*l_m; l_i++) {
    l_a_dense[l_i] = (REALTYPE)0.0;
  }
  /* init dense A using sparse A */
  for ( l_i = 0; l_i < l_m; l_i++ ) {
    l_elems = l_rowptr[l_i+1] - l_rowptr[l_i];
    for ( l_z = 0; l_z < l_elems; l_z++ ) {
      l_a_dense[(l_i*l_k)+l_colidx[l_rowptr[l_i]+l_z]] = l_a_sp[l_rowptr[l_i]+l_z];
    }
  }

  /* touch C */
  for ( l_i = 0; l_i < l_m*l_n; l_i++) {
    l_c[l_i] = (REALTYPE)0.0;
    l_c_gold[l_i] = (REALTYPE)0.0;
    l_c_dense[l_i] = (REALTYPE)0.0;
  }
  
  /* compute golden results */
  printf("computing golden solution...\n");
  for ( l_j = 0; l_j < l_n; l_j++ ) {
    for (l_i = 0; l_i < l_m; l_i++ ) {
      l_elems = l_rowptr[l_i+1] - l_rowptr[l_i];
      for (l_z = 0; l_z < l_elems; l_z++) {
        l_c_gold[(l_n*l_i) + l_j] +=  l_a_sp[l_rowptr[l_i]+l_z] * l_b[(l_n*l_colidx[l_rowptr[l_i]+l_z])+l_j];
      }
    }
  }
  printf("...done!\n");

  /* libxsmm generated code */
  printf("computing libxsmm (A sparse) solution...\n");
  libxsmm_code(l_b, l_c);
  printf("...done!\n");

  /* BLAS code */
  printf("computing BLAS (A dense) solution...\n");
  double alpha = 1.0;
  double beta = 1.0;
//.........这里部分代码省略.........

/*determines if star with prob p should be separrated into stream*/
int prob_ok(int n, double* p)
{
    int ok;
    double r;
    double step1, step2, step3;

    r = drand48();

    switch (n)
    {
    case 1:
        if ( r > p[0] )
        {
            ok = 0;
        }
        else
        {
            ok = 1;
        }
        break;
    case 2:
        step1 = p[0] + p[1];
        if ( r > step1 )
        {
            ok = 0;
        }
        else if ( (r < p[0]) )
        {
            ok = 1;
        }
        else if ( (r > p[0]) && (r <= step1) )
        {
            ok = 2;
        }
        break;
    case 3:
        step1 = p[0] + p[1];
        step2 = p[0] + p[1] + p[2];
        if ( r > step2 )
        {
            ok = 0;
        }
        else if ( (r < p[0]) )
        {
            ok = 1;
        }
        else if ( (r > p[0]) && (r <= step1) )
        {
            ok = 2;
        }
        else if ( (r > step1) && (r <= step2) )
        {
            ok = 3;
        }
        break;
    case 4:
        step1 = p[0] + p[1];
        step2 = p[0] + p[1] + p[2];
        step3 = p[0] + p[1] + p[2] + p[3];
        if ( r > step3 )
        {
            ok = 0;
        }
        else if ( (r <= p[0]) )
        {
            ok = 1;
        }
        else if ( (r > p[0]) && (r <= step1) )
        {
            ok = 2;
        }
        else if ( (r > step1) && (r <= step2) )
        {
            ok = 3;
        }
        else if ( (r > step2) && (r <= step3) )
        {
            ok = 4;
        }
        break;
    default:
        fprintf(stderr,
                "ERROR:  Too many streams to separate using current code; "
                "please update the switch statement in probability.c->prob_ok to handle %d streams", n);
        exit(EXIT_SUCCESS);
    }
    return ok;
}

void *
fail_prone_malloc(size_t size)
{
   return drand48() < ALLOC_ERR_PROB ? NULL : __libc_malloc(size);
}

void *
fail_prone_calloc(size_t nitems, size_t size)
{
   return drand48() < ALLOC_ERR_PROB ? NULL : __libc_calloc(nitems, size);
}

int main()
{
  int i, j, k, d = 10,d2=5;
  float * a = fvec_new (d * d);
  float * b = fvec_new (d * d);
  float * b0 = fvec_new (d * d);
#define B0(i,j) b0[(i)+(j)*d]
#define A(i,j) a[(i)+(j)*d]
#define B(i,j) b[(i)+(j)*d] 
  float * lambda = fvec_new (d);
  float * v = fvec_new (d * d);
  float *v_part=fvec_new (d * d2);

  for (i = 0 ; i < d ; i++)
    for (j = 0 ; j  <= i ; j++) {
      A(i,j) = A(j,i) = drand48(); 
      B0(i,j)=drand48();
      B0(j,i)=drand48();
      
/*      B(i,j) = B(j,i) = drand48(); */
    }
  /* make a positive definite b (with b=b0*b0') */
  for (i = 0 ; i < d ; i++)
    for (j = 0 ; j  < d ; j++) {
      double accu=0;
      for(k=0;k<d;k++) 
        accu+=B0(i,k)*B0(j,k);
      B(i,j)=accu;
    }

  printf ("a = ");
  fmat_print(a,d,d);

  printf ("\nb = "); 
  fmat_print(b,d,d);

  printf ("Solution of the eigenproblem Av=lambda v\n");
  
  printf ("\n");
  int ret=eigs_sym (d, a, lambda, v);
  assert(ret==0);
  printf ("\n");


  printf("Eigenvectors:\n");
  fmat_print(v,d,d);

  fprintf(stdout, "lambda = ");
  fvec_print (lambda, d);
  printf ("\n");

  printf("Partial eigenvalues/vectors:\n");

  printf ("\n");
  ret=eigs_sym_part (d, a, d2, lambda, v_part);
  assert(ret>0);
  if(ret<d2) 
    printf("!!! only %d / %d eigenvalues converged\n",ret,d2);

  printf ("\n");

  printf("Eigenvectors:\n");
  fmat_print(v_part,d,d2);
  
  fprintf(stdout, "lambda = ");
  fvec_print (lambda, d2);
  printf ("\n");

   
  printf ("Solution of the generalized eigenproblem Av=lambda B v\n");

  printf ("\n");
  ret=geigs_sym (d, a, b, lambda, v);
  assert(ret==0);
  printf ("\n");

  fmat_print(v,d,d);

  fprintf(stdout, "lambda = ");
  fvec_print (lambda, d);
  printf ("\n");

  free (a);
  free (lambda);
  free (v);

  return 0;
}

int doAllFeatures()
{
	/* Initial biases */
	{
		int u,m,f;
		
		for(u=0;u<NUSERS;u++) {
			for(f=0;f<NFEATURES;f++)
			    bU[u][f]=drand48()*0.01-0.005;
		}
		for(m=0;m<NMOVIES;m++) {
			for(f=0;f<NFEATURES;f++)
			    bV[m][f]=drand48()*0.01-0.005;
		}
	}
	
	
	/* Initial estimation for current feature */
	{
		int u,m,f;
		
		for(u=0;u<NUSERS;u++) {
			for(f=0;f<NFEATURES;f++)
			    sU[u][f]=drand48()*0.1-0.04;
		}
		for(m=0;m<NMOVIES;m++) {
			for(f=0;f<NFEATURES;f++) {
			    sV[m][f]=drand48()*0.05-0.025;
			    sY[m][f]=drand48()*0.02-0.01;
			}
		}
	}
	
	/* Optimize current feature */
	double nrmse=2., last_rmse=10.;
	double thr=sqrt(1.-E);
	int loopcount=0;
	    //thr=sqrt(1.-E2);
	double Gamma2 = G2;
	double Gamma0 = G0;
	while( ( nrmse < (last_rmse-E) ) || loopcount++ < 20) {
		last_rmse=nrmse;
		clock_t t0=clock();

		int u,m, f;
		for(u=0;u<NUSERS;u++) {

			// Calculate sumY and NuSY for each factor
			double sumY[NFEATURES];
			ZERO(sumY);
			double lNuSY[NFEATURES];
			ZERO(lNuSY);
			int base0=useridx[u][0];
			int d0=UNTRAIN(u);
			int j;
			int f;
			int dall=UNALL(u);
			double NuS = 1.0/sqrt(dall);
			for(j=0;j<d0;j++) {
				int mm=userent[base0+j]&USER_MOVIEMASK;
				for(f=0;f<NFEATURES;f++)
					sumY[f]+=sY[mm][f];
			}
//if ( loopcount > 1 ) {
//printf("sumY: %f\n", sumY);
//fflush(stdout);
//}
			for(j=0;j<d0;j++) {
				int mm=userent[base0+j]&USER_MOVIEMASK;
				for(f=0;f<NFEATURES;f++) {
					lNuSY[f] = NuS * sumY[f]; 
//if ( loopcount > 1 ) {
//printf("lNuSY: %f\n", lNuSY[f]);
//fflush(stdout);
//}
				}
			}

			double ycontrib[NFEATURES];
			ZERO(ycontrib);

			// For all rated movies
			double bdampen = d0/1.1;
			for(j=0;j<d0;j++) {
				int m=userent[base0+j]&USER_MOVIEMASK;

				// Figure out the current error
				double ee=err[base0+j];
				double e2 = ee;
				for (f=0; f<NFEATURES; f++) {
				    e2 -= (bU[u][f] + bV[m][f]);
					e2 -= ((sU[u][f]+lNuSY[f])*sV[m][f]);
				}


				// update U V and slope component of Y
				//double yfactor = NuS/sqrt(moviecount[m]); 
				//double yfactor = NuS;
				double yfactor = NuS/d0;
				for (f=0; f<NFEATURES; f++) {
//.........这里部分代码省略.........

void *elCamino(float *x_obs, float *x_aComparar,float *y_ayudaAx, float *xDelMomento,float *xCandidato, float *a, float *b, float *c, float *d,float *l,float paso, int t, int iteraciones, int burn)
{
  int i;
  float aPrima;
  float bPrima;
  float cPrima;
  float dPrima;
  float lPresente;
  float lCandidato;
  float gamma;
  float beta;
  float alpha;
  float x_0 = 15 ;
  float y_0 = 13; 

  const gsl_rng_type * T;
  gsl_rng * r;
  gsl_rng_env_setup();
  T = gsl_rng_default;
  r = gsl_rng_alloc (T);

  for(i=0; i < (iteraciones-1) ; i++)
    {
      aPrima = gsl_ran_gaussian(r, 0.1)+a[i];
      //bPrima = gsl_ran_gaussian(r, 0.1)+b[i];
      //cPrima = gsl_ran_gaussian(r, 0.1)+c[i];
      //dPrima = gsl_ran_gaussian(r, 0.1)+d[i];

      xDelMomento = RungeKutta(x_aComparar, y_ayudaAx, x_0, y_0, a[i], b[i], c[i], d[i], t, paso);
      xCandidato = RungeKutta(x_aComparar, y_ayudaAx, x_0, y_0, aPrima, b[i], c[i], d[i], t, paso);

      lPresente = likelyhood(x_obs, xDelMomento, t);
      lCandidato = likelyhood(x_obs, xCandidato , t);
 
      gamma = (lCandidato - lPresente);
      if(gamma>=0.00)
	{
	  a[i+1]=aPrima;
	}
      else
	{
	  beta = drand48();
	  alpha = exp(gamma);
	  if(beta<=alpha)
	    {
	      a[i+1]=aPrima;
	    }
	  else
	    {
	      a[i+1]=a[i];
	    }
	}
      //aPrima = gsl_ran_gaussian(r, 0.1)+a[i];
      bPrima = gsl_ran_gaussian(r, 0.1)+b[i];
      //cPrima = gsl_ran_gaussian(r, 0.1)+c[i];
      //dPrima = gsl_ran_gaussian(r, 0.1)+d[i];

      xDelMomento = RungeKutta(x_aComparar, y_ayudaAx, x_0, y_0, a[i], b[i], c[i], d[i], t, paso);
      xCandidato = RungeKutta(x_aComparar, y_ayudaAx, x_0, y_0, a[i], bPrima, c[i], d[i], t, paso);

      lPresente = likelyhood(x_obs, xDelMomento, t);
      lCandidato = likelyhood(x_obs, xCandidato , t);

      gamma = (lCandidato - lPresente);
      if(gamma>=0.00)
	{
	  b[i+1]=bPrima;
	}
      else
	{
	  beta = drand48();
	  alpha = exp(gamma);
	  if(beta<=alpha)
	    {
	      b[i+1]=bPrima;
	    }
	  else
	    {
	      b[i+1]=b[i];
	    }
	}

      //aPrima = gsl_ran_gaussian(r, 0.1)+a[i];
      //bPrima = gsl_ran_gaussian(r, 0.1)+b[i];
      cPrima = gsl_ran_gaussian(r, 0.1)+c[i];
      //dPrima = gsl_ran_gaussian(r, 0.1)+d[i];


      xDelMomento = RungeKutta(x_aComparar, y_ayudaAx, x_0, y_0, a[i], b[i], c[i], d[i], t, paso);
      xCandidato = RungeKutta(x_aComparar, y_ayudaAx, x_0, y_0, a[i], b[i], cPrima, d[i], t, paso);

      lPresente = likelyhood(x_obs, xDelMomento, t);
      lCandidato = likelyhood(x_obs, xCandidato , t);

      gamma = (lCandidato - lPresente);
      if(gamma>=0.00)
	{
	  c[i+1]=cPrima;
	}
      else
//.........这里部分代码省略.........

//.........这里部分代码省略.........
		g_latest = g;
		g_size_latest = gsize;
		g_count_latest = count;
		fflush(stdout);
	}

	while(gsize < 206)
	{
		/*
		 * if we have a counter example
		 */
		if(count == 0)
		{
			printf("Eureka!  Counter-example found!\n");
			PrintGraph(g,gsize);
			fflush(stdout);
			/*
			 * make a new graph one size bigger
			 */
			new_g = (int *)malloc((gsize+1)*(gsize+1)*sizeof(int));
			if(new_g == NULL)
				exit(1);
			/*
			 * copy the old graph into the new graph leaving the
			 * last row and last column alone
			 */
			CopyGraph(g,gsize,new_g,gsize+1);

			/*
			 * zero out the last column and last row
			 */
			for(i=0; i < (gsize+1); i++)
			{
				if(drand48() > 0.5) {
					new_g[i*(gsize+1) + gsize] = 0; // last column
					new_g[gsize*(gsize+1) + i] = 0; // last row
				}
				else
				{
					new_g[i*(gsize+1) + gsize] = 1; // last column
					new_g[gsize*(gsize+1) + i] = 1; // last row
				}
			}

			/*
			 * throw away the old graph and make new one the
			 * graph
			 */
			free(g);
			g = new_g;
			gsize = gsize+1;

			/*
			 * reset the taboo list for the new graph
			 */
			taboo_list = FIFOResetEdge(taboo_list);

			/*
			 * keep going
			 */
			//continue;
			break; //Go to first while loop
		}

		/*
		 * otherwise, we need to consider flipping an edge

// Test node data sorting
void SortNodeData() {
  int NNodes = 10000;
  int NEdges = 100000;

  TPt <TNodeEDatNet<TInt, TInt> > Net;
  TPt <TNodeEDatNet<TInt, TInt> > Net1;
  TPt <TNodeEDatNet<TInt, TInt> > Net2;
  int i;
  int n;
  int NCount;
  int x,y;
  bool t;
  int NodeId;
  int NodeDat;
  bool ok;
  bool Sorted;
  int Min;
  int Value;

  Net = TNodeEDatNet<TInt, TInt>::New();
  t = Net->Empty();

  // create the nodes
  for (i = 0; i < NNodes; i++) {
    Net->AddNode((i*13) % NNodes);
  }
  t = Net->Empty();
  n = Net->GetNodes();

  // create random edges
  NCount = NEdges;
  while (NCount > 0) {
    x = (long) (drand48() * NNodes);
    y = (long) (drand48() * NNodes);
    // Net->GetEdges() is not correct for the loops (x == y),
    // skip the loops in this test
    if (x != y  &&  !Net->IsEdge(x,y)) {
      n = Net->AddEdge(x, y);
      NCount--;
    }
  }
  PrintNStats("SortNodeData:Net", Net);

  // add data to nodes, square of node ID % NNodes
  for (TNodeEDatNet<TInt, TInt>::TNodeI NI = Net->BegNI(); NI < Net->EndNI(); NI++) {
    NodeId = NI.GetId();
    NodeDat = (NI.GetId()*NI.GetId()) % NNodes;
    Net->SetNDat(NodeId, NodeDat);
  }

  // test node data
  ok = true;
  for (TNodeEDatNet<TInt, TInt>::TNodeI NI = Net->BegNI(); NI < Net->EndNI(); NI++) {
    NodeDat = Net->GetNDat(NI.GetId());
    Value = (NI.GetId()*NI.GetId()) % NNodes;
    if (NodeDat != Value) {
      ok = false;
    }
  }
  printf("network SortNodeData:Net, status1 %s\n", (ok == true) ? "ok" : "ERROR");

  // test sorting of node IDs (unsorted)
  Min = -1;
  Sorted = true;
  for (TNodeEDatNet<TInt, TInt>::TNodeI NI = Net->BegNI(); NI < Net->EndNI(); NI++) {
    Value = NI.GetId();
    if (Min > Value) {
      Sorted = false;
    }
    Min = Value;
  }
  printf("network SortNodeData:Net, status2 %s\n", (Sorted == false) ? "ok" : "ERROR");

  // sort the nodes by node IDs (sorted)
  Net->SortNIdById();

  // test sorting of node IDs
  Min = -1;
  Sorted = true;
  for (TNodeEDatNet<TInt, TInt>::TNodeI NI = Net->BegNI(); NI < Net->EndNI(); NI++) {
    Value = NI.GetId();
    if (Min > Value) {
      Sorted = false;
    }
    Min = Value;
  }
  printf("network SortNodeData:Net, status3 %s\n", (Sorted == true) ? "ok" : "ERROR");

  // test sorting of node data (unsorted)
  Min = -1;
  Sorted = true;
  for (TNodeEDatNet<TInt, TInt>::TNodeI NI = Net->BegNI(); NI < Net->EndNI(); NI++) {
    Value = Net->GetNDat(NI.GetId());
    if (Min > Value) {
      Sorted = false;
    }
    Min = Value;
  }
  printf("network SortNodeData:Net, status4 %s\n", (Sorted == false) ? "ok" : "ERROR");
//.........这里部分代码省略.........

int main(void)
{
    // seed pseudorandom number generator
    srand48(time(NULL));

    // instantiate window
    GWindow window = newGWindow(WIDTH, HEIGHT);

    // instantiate bricks
    initBricks(window);

    // instantiate ball, centered in middle of window
    GOval ball = initBall(window);

    // instantiate paddle, centered at bottom of window
    GRect paddle = initPaddle(window);

    // instantiate scoreboard, centered in middle of window, just above ball
    GLabel label = initScoreboard(window);

    // number of bricks initially
    int bricks = COLS * ROWS;

    // number of lives initially
    int lives = LIVES;

    // number of points initially
    int points = 0;
    
    // velocity of ball
    // TODO use drand for velocityX
    double velocityX = drand48() * 3.0;
    double velocityY = 3.0;

    // keep playing until game over
    waitForClick();    
    while (lives > 0 && bricks > 0)
    {
        // check for mouse event
        GEvent event = getNextEvent(MOUSE_EVENT);


        updateScoreboard(window, label, points);
        // set velocity of ball
        move(ball, velocityX, velocityY);

        // detect collision
        GObject collision = detectCollision(window, ball);

        if (getX(ball) + getWidth(ball) >= getWidth(window))
        {
            velocityX = -velocityX;
        }
        // bounce off left edge of window
        else if (getX(ball) <= 0)
        {
            velocityX = -velocityX;
        }
        else if (getY(ball) + getHeight(ball) >= getHeight(window))
        {
            lives -= 1;
            waitForClick();
            setLocation(ball, 190, 290);
            move(ball, velocityX, -velocityY);
        }
        else if (getY(ball) <= 0)
        {
            velocityY = -velocityY;
        }
        else if (collision != NULL)
        {
            if (collision == paddle)
            {
                velocityY = -velocityY;
            }
            else if (strcmp(getType(collision), "GRect") == 0)
            {
                // TODO
                velocityY = -velocityY;
                removeGWindow(window, collision);
                points += 1;
                bricks -= 1;
            }
        }
        pause(10);
        


        if (event != NULL)
        {
                if (getEventType(event) == MOUSE_MOVED)
                {
                    // ensure circle follows top cursor
                    double x = getX(event) - getWidth(paddle) / 2;
                    double y = 525;
                    setLocation(paddle, x, y);
                }
        }

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

// Test update edge data
void UpdateEdgeData() {
  int NNodes = 10000;
  int NEdges = 100000;

  TPt <TNodeEDatNet<TInt, TInt> > Net;
  TPt <TNodeEDatNet<TInt, TInt> > Net1;
  TPt <TNodeEDatNet<TInt, TInt> > Net2;
  int i;
  int n;
  int NCount;
  int x,y;
  bool t;
  int SrcNId;
  int DstNId;
  int EdgeDat;
  int Value;
  bool ok;

  Net = TNodeEDatNet<TInt, TInt>::New();
  t = Net->Empty();

  // create the nodes
  for (i = 0; i < NNodes; i++) {
    Net->AddNode(i);
  }
  t = Net->Empty();
  n = Net->GetNodes();

  // create random edges and edge data x+y+10
  NCount = NEdges;
  while (NCount > 0) {
    x = (long) (drand48() * NNodes);
    y = (long) (drand48() * NNodes);
    // Net->GetEdges() is not correct for the loops (x == y),
    // skip the loops in this test
    if (x != y  &&  !Net->IsEdge(x,y)) {
      n = Net->AddEdge(x, y, x+y+10);
      NCount--;
    }
  }
  PrintNStats("UpdateEdgeData:Net", Net);

  // verify edge data, x+y+10
  ok = true;
  for (TNodeEDatNet<TInt, TInt>::TEdgeI EI = Net->BegEI(); EI < Net->EndEI(); EI++) {
    SrcNId = EI.GetSrcNId();
    DstNId = EI.GetDstNId();
    EdgeDat = Net->GetEDat(SrcNId, DstNId);
    Value = SrcNId+DstNId+10;
    if (EdgeDat != Value) {
      ok = false;
    }
  }
  printf("network UpdateEdgeData:Net, status1 %s\n", (ok == true) ? "ok" : "ERROR");

  // update edge data, x+y+5
  for (TNodeEDatNet<TInt, TInt>::TEdgeI EI = Net->BegEI(); EI < Net->EndEI(); EI++) {
    Net->SetEDat(EI.GetSrcNId(),EI.GetDstNId(),EI.GetSrcNId()+EI.GetDstNId()+5);
  }

  // verify edge data, x+y+5
  ok = true;
  for (TNodeEDatNet<TInt, TInt>::TEdgeI EI = Net->BegEI(); EI < Net->EndEI(); EI++) {
    SrcNId = EI.GetSrcNId();
    DstNId = EI.GetDstNId();
    EdgeDat = Net->GetEDat(SrcNId, DstNId);
    Value = SrcNId+DstNId+5;
    if (EdgeDat != Value) {
      ok = false;
    }
  }
  printf("network UpdateEdgeData:Net, status2 %s\n", (ok == true) ? "ok" : "ERROR");
}

// Test node, edge creation
void ManipulateNodesEdges() {
  int NNodes = 10000;
  int NEdges = 100000;
  const char *FName = "demo.net.dat";

  TPt <TNodeEDatNet<TInt, TInt> > Net;
  TPt <TNodeEDatNet<TInt, TInt> > Net1;
  TPt <TNodeEDatNet<TInt, TInt> > Net2;
  int i;
  int n;
  int NCount;
  int ECount1;
  int ECount2;
  int x,y;
  bool t;

  Net = TNodeEDatNet<TInt, TInt>::New();
  t = Net->Empty();

  // create the nodes
  for (i = 0; i < NNodes; i++) {
    Net->AddNode(i);
  }
  t = Net->Empty();
  n = Net->GetNodes();

  // create random edges
  NCount = NEdges;
  while (NCount > 0) {
    x = (long) (drand48() * NNodes);
    y = (long) (drand48() * NNodes);
    // Net->GetEdges() is not correct for the loops (x == y),
    // skip the loops in this test
    if (x != y  &&  !Net->IsEdge(x,y)) {
      n = Net->AddEdge(x, y);
      NCount--;
    }
  }
  PrintNStats("ManipulateNodesEdges:Net", Net);

  // get all the nodes
  NCount = 0;
  for (TNodeEDatNet<TInt, TInt>::TNodeI NI = Net->BegNI(); NI < Net->EndNI(); NI++) {
    NCount++;
  }

  // get all the edges for all the nodes
  ECount1 = 0;
  for (TNodeEDatNet<TInt, TInt>::TNodeI NI = Net->BegNI(); NI < Net->EndNI(); NI++) {
    for (int e = 0; e < NI.GetOutDeg(); e++) {
      ECount1++;
    }
  }

  // get all the edges directly
  ECount2 = 0;
  for (TNodeEDatNet<TInt, TInt>::TEdgeI EI = Net->BegEI(); EI < Net->EndEI(); EI++) {
    ECount2++;
  }
  printf("network ManipulateNodesEdges:Net, nodes %d, edges1 %d, edges2 %d\n",
      NCount, ECount1, ECount2);

  // assignment
  Net1 = TNodeEDatNet<TInt, TInt>::New();
  *Net1 = *Net;
  PrintNStats("ManipulateNodesEdges:Net1",Net1);

  // save the network
  {
    TFOut FOut(FName);
    Net->Save(FOut);
    FOut.Flush();
  }

  // load the network
  {
    TFIn FIn(FName);
    Net2 = TNodeEDatNet<TInt, TInt>::Load(FIn);
  }
  PrintNStats("ManipulateNodesEdges:Net2",Net2);

  // remove all the nodes and edges
  for (i = 0; i < NNodes; i++) {
    n = Net->GetRndNId();
    Net->DelNode(n);
  }
  PrintNStats("ManipulateNodesEdges:Net",Net);

  Net1->Clr();
  PrintNStats("ManipulateNodesEdges:Net1",Net1);
}

// Test update node data
void UpdateNodeData() {
  int NNodes = 10000;
  int NEdges = 100000;

  TPt <TNodeEDatNet<TInt, TInt> > Net;
  TPt <TNodeEDatNet<TInt, TInt> > Net1;
  TPt <TNodeEDatNet<TInt, TInt> > Net2;
  int i;
  int n;
  int NCount;
  int x,y;
  bool t;
  int NodeDat;
  int Value;
  bool ok;

  Net = TNodeEDatNet<TInt, TInt>::New();
  t = Net->Empty();

  // create the nodes
  for (i = 0; i < NNodes; i++) {
    Net->AddNode(i,i+5);
  }
  t = Net->Empty();
  n = Net->GetNodes();

  // create random edges
  NCount = NEdges;
  while (NCount > 0) {
    x = (long) (drand48() * NNodes);
    y = (long) (drand48() * NNodes);
    // Net->GetEdges() is not correct for the loops (x == y),
    // skip the loops in this test
    if (x != y  &&  !Net->IsEdge(x,y)) {
      n = Net->AddEdge(x, y);
      NCount--;
    }
  }
  PrintNStats("UpdateNodeData:Net", Net);

  // read and test node data
  ok = true;
  for (TNodeEDatNet<TInt, TInt>::TNodeI NI = Net->BegNI(); NI < Net->EndNI(); NI++) {
    NodeDat = Net->GetNDat(NI.GetId());
    Value = NI.GetId()+5;
    if (NodeDat != Value) {
      ok = false;
    }
  }
  printf("network UpdateNodeData:Net, status1 %s\n", (ok == true) ? "ok" : "ERROR");

  // update node data, node ID + 10
  for (TNodeEDatNet<TInt, TInt>::TNodeI NI = Net->BegNI(); NI < Net->EndNI(); NI++) {
    Net->SetNDat(NI.GetId(), NI.GetId()+10);
  }

  // read and test node data
  ok = true;
  for (TNodeEDatNet<TInt, TInt>::TNodeI NI = Net->BegNI(); NI < Net->EndNI(); NI++) {
    NodeDat = Net->GetNDat(NI.GetId());
    Value = NI.GetId()+10;
    if (NodeDat != Value) {
      ok = false;
    }
  }
  printf("network UpdateNodeData:Net, status2 %s\n", (ok == true) ? "ok" : "ERROR");
}

float Lowpass::main_body_quality(
  Streamline *st,
  float radius,
  int nsamples,
  float dist,
  Window2d *win,
  int &which_side
)
{
  int i;
  float rad;
  float delta;
  float x,y;
  float sum = 0;
  int count = 0;
  VectorField *vf = st->vf;

  /* kluge to work around problems at the image borders */

  float xmin = 2.0 / xsize;
  float ymin = 2.0 / ysize;
  float xmax = 1 - 2.0 / xsize;
  float ymax = image->getaspect() - 2.0 / ysize;

  /* measure the quality at several random positions along the streamline */

  for (i = 0; i < nsamples; i++) {

    int index = (int) (st->samples * drand48());
    SamplePoint *sample = &st->pts[index];

    x = sample->x;
    y = sample->y;

    /* don't measure near the borders */
    if (x < xmin || x > xmax || y < ymin || y > ymax)
      continue;

    rad = radius * rad_image->get_value(x,y) / xsize;
    delta = 0.3 * rad / nsamples;

    /* find out how to move perpendicular to the vector field */
    float dx = - delta * vf->yval(x,y);
    float dy =   delta * vf->xval(x,y);

    /* sample along either side of the streamline */
    float x1 = x + dx;
    float y1 = y + dy;
    float x2 = x - dx;
    float y2 = y - dy;

    /* measure the quality at the two points */
    float value1 = image->get_value(x1,y1);
    float value2 = image->get_value(x2,y2);
   // printf("value1: %f, value2: %f\n",value1,value2);

    sum += value1 - value2;

    count++;

  }

  /* return quality */

  if (count == 0)
    return (0.0);
  else {

    if (sum > 0)
      which_side = LEFT;
    else
      which_side = RIGHT;

    return (fabs (sum / count));
  }
}