program Bp_Neutral_Net 
	character*1 cc
	character*15 wfilename
	integer m_samplenum
	common/bp/m_inputnum,m_hidenum,m_outnum,m_samplenum
	common/nn/n2
	common/ee/ee,es1,es
cc**********************************************************
	double precision x(30),x1(30),y(30),dy
	double precision w_ih(30,30),w_ho(30,30)
	double precision sample_in(200,30),sample_t(200,30),
	&sample_out(200,30),sample_x1(200,30)
	common/weight/w_ih,w_ho
	common/sample/sample_in,sample_t,sample_out,sample_x1
	ee=1.0/4000.0
	es1=0.01
	es=0.01
	write(*,*)'Do you want to learn(l) or calculate(c) or derive(d)?'
	read(*,*)cc
	if(cc.eq.'l'.or.cc.eq.'L') then
	call bp_ini	
	call readsampledata
	call Learning
	write(*,*)'Enter the name of bp file!'
	read(*,*)wfilename
	open(2,file=wfilename)
	write(2,*)n2
	write(2,*)m_inputnum,m_hidenum,m_outnum
	write(2,*)((w_ih(i,j),j=1,m_inputnum+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='bp.dat')
	read(2,*)n2
	read(2,*)m_inputnum,m_hidenum,m_outnum
	read(2,*)((w_ih(i,j),j=1,m_inputnum+1),i=1,m_hidenum)
	read(2,*)((w_ho(i,j),j=1,m_hidenum+1),i=1,m_outnum)
	close(2)
cc	write(*,*)((w_ih(i,j),j=1,m_inputnum+1),i=1,m_hidenum)
	write(*,*)'enter x:'
	read(*,*)(x(i),i=1,m_inputnum)
	call calculate_output(x,x1,y)
	write(*,*)(y(i),i=1,m_outnum)
	elseif(cc.eq.'d'.or.cc.eq.'D') then
	open(2,file='bp.dat')
	read(2,*)m_inputnum,m_hidenum,m_outnum
	read(2,*)((w_ih(i,j),j=1,m_inputnum+1),i=1,m_hidenum)
	read(2,*)((w_ho(i,j),j=1,m_hidenum+1),i=1,m_outnum)
	close(2)
	read(*,*)(x(i),i=1,m_inputnum)
	read(*,*)i,j
	call derivevalue(x,i,j,dy)
	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-2.0
240	continue
	endif
	end
			 
	subroutine bp_ini
	double precision w_ih(30,30),w_ho(30,30)
	 double precision rn
	common/bp/m_inputnum,m_hidenum,m_outnum,m_samplenum
	common/weight/w_ih,w_ho
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
	r1=10
	r2=100
	do 20 i=1,m_hidenum
		do 20 j=1,m_inputnum+1
20			w_ih(i,j)=rn(r1)
	do 40 i=1,m_outnum
			do 40 j=1,m_hidenum+1
40			w_ho(i,j)=rn(r2)
	end
	
	function rn(R)
      double precision rn
	s=65536
	u=2053
	v=13849
	r=u*r+v
	m=r/s
	r=r-m*s
	rn=r/s
	return 
	end

	subroutine readsampledata
	character*15 filename
	common/nn/n2
	common/ee/ee,es1,es
	common/bp/m_inputnum,m_hidenum,m_outnum,m_samplenum
	double precision sample_in(200,30),sample_t(200,30),
     &sample_out(200,30),sample_x1(200,30)
	common/sample/sample_in,sample_t,sample_out,sample_x1
cc     	write(*,*)m_inputnum,m_hidenum,m_outnum
	write(*,*)'Enter the name of sample file:'
	read(*,*)filename
	open(1,file=filename)
	read(1,*)m_samplenum
CC	WRITE(*,*)M_SAMPLENUM
	do 50 i=1,m_samplenum
	read(1,*)(sample_in(i,j),j=1,m_inputnum)
	read(1,*)(sample_out(i,j),j=1,m_outnum)
cc	WRITE(*,*)(sample_in(i,j),j=1,m_inputnum)
cc	WRITE(*,*)(sample_out(i,j),j=1,m_outnum)
50	continue
	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)=(2+sample_out(i,j))*es
53	continue
	endif	
49	continue
	close(1)
cc	WRITE(*,*)(sample_in(i,j),j=1,m_inputnum)
cc	WRITE(*,*)(sample_out(i,j),j=1,m_outnum)
	end
	
	subroutine outputgenerate(x,x1,y)
	double precision x(30),x1(30),y(30),s,sig
	double precision w_ih(30,30),w_ho(30,30)
	common/bp/m_inputnum,m_hidenum,m_outnum
	common/weight/w_ih,w_ho
	x1(m_hidenum+1)=-1.
	x(m_inputnum+1)=-1.
	do 80 i=1,m_hidenum
		s=0.
		do 90 j=1,m_inputnum+1
		s=s+w_ih(i,j)*x(j)
90	continue	
		x1(i)=sig(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)
	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
	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)
	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
			do 110	j=1,m_inputnum+1
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
	err_whole=err_whole+err
130	continue
		uu=0.5
	write(*,*)'err_whole=',err_whole,'err_whole1=',err_whole1
	if(err_whole-err_whole1.lt.0) then
		uu=1.20*uu
		alfa=0.9
	else
		uu=0.60*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+1
	         ss=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,l))*
     &sample_t(i,l)*w_ho(l,j)
210	continue
		ss=ss+s*(1-sample_x1(i,j))*sample_x1(i,j)*sample_in(i,k)
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)
190	continue
180	continue
cc	write(*,*)((w_ho(i,j),j=1,hide_num+1),i=1,out_num)
cc	write(*,*)'err_whole=',err_whole
	err_whole1=err_whole
	if(err_whole.gt.err_xy) goto 120
	write(*,*)'Training is over!'
	write(*,*)'err_whole=',err_whole
	write(*,*)'The training number is',training_num
	end
	
	subroutine derivevalue(x,iy,jx,dy)
	double precision x(30),dy
	integer iy,jx
	double precision y(30),x1(30)
	double precision w_ih(30,30),w_ho(30,30)
	common/bp/m_inputnum,m_hidenum,m_outnum,m_samplenum
	common/weight/w_ih,w_ho
	call outputgenerate(x,x1,y)
	s=0.
	do 220	k=1,m_hidenum+1
	s=s+(1-y(iy))*y(iy)*w_ho(iy,k)*(1-x1(k))*x1(k)*w_ih(k,jx)
220	continue      
	dy=s
	end			 
	
	 function sig_d(x)
	double precision sig_d,sig
	double precision x
	sig_d=sig(x)*(1.-sig(x))
	end
	 
	 function sig(x)
	double precision sig
	double precision x
	sig=1./(1.+exp(-x))
	end