Subj : Re: Compiler and an interpreter
To   : comp.programming
From : Gerry Quinn
Date : Sat Aug 06 2005 10:45 pm

In article <42f4f604$0$91541$ed2e19e4@ptn-nntp-reader04.plus.net>, 
usenet@jdh30.plus.com says...
> Gerry Quinn wrote:

> The set of "n"th-nearest neighbours is defined as:
> 
> nth(i, 0) = {i}
> nth(i, 1) = neighbours(i)
> nth(i, n) = union(j in nth(i,n-1), nth(i, 1)) - nth(i, n-1) - nth(i, n-2)
> 
> where neighbours(i) gives the set of neighbours of "i".

I don't understand the last line, so I may have solved a different 
problem!  It is, at least, a similar problem.  I calculate a layer of 
neighbours expanding outward from a single point - in a simple cubic 
lattice they would be a series of octahedral layers of increasing size.  

In other words, I am saying:
nth( i, n ) = all neighbours of nth( i, n-1 ), that don't appear in 
nth( i, n-1 ) or nth( i, n-2 )

The following compiles, but I haven't tested it so it could be bugged.  
I haven't used sets before but they seem straightforward.  I use 
vectors where possible and don't use complicated nested templates. I 
think it's fairly comprehensible what I'm doing.  At least it's 
comprehensible to me whereas your code isn't, but perhaps that 
situation is symmetric!  

I did hardcode 3D, for efficiency and because it's faster to type - it 
would be easy to use a vector.

I'd be surprised if it was particularly slow.  However, it just occurs 
to me that one could modify the FindNeighbours() function to check the 
m_Prev and m_Current layers before inserting a neighbour into the 
m_Next.  This would obviate the need for set_difference() functions 
because no invalid neighbour would ever be inserted.  I don't know if 
it's a win or not.

I've a feeling your nested templates are the cause of slowness 
problems.


// A 'neighbour' is a real atom with an offset
class Neighbour
{
public:
	int m_Cell[ 3 ];		
	int m_nAtomID;
public:
	bool operator < ( const Neighbour & other ) const;
};


// An atom is a defined virtual atom in the unit cell
class Atom
{
public:
	int m_nID;
	vector< Neighbour > m_Neighbours;
};


typedef set< Neighbour > nbset;


// Definition of lattice unit cell
class Cell  
{
public:
	vector< Atom > m_Atoms;		
public:
	void FindNeighbours( const nbset & current, 
			nbset & neighbours ) const;
	void FindNthNeighbours( int nAtomID, int nth, 
			set< Neighbour > & result ) const;
}; 


// Required for set operations that require sorting
// You don't need an AtomCompare struct!
bool Neighbour::operator < ( const Neighbour & other ) const 
{
	if ( m_nAtomID < other.m_nAtomID ) return true;
	if ( m_nAtomID != other.m_nAtomID ) return false;
	for ( int iDim = 0; iDim < 3; iDim++ )
	{		
		if ( m_Cell[ iDim ] < other.m_Cell[ iDim ] ) return true;
		if ( m_Cell[ iDim ] != other.m_Cell[ iDim ] ) return false;
	}
	return false;
}


// Find the unique neighbours of a given set of neighbours
void Cell::FindNeighbours( const nbset & current, 
		nbset & neighbours ) const
{
	for ( nbset::const_iterator it = current.begin(); 
			it != current.end(); it++ )
	{
		const Neighbour & nbwork = *it;
		const Atom & atom = m_Atoms[ nbwork.m_nAtomID ];
		for ( int iNg = 0; iNg < atom.m_Neighbours.size(); iNg++ )
		{
			const Neighbour & nb = atom.m_Neighbours[ iNg ];
			Neighbour newNb;
			for ( int iDim = 0; iDim < 3; iDim++ )
			{
newNb.m_Cell[ iDim ] = nbwork.m_Cell[ iDim ] + nb.m_Cell[ iDim ];
				newNb.m_nAtomID = nb.m_nAtomID;
			}			
			if ( neighbours.find( newNb ) != neighbours.end() )
			{
				neighbours.insert( newNb );
			}
		}
	}
}


// Find the nth layer of neighbours
void Cell::FindNthNeighbours( int nAtomID, int nth, 
		set< Neighbour > & result ) const
{
	set< Neighbour > m_Prev;		// n-2 layer
	set< Neighbour > m_Current;		// n-1 layer
	set< Neighbour > m_Next;		// n layer
	
	// Starting atom, converted to neighbour at offset 0,0,0
	Neighbour start;
	for ( int iDim = 0; iDim < 3; iDim++ )
	{
		start.m_Cell[ iDim ] = 0;
	}
	start.m_nAtomID = nAtomID;
	m_Prev.insert( start );
	if ( nth == 0 ) 
	{
		swap( m_Prev, result );
		return;
	}

	// Current layer is layer 1
	FindNeighbours( m_Prev, m_Current );
	if ( nth == 1 ) 
	{
		swap( m_Current, result );
		return;
	}

	// Do all the layers
	for ( int iLayers = 2; iLayers <= nth; iLayers++ )
	{
		set< Neighbour > work;
		set< Neighbour > work2;
		FindNeighbours( m_Current, work );
		set_difference( work.begin(), work.end(), 
m_Current.begin(), m_Current.end(), work2.begin() );		
		set_difference( work2.begin(), work2.end(), 
m_Prev.begin(), m_Prev.end(), m_Next.begin() );		
		
		if ( iLayers == nth )
		{
			swap( m_Next, result );
			return;
		}

		swap( m_Prev, m_Current );
		swap( m_Current, m_Next );
		m_Next.clear();
	}
}


- Gerry Quinn

.