File Coverage

blib/lib/AI/NeuralNet/SOM.pm

Criterion	Covered	Total	%
statement	279	582	47.9
branch	46	176	26.1
condition	15	144	10.4
subroutine	25	34	73.5
pod	17	17	100.0
total	382	953	40.0

#!/usr/bin/perl

# $Id: SOM.pm, v0.01 2000/09/28
#
# Copyright (c) 2000  Alexander Voischev, Russian Federation
#
# See AUTHOR section in pod text below for usage and distribution rights.
#

BEGIN {
	$AI::NeuralNet::SOM::VERSION = "0.02";
}

package AI::NeuralNet::SOM;

require 5.005_62;
use strict;
use warnings;
use Carp;
use POSIX;
use locale;

require Exporter;
use AutoLoader qw(AUTOLOAD);

our @ISA = qw(Exporter);

our $VERSION = '0.1';

$AI::NeuralNet::SOM::INV_ALPHA_CONSTANT  = 100.0;

#####################################################
#
# "Public" methods.
#
#####################################################

	##################################################
	#
	#
	#
	##################################################
	sub new {
		my $class = shift;
		my $self = {};
		bless($self, $class);
		$self->{XDIM} = 0;
		$self->{YDIM} = 0;
		$self->{IDIM} = 0;
		$self->{MAP} = ();
		$self->clear_all_labels;
		return $self;
	}

	##################################################
	#
	#
	#
	##################################################
	sub initialize {
	   my $self = shift;
		my $xdim = shift; # X Dimension of map
		my $ydim = shift; # Y Dimension of map
		my $idim = shift; # Dimension of input data
		my $topology = shift; # Map topolgy
		my $neighborhood = shift; # neighborhood function type
		my $init_type = shift; # Initialization type
		my $random_seed = shift; # Random seed
		my $data = shift; # Data values for initialize

		# Checking error in arguments
		croak "Invalid XDim parameter" if ($xdim < 1);
		croak "Invalid YDim parameter" if ($ydim < 1);
		croak "Invalid IDim parameter" if ($idim < 1);
		croak "Invalid topology parameter" if ($topology ne 'hexa') and ($topology ne 'rect');
		croak "Invalid neighborhood parameter" if ($neighborhood ne 'bubble') and ($neighborhood ne 'gaussian');
		croak "Invalid initialization type parameter" if ($init_type ne 'random') and ($init_type ne 'linear');
		croak "Invalid Random seed parameter" if ($random_seed < 0);
		croak "Invalid Data parameter" if ((ref($data) eq 'ARRAY') and (@$data % $idim));

		$self->{XDIM} = $xdim;
		$self->{YDIM} = $ydim;
		$self->{IDIM} = $idim;
		$self->{TOPOLOGY} = $topology;
		$self->{NEIGHBORHOOD} = $neighborhood;
		$self->{RANDOMSEED} = $random_seed;

		my $patterns_count = @$data / $idim;

		srand($random_seed);
		if ($init_type eq 'random') {
			my @max_vals = ();
			my @min_vals = ();
			for my $i (0..$idim-1) {
				my $max = 0.0 - DBL_MAX;
				my $min = DBL_MAX;
				for my $p (0..$patterns_count-1) {
					if ($data->[$i + $p * $idim] =~ /^-?\d+\.?\d*$/ ) { # Is a real number
						$max = $data->[$i + $p * $idim] if ($max < $data->[$i + $p * $idim]);
						$min = $data->[$i + $p * $idim] if ($min > $data->[$i + $p * $idim]);
					}
				}
		   	$max_vals[$i] = $max;
				$min_vals[$i] = $min;
			}

			my $c = 0;
			for my $z (0..$xdim * $ydim-1) {
				for my $i (0..$idim-1) {
					$self->{MAP}->[$c++] = $min_vals[$i] + rand($max_vals[$i]-$min_vals[$i]);
				}
			}
		}
		if ($init_type eq 'linear') {
			my $ind = 0;
			my $eigen1;
			my $eigen2;
			my $mean;
			my $xf;
			my $yf;
			($eigen1, $eigen2, $mean) = $self->_find_two_eigenvectors_and_mean($data);

			for my $z (0..$self->{XDIM}*$self->{YDIM}-1) {
				$xf = 4.0 * int($ind % $self->{XDIM}) / ($self->{XDIM} - 1.0) - 2.0;
				$yf = 4.0 * int($ind / $self->{XDIM}) / ($self->{YDIM} - 1.0) - 2.0;
				for my $i (0..$self->{IDIM}-1) {
					$self->{MAP}->[$z * $self->{IDIM} + $i] = $mean->[$i] + $xf * $eigen1->[$i] + $yf * $eigen2->[$i];
				}
				$ind++;
			}
		}
	}

	##################################################
	#
	#
	#
	##################################################
	sub train {
		my $self = shift;
		my ($train_length, $alpha, $radius, $alpha_type, $data) = @_;

	  	croak "Invalid Data parameter" if ((ref($data) eq 'ARRAY') and (@$data % $self->{IDIM}));

		my ($adapt_func, $alpha_func);
		if ($self->{NEIGHBORHOOD} eq 'bubble') {
			$adapt_func = '_adapt_bubble';
		}
		else {
			$adapt_func = '_adapt_gaussian';
		}

		if ($alpha_type eq 'linear') {
			$alpha_func = '_alpha_linear';
		}
		else {
			$alpha_func = '_alpha_inverse_t';
		}

		my $p = 0;
		my $patterns_count = @$data / $self->{IDIM};
		my ($trad, $talp);
		my ($xwin, $ywin, $min_diff);

		my @data_slice;
		for my $len (0..$train_length-1) {

			if (++$p == $patterns_count) {$p = 0;}

			$trad = 1.0 + ($radius - 1.0) * ($train_length - $len) / $train_length;
			$talp = $self->$alpha_func($len, $train_length, $alpha);

			@data_slice = @$data[$p*$self->{IDIM}..($p+1)*$self->{IDIM}-1];
			($xwin, $ywin, $min_diff) = $self->winner(\@data_slice);

			$self->$adapt_func($xwin, $ywin, $trad, $talp, \@data_slice);
		}
	}

	##################################################
	#
	#
	#
	##################################################
	sub qerror {
		my $self = shift;
		my $data = shift;

	  	croak "Invalid Data parameter" if ((ref($data) eq 'ARRAY') and (@$data % $self->{IDIM}));

		my $patterns_count = @$data / $self->{IDIM};
		my ($xwin, $ywin, $min_diff);
		my $qerror=0;

		my @data_slice;
		for my $p (0..$patterns_count-1) {
			@data_slice = @$data[$p*$self->{IDIM}..($p+1)*$self->{IDIM}-1];
			($xwin, $ywin, $min_diff) = $self->winner(\@data_slice);

			$qerror += $min_diff;
		}
		return ($qerror/$patterns_count);
	}

	##################################################
	#
	#
	#
	##################################################
	sub winner {
		my $self = shift;
		my $data = shift;

		my ($xwin, $ywin);
		my ($diff, $difference, $masked);
		my $min_diff = DBL_MAX;

		for my $y (0..$self->{YDIM}-1) {
			for my $x (0..$self->{XDIM}-1) {
				$masked = 0;
				$difference = 0.0;
				for my $i (0..$self->{IDIM}-1) {
					if ($data->[$i] =~ /^-?\d+\.?\d*$/) { # Is a real number
						$diff = $self->{MAP}->[($y * $self->{XDIM} + $x) * $self->{IDIM} + $i] - $data->[$i];
						$difference += $diff ** 2;
					}
					else {
						$masked++;
					}
				}
				# If data pattern is empty
				croak "Empty pattern" if ($masked == $self->{IDIM});

				# If distance is smaller than previous distances
				if ($difference < $min_diff) {
					$xwin = $x;
					$ywin = $y;
					$min_diff = $difference;
				}
			}
		}
		return ($xwin, $ywin, sqrt($min_diff));
	}

	##################################################
	#
	#
	#
	##################################################
	sub set_label {
		my $self = shift;
		my ($x, $y, $label) = @_;
		croak "Invalid argument" if (($x<0) or ($x>=$self->{XDIM}) or ($y<0) or ($y>=$self->{YDIM}) or not defined($label));
		$self->{LABELS}->[$y * $self->{XDIM} + $x] =$label;
	}

	##################################################
	#
	#
	#
	##################################################
	sub clear_all_labels {
		my $self = shift;
		$self->{LABELS} = ();
	}

	##################################################
	#
	#
	#
	##################################################
	sub save {
		my $self = shift;
		my $file = shift;

		my $xdim = $self->{XDIM};
		my $ydim = $self->{YDIM};
		my $idim = $self->{IDIM};
		my $topol = $self->{TOPOLOGY};
		my $neigh = $self->{NEIGHBORHOOD};
		print $file "# Created by SOM.pm v0.01 by Voischev Alexander e-mail: voischev\@mail.ru\n";
		print $file "$idim $topol $xdim $ydim $neigh\n";

		for my $z (0..$self->{XDIM}*$self->{YDIM}-1) {
			for my $i (0..$self->{IDIM}-1) {
				my $value = $self->{MAP}->[$z * $self->{IDIM} + $i];
				print $file "$value ";
			}
			my $label = $self->{LABELS}->[$z];
			print $file "$label" if (defined($label));
			print $file "\n";
		}
	}

	##################################################
	#
	#
	#
	##################################################
	sub load {
		my $self = shift;
		my $file = shift;
		my $header;

		while (<$file>) {
			if (!/^ *#/) {
				s/^ *//;
				$header = $_;
				last;
			}
		}
		chomp($header);
		my ($idim, $topology, $xdim, $ydim, $neighborhood) = split(/ /, $header);

		croak "Invalid XDim parameter" if ($xdim < 1);
		croak "Invalid YDim parameter" if ($ydim < 1);
		croak "Invalid IDim parameter" if ($idim < 1);
		croak "Invalid topology parameter" if ($topology ne 'hexa') and ($topology ne 'rect');
		croak "Invalid neighborhood parameter" if ($neighborhood ne 'bubble') and ($neighborhood ne 'gaussian');

		$self->{XDIM} = $xdim;
		$self->{YDIM} = $ydim;
		$self->{IDIM} = $idim;
		$self->{TOPOLOGY} = $topology;
		$self->{NEIGHBORHOOD} = $neighborhood;

		my @pattern;
		my @line;
		my @data;
		my $z = 0;
		while (<$file>) {
			if (!/^ *#/) {
				chomp();
				if (/#/) {
				    @line = split(/ /,$`);
				}
				else {
				    @line = split(/ /);
				}
				@pattern = splice(@line, 0, $self->{IDIM});
				push (@data, @pattern);
				if (defined($line[0])) {
					$self->{LABELS}->[$z] = $line[0];
				}
				$z++;
			}
		}
		$self->{MAP} = \@data;
	}

	##################################################
	#
	#
	#
	##################################################
	sub umatrix {
		my $self = shift;
		my ($i,$j,$k,$count,$bx,$by,$bz);
		my ($dx,$dy,$dz1,$dz2,$dz,$temp,$max,$min,$bw);
		my @medtable;
		my $tmp;

		my @umat = ();

		if ($self->{XDIM}<=0 or $self->{YDIM}<=0 or $self->{IDIM}<=0) {
			return undef;
		}
		$max = 0;
		$min = 0;

		if ($self->{TOPOLOGY} eq 'rect') {
		# Rectangular topology
			for $j (0..$self->{YDIM}-1) {
				for $i (0..$self->{XDIM}-1) {
					$dx=0; $dy=0; $dz1=0; $dz2=0; $count=0;
					$bx=0; $by=0; $bz=0;
					for $k (0..$self->{IDIM}-1) {
						if ($i < $self->{XDIM}-1) {
							$temp = $self->map($i,$j,$k) - $self->map($i+1,$j,$k);
							$dx += $temp ** 2;
							$bx = 1;
						}
						if ($j < $self->{YDIM}-1) {
							$temp = $self->map($i,$j,$k) - $self->map($i,$j+1,$k);
							$dy += $temp ** 2;
							$by = 1;
						}
						if ($j < $self->{YDIM}-1 and $i < $self->{XDIM}-1) {
							$temp = $self->map($i,$j,$k) - $self->map($i+1,$j+1,$k);
							$dz1 += $temp ** 2;
							$temp = $self->map($i,$j+1,$k) - $self->map($i+1,$j,$k);
							$dz2 += $temp ** 2;
							$bz=1;
						}
					}
					$dz = (sqrt($dz1)/sqrt(2.0)+sqrt($dz2)/sqrt(2.0))/2;

					if ($bx) {
						$umat[2*$i+1+(2*$j)*($self->{XDIM}*2-1)] = sqrt($dx);
					}
					if ($by) {
						$umat[2*$i+(2*$j+1)*($self->{XDIM}*2-1)] = sqrt($dy);
					}
					if ($bz) {
						$umat[2*$i+1+(2*$j+1)*($self->{XDIM}*2-1)] = $dz;
					}
				}
			}
		}
		else {
		# Hexagonal topology
			for $j (0..$self->{YDIM}-1) {
				for $i (0..$self->{XDIM}-1) {
					$dx=0; $dy=0; $dz=0; $count=0;
					$bx=0; $by=0; $bz=0;
					$temp=0;
					if ($i<$self->{XDIM}-1)
					{
						for $k (0..$self->{IDIM}-1) {
							$temp = $self->map($i,$j,$k) - $self->map($i+1,$j,$k);
							$dx += $temp ** 2;
							$bx = 1;
						}
					}
					$temp=0;
					if ($j < $self->{YDIM}-1) {
						if ($j%2) {
							for $k (0..$self->{IDIM}-1)
							{
								$temp = $self->map($i,$j,$k) - $self->map($i,$j+1,$k);
								$dy += $temp ** 2;
								$by=1;
							}
						}
						else {
							if ($i>0) {
								for $k (0..$self->{IDIM}-1) {
									$temp = $self->map($i,$j,$k) - $self->map($i-1,$j+1,$k);
									$dy += $temp ** 2;
									$by=1;
								}
							}
							else {
								$temp=0;
							}
						}
					}
					$temp=0;
					if ($j < $self->{YDIM}-1) {
						if (!($j%2)) {
							for $k (0..$self->{IDIM}-1) {
								$temp = $self->map($i,$j,$k) - $self->map($i,$j+1,$k);
								$dz += $temp ** 2;
							}
							$bz=1;
						}
						else {
							if ($i < $self->{XDIM}-1) {
								for $k (0..$self->{IDIM}-1) {
									$temp = $self->map($i,$j,$k) - $self->map($i+1,$j+1,$k);
									$dz += $temp ** 2;
								}
								$bz=1;
							}
						}
					}
					else {
						$temp=0;
					}

					if ($bx) {
						$umat[2*$i+1+(2*$j)*($self->{XDIM}*2-1)] = sqrt($dx);
					}
	   	 	   if ($by) {
						if ($j%2) {
							$umat[2*$i+(2*$j+1)*($self->{XDIM}*2-1)] = sqrt($dy);
						}
						else {
							$umat[2*$i-1+(2*$j+1)*($self->{XDIM}*2-1)] = sqrt($dy);
						}
					}
					if ($bz) {
						if ($j%2) {
							$umat[2*$i+1+(2*$j+1)*($self->{XDIM}*2-1)] = sqrt($dz);
						}
						else {
							$umat[2*$i+(2*$j+1)*($self->{XDIM}*2-1)] = sqrt($dz);
						}
					}
				}
			}
		}

		# Set the values corresponding to the model vectors themselves
		# to medians of the surrounding values
		if ($self->{TOPOLOGY} eq 'rect') {
		# Rectangular topology
		# medians of the 4-neighborhood
			for ($j=0; $j<$self->{YDIM} * 2 - 1; $j+=2) {
				for ($i=0;$i<$self->{XDIM} * 2 - 1; $i+=2) {
					if($i>0 and $j>0 and $i<$self->{XDIM} * 2 - 2 and $j<$self->{YDIM} * 2 - 2) {
					# in the middle of the map
						$medtable[0] = $umat[$i-1+$j*($self->{XDIM}*2-1)];
						$medtable[1] = $umat[$i+1+$j*($self->{XDIM}*2-1)];
						$medtable[2] = $umat[$i+($j-1)*($self->{XDIM}*2-1)];
						$medtable[3] = $umat[$i+($j+1)*($self->{XDIM}*2-1)];
						$#medtable = 3;
						@medtable = sort{$a <=> $b} @medtable;
						# Actually mean of two median values
						$umat[$i+$j*($self->{XDIM}*2-1)]=($medtable[1]+$medtable[2])/2.0;
					}
					elsif($j==0 and $i>0 and $i<$self->{XDIM} * 2 - 2) {
					# in the upper edge
						$medtable[0]=$umat[$i-1+$j*($self->{XDIM}*2-1)];
						$medtable[1]=$umat[$i+1+$j*($self->{XDIM}*2-1)];
						$medtable[2]=$umat[$i+($j+1)*($self->{XDIM}*2-1)];
						$#medtable = 2;
						@medtable = sort{$a<=>$b} @medtable;
						$umat[$i+$j*($self->{XDIM}*2-1)]=$medtable[1];
					}
					elsif($j==$self->{YDIM} * 2 - 2 and $i>0 and $i<$self->{XDIM} * 2 - 2) {
					# in the lower edge
						$medtable[0]=$umat[$i-1+$j*($self->{XDIM}*2-1)];
						$medtable[1]=$umat[$i+1+$j*($self->{XDIM}*2-1)];
						$medtable[2]=$umat[$i+($j-1)*($self->{XDIM}*2-1)];
						$#medtable = 2;
						@medtable = sort{$a<=>$b} @medtable;
						$umat[$i+$j*($self->{XDIM}*2-1)]=$medtable[1];
					}
					elsif($i==0 and $j>0 and $j<$self->{YDIM} * 2 - 2) {
					# in the left edge
						$medtable[0]=$umat[$i+1+$j*($self->{XDIM}*2-1)];
						$medtable[1]=$umat[$i+($j-1)*($self->{XDIM}*2-1)];
						$medtable[2]=$umat[$i+($j+1)*($self->{XDIM}*2-1)];
						$#medtable = 2;
						@medtable = sort{$a<=>$b} @medtable;
						$umat[$i+$j*($self->{XDIM}*2-1)]=$medtable[1];
					}
					elsif($i==$self->{XDIM} * 2 - 2 and $j>0 and $j<$self->{YDIM} * 2 - 2) {
					# in the right edge
						$medtable[0]=$umat[$i-1+$j*($self->{XDIM}*2-1)];
						$medtable[1]=$umat[$i+($j-1)*($self->{XDIM}*2-1)];
						$medtable[2]=$umat[$i+($j+1)*($self->{XDIM}*2-1)];
						$#medtable = 2;
						@medtable = sort{$a<=>$b} @medtable;
						$umat[$i+$j*($self->{XDIM}*2-1)]=$medtable[1];
					}
					elsif($i==0 && $j==0) {
					# the upper left-hand corner
						$umat[$i+$j*($self->{XDIM}*2-1)]=($umat[$i+1+$j*($self->{XDIM}*2-1)]+$umat[$i+($j+1)*($self->{XDIM}*2-1)])/2.0;
					}
					elsif($i==$self->{XDIM} * 2 - 2 and $j==0) {
					# the upper right-hand corner
						$umat[$i+$j*($self->{XDIM}*2-1)]=($umat[$i-1+$j*($self->{XDIM}*2-1)]+$umat[$i+($j+1)*($self->{XDIM}*2-1)])/2.0;
					}
					elsif($i==0 and $j==$self->{YDIM} * 2 - 2) {
					# the lower left-hand corner
						$umat[$i+$j*($self->{XDIM}*2-1)]=($umat[$i+1+$j*($self->{XDIM}*2-1)]+$umat[$i+($j-1)*($self->{XDIM}*2-1)])/2.0;
					}
					elsif($i==$self->{XDIM} * 2 - 2 and $j==$self->{YDIM} * 2 - 2) {
					# the lower right-hand corner
						$umat[$i+$j*($self->{XDIM}*2-1)]=($umat[$i-1+$j*($self->{XDIM}*2-1)]+$umat[$i+($j-1)*($self->{XDIM}*2-1)])/2.0;
					}
				}
			}
		}
		else {
		# Hexagonal topology
			for ($j=0; $j<$self->{YDIM}*2-1; $j+=2) {
				for ($i=0; $i<$self->{XDIM}*2-1; $i+=2) {
					if($i>0 and $j>0 and $i<$self->{XDIM} * 2 - 2 and $j<$self->{YDIM} * 2 - 2) {
					# in the middle of the map
						$medtable[0]=$umat[$i-1+$j*($self->{XDIM}*2-1)];
						$medtable[1]=$umat[$i+1+$j*($self->{XDIM}*2-1)];
						if(!($j%4)) {
							$medtable[2]=$umat[$i-1+($j-1)*($self->{XDIM}*2-1)];
							$medtable[3]=$umat[$i+($j-1)*($self->{XDIM}*2-1)];
							$medtable[4]=$umat[$i-1+($j+1)*($self->{XDIM}*2-1)];
							$medtable[5]=$umat[$i+($j+1)*($self->{XDIM}*2-1)];
						}
						else {
							$medtable[2]=$umat[$i+($j-1)*($self->{XDIM}*2-1)];
							$medtable[3]=$umat[$i+1+($j-1)*($self->{XDIM}*2-1)];
							$medtable[4]=$umat[$i+($j+1)*($self->{XDIM}*2-1)];
							$medtable[5]=$umat[$i+1+($j+1)*($self->{XDIM}*2-1)];
						}
						$#medtable = 5;
						@medtable = sort{$a<=>$b} @medtable;
						# Actually mean of two median values
						$umat[$i+$j*($self->{XDIM}*2-1)]=($medtable[2]+$medtable[3])/2;
					}
					elsif($j==0 and $i>0 and $i<$self->{XDIM} * 2 - 2) {
					# in the upper edge
						$medtable[0]=$umat[$i-1+$j*($self->{XDIM}*2-1)];
						$medtable[1]=$umat[$i+1+$j*($self->{XDIM}*2-1)];
						$medtable[2]=$umat[$i+($j+1)*($self->{XDIM}*2-1)];
						$medtable[3]=$umat[$i-1+($j+1)*($self->{XDIM}*2-1)];
						$#medtable = 3;
						@medtable = sort{$a<=>$b} @medtable;
						# Actually mean of two median values
						$umat[$i+$j*($self->{XDIM}*2-1)]=($medtable[1]+$medtable[2])/2;
					}
					elsif($j==$self->{YDIM} * 2 - 2 and $i>0 and $i<$self->{XDIM} * 2 - 2) {
					# in the lower edge
						$medtable[0]=$umat[$i-1+$j*($self->{XDIM}*2-1)];
						$medtable[1]=$umat[$i+1+$j*($self->{XDIM}*2-1)];
						if(!($j%4)) {
							$medtable[2]=$umat[$i-1+($j-1)*($self->{XDIM}*2-1)];
							$medtable[3]=$umat[$i+($j-1)*($self->{XDIM}*2-1)];
						}
						else {
							$medtable[2]=$umat[$i+($j-1)*($self->{XDIM}*2-1)];
							$medtable[3]=$umat[$i+1+($j-1)*($self->{XDIM}*2-1)];
						}
						$#medtable = 3;
						@medtable = sort{$a<=>$b} @medtable;
						# Actually mean of two median values
						$umat[$i+$j*($self->{XDIM}*2-1)]=($medtable[1]+$medtable[2])/2;
					}
					elsif($i==0 and $j>0 and $j<$self->{YDIM} * 2 - 2) {
					# in the left edge
						$medtable[0]=$umat[$i+1+$j*($self->{XDIM}*2-1)];
						if(!($j%4)) {
							$medtable[1]=$umat[$i+($j-1)*($self->{XDIM}*2-1)];
							$medtable[2]=$umat[$i+($j+1)*($self->{XDIM}*2-1)];
							$#medtable = 2;
							@medtable = sort{$a<=>$b} @medtable;
							$umat[$i+$j*($self->{XDIM}*2-1)]=$medtable[1];
						}
						else {
							$medtable[1]=$umat[$i+($j-1)*($self->{XDIM}*2-1)];
							$medtable[2]=$umat[$i+1+($j-1)*($self->{XDIM}*2-1)];
							$medtable[3]=$umat[$i+($j+1)*($self->{XDIM}*2-1)];
							$medtable[4]=$umat[$i+1+($j+1)*($self->{XDIM}*2-1)];
							$#medtable = 4;
							@medtable = sort{$a<=>$b} @medtable;
							$umat[$i+$j*($self->{XDIM}*2-1)]=$medtable[2];
						}
					}
					elsif($i==$self->{XDIM} * 2 - 2 and $j>0 and $j<$self->{YDIM} * 2 - 2) {
					# in the right edge
						$medtable[0]=$umat[$i-1+$j*($self->{XDIM}*2-1)];
						if($j%4) {
							$medtable[1]=$umat[$i+($j-1)*($self->{XDIM}*2-1)];
							$medtable[2]=$umat[$i+($j+1)*($self->{XDIM}*2-1)];
							$#medtable = 2;
							@medtable = sort{$a<=>$b} @medtable;
							$umat[$i+$j*($self->{XDIM}*2-1)]=$medtable[1];
						}
						else {
							$medtable[1]=$umat[$i+($j-1)*($self->{XDIM}*2-1)];
							$medtable[2]=$umat[$i-1+($j-1)*($self->{XDIM}*2-1)];
							$medtable[3]=$umat[$i+($j+1)*($self->{XDIM}*2-1)];
							$medtable[4]=$umat[$i-1+($j+1)*($self->{XDIM}*2-1)];
							$#medtable = 4;
							@medtable = sort{$a<=>$b} @medtable;
							$umat[$i+$j*($self->{XDIM}*2-1)]=$medtable[2];
						}
					}
					elsif($i==0 and $j==0) {
					# the upper left-hand corner
						$umat[$i+$j*($self->{XDIM}*2-1)]=($umat[$i+1+$j*($self->{XDIM}*2-1)]+$umat[$i+($j+1)*($self->{XDIM}*2-1)])/2.0;
					}
					elsif($i==$self->{XDIM} * 2 - 2 and $j==0) {
					# the upper right-hand corner
						$medtable[0]=$umat[$i-1+$j*($self->{XDIM}*2-1)];
						$medtable[1]=$umat[$i-1+($j+1)*($self->{XDIM}*2-1)];
						$medtable[2]=$umat[$i+($j+1)*($self->{XDIM}*2-1)];
						$#medtable = 2;
						@medtable = sort{$a<=>$b} @medtable;
						$umat[$i+$j*($self->{XDIM}*2-1)]=$medtable[1];
					}
					elsif($i==0 and $j==$self->{YDIM} * 2 - 2) {
					# the lower left-hand corner
						if(!($j%4)) {
							$umat[$i+$j*($self->{XDIM}*2-1)]=($umat[$i+1+$j*($self->{XDIM}*2-1)]+$umat[$i+($j-1)*($self->{XDIM}*2-1)])/2.0;
						}
						else {
							$medtable[0]=$umat[$i+1+$j*($self->{XDIM}*2-1)];
							$medtable[1]=$umat[$i+($j-1)*($self->{XDIM}*2-1)];
							$medtable[2]=$umat[$i+1+($j-1)*($self->{XDIM}*2-1)];
							$#medtable = 2;
							@medtable = sort{$a<=>$b} @medtable;
							$umat[$i+$j*($self->{XDIM}*2-1)]=$medtable[1];
						}
					}
					elsif($i==$self->{XDIM} * 2 - 2 and $j==$self->{YDIM} * 2 - 2) {
					# the lower right-hand corner
						if($j%4) {
							$umat[$i+$j*($self->{XDIM}*2-1)]=($umat[$i-1+$j*($self->{XDIM}*2-1)]+$umat[$i+($j-1)*($self->{XDIM}*2-1)])/2.0;
						}
						else {
							$medtable[0]=$umat[$i-1+$j*($self->{XDIM}*2-1)];
							$medtable[1]=$umat[$i+($j-1)*($self->{XDIM}*2-1)];
							$medtable[2]=$umat[$i-1+($j-1)*($self->{XDIM}*2-1)];
							$#medtable = 2;
							@medtable = sort{$a<=>$b} @medtable;
							$umat[$i+$j*($self->{XDIM}*2-1)]=$medtable[1];
						}
					}
				}
			}
		}

		# scale values to (0..1)

		my @umat_sort = sort{$a<=>$b} @umat;
		$bw = $umat_sort[$#umat_sort] - $umat_sort[0];
		$min = $umat_sort[0];
		for $i (0..$self->{XDIM} * 2 - 2) {
			for $j (0..$self->{YDIM} * 2 - 2) {
				$umat[$i+$j*($self->{XDIM}*2-1)] = ($umat[$i+$j*($self->{XDIM}*2-1)]-$min)/$bw;
			}
		}

		return \@umat;
	}

	##################################################
	#
	#
	#
	##################################################

	sub x_dim {
		my $self = shift;
		return $self->{XDIM};
	}

	sub y_dim {
		my $self = shift;
		return $self->{YDIM};
	}

	sub i_dim {
		my $self = shift;
		return $self->{IDIM};
	}

	sub topology {
		my $self = shift;
		return $self->{TOPOLOGY};
	}

	sub neighborhood {
		my $self = shift;
		return $self->{NEIGHBORHOOD};
	}

	sub map {
		my $self = shift;
		my ($x, $y, $z) = @_;
		return $self->{MAP}->[($y*$self->{XDIM} + $x) * $self->{IDIM} + $z];
	}

	sub label {
		my $self = shift;
		my ($x, $y) = @_;
		croak "Invalid argument" if (($x<0) or ($x>=$self->{XDIM}) or ($y<0) or ($y>=$self->{YDIM}));
		return $self->{LABELS}->[$y * $self->{XDIM} + $x];
	}

################################################################################
#
# "Private" methods.
#
################################################################################

	sub _find_two_eigenvectors_and_mean {
		my $self = shift;
		my $data = shift;

		my $k = 0;
		my $i = 0;
		my $j = 0;

		my $n = $self->{IDIM};

		my @r = ();
		my @m = ();
		my @u = ();
		my @v = ();
		my @k2 = ();

		my @mu = ();
		my $patterns_count = @$data / $n;

		for ($k=0; $k<$patterns_count; $k++) {
			for ($i=0; $i<$n; $i++) {
				if ($data->[$i + $k * $n] =~ /^-?\d+\.?\d*$/) { # Is a real number
					$m[$i] += $data->[$i + $k * $n];
					$k2[$i]++;
				}
			}
		}

		$i = 0;
		foreach my $k2item (@k2) {
			$m[$i++] /= $k2item;
		}

		for ($k=0; $k < $patterns_count; $k++) {
			for ($i=0; $i<$n; $i++) {
				if ($data->[$k * $n + $i] =~ /^-?\d+\.?\d*$/) { # Is a real number
					for ($j=$i; $j<$n; $j++) {
						if ($data->[$k * $n + $j] =~ /^-?\d+\.?\d*$/) { # Is a real number
							$r[$i * $n + $j] += ($data->[$k * $n + $i] - $m[$i]) * ($data->[$k * $n + $j] - $m[$j]);
						}
					}
				}
			}
		}

		for ($i=0; $i<$n; $i++) {
			for ($j=$i; $j<$n; $j++) {
				$r[$j * $n + $i]=$r[$i * $n + $j] /= $patterns_count;
			}
		}

		for ($i=0; $i<2; $i++) {
			for ($j=0; $j<$n; $j++) {
				$u[$i*$n+$j] = rand(2) - 1.0;
			}

			$self->_normalize(\@u, $i * $n);
			$mu[$i] = 1.0;
		}

		for ($k=0; $k<10; $k++) {
			for ($i=0; $i<2; $i++) {
				for ($j=0; $j<$n; $j++) {
					$v[$i * $n + $j] = $mu[$i] * $self->_inner_product(\@r, $j*$n, \@u, $i*$n, $n) + $u[$i * $n + $j];
				}
			}

			$self->_gram_schmidt(\@v, $n, 2);
			my $sum = 0.0;

			for ($i=0; $i<2; $i++) {
				for ($j=0; $j<$n; $j++) {
					$sum += abs($v[$i * $n + $j] / $self->_inner_product(\@r, $j*$n, \@v, $i*$n, $n));
				}
				$mu[$i] = $sum / $n;
			}
			@u = @v;
		}

		my @eigen1 = ();
		my @eigen2 = ();

    	for ($j=0; $j<$n; $j++) {
			$eigen1[$j] = $u[$j] / sqrt($mu[0]);
		}
	   for ($j=0; $j<$n; $j++) {
			$eigen2[$j] = $u[$j+$n] / sqrt($mu[1]);
		}

		return \(@eigen1, @eigen2, @m)
	}

	sub _inner_product {
		my $self = shift;
		my ($data1, $start_index1, $data2, $start_index2, $count) = @_;

		my $sum = 0;
		for (my $i=0; $i<$count; $i++) {
			$sum += $data1->[$start_index1++] * $data2->[$start_index2++]
		}
		return $sum;
	}

	sub _normalize {
		my $self = shift;
		my $data = shift;
		my $start_index = shift;

		my $sum = sqrt($self->_inner_product($data, $start_index, $data, $start_index, $self->{IDIM}));

		for (my $i=$start_index; $i<$start_index + $self->{IDIM}; $i++) {
			$data->[$i] /= $sum;
		}
	}

	sub _gram_schmidt {
		my $self = shift;
		my $data = shift;
		my $n = shift;
		my $e = shift;

		my @w = ();
		my $sum = 0;

		for (my $i=0; $i<$e; $i++) {
			for (my $t=0; $t<$n; $t++) {
				$sum = $data->[$i * $n + $t];
					for (my $j=0; $j<$i; $j++) {
						for (my $p=0; $p<$n; $p++) {
							$sum -= $w[$j * $n + $t] * $w[$j * $n + $p] * $data->[$i * $n + $p];
						}
					}
				$w[$i * $n + $t] = $sum;
			}
			$self->_normalize(\@w, $i * $n);
		}
		@$data = @w;
	}

	sub _get_hexa_dist {
		my $self = shift;
		my ($bx, $by, $tx, $ty) = @_;

		my $diff = $bx - $tx;

		if ((($by - $ty) % 2) != 0) {
			if (($by % 2) == 0) {
				$diff -= 0.5;
			}
			else {
				$diff += 0.5;
			}
		}

		my $temp = $diff ** 2;
		$diff = $by - $ty;
		$temp += 0.75 * $diff ** 2;
		return sqrt($temp);
	}

	sub _get_rect_dist {
		my $self = shift;
		my ($bx, $by, $tx, $ty) = @_;

		my $diff = $bx - $tx;
		my $temp = $diff ** 2;
		$diff = $by - $ty;
		$temp += $diff ** 2;
		return sqrt($temp);
	}

	sub _adapt_bubble {
		my $self = shift;
		my ($bx, $by, $radius, $alpha, $data) = @_;

		my $dist_func;
		if ($self->{TOPOLOGY} eq 'rect') {
			$dist_func = "_get_rect_dist";
		}
		else {
			$dist_func = "_get_hexa_dist";
		}

		for(my $x=0; $x<$self->{XDIM}; $x++) {
			for(my $y=0; $y<$self->{YDIM}; $y++) {
				if ($self->$dist_func($bx, $by, $x, $y) <= $radius) {
					for (my $i=0; $i<$self->{IDIM}; $i++) {
						if ($data->[$i] =~ /^-?\d+\.?\d*$/) { # Is a real number
							$self->{MAP}->[($y * $self->{XDIM} + $x) * $self->{IDIM} + $i] +=
								$alpha * ($data->[$i] - $self->{MAP}->[($y * $self->{XDIM} + $x) * $self->{IDIM} + $i]);
						}
					}
				}
			}
		}
	}

	sub _adapt_gaussian {
		my $self = shift;
		my ($bx, $by, $radius, $alpha, $data) = @_;

		my $dist_func;
		if ($self->{TOPOLOGY} eq 'rect') {
			$dist_func = "_get_rect_dist";
		}
		else {
			$dist_func = "_get_hexa_dist";
		}

		for(my $x=0; $x<$self->{XDIM}; $x++) {
			for(my $y=0; $y<$self->{YDIM}; $y++) {
				my $dd = $self->$dist_func($bx, $by, $x, $y);
				my $alp = $alpha * exp(-$dd ** 2 / (2.0 * $radius ** 2)); # -$dd**2 - ????
				for (my $i=0; $i<$self->{IDIM}; $i++) {
					if ($data->[$i] =~ /^-?\d+\.?\d*$/) { # Is a real number
						$self->{MAP}->[($y * $self->{XDIM} + $x) * $self->{IDIM} + $i] +=
							$alp * ($data->[$i] - $self->{MAP}->[($y * $self->{XDIM} + $x) * $self->{IDIM} + $i]);
					}
				}
			}
		}
	}

	sub _alpha_linear {
		my $self = shift;
		my ($iter, $epoches, $alpha) = @_;
		return $alpha * ($epoches - $iter) / $epoches;
	}

	sub _alpha_inverse_t {
		my $self = shift;
		my ($iter, $epoches, $alpha) = @_;
		my $c = $epoches / $AI::NeuralNet::SOM::INV_ALPHA_CONSTANT;
		return $alpha * $c /($c + $iter) / $epoches;
	}

1;

__END__