text2.f90

函数逼近和预测

program Bp_Neutral_Net 
	character*1 cc
	integer m_samplenum,m_testnum
	common/bp/m_inputnum,m_hidenum,m_outnum,m_samplenum
	common/nn/n2
	common/ee/ee,es1,es		
	common/weight/w_ih,w_ho,b,a
	common/sample/sample_in,sample_t,sample_out,sample_x1
    double precision x(30),x1(30),y(30)
	double precision w_ih(30,30),w_ho(30,30),b(30,1),a(30,1)
	double precision sample_in(200,30),sample_t(200,30),sample_out(200,30),sample_x1(200,30)
    ee=1.0/2000.0
	es1=0.01
	es=-0.01
	write(*,*)'Do you want to learn(l) or calculate(c)?'
	read(*,*)cc
	if(cc.eq.'l'.or.cc.eq.'L') then
	call bp_ini	
	call readsampledata
	call Learning
	open(2,file="quanzhijieguo.txt")
	write(2,*)n2
	write(2,*)m_inputnum,m_hidenum,m_outnum
	write(2,*)((w_ih(i,j),j=1,m_inputnum),b(i,1),a(i,1),i=1,m_hidenum)
	write(2,*)((w_ho(i,j),j=1,m_hidenum+1),i=1,m_outnum)
	close(2)
	elseif(cc.eq.'c'.or.cc.eq.'C') then
	open(2,file="quanzhijieguo.txt")
	read(2,*)n2
	read(2,*)m_inputnum,m_hidenum,m_outnum
	read(2,*)((w_ih(i,j),j=1,m_inputnum),b(i,1),a(i,1),i=1,m_hidenum)
	read(2,*)((w_ho(i,j),j=1,m_hidenum+1),i=1,m_outnum)
	close(2)
	open(3,file='testdata.txt')
	open(4,file='testresult.txt')
	read(3,*)m_testnum
	do i=1,m_testnum
	   read(3,*)(x(i),i=1,m_inputnum)
       call calculate_output(x,x1,y)
       write(4,*)(y(i),i=1,m_outnum)
	enddo
	close(3)
	endif
	end
	


	subroutine calculate_output(x,x1,y)
	double precision x(30),x1(30),y(30)
	common/bp/m_inputnum,m_hidenum,m_outnum,m_samplenum
	common/ee/ee,es1,es
	common/nn/n2
	do 230	i=1,m_inputnum
	x(i)=x(i)*ee
230	continue
  ! write(*,*)ee,es1,es,n2
	call outputgenerate(x,x1,y)
  !	write(*,*)(y(i),i=1,m_outnum)
	if(n2.eq.0) then
		do 235 i=1,m_outnum
		y(i)=y(i)/es1
235		continue
	elseif(n2.gt.0) then
		do 240 i=1,m_outnum
		y(i)=y(i)/es
240	continue
	endif
	end
			 
	
	
	subroutine bp_ini
	double precision w_ih(30,30),w_ho(30,30),b(30,1),a(30,1)
	integer iseed
	common/bp/m_inputnum,m_hidenum,m_outnum,m_samplenum
	common/weight/w_ih,w_ho,b,a
