c*******************************************************
c*******************************************************
C      piecewise linear CSE with adaptive knots
c*******************************************************
c*******************************************************

      real*8 xg(16),aux,xres(16),oldknot(17),diff
      integer i,j,l,nn,loc,n
      character rep

      real*8 x(16),fvec(16),diag(16),fjac(16,16),r(136),
     *                 qtf(16),wa1(16),wa2(16),wa3(16),wa4(16)
      real*8 yknot(18),alph(17),bet(17),ti
      real*8 xsc(16),fsc(16),rp(5),fvec(16)
      integer ip(6)
      logical check


      real*8 knot(17),p(2),u(10000,4)
      integeR mc,np
      logical cvrg
      common knot,p,u,mc,np,cvrg

      real*8 moy1,moy2,std1,std2,era,fn
      external dneqnf,funcv,erset,dn4qbj,dneqbf
      intrinsic dlog,dexp
      open(unit=60,file='data_sym_10_4.dat',status='old')

c n=# of parameter, mc= size of Monte carlo
c np= number of players, p(.)=parameters of the distribution

      call erset(0,1,0)
      nn=16
      mc=10000
      np=4
      p(1)=3.0d0
      p(2)=3.0d0

c read mc private signals for each of the 4 players

      do 310 m=1,mc
      do 310 i=1,np
         read(60,100)u(m,i)
  310 continue

      moy1=0.0d0
      std1=0.0d0
      moy2=0.0d0
      std2=0.0d0

      do 313 i=1,np
      do 313 m=1,mc
             if (i.lt.3) then
                 moy1=moy1+u(m,i)
             else
                 moy2=moy2+u(m,i)
             endif
  313 continue

      moy1=moy1/(m*2.0d0)
      moy2=moy2/(m*2.0d0)
      do 315 m=1,mc
      do 315 i=1,np
         if (i.lt.3) then
             std1=std1+(u(m,i)-moy1)**2
         else
             std2=std2+(u(m,i)-moy2)**2
         endif
  315 continue
      std1=std1/(m*2.0d0)
      std2=std2/(m*2.0d0)

      write(6,*)"moy et std de player 1",moy1,std1**0.5
      write(6,*)"moy et std de player 2",moy2,std2**0.5
      write(6,*)' '
      read *,rep



c parameters initialization

      oldknot(1)=0.0d0
      knot(1)=0.0d0
      knot(2)=2.0d0/3.0d0
      knot(3)=1.0d0
      x(1)=dlog(0.25d0)
      x(2)=dlog(0.25d0)
      n=2

c enter a loop that calculates the CSE for different k

  700 if (n.le.nn) then

c         cvrg=.true.
         cvrg=.false.

         xsc=1.0d0
         fsc=1.0d0
         call dn4qbj(ip,rp)
         ip(3)=500
         ip(4)=1000
c         ip(6)=1

         max=200
         era=1.e-12
         xg=x
         call dneqnf(funcv,era,n,max,xg,x,fn)
c         call dneqbf(funcv,n,xg,xsc,fsc,ip,rp,x,fvec)

         cvrg=.true.
         call funcv(x,fvec,n)
         write(6,*)'The agorithm has converged. n=',n
         read *,rep

c update the parameters before moving to a CSE with a higher k

       yknot(1)=0.0d0
       do 456 l=1,n
          yknot(l+1)=yknot(l)+dexp(x(l))
          alph(l)=yknot(l)
          bet(l)=(yknot(l+1)-yknot(l))/(knot(l+1)-knot(l))
  456  continue

          diff=0.0d0
          xres=x
          do 201 l=2,n
             oldknot(l)=knot(l)
             aux=abs(bet(l)-bet(l-1))
             if (aux.gt.diff) then
                diff=aux
                loc=l
             endif
  201     continue
          oldknot(n+1)=1.0d0

          n=n+2
          knot(loc+1)=(oldknot(loc-1)+oldknot(loc)+oldknot(loc+1))/3.0d0
          knot(loc)=(oldknot(loc-1)+2.0d0*knot(loc+1))/3.0d0
          knot(n+1)=1.0d0

          do 202 l=1,loc-1
             knot(loc-l+1)=(oldknot(loc-l)+2.0d0*knot(loc+2-l))/3.0d0
  202     continue

          do 203 l=loc+2,n
             knot(l)=(oldknot(l-1)+2.0d0*knot(l-1))/3.0d0
  203     continue


        yknot(1)=0.0d0
        do 777 ll=1,n
           check=.true.
           l=0
           ti=knot(ll+1)
  608      if (check) then
              l=l+1
              if (ti.ge.oldknot(l).and.ti.le.oldknot(l+1)) then
                 yknot(ll+1)=alph(l)+bet(l)*(ti-oldknot(l))
                 check=.false.
              endif
              goto 608
           endif
           x(ll)=dlog(yknot(ll+1)-yknot(ll)) 
