请教算法：“树“的构造

[xiongzm]

就是构造一棵“完全2叉树”， 只是在父节点与子节点之间还有sp-1个节点，只是这些节点只有
一个子节点，比如说构造一个深度为3（或者说2代）
，“父”节点与“子”节点之间还有2（sp=3,sp-1=2)个节点的一棵这样的树，共有14×3＋1
＝43个节点。这些节点都是一个一个的在三位空间中随机生成的，当然这些节点是不能重叠的， 因此要得到生成的每个节点的三维坐标， 同时要设计好每个节点的
编号，这样以便打印出这棵树中各个节点的连接性。怎么设计数据结构，线性的还是链式结构了？想模拟一个小问题，如果数据结构设计的好，以后的数据分析就很方便了。

[xiongzm]

我原来写过一个程序，只生成了一棵树，现在想生成更多的，比如30棵，为了补重叠，密度可以非常小。原来的程序很乱的，先帖出来。


#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define MAXATOMS 10000


int gen, nseg, natoms, nbonds, nodes, gnodes;

double x[MAXATOMS], y[MAXATOMS], z[MAXATOMS];
int imagex[MAXATOMS], imagey[MAXATOMS], imagez[MAXATOMS];
double bondlength;
double xlo, xhi, ylo, yhi, zlo, zhi, xprd, yprd, zprd;
double log2;

unsigned int seed;
struct tms *tt;

void getcos( double *, double *, double * );
void pbc( double *, double *, double * );
void image( double tx, double ty, double tz, int n );
int check( double x, double y, double z, int inode, int offset );
int getnode();
void init();
void writeheader();
void writeatoms();
void writebonds();





int main()
{
int height, repeat=0, saw;
int inode, i, j, nextnode;
int parent, child;
int goffset1, goffset2, goffset, offset;
double xsurf, ysurf, zsurf;
double xx, yy, zz, tx, ty, tz, xij, yij, zij, rij;
/*
unsigned int seed;
struct tms *tt;
*/
seed = times( tt );
srand( seed );

init();

x[0] = 0.0;
y[0] = 0.0;
z[0] = 0.0;
imagex[0] = 0;
imagey[0] = 0;
imagez[0] = 0;

// node information store in the front of array, while their children
// store after all the nodes
for( inode = 0; inode < gnodes; inode++ ) // inode: from 0
{
printf( "================inode = %d=================\n", inode );

height = (int) ( log( inode + 2 ) / log2 );
if( height < log( inode + 2 )/log2 )
height += 1; // height = gen + 1

goffset1 = (int) ( pow( 2, height ) - 2 ); // group offset1: the number of nodes until the height of the inode
// that is the number of branches
// via each node to generate its two branches

goffset2 = 2 * (int) ( inode - pow( 2, height - 1 ) + 1 ); //

goffset = goffset1 + goffset2; // the generated branches

offset = ( nseg - 1 ) * goffset; // the offset number of atoms related to the nodes

// generate left branch
for( i = 0; i < nseg; i++ ) // for every node, generate the atom until the next generation's node
{
do
{
getcos( &xsurf, &ysurf, &zsurf );

if( i == 0 )
{
parent = inode;
child = nodes + offset;
}
else
{
parent = nodes + offset + i - 1;
child = parent + 1;
}

xx = x[parent] + xsurf * bondlength;
yy = y[parent] + ysurf * bondlength;
zz = z[parent] + zsurf * bondlength;
tx = xx;
ty = yy;
tz = zz;

pbc( &xx, &yy, &zz );
repeat = check( xx, yy, zz, inode, child );

/* if( i == nseg - 1 )
{
nextnode = getnode( inode );
repeat = check( xx, yy, zz, nextnode - 1, child );
}
else
repeat = check( xx, yy, zz, inode, child );
*/
}
while( repeat == 1 );

x[child] = xx;
y[child] = yy;
z[child] = zz;

// printf( "x=%f y=%f z=%f\n", xx, yy, zz );

image( tx, ty, tz, child );
}

nextnode = getnode( inode );

x[nextnode] = xx;
y[nextnode] = yy;
z[nextnode] = zz;


image( tx, ty, tz, nextnode );

// generate right branch
for( i = 0; i < nseg; i++ )
{
do
{
getcos( &xsurf, &ysurf, &zsurf );

if( i == 0 )
{
parent = inode;
child = nodes + offset + nseg - 1;
}
else
{
parent = nodes + offset + nseg + i - 2;
child = parent + 1;
}

xx = x[parent] + xsurf * bondlength;
yy = y[parent] + ysurf * bondlength;
zz = z[parent] + zsurf * bondlength;
tx = xx;
ty = yy;
tz = zz;

pbc( &xx, &yy, &zz );
repeat = check( xx, yy, zz, inode, child );

// check function: you should make sure how many nodes should test?
// so, we take search method in check function, not judge here

/* if( i == nseg - 1)
repeat = check( xx, yy, zz, nextnode, child );
else
repeat = check( xx, yy, zz, inode, child );
*/
}
while( repeat == 1 );



x[child] = xx;
y[child] = yy;
z[child] = zz;

image( tx, ty, tz, child );
}

nextnode += 1;
x[nextnode] = xx;
y[nextnode] = yy;
z[nextnode] = zz;
image( tx, ty, tz, nextnode );

}

saw = 0;
for( i = 0; i < natoms - 1; i++ )
for( j = i + 1; j < natoms; j++ )
{
xij = x[i] - x[j];
yij = y[i] - y[j];
zij = z[i] - z[j];
rij = xij * xij + yij * yij + zij * zij;
if( rij < 0.98 * bondlength * bondlength )
{
saw += 1;
printf("%d---%d %f\n", i, j, rij );

}
}

printf("\n\n ------- saw = %d \n\n", saw );

writeatoms();
writebonds();

}