5	write(*,*)'The element number of input layer:'
	read(*,*)m_inputnum
	write(*,*)'The element number of hide layer:'
	read(*,*)m_hidenum
	write(*,*)'The element number of output layer:'
	read(*,*)m_outnum
	if(m_inputnum.gt.30.or.m_hidenum.gt.30.or.m_outnum.gt.30) then
	write(*,*)'The net is too large!'
	goto 5
	endif
	iseed=563800
	do i=1,m_hidenum
	    b(i,1)=i*m_samplenum/m_hidenum
		a(i,1)=10
		do j=1,m_inputnum
			w_ih(i,j)=-ran(iseed)
			enddo
	enddo
	do i=1,m_outnum
		do j=1,m_hidenum+1
		    w_ho(i,j)=-ran(iseed)
		enddo
	enddo
	end

	
	subroutine readsampledata
	common/nn/n2
	common/ee/ee,es1,es
	common/bp/m_inputnum,m_hidenum,m_outnum,m_samplenum	 
	common/sample/sample_in,sample_t,sample_out,sample_x1
    double precision sample_in(200,30),sample_t(200,30),sample_out(200,30),sample_x1(200,30)
	write(*,*)'sample number is:'
	read(*,*)m_samplenum
	open(1,file="sampledata.txt")
	do  i=1,m_samplenum
	 	do j=1,m_inputnum
	      read(1,*)sample_in(i,j)
	    enddo
	enddo
	do i=1,m_samplenum
	    do j=1,m_outnum
            read(1,*)sample_out(i,j)
	    enddo
	enddo	
	n2=0
	do 55 i=1,m_samplenum
			do 54 j=1,m_outnum
		if(sample_out(i,j).lt.0) n2=n2+1
54		continue
55	continue
	do 49 i=1,m_samplenum 
	 		do 56 j=1,m_inputnum
		sample_in(i,j)=sample_in(i,j)*ee
56	continue
          if(n2.eq.0) then	
		do 57 j=1,m_outnum
	sample_out(i,j)=sample_out(i,j)*es1
57	continue
	elseif(n2.gt.0) then
		do 53 j=1,m_outnum
	sample_out(i,j)=(sample_out(i,j))*es
53	continue
	endif	
49	continue
	close(1)
	end
	
	subroutine outputgenerate(x,x1,y)
	double precision x(30),x1(30),y(30),s,morlet,sig
	double precision w_ih(30,30),w_ho(30,30),b(30,1),a(30,1)
	common/bp/m_inputnum,m_hidenum,m_outnum
	common/weight/w_ih,w_ho,b,a
	x1(m_hidenum+1)=-1.
	do 80 i=1,m_hidenum
		s=0.
		do 90 j=1,m_inputnum
		s=s+(w_ih(i,j)*x(j)-b(i,1))/a(i,1)
90	continue	
		x1(i)=morlet(s)
80	continue
	do 100	i=1,m_outnum
		s=0.
		do 110 j=1,m_hidenum+1
		s=s+w_ho(i,j)*x1(j)
110	continue
		y(i)=sig(s)
100	continue
	end		
	
	subroutine Learning
	double precision w_ih(30,30),w_ho(30,30),b(30,1),a(30,1),u(30,1)
	double precision sample_in(200,30),sample_t(200,30),sample_out(200,30),sample_x1(200,30)
	common/bp/m_inputnum,m_hidenum,m_outnum,m_samplenum
	common/weight/w_ih,w_ho,b,a
	common/sample/sample_in,sample_t,sample_out,sample_x1
	integer training_num
	double precision  x(30),x1(30),y(30)
	double precision w1(30,30),w2(30,30),w11(30,30),w22(30,30),b1(30,1),a1(30,1),b11(30,1),a11(30,1)
	double precision err,err_whole,err_whole1
	write(*,*)'Enter the error:'
	read(*,*)err_xy
	training_num=0
	do 105	i=1,m_outnum
			do 105	j=1,m_hidenum+1
105	w1(i,j)=w_ho(i,j)
	do 110	i=1,m_hidenum
	        b1(i,1)=b(i,1)
			a1(i,1)=a(i,1)
			do 110	j=1,m_inputnum
110	w2(i,j)=w_ih(i,j)
	err_whole1=0.
120	training_num=training_num+1
	err_whole=0.
	do 130 i=1,m_samplenum
			do 135 k=1,m_inputnum
			x(k)=sample_in(i,k)
135			continue
		call outputgenerate(x,x1,y)
			do 136 k=1,m_hidenum+1
			sample_x1(i,k)=x1(k)
136			continue
			do 137 k=1,m_outnum
			sample_t(i,k)=y(k)
137			continue
			err=0.
		do 140 j=1,m_outnum
			err=err+(sample_t(i,j)-sample_out(i,j))**2