c           write(6,*)ll,knot(ll),yknot(ll),x(ll)
  777   continue
c        write(6,*)yknot(n+1),yknot(n)

         goto 700
      endif
  100 format (2d20.12)
      end


c********************************************:wq
*********************
c*****************************************************************
c      This subroutine calculates the objective function
c*****************************************************************
c*****************************************************************

      subroutine funcv(x,fvec,n)
c      subroutine funcv(n,x,fvec)

      integer n,iflag
      real*8 x(n),fvec(n)

      real*8 alph(n),bet(n),da(n),aux,bi,ti,erre
      real*8 dbdt,dbdp,cumu,dcumu,cumu1,dcumu1
      real*8 yknot(17)
      integer m,i,l
      character rep
      logical check
      real*8 dbetdf,dbeta,const
      external dbetdf,dbeta

      real*8 knot(17),p(2),u(10000,4)
      integer mc,np
      logical cvrg
      common knot,p,u,mc,np,cvrg

c  set the value of the constrained strategy parameters
c  such that the strategy is continuous
c alph= constant, bet=slope, da(l)= derivative wrt l parameter


       const=dbeta(p(1),p(2))

       da=0.0d0
       yknot(1)=0.0d0
c       yknot(2)=dexp(x(1))
c       alph(1)=0.0d0
c       bet(1)=yknot(2)/knot(2)
       do 456 l=1,n
c          bet(l)=dexp(x(l))
c          alph(l)=alph(l-1)+bet(l-1)*(knot(l)-knot(l-1))
          yknot(l+1)=yknot(l)+dexp(x(l)) 
          alph(l)=yknot(l)
          bet(l)=(yknot(l+1)-yknot(l))/(knot(l+1)-knot(l))
  456  continue

      do 898 m=1,mc

         ti=u(m,1)
         check=.true.
         l=0
  600    if (check) then
            l=l+1
            if (ti.ge.knot(l).and.ti.le.knot(l+1)) then
                bi=alph(l)+bet(l)*(ti-knot(l))
                dbdt=bet(l)
c                dbdp=ti-knot(l)
                dbdp=(ti-knot(l))/(knot(l+1)-knot(l))
                check=.false.

                cumu1=dbetdf(ti,p(1),p(2))
                dcumu1=(ti**(p(1)-1))*((1-ti)**(p(2)-1))/const
                cumu=cumu1**(np-1.0d0)
                dcumu=(np-1.0d0)*dcumu1*cumu1**(np-2.0d0)

                da(l)=da(l)+dbdp*((ti-bi)*dcumu/dbdt-cumu)
            endif
            goto 600
         endif

  898 continue

       do 763 l=1,n
c  763     fvec(l)=da(l)*bet(l)/mc
  763     fvec(l)=da(l)*yknot(l+1)/mc

      if (cvrg) then
         write(6,*)'       t              CSE(t)'
         do 610 m=1,101
            ti=1.0d0*(m*1.0d0-1.0d0)/100.0d0
            check=.true.
            l=0
  700       if (check) then
               l=l+1
               if (ti.ge.knot(l).and.ti.le.knot(l+1)) then
                   bi=alph(l)+bet(l)*(ti-knot(l))
                   check=.false.
               endif
               goto 700
            endif
            write(6,151)ti,bi
  610    continue
      write(6,*)' '

      write(6,*)
     &'       l             knot(l-1)       knot(l)      CSE[knot(l)]'
      do 608 l=1,n-1
         write(6,151)l*1.0d0,knot(l),knot(l+1),alph(l+1)
  608 continue
      write(6,151)l*1.0d0,knot(n),knot(n+1),alph(n)+bet(n)*(1-knot(n))
      write(6,*)' '
      write(6,*)'       l             alph(l)         bet(l)'
      do 609 l=1,n
         write(6,151)l*1.0d0,alph(l),bet(l)
  609 continue
      write(6,*)' '
      erre=0.0d0
      write(6,*)'       l              x(l)            f(l)'
      do 607 l=1,n
         write(6,151)l*1.0d0,x(l),fvec(l)
         erre=erre+(fvec(l))**2
  607 continue
      write(6,*)' '
      write(6,*)'Norm of derivatives:',erre**0.5

      endif
c      read *,rep

  151 format (4d16.8)
      return
      end