[xiongzm]

// ----------------------------------------------------------------
// init() get the parameter
//
void init()
{
double rhostar, volume;
int i;

// get input parameter

// printf( "Please input how many generations of your dendrimer: " );
scanf( "%d", &gen );
//printf( "Please input the spacer segments: " );
scanf( "%d", &nseg );
//printf( "Please input the bond length (float):" );
scanf( "%lf", &bondlength );
//printf( "%f", bondlength );
//printf( "Please input the number density of your dendrimer monomers: " );
scanf( "%lf", &rhostar );

// calculate the basic parameters of the dendrimer
natoms = nseg * pow( 2, gen + 2 ) - 2 * nseg + 1;
// natoms = ( 1 + nseg ) * pow( 2, gen + 1 ) + pow( 2, gen ) - 1;
//printf( "natoms = %d\n", natoms );
nodes = pow( 2, gen + 2 ) - 1; // # of nodes, for example: gen=1, nodes=7
// but index from
//printf( "nodes = %d\n", nodes );
gnodes= pow(2, gen + 1 ) -1; // # ofmn nodes, to generation
//printf( "gnodes =%d\n", gnodes );

// calculate the cell information
volume = 1.0 * natoms/rhostar;
xprd = pow( volume, 1.0/3.0 );
yprd = xprd;
zprd = xprd;
xlo = -xprd/2.0;
xhi = -xlo;
ylo = xlo;
yhi = xhi;
zlo = xlo;
zhi = xhi;

for( i = 0; i < natoms; i++ )
{
imagex[i] = 0;
imagey[i] = 0;
imagez[i] = 0;
x[i] = 9999.0;
}

//printf( "volume=%f\n", volume );
//printf( "----------xlo=%f xhi=%f------------------------\n", xlo, xhi );
log2 = log(2);


}



// ----------------------------------------------------------------
// image convention
// it's a reference, child's image box refered to the parent
//
/* the absolute version of image */

void image( double tx, double ty, double tz, int n )
{
if( tx > xhi )
imagex[n] = 1 + (int)( ( tx - xhi )/xprd );
else if( tx < xlo )
imagex[n] = (int)( ( tx + xlo )/xprd );

if( ty > yhi )
imagey[n] = 1 + (int)( ( ty - yhi )/yprd );
else if( ty < ylo )
imagey[n] = (int)( ( ty + ylo )/yprd );

if( tz > zhi )
imagez[n] = 1 + (int)( ( tz - zhi )/zprd );
else if( tz < zlo )
imagez[n] = (int)( ( tz + zlo )/zprd );

}



