main.c

一类球型分布问题的实现

/***************************************************************************
 Hyperspherical distribution
                          main.c  -  description
                             -------------------
    begin                : 2003-07-04
    last update      : 2003-07-04
    copyright            : (C) 2003 by 
    email                : [e-mail]
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
 /*
 It generates an uniform distribution inside a D-sphere (radius 1, center: origin).

 The method itself is proved. However, there may be some bugs.
 
*/
                                                 
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <syscall.h>

#define Max_D 			100 // Max number of dimensions

// Structures
struct vector		{int size;double x[Max_D];};

// Subroutines
double 			alea_float(double min, double max);
double  alea_normal(double mean,double std_dev);
void display_vector(struct vector x);
void  save_vector(FILE *f_run,struct vector x);

// Global variables
double   pi;

//=============================================================================== 
int main(int argc, char *argv[])
{
 FILE *f_run;

 int                 D;
 int                 i;
 int                 j;
 double           length;
 int                 N;
 double           pw;
 double           r;
 struct  vector   x;

 pi=2*acos(0);

  f_run =fopen("run.txt","w");  // File to save result
   
 
 //---------------------------------------   Interactive input                           
 input:
  printf("\n Dimension?  (0=stop):  ");
  scanf("%i",&D);
  if(D==0) goto end;
  if (D>Max_D )
  {
     printf("\n No more than %i dimensions, please.",Max_D);
     goto input;
  }
  
 printf("\n Number of points to generate? (0=stop): ");
 scanf("%i",&N);
 if (N==0) goto end;

 x.size=D;
 pw=1/(double)D;
 //------------------------------------------------------------------------------------------- LOOP on points to generate
 // This just to see how "looks like" the distribution.  Normally, in real applications,
 // you generate just _one_  point
for (i=0;i<N;i++)
{

// ----------------------------------- Step 1.  Direction
    length=0;
   for (j=0;j<D;j++)
   {
          x.x[j]=alea_normal(0,1);
          length=length+  x.x[j]*x.x[j];
   }

   length=sqrt(length);

 //----------------------------------- Step 2.   Random radius

  r=alea_float(0,1);
  r=pow(r,pw);

       for (j=0;j<D;j++)
   {
          x.x[j]=r*x.x[j]/length;
   }

 
 //--------------------------------------------------------
 display_vector(x);
 save_vector(f_run,x);
 } //------------------------------------------------------------------------ End of loop on points
 
 goto input;   // Ask for a new set of points

 end:
 printf("\n %i generated points saved on run.txt" ,N);
 return EXIT_SUCCESS;
}

// ========================================================================= ALEA_FLOAT
double alea_float(double min, double max) /* Random  number between min and max */
{
double 	r;

// Normally, RAND_MAX = 32767 = 2^31-1

	r=(double)rand()/RAND_MAX; // r in [0,1]
	r=min+r*(max-min);
	return r;
}
//==========================================================================  ALEA_NORMAL
double alea_normal(double mean, double std_dev)
 {
 /*
     Use the polar form of the Box-Muller transformation to obtain
     a pseudo random number from a Gaussian distribution
 */
        double x1, x2, w, y1;
        //double  y2;

         do {
                 x1 = 2.0 * alea_float(0,1) - 1.0;
                 x2 = 2.0 * alea_float(0,1) - 1.0;
                 w = x1 * x1 + x2 * x2;
         } while ( w >= 1.0);

         w = sqrt( (-2.0 * log( w ) ) / w );
         y1 = x1 * w;
        // y2 = x2 * w;
          y1=y1*std_dev+mean;
         return y1;
  }

//========================================================================== DISPLAY_VECTOR
void display_vector(struct vector x)
{
int				d;
int				ncar=0;

//printf("\n Vector size %i\n ",x.size);
printf("\n");
for (d=0;d<x.size;d++)
{
	printf("%6.6f ",x.x[d]);
	ncar=ncar+10;
	if (ncar>75)
	{
		ncar=0;
		printf("\n");
	}
}
}



 //========================================================================== SAVE_VECTOR
void  save_vector(FILE *f_run,struct vector x)
{
int				d;
int				ncar=0;

fprintf(f_run,"\n");
for (d=0;d<x.size;d++)
{
	fprintf(f_run,"%6.6f ",x.x[d]);
	ncar=ncar+10;
	if (ncar>75)
	{
		ncar=0;
		fprintf(f_run,"\n");
	}
}
}