function [Best_tour_length= tspsiman(EUC_2D) 


flops(0) ;
t0= clock ;
xy= EUC_2D ; 
n_cities= length(xy) ;
rand('state',sum(100*clock));%RAND('state',sum(100*clock)) resets the generator to a different state each time.
 

% Create the distance matrix for all of the cities given x-y coordinates
distance_matrix = zeros(n_cities) ;
for n_cities_x = 1: n_cities,
     for n_cities_y = 1:n_cities_x
          x = xy(n_cities_x, 1) ;
          y = xy(n_cities_x, 2) ;
          xx = xy(n_cities_y, 1) ;
          yy = xy(n_cities_y, 2) ;
          distance_matrix(n_cities_x, n_cities_y)= ceil(sqrt((x - xx)^2 + (y - yy)^2)) ;
          distance_matrix(n_cities_y, n_cities_x)= distance_matrix(n_cities_x, n_cities_y) ;
 	end
end 
% End of matrix construction

% Construct an initial tour using Nearest Neighbor heuristic 
lenbestNN= inf ;
pbestNN= [] ;

   prand= randperm(n_cities) ; % A random selection of the starting city
		f=find(prand==1) ;
		prand(f)= prand(1) ;
 	p= [1 prand(2)] ;
	i= prand(3) ;
	count= 3 ;
	 
	while count <= n_cities
     	NNdist= inf ;
     	pp= i ;
     	for j= 1: n_cities
          	if (distance_matrix(i, j) < NNdist) & (j~=i) & ((j~=p) == ones(1,length(p)))%333??????????????????
                NNdist= distance_matrix(i, j) ; 
                pp= j ;
          	end           
     	end
     	p= [p pp] ; 
      i= pp ;
     	count= count + 1 ;
	end
	% Computing tour cost or length using Tourdist.m function
	len= tourdist(p, distance_matrix) ;
 
	if len < lenbestNN
		lenbestNN= len ; 
		pbestNN= p ; 
	end 
% End of initial tour construction

solnn= [] ;
lenn= [] ; temp= [] ;
soln= 1 ;
% ========================
% A 2-Opt local search
% ========================
lencurr= lenbestNN; 
Best_tour_length= lenbestNN 
pcurr= pbestNN ; 
pbest= pbestNN ; 
 
% ========================
% Temperature control
% ========================
restart= 1 ;
Tstart= 30 ; % Start temperature
Tend= 1 ; % Stop temperature
Tred= 0.97 ;
T= Tstart ;
Nochange= 2 ; % If after Nochange neighborhood searches, no improvements or 
              %  changes in tour search, annealing complete, break search.
% ========================

lenn= [lenn lencurr] ;
temp= [temp T] ;
solnn= [solnn soln] ;

bb= 0 ;  

while T >= Tend  
	big= n_cities - 1 ; 
	while big >= 3   
		small= big - 2 ;
		while small >= 1
		 
			curropt= distance_matrix(pcurr(big),pcurr(big+1)) + distance_matrix(pcurr(small),pcurr(small+1)) ;  
	 		swap2= distance_matrix(pcurr(small),pcurr(big)) + distance_matrix(pcurr(small+1),pcurr(big+1)) ; 
  
 			soln= soln + 1 ;
	 	
			if swap2 < curropt
 				order2= 1: n_cities ;  
 				order2=[1:small big:-1:small+1 big+1:n_cities] ; 
 
	 			pcurr= pcurr(order2) ;
				lencurr= tourdist(pcurr, distance_matrix) ;
				 
				lenn= [lenn lencurr] ; 
				temp= [temp T] ; 
				solnn= [solnn soln] ; 
 
 		 		if lencurr < Best_tour_length
				  	Best_tour_length= lencurr       
					pbest= pcurr ;  
					Temperature_of_best_tour_length= T
					Solution_count= soln  
					T= Tred * T ; 
					if T <= 3
						T= 50 ;
					end
 				end
				Tcurr= T ;      
   			bb= 0 ;
				big= n_cities - 1 ; 
				small= big - 1 ;
				if T <= 3
					T= 10 ;
				end
 
				if T <= Tend ;
					big= 2.9 ; 
					break
 				end	
		
			elseif swap2 > curropt
            %r= abs(randn) ;
            r= rand; % where r ranges from 0.0 to 1.0 
            diff= swap2 - curropt ;
            %if r < exp(-(diff) / T)
        		if r <= exp(-(diff) / T)
					order2= 1: n_cities ;
					order2=[1:small big:-1:small+1 big+1:n_cities] ; 
 					pcurr= pcurr(order2) ;
           			lencurr= tourdist(pcurr, distance_matrix) ;
					
					T= Tred * T ; 
  					bb= 0 ;  
 	
				end
		 	end
     		small= small - 1 ;
		end
		big= big - 1 ;   
	end
	bb= bb + 1 ;   
 	if T <= Tend | bb > Nochange ;
 	
		clc
		Best_tour_length
		besttour= [pbest -pbest(1)]  
		Temperature_of_best_tour_length
		Solution_count
		Search_stop_temperature= T  
		Elapsed_time= etime(clock, t0) % In seconds
      Solutions_generated= soln 
      Floating_point_operations = flops
	 	if bb > Nochange
			No_change= bb
		end
		disp('Press ENTER to display plot (Temperature vs. Tour Length) or Ctrl^C to end search.') 
 		pause
		
		clc
		plot(temp, lenn)
		title('Simulated Annealing w/ 2-Opt local search')
		xlabel('Temperature (not scaled)')
		ylabel('Tour Lengths/Costs')
		grid
 	 	
		disp('Press ENTER to display plot (Number of Solutions vs. Tour Length) or Ctrl^C to end search.') 
 		pause

		clc
		plot(solnn, lenn)
		title('Simulated Annealing w/ 2-Opt local search')
		xlabel('Number of Solutions')
		ylabel('Tour Lengths/Costs')
		grid

		disp('Press ENTER to restart search (if var restart > 0) or Ctrl^C to end search.')
		pause
	 	if restart > 0
 
			clc
			T= Tstart ; bb= 0 ;  
			solnn= []; lenn= []; temp= [] ;
			
		 	% =======================================================
			% This time randomly generate tours and restart annealing
			% =======================================================
			prand= randperm(n_cities) ;
				f=find(prand==1) ;
				prand(f)= prand(1) ; prand(1)= 1 ;

		   lencurr= Best_tour_length 
         pcurr= pbest ;
         % =======================================================

 		end 
 	end
end	
% End of local search