program svd
! find best matrix w and vector b in equation w*p +b = t
! where p is matrix of column vector of inputs
! t is matrix of column vector of targets
! (i suppose b should be a matrix of column vector of b)
! uses numerical recipes subroutines
! based on pseudo-inverse obtained form singular value decomposition
! pages 431-436 in book "learning from data"
use precis
use nrtype; use nrutil; use nr
use utilmod

implicit none

integer km,im,jm,iseed,irec
real(mrl), allocatable, dimension(:,:) :: wb,wbt,ps,p,t,pt,pn,u,v,w,s
real(mrl), allocatable, dimension(:) :: b,e,p1,t1,sig
integer, allocatable, dimension(:) :: irr
real(mrl), allocatable, dimension(:) :: arr
integer nread,nbatch,nalpha,verbose,init,ima(0:3),ifuna(3),idr
real(mrl) alpha0,alphai,avgerror,funcoa(3),dummy(1),rmstarget
character(80) filea,fileb,filec,netin,netout
character(1000) aline
namelist/params/ima,ifuna,funcoa,nread,nbatch,iseed,&
alpha0,nalpha,alphai,filea,fileb,filec,verbose,init,netin,netout

print *,' usage: svd < namelist_input_file > summary_output_file'
read (*,nml=params)
jm=ima(0)
im=ima(1)
if (ima(2)>0) im=ima(2)
if (ima(3)>0) im=ima(3)
km=jm+1

print *,'will "train" svd with ',trim(filea),','
open (unit=89,file=trim(filea),status='old')

irec=0
5 call read_input_target_file(89,dummy,dummy,aline,idr)
if (idr>2) goto 10
if (idr>1) goto 5
irec=irec+1
goto 5
10 print *,' counted irec=',irec

rewind (89)

allocate(ps(jm,irec),t(im,irec),pt(irec,km),pn(irec,km),p(km,irec))
allocate(w(im,jm),b(im),e(im),p1(jm),t1(im),wb(im,km),wbt(km,im))
allocate(s(km,km),u(irec,km),v(km,km),sig(km))

i=1
15 call read_input_target_file(89,ps(:,i),t(:,i),aline,idr)
if (idr>2) goto 20
if (idr>1) goto 15
p(1:jm,i)=ps(1:jm,i)
p(km,i)=1._mrl
i=i+1
goto 15
20 print *,' finished reading in p and t'

print *,'p:'
if (verbose >= 1) call printmat(p)

pt=transpose(p)
print *,'pt:'
if (verbose >1) call printmat(pt)

u=pt

!print *, u(i-1,:)
call svdcmp(u,sig,v)
!print *, u(i-1,:)
!stop

print *,'singular value decomp:'
if (verbose >= 1) call printmat(u)
if (verbose >= 1) call printcol(sig)
if (verbose >= 1) call printmat(v)

s=0._mrl
do i=1,km
s(i,i)=sig(i)
end do
if (verbose >= 1) call printmat(s)

pn=matmul(u,matmul(s,transpose(v)))
print *,'test svd reconstruction:'
print *, pn(1,:)-p(:,1)
if (verbose >= 1) call printmat(pn)

s=0._mrl
do i=1,km
s(i,i)=1./sig(i)
end do

print *, 'sig+:'
if (verbose >= 1) call printmat(s)

p=matmul(v,matmul(s,transpose(u)))
print *,'x+:'
if (verbose >= 1) call printmat(p)
wb=transpose(matmul(p,transpose(t)))

print *, 'wb'
if (verbose >= 1) call printmat(wb)

w=wb(:,1:jm)
b=wb(:,km)
print *, 'w:'
if (verbose >= 0) call printmat(w)
print *, 'b:'
if (verbose >= 0) call printcol(b)

close (unit=89)

call make_pred_file(filea,'svd_a.dat')
call make_pred_file(fileb,'svd_b.dat')
call make_pred_file(filec,'svd_c.dat')

print *,'sorted absolute values of first row of weights:'
allocate(arr(jm),irr(jm))
do i=1,jm
arr(i)=abs(w(1,i))
irr(i)=i
end do

call sort_down2(arr,irr)
do i=1,jm
print *,irr(i),arr(i)
end do

contains
subroutine make_pred_file(infile,outfile)
character(*) infile, outfile
print *,'will verify performance of svd with ',trim(infile),','
open (unit=89,file=trim(infile),status='old',err=101)
print *,'prediction vs. target is output in file ',trim(outfile)
open (unit=91,file=trim(outfile),status='unknown')
write (91,*) '# svd predicted vector, target vector for ',trim(infile)
write (91,*) '# ',trim(infile),' had these comments:'
avgerror=0.
rmstarget=0.
irec=0
25 call read_input_target_file(89,p1,t1,aline,idr)
if (idr>2) goto 30
if (idr == 2) then



write(91,'(a)') trim(aline)
else



irec=irec+1



e=matmul(w,p1)+b



write (91,'(30e17.8)') e,t1



e=e-t1



avgerror=avgerror+dot_product(e,e)



rmstarget=rmstarget+dot_product(t1,t1)
end if
goto 25
30 close (unit=89)
avgerror=avgerror/irec
rmstarget=rmstarget/irec
write (91,*) '#',trim(infile),' mse=',avgerror,' mst=',rmstarget
close (unit=91)

print *, ' average square error with ',trim(infile),' is mse=',avgerror
print *, ' average square of target in ',trim(infile),' is mst=',rmstarget
101 return
end subroutine make_pred_file

end program svd