/*
void image( int n )
{
if( abs( x[n] - x[n-1] ) < 2.0 * bondlength )
imagex[n] = imagex[n-1];
else if( ( x[n] - x[n-1] ) < 0.0 )
imagex[n] = imagex[n-1] + 1;
else if( ( x[n] - x[n-1] ) > 0.0 )
imagex[n] = imagex[n-1] - 1;

if( abs( y[n] - y[n-1] ) < 2.0 * bondlength )
imagey[n] = imagey[n-1];
else if( ( y[n] - y[n-1] ) < 0.0 )
imagey[n] = imagey[n-1] + 1;
else if( ( y[n] - y[n-1] ) > 0.0 )
imagey[n] = imagey[n-1] - 1;

if( abs( z[n] - z[n-1] ) < 2.0 * bondlength )
imagez[n] = imagez[n-1];
else if( ( z[n] - z[n-1] ) < 0.0 )
imagez[n] = imagez[n-1] + 1;
else if( ( z[n] - z[n-1] ) > 0.0 )
imagez[n] = imagez[n-1] - 1;

}
*/

// ----------------------------------------------------------------
// periodic boundary conditions
// map atom back into periodic box
//
void pbc( double *x, double *y, double *z )
{
if( *x < xlo )
*x += xprd;
if( *x >= xhi )
*x -= xprd;
if( *y < ylo )
*y += yprd;
if( *y >= yhi )
*y -= yprd;
if( *z < zlo )
*z += zprd;
if( *z >= zhi )
*z -= zprd;
}


// ----------------------------------------------------------------
// get random cosine value
// first get random point inside sphere of unit radius
// then project point to surface of sphere of unit radius
// which uses analogue triangle
// that is, get three random cosine value
//
void getcos( double *xsurf, double *ysurf, double *zsurf )
{
double xinner, yinner, zinner, rsq, r;

do
{
xinner = 2.0 * (double) rand() / (double)(RAND_MAX) - 1.0;
yinner = 2.0 * (double) rand() / (double)(RAND_MAX) - 1.0;
zinner = 2.0 * (double) rand() / (double)(RAND_MAX) - 1.0;

rsq = xinner * xinner + yinner * yinner + zinner * zinner;
}
while( rsq > 1.0 );

r = sqrt( rsq );
*xsurf = xinner/r;
*ysurf = yinner/r;
*zsurf = zinner/r;
}


// ---------------------------------------------------------------------
// getnode()
// from inode to get nextnode, its left node
//
int getnode( int inode )
{
int height, goffset1, goffset2;
height = (int) ( log( inode + 2 ) / log2 );
if( height < log( inode + 2 )/log2 )
height += 1;
goffset1 = (int) ( pow( 2, height ) - 2 ); // group offset1: the number of nodes until the height of the inode
// that is the number of branches
// via each node to generate its two branches

goffset2 = 2 * (int) ( inode - pow( 2, height - 1 ) + 1 ); //

return goffset1 + goffset2 + 1;

}


