evaluate_2lsom Subroutine

   subroutine evaluate_2lsom(kohonen_map,input_data,results)
   !========================================================================================
!!   Subroutine to calculate some clustering statistics of a two-level self_organized_map 
      class(two_level_self_organizing_map) :: kohonen_map
!! A `two_level_self_organizing_map` object
      type(kohonen_pattern),dimension(:),intent(inout) :: input_data
!! A `kohonen_pattern` array
      real(kind=wp),dimension(:),optional :: results
!! A real array
      integer :: ipattern1,ipattern2,current_cluster,ic,is,pos
      integer,dimension(size(input_data)) :: indicator,positions
      integer,dimension(kohonen_map%number_clusters,size(input_data)) :: sample_positions
      integer,dimension(kohonen_map%number_clusters) :: number_samples_cluster
      type(kohonen_pattern) :: current_pattern1,current_pattern2
      type(kohonen_prototype) :: current_prototype1,current_prototype2
      real(kind=wp),dimension(kohonen_map%number_variables1,kohonen_map%number_variables2) :: &
                           current_values1,current_values2
      real(kind=wp) :: current_dissimilarity
      real(kind=wp),dimension(size(input_data)) :: mean_dissimilarity_a,silhouette,min_b,W,B
      real(kind=wp),dimension(size(input_data),kohonen_map%number_clusters) :: mean_dissimilarity_b
   !  Silouette or whatever
   !  find samples in each cluster
      min_b=10.0d8;
      B=0.0_wp;W=0.0_wp;
      do is=1,size(input_data) 
         positions(is)=is;
      enddo
   !   
      do ic=1,kohonen_map%number_clusters
         indicator=0;
         pos=0;
         where(kohonen_map%cluster_samples .eq. ic)
         indicator=1;        
         end where
         number_samples_cluster(ic)=sum(indicator);
         if(number_samples_cluster(ic) .eq. 0) then
         write(6,*) 'WARNING: Empty cluster. No cluster evaluation is done'
         if(present(results)) then
            results=0.0_wp;
         endif
         return
         endif
         do ipattern1=1,size(indicator)
            if(indicator(ipattern1) .eq. 1) then
            pos=pos+1;
            sample_positions(ic,pos)=ipattern1;
            endif
         enddo
      enddo!ic
   !
      do ipattern1=1,size(input_data)
         current_cluster=kohonen_map%cluster_samples(ipattern1);
   ! get current prototype
         current_pattern1=input_data(ipattern1);
         call current_pattern1%get(current_prototype1);
         call current_prototype1%get_prototype(current_values1);
         current_dissimilarity=0.0_wp;
         do ipattern2=1,number_samples_cluster(current_cluster);
            if (sample_positions(current_cluster,ipattern2) .ne. ipattern1) then
               current_pattern2=input_data(ipattern2);
               call current_pattern2%get(current_prototype2);
               call current_prototype2%get_prototype(current_values2);
   !	    if(kohonen_map%number_variables1 .eq. 1 .or. kohonen_map%number_variables2 .eq. 1 ) then
               current_dissimilarity=current_dissimilarity+sum((current_values1-current_values2)**2);
   !            else 
   !              current_dissimilarity=current_dissimilarity+sum(sum((current_values1-current_values2)**2));
   !            endif
            endif
         enddo!ipattern2
         mean_dissimilarity_a(ipattern1)=current_dissimilarity/max(1.0,real(number_samples_cluster(current_cluster)-1));
         W(ipattern1)=current_dissimilarity;
         !
         do ic=1,kohonen_map%number_clusters
            current_dissimilarity=0.0_wp;
            if(ic .ne. current_cluster) then
               do ipattern2=1,number_samples_cluster(ic);
                  current_pattern2=input_data(ipattern2);
                  call current_pattern2%get(current_prototype2);
                  call current_prototype2%get_prototype(current_values2);
   !               if(kohonen_map%number_variables1 .eq. 1 .or. kohonen_map%number_variables2 .eq. 1 ) then
                  current_dissimilarity=current_dissimilarity+sum((current_values1-current_values2)**2);
   !               else
   !                 current_dissimilarity=current_dissimilarity+sum(sum((current_values1-current_values2)**2));
   !               endif
               enddo!ipattern2
            endif
            mean_dissimilarity_b(ipattern1,ic)=current_dissimilarity/max(1.0,real(number_samples_cluster(ic)-1));       
            B(ipattern1)=B(ipattern1)+current_dissimilarity;
         enddo!ic
         
         !
         do ic=1,kohonen_map%number_clusters
            if(ic .ne. current_cluster .and. mean_dissimilarity_b(ipattern1,ic) .lt. min_b(ipattern1) ) then
               min_b(ipattern1)=mean_dissimilarity_b(ipattern1,ic)
            endif
         enddo
         !
         if(mean_dissimilarity_a(ipattern1) .lt. min_b(ipattern1)) then
         silhouette(ipattern1)=1.0_wp-(mean_dissimilarity_a(ipattern1)/min_b(ipattern1));
         else
         silhouette(ipattern1)=(min_b(ipattern1)/mean_dissimilarity_a(ipattern1))-1.0_wp;
         endif
         !write(6,*) ipattern1,mean_dissimilarity_a(ipattern1),min_b(ipattern1),silhouette(ipattern1) !(mean_dissimilarity_b(ipattern1,ic),ic=1,3);
      enddo!ipattern1
      write(6,*) ''
      write(6,*) 'CH= ',(sum(B)/dble(kohonen_map%number_clusters-1))&
                  /(sum(W)/(dble(size(input_data)-kohonen_map%number_clusters))); !,sum(B),sum(W),
      write(6,*) 'Silhouette= ',sum(silhouette)/real(size(silhouette))
      write(6,*) ''
   !
      if(present(results)) then
      results(1)=sum(B);
      results(2)=sum(W);
      results(3)=sum(silhouette)/real(size(silhouette));
      endif
   !
   end subroutine evaluate_2lsom