140		continue
		err=err*0.5
!		write(*,*)err
	err_whole=err_whole+err
!	write(*,*)err_whole
130	continue
		uu=0.2
	write(*,*)'err_whole=',err_whole,'err_whole1=',err_whole1
	if(err_whole-err_whole1.le.0) then
		uu=1.20*uu
		alfa=0.9
	else
		uu=0.6*uu
		alfa=0.
	endif
	do 150	j=1,m_outnum
		do 160	k=1,m_hidenum+1
			s=0.
			do 170 i=1,m_samplenum
		s=s+(sample_t(i,j)-sample_out(i,j))*(1-sample_t(i,j))*sample_t(i,j)*sample_x1(i,k)
170	continue
	w11(j,k)=w_ho(j,k)-uu*s+alfa*(w_ho(j,k)-w1(j,k))
	w1(j,k)=w_ho(j,k)
	w_ho(j,k)=w11(j,k)
160	continue
150	continue
	do 180	j=1,m_hidenum
		do 190	k=1,m_inputnum
		     u(j,k)=(w_ih(j,k)*x(k)-b(j,1))/a(j,1)
	         ss=0.
			 sss=0
			 ssss=0
			do 200	i=1,m_samplenum
					s=0.
				do 210 l=1,m_outnum
			s=s+(sample_t(i,l)-sample_out(i,l))*(1-sample_t(i,j))*sample_t(i,j)*w_ho(l,j)
210	continue
		ss=ss+s*(-1.75*tan(1.75*u(j,k))-u(j,k))*sample_in(i,j)*sample_x1(i,k)/a(j,1)!*(1-sample_x1(i,k))	
		sss=sss+s*(-1.75*tan(1.75*u(j,k))-u(j,k))*sample_x1(i,k)*(-1.)/a(j,1)!*(1-sample_x1(i,k))
		ssss=ssss+s*(-1.75*tan(1.75*u(j,k))-u(j,k))*sample_x1(i,k)*u(j,k)/a(j,1)*(-1.)!*(1-sample_x1(i,k))	
	!	ss=ss+s*w_ho(j,k)*(cos(1.75*u(j,k)))*exp(u(j,k)**2/(-2))*sample_in(i,j)/a(j,1)*sample_t(i,j)*(1-sample_t(i,j))	
	!	sss=sss+s*w_ho(j,k)*(cos(1.75*u(j,k)))*exp(u(j,k)**2/(-2))/a(j,1)*sample_t(i,j)*(1-sample_t(i,j))
	!	ssss=ssss+s*w_ho(j,k)*(cos(1.75*u(j,k)))*exp(u(j,k)**2/(-2))*sample_in(i,j)*w_ih(j,k)/a(j,1)*sample_t(i,j)*(1-sample_t(i,j))				  
200	continue
	w22(j,k)=w_ih(j,k)-uu*ss+alfa*(w_ih(j,k)-w2(j,k))
	w2(j,k)=w_ih(j,k)
	w_ih(j,k)=w22(j,k)
	b11(j,1)=b(j,1)-uu*sss+alfa*(b(j,1)-b1(j,1))
	b1(j,1)=b(j,1)
	b(j,1)=b11(j,1)
	a11(j,1)=a(j,1)-uu*ssss+alfa*(a(j,1)-a1(j,1))
	a1(j,1)=a(j,1)
	a(j,1)=a11(j,1)
190	continue
180	continue
	err_whole1=err_whole
	if(training_num.gt.500000) goto 250
	if(err_whole.gt.err_xy) goto 120
250	write(*,*)'Training is over!'
!	write(*,*)'err_whole=',err_whole
	write(*,*)'The training number is',training_num
	end
	 
	
	function morlet(x)
	double precision morlet
	double precision x
	morlet=cos(1.75*x)*exp(-x**2/2.)
	end


	function sig(x)
	double precision sig
	double precision x
	sig=1./(1.+exp(-x))
	end