program f90_MDS !******************************************** ! Multi Dimension Scaling !******************************************** implicit none integer::knor ! Index for normalization (0: distance matrix, 1: normalization) character(len=100)::strcom ! Comment integer::mcol ! Number of variables integer::ndata ! Number of data sets real(8),allocatable::xd(:,:) ! Array for input data real(8),allocatable::mat(:,:) ! Array for Young-Householder transformation real(8),allocatable::ary(:,:) ! Array for eigen value analysis real(8),allocatable::ev(:) ! Eigenvalue real(8),allocatable::vec(:,:) ! Eigenvector real(8),allocatable::zd(:,:) ! Score real(8),allocatable::tr ! Sum total of eigenvalues real(8),parameter::p=2 ! Dimension of distance (p=2: Euclidean distance) character(len=50),allocatable::dname(:) integer::i,j,k character(len=50)::fnameR,fnameW character(len=100)::dummy real(8)::sum call getarg(1,dummy);read(dummy,*) knor call getarg(2,fnameR) call getarg(3,fnameW) !Data input if(knor==0)then !Distance matrix open(11,file=fnameR,status='old') read(11,'(a)') strcom read(11,*) ndata mcol=ndata allocate(dname(1:ndata)) allocate(xd(1:ndata,1:mcol),zd(1:ndata,1:ndata)) allocate(ary(1:ndata,1:ndata),mat(1:ndata,1:ndata)) allocate(ev(1:ndata),vec(1:ndata,1:ndata)) do i=1,ndata do j=1,ndata ary(i,j)=0.0D0 end do end do do k=1,ndata read(11,*) dname(k),(ary(k,j),j=1,k) end do close(11) do k=1,ndata do j=1,i ary(k,j)=ary(j,k) end do end do end if if(1<=knor)then !Vector data open(11,file=fnameR,status='old') read(11,'(a)') strcom read(11,*) mcol,ndata allocate(xd(1:ndata,1:mcol),zd(1:ndata,1:ndata)) allocate(ary(1:ndata,1:ndata),mat(1:ndata,1:ndata)) allocate(ev(1:ndata),vec(1:ndata,1:ndata)) allocate(dname(1:ndata)) do k=1,ndata read(11,*) (xd(k,j),j=1,mcol) write(dname(k),*) 'Score(',k,')' end do close(11) !Normalization if(knor==1)call STAND(ndata,mcol,xd) !Calculation of distance do i=1,ndata do j=1,ndata ary(i,j)=0.0D0 end do end do !Minkowski model do i=1,ndata do j=i,ndata sum=0.0D0 do k=1,mcol sum=sum+(abs(xd(i,k)-xd(j,k)))**p end do sum=sum**(1.0D0/p) ary(i,j)=sum end do end do do i=1,ndata do j=1,i ary(i,j)=ary(j,i) end do end do end if !Squared distance do k=1,ndata do j=1,i ary(k,j)=ary(k,j)**2.0D0 end do end do !Centering matrix do i=1,ndata do j=1,ndata mat(i,j)=-1.0D0/dble(ndata) end do mat(i,i)=mat(i,i)+1.0D0 end do !Young-Housholder transformation do i=1,ndata do j=1,ndata zd(i,j)=0.0D0 do k=1,ndata zd(i,j)=zd(i,j)+ary(i,k)*mat(k,j) end do end do end do tr=0.0D0 do i=1,ndata do j=1,ndata ary(i,j)=0.0D0 do k=1,ndata ary(i,j)=ary(i,j)+mat(i,k)*zd(k,j) end do ary(i,j)=-0.5D0*ary(i,j) end do tr=tr+ary(i,i) end do !Eigenvalue, eigenvector call JACOBI(ndata,ary,ev,vec,ndata) do i=1,mcol write(6,'("ev(",i0,")=",e15.7)') i,ev(i) end do !Calculation of Coordinates do k=1,ndata do i=1,mcol if(0.0D0<=ev(i))then zd(k,i)=sqrt(ev(i))*vec(k,i) else zd(k,i)=0.0D0 end if end do end do !Print out call PROUT(fnameW,strcom,mcol,ndata,ev,tr,vec,dname,zd) stop contains subroutine PROUT(fnameW,strcom,mcol,ndata,ev,tr,vec,dname,zd) !Print out character(len=50),intent(in)::fnameW character(len=100),intent(in)::strcom integer,intent(in)::mcol integer,intent(in)::ndata real(8),intent(in)::ev(1:mcol) real(8),intent(in)::vec(1:mcol,1:mcol) real(8),intent(in)::zd(1:ndata,1:mcol) real(8),intent(in)::tr character(len=50),intent(in)::dname(1:ndata) character(len=100)::dummy character(len=1000)::linebuf open(12, file=fnameW, status='replace') write(12,'(a)') trim(adjustl(strcom)) write(linebuf,*) 'mcol=,',mcol call del_spaces(linebuf) write (12,'(a)') trim(linebuf) write(linebuf,*) 'ndata=,',ndata call del_spaces(linebuf) write (12,'(a)') trim(linebuf) write(dummy,*) 'No.' linebuf=trim(dummy) do j=1,mcol write(dummy,*) j linebuf=trim(linebuf)//','//trim(dummy) end do call del_spaces(linebuf) write(12,*) trim(linebuf) write(dummy,*) 'Eigen_value' linebuf=trim(dummy) do j=1,mcol write(dummy,*) ev(j) linebuf=trim(linebuf)//','//trim(dummy) end do call del_spaces(linebuf) write(12,*) trim(linebuf) write(dummy,*) 'Proportion' linebuf=trim(dummy) do j=1,mcol write(dummy,*) ev(j)/tr linebuf=trim(linebuf)//','//trim(dummy) end do call del_spaces(linebuf) write(12,*) trim(linebuf) do i=1,mcol write(dummy,*) 'Eigen_vec(',i,')' linebuf=trim(dummy) do j=1,mcol write(dummy,*) vec(i,j) linebuf=trim(linebuf)//','//trim(dummy) end do call del_spaces(linebuf) write(12,*) trim(linebuf) end do do i=1,ndata write(dummy,*) dname(i) linebuf=trim(dummy) do j=1,mcol write(dummy,*) zd(i,j) linebuf=trim(linebuf)//','//trim(dummy) end do call del_spaces(linebuf) write(12,*) trim(linebuf) end do close(12) end subroutine PROUT !************************************************** ! Normalization !************************************************** subroutine STAND(ndata,mcol,xd) integer,intent(in)::ndata,mcol real(8),intent(inout)::xd(1:ndata,1:mcol) integer::i,k real(8)::xm(1:mcol),sd(1:mcol) !Average do i=1,mcol xm(i)=0.0D0 do k=1,ndata xm(i)=xm(i)+xd(k,i) end do xm(i)=xm(i)/dble(ndata) end do !Standard deviation do i=1,mcol sd(i)=0.0D0 do k=1,ndata sd(i)=sd(i)+(xd(k,i)-xm(i))*(xd(k,i)-xm(i)) end do sd(i)=sqrt(sd(i)/dble(ndata-1)) end do !Normalization do i=1,mcol do k=1,ndata xd(k,i)=(xd(k,i)-xm(i))/sd(i) end do end do end subroutine STAND !************************************************** ! Jacobi method for eigen value analysis !************************************************** subroutine JACOBI(n,array,ev,vec,nt) !Jacobi method by Dr.Okumura !Eigenvalue: ev(k), eigenvector: vec(k,i) integer,intent(in)::n,nt real(8),intent(inout)::array(1:nt,1:nt) real(8),intent(out)::ev(1:nt),vec(1:nt,1:nt) integer::i,j,k,im real(8)::s,c,aoff,eps,t real(8)::aa,aj,ak,vj,vk,vmax,bb do i=1,n do j=1,n vec(i,j)=0.0 end do vec(i,i)=1.0 end do aoff=0.0 do j=1,n-1 do k=j+1,n aoff=aoff+array(j,k)*array(j,k) end do end do eps=aoff*2.0 do j=1,n eps=eps+array(j,j)*array(j,j) end do eps=eps*0.00000000005 do while(aoff>eps) do j=1,n-1 do k=j+1,n if(array(j,k)==0.0)then t=0.0 else aa=(array(k,k)-array(j,j))/(2.0*array(j,k)) end if if(aa>=0)then t=1.0/(aa+sqrt(1.0+aa*aa)) else t=1.0/(aa-sqrt(1.0+aa*aa)) end if c=1.0/sqrt(1.0+t*t) s=t*c array(j,j)=array(j,j)-t*array(j,k) array(k,k)=array(k,k)+t*array(j,k) array(j,k)=0.0 do i=1,j-1 aj=array(i,j) ak=array(i,k) array(i,j)=aj*c-ak*s array(i,k)=aj*s+ak*c end do do i=j+1,k-1 aj=array(j,i) ak=array(i,k) array(j,i)=aj*c-ak*s array(i,k)=aj*s+ak*c end do do i=k+1,n aj=array(j,i) ak=array(k,i) array(j,i)=aj*c-ak*s array(k,i)=aj*s+ak*c end do do i=1,n vj=vec(i,j) vk=vec(i,k) vec(i,j)=vj*c-vk*s vec(i,k)=vj*s+vk*c end do end do end do aoff=0.0 do j=1,n-1 do k=j+1,n aoff=aoff+array(j,k)*array(j,k) end do end do end do do k=1,n ev(k)=array(k,k) end do do k=1,n-1 vmax=-1.0E+30 do i=k,n if(ev(i)>=vmax)then im=i vmax=ev(i) end if end do aa=ev(k) ev(k)=ev(im) ev(im)=aa do i=1,n bb=vec(i,k) vec(i,k)=vec(i,im) vec(i,im)=bb end do end do endsubroutine JACOBI !************************************************** ! Deleate space !************************************************** subroutine del_spaces(s) character (*), intent (inout) :: s character (len=len(s)) tmp integer i, j j = 1 do i = 1, len(s) if (s(i:i)==' ') cycle tmp(j:j) = s(i:i) j = j + 1 end do s = tmp(1:j-1) end subroutine del_spaces end program f90_MDS