calculate_sum2_clusters_grid Subroutine

   subroutine calculate_sum2_clusters_grid(kohonen_map,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
   real(kind=wp),dimension(:),optional :: results
!! A real array
   integer :: ix1,iy1,iz1,ic,current_cluster,j,ipos,pos,isample
   integer :: nx,ny,nz,cix,ciy,ciz,current_pos,sample_pos
   real(kind=wp),dimension(kohonen_map%number_variables1,&
                  kohonen_map%number_variables2) :: mean_unit
   real(kind=wp),dimension(kohonen_map%number_clusters) :: W,B
   real(kind=wp) :: current_dissimilarity
   real(kind=wp),dimension(kohonen_map%number_variables1*&
                  kohonen_map%number_variables2,1) :: current_values1,current_values2
   type(kohonen_prototype) :: current_prototype1,current_prototype2
   integer,dimension(kohonen_map%number_nodes) :: indicator
   integer,dimension(kohonen_map%number_clusters) :: number_samples_cluster
   integer, dimension(kohonen_map%number_nodes) ::  positions
   integer,dimension(kohonen_map%number_clusters,kohonen_map%number_nodes) :: sample_positions
   real(kind=wp),dimension(kohonen_map%number_nodes) :: mean_dissimilarity_a,min_b,silhouette
   real(kind=wp),dimension(kohonen_map%number_nodes,kohonen_map%number_clusters) :: mean_dissimilarity_b
   
   !  
   positions(1:kohonen_map%number_nodes)=(/(ipos,ipos=1,kohonen_map%number_nodes)/)
   !
   nx=kohonen_map%parameters(1)%number_nodes_nx;
   ny=kohonen_map%parameters(1)%number_nodes_ny;
   nz=kohonen_map%parameters(1)%number_nodes_nz;
   !  find samples in each cluster
      min_b=10.0d8;
   !   
      B=0.0_wp;W=0.0_wp;mean_unit=0.0_wp
   !
   do iz1=1,size(kohonen_map%grid,3);
      do iy1=1,size(kohonen_map%grid,2);
         do ix1=1,size(kohonen_map%grid,1);
            current_cluster=kohonen_map%grid_cluster(ix1,iy1,iz1);
            current_prototype1=kohonen_map%grid(ix1,iy1,iz1);
            call current_prototype1%get_prototype(current_values1);
            current_prototype2=kohonen_map%cluster_layer(current_cluster);
            call current_prototype2%get_prototype(current_values2)              
            W(current_cluster)=W(current_cluster)+sum((current_values1-current_values2)**2);
            mean_unit=mean_unit+current_values1;
         enddo!ix
      enddo!iy
   enddo!iz
   mean_unit=mean_unit/dble(nx*ny*nz);
   !
   do ic=1,kohonen_map%number_clusters
      current_prototype1=kohonen_map%cluster_layer(ic);
      call current_prototype1%get_prototype(current_values1);
      B(ic)=B(ic)+sum((current_values1-mean_unit)**2);
   enddo!ic
   !
   if(present(results)) then
      results(1)=sum(W);results(2)=sum(B);
   !  else 
   !     write(6,*) 'Number clusters,W,B= ',kohonen_map%number_clusters,sum(W),sum(B)
   endif
   !
   ! Silhouette for grid layer
   !
   do ic=1,kohonen_map%number_clusters
      indicator=0;
      where(kohonen_map%cluster_cells_index(:,4) == ic)
         indicator=1;
      end where
      number_samples_cluster(ic)=sum(indicator);
      pos=0;
      do isample=1,size(indicator)
         if(indicator(isample) .eq. 1) then
            pos=pos+1;
            sample_positions(ic,pos)=isample;
         endif
      enddo
      !write(6,*) (sample_positions(ic,isample),isample=1,number_samples_cluster(ic));
   enddo
   !
   current_pos=0;
   do iz1=1,size(kohonen_map%grid,3);
      do iy1=1,size(kohonen_map%grid,2);
         do ix1=1,size(kohonen_map%grid,1);
            current_pos=current_pos+1;
            current_cluster=kohonen_map%grid_cluster(ix1,iy1,iz1);
            current_prototype1=kohonen_map%grid(ix1,iy1,iz1);
            call current_prototype1%get_prototype(current_values1);
            current_dissimilarity=0.0_wp;
            do ic=1,number_samples_cluster(current_cluster);
               if(sample_positions(current_cluster,ic) .ne. current_pos ) then
                  sample_pos=sample_positions(current_cluster,ic);
                  cix=kohonen_map%cluster_cells_index(sample_pos,1);
                  ciy=kohonen_map%cluster_cells_index(sample_pos,2);
                  ciz=kohonen_map%cluster_cells_index(sample_pos,3);
                  current_prototype2=kohonen_map%grid(cix,ciy,ciz);
                  call current_prototype2%get_prototype(current_values2);
                  current_dissimilarity=current_dissimilarity+sum((current_values1-current_values2)**2);
               endif
            enddo
            mean_dissimilarity_a(current_pos)=current_dissimilarity/dble(number_samples_cluster(current_cluster));
   ! 
            do ic=1,kohonen_map%number_clusters
               if(ic /= current_cluster) then
                  current_dissimilarity=0.0_wp;
                  do isample=1,number_samples_cluster(ic);
                     sample_pos=sample_positions(ic,isample);
                     cix=kohonen_map%cluster_cells_index(sample_pos,1);
                     ciy=kohonen_map%cluster_cells_index(sample_pos,2);
                     ciz=kohonen_map%cluster_cells_index(sample_pos,3);
                     current_prototype2=kohonen_map%grid(cix,ciy,ciz);
                     call current_prototype2%get_prototype(current_values2);
                     current_dissimilarity=current_dissimilarity+sum((current_values1-current_values2)**2);
                  enddo!isample
                  mean_dissimilarity_b(current_pos,ic)=current_dissimilarity/dble(number_samples_cluster(ic));
               endif              
            enddo!ic
            !
            do ic=1,kohonen_map%number_clusters
               if(ic .ne. current_cluster .and. mean_dissimilarity_b(current_pos,ic) .lt. min_b(current_pos) ) then
                  min_b(current_pos)=mean_dissimilarity_b(current_pos,ic)
               endif
            enddo
   !
            if(mean_dissimilarity_a(current_pos) .lt. min_b(current_pos)) then
               silhouette(current_pos)=1.0_wp-(mean_dissimilarity_a(current_pos)/min_b(current_pos));
            else
               silhouette(current_pos)=(min_b(current_pos)/mean_dissimilarity_a(current_pos))-1.0_wp;
            endif
   !
         enddo!ix1
      enddo!iy1
   enddo!iz1
   write(*,*) 'Number clusters,W,B,sil= ',kohonen_map%number_clusters,sum(W),sum(B),sum(silhouette)/dble(kohonen_map%number_nodes);
   if(present(results)) then
      results(3)=sum(silhouette)/dble(kohonen_map%number_nodes);
   endif
   !  write(*,*) 'b',min_b !mean_dissimilarity_b
   !  write(*,*) 'a',mean_dissimilarity_a
   !  write(*,*) 'a,b',sum(mean_dissimilarity_a)/dble(kohonen_map%number_clusters),&
   !                   sum(mean_dissimilarity_b)/dble(kohonen_map%number_clusters)
   !
   !
   end subroutine calculate_sum2_clusters_grid