// ---------------------------------------------------------------------
// check function, return whether the given point overlaps with
// other atoms. explicit parameter inode and offset avoid to
// calculate it again
//
int check( double xx, double yy, double zz, int inode, int child )
{
int i, height, thenode;
double dx, dy, dz, rsq;

height = (int) ( log( inode + 2 ) / log2 );
thenode = pow( 2, height + 1 ) - 2;

for( i = 0; i < nodes; i++ )
{
if( x[i] != 9999.0 && x[i] != xx )
{

dx = ( x[i] - xx ) - (int)( ( x[i] - xx )/xprd ) * xprd;
dy = ( y[i] - yy ) - (int)( ( y[i] - yy )/yprd ) * yprd;
dz = ( z[i] - zz ) - (int)( ( z[i] - zz )/zprd ) * zprd;

if( dx > xhi )
dx = xprd - dx;
else if( dx < xlo )
dx = xprd + dx;
if( dy > yhi )
dy = yprd - dy;
else if( dy < ylo )
dy = yprd + dy;
if( dz > zhi )
dz = zprd - dz;
else if( dz < zlo )
dz = zprd + dz;

// which method is quick?
// r by sqrt function or rsq use bondlength times bondlength
rsq = dx * dx + dy * dy + dz * dz;
if( rsq < 0.98 * bondlength * bondlength )
return 1;

}
}

for( i = nodes; i < child; i++ )
{
dx = ( x[i] - xx ) - (int)( ( x[i] - xx )/xprd ) * xprd;
dy = ( y[i] - yy ) - (int)( ( y[i] - yy )/yprd ) * yprd;
dz = ( z[i] - zz ) - (int)( ( z[i] - zz )/zprd ) * zprd;

if( dx > xhi )
dx = xprd - dx;
else if( dx < xlo )
dx = xprd + dx;
if( dy > yhi )
dy = yprd - dy;
else if( dy < ylo )
dy = yprd + dy;
if( dz > zhi )
dz = zprd - dz;
else if( dz < zlo )
dz = zprd + dz;

rsq = dx * dx + dy * dy + dz * dz;

if( rsq < 0.98* bondlength * bondlength )
return 1;
}

return 0;
}

[xiongzm]

// ----------------------------------------------------------------
// data file begins here:
// simple printf is just used, if you don't specify pipe file,
// the data file output to the stdout, so you should run the program
// as: dendrimer <input-file >output-file
//
// ----------------------------------------------------------------
//



// ---------------------------------------------------------------------
// writeatoms()
//
void writeatoms()
{
int i;
printf( "%s\n", "Atoms" );
printf( "\n" );

for( i = 0; i < natoms; i++ )
printf( "%d %d %d %f %f %f %d %d %d\n", i, 1, 1, x[i], y[i], z[i], imagex[i], imagey[i], imagez[i] );
printf( "\n" );

}

// ----------------------------------------------------------------------
// writebonds()

void writebonds()
{
int inode, i;
int height, goffset1, goffset2, goffset, offset;
int batom1, batom2;
float xij, yij, zij;


printf( "%s\n", "Bonds" );
printf( "\n" );



nbonds = 0; // count the number of bonds

for( inode = 0; inode < gnodes; inode++ )
{
height = (int) ( log( inode + 2 ) / log2 );
if( height < log( inode + 2 )/log2 )
height += 1; // height = gen + 1

goffset1 = (int) ( pow( 2, height ) - 2 ); // group offset1: the number of nodes until the height of the inode
// that is the number of branches
// via each node to generate its two branches
goffset2 = 2 * (int) ( inode - pow( 2, height - 1 ) + 1 ); //
goffset = goffset1 + goffset2; // the generated branches
offset = ( nseg - 1 ) * goffset;

// write bonds on left branch

for( i = 0; i < nseg; i++ )
{
if( i == 0 )
{
batom1 = inode;
batom2 = nodes + offset;
}
else if( i == nseg - 1 )
{
batom1 = nodes + offset + i - 1;
batom2 = getnode( inode );
}
else
{
batom1 = nodes + offset + i - 1;
batom2 = batom1 + 1;
}

nbonds += 1;
printf( "%d %d %d %d\n", nbonds, 1, batom1, batom2 );

xij = x[batom1] - x[batom2];
yij = y[batom1] - y[batom2];
zij = z[batom1] - z[batom2];

printf( "%d---%d %f\n", batom1, batom2, xij * xij + yij * yij + zij * zij );
}

// write bonds on right branch
for( i = 0; i < nseg; i++ )
{
if( i == 0 )
{
batom1 = inode;
batom2 = nodes + offset + nseg - 1;
}
else if( i == nseg - 1 )
{
batom1 = nodes + offset + nseg + i - 2;
batom2 = getnode( inode ) + 1;
}
else
{
batom1 = nodes + offset + nseg + i - 2;
batom2 = batom1 + 1;
}

nbonds += 1;
printf( "%d %d %d %d\n", nbonds, 1, batom1, batom2 );

xij = x[batom1] - x[batom2];
yij = y[batom1] - y[batom2];
zij = z[batom1] - z[batom2];

printf( "%d---%d %f\n", batom1, batom2, xij * xij + yij * yij + zij * zij );
}
}

}