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

      real*8 aux,oldknot(2,17),diff
      integer robo,i,j,l,nn,loc,m,n
      character rep

      real*8 x(32),fvec(32)
      real*8 xh(32),xg(32),yknot(2,18),alph(2,17),bet(2,17),ti
      logical check
      real*8 xsc(32),fsc(32),rp(5),fvec(32)
      integer ip(6)

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

      real*8 moy1,moy2,std1,std2,era,fn,xres(32)
      external dneqnf,funcv,erset,drnun,rnset,dn4qbj,dneqbf

      intrinsic dlog,dexp
      open(unit=60,file='data_asym_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,0,0)
      call rnset(778354)
      nn=16
      mc=10000
      np=4
      p(1)=3.0d0
      p(2)=3.0d0
      p(3)=5.0d0
      p(4)=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

        x(           1 )=  -2.84827109173688     
        x(           2 )=  -2.85209052797506     
        x(           3 )=  -2.85688597597939     
        x(           4 )=  -2.86243170560486     
        x(           5 )=  -2.87005433886236     
        x(           6 )=  -2.87972857171734     
        x(           7 )=  -2.89195095174828     
        x(           8 )=  -2.90986034638539     
        x(           9 )=  -2.93411212816428     
        x(          10 )=  -2.96990679761558     
        x(          11 )=  -3.03212449908756     
        x(          12 )=  -3.14051739544213     
        x(          13 )=  -3.33100693581025     
        x(          14 )=  -3.69184084490141     
        x(          15 )=  -5.61683383184477     
        x(          16 )=  -7.41860036191254     
        x(          17 )=   4.37486841314759     
        x(          18 )=   3.55315379806175     
        x(          19 )=   3.09860286531008     
        x(          20 )=   2.74440338704673     
        x(          21 )=   2.45415125926837     
        x(          22 )=   2.19413505997764     
        x(          23 )=   1.94448901027584     
        x(          24 )=   1.71125414593309     
        x(          25 )=   1.47745377147087     
        x(          26 )=   1.23478257742541     
        x(          27 )=   1.01951187118377     
        x(          28 )=  0.846654741348483     
        x(          29 )=  0.726435662568677     
        x(          30 )=  0.619443991395723     
        x(          31 )=   1.12697386952317     
        knot(1,           1 )=  0.000000000000000E+000
        knot(1,           2 )=  6.250000000000000E-002
        knot(1,           3 )=  0.125000000000000     
        knot(1,           4 )=  0.187500000000000     
        knot(1,           5 )=  0.250000000000000     
        knot(1,           6 )=  0.312500000000000     
        knot(1,           7 )=  0.375000000000000     
        knot(1,           8 )=  0.437500000000000     
        knot(1,           9 )=  0.500000000000000     
        knot(1,          10 )=  0.562500000000000     
        knot(1,          11 )=  0.625000000000000     
        knot(1,          12 )=  0.687500000000000     
        knot(1,          13 )=  0.750000000000000     
        knot(1,          14 )=  0.812500000000000     
        knot(1,          15 )=  0.875000000000000     
        knot(1,          16 )=  0.937500000000000     
        knot(1,          17 )=   1.00000000000000  

       n=16
       robo=15
       do 784 l=1,n+1
          knot(2,l)=knot(1,l)
  784  continue

       xh=x

c enter a loop that calculates the CSE for different k

  700 if (n.le.nn) then

c        do 688 l=1,2*n-1
c           call drnun(1,aux)
c           x(l)=xh(l)*(0.995d0+aux/100.0d0)
c  688   continue

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

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

c         cvrg=.true.
c         call funcv(x,fvec,2*n,1)

      
      write(6,*)'fn',fn
      if (fn.gt.1.0e-20) then
         goto 700
      else
         read *,rep
      endif

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

       xres=x

       yknot(1,1)=0.0d0
       yknot(2,1)=0.0d0
       do 458 l=1,n-1
          yknot(1,l+1)=yknot(1,l)+dexp(x(l))
          alph(1,l)=yknot(1,l)
          bet(1,l)=(yknot(1,l+1)-yknot(1,l))/(knot(1,l+1)-knot(1,l))
          aux=dexp(x(n+l))
          aux=aux/(1.0d0+aux)
          yknot(2,l+1)=yknot(2,l)+(yknot(1,l+1)-yknot(2,l))*aux
          alph(2,l)=yknot(2,l)
          bet(2,l)=(yknot(2,l+1)-yknot(2,l))/(knot(2,l+1)-knot(2,l))
  458  continue
       yknot(1,n+1)=yknot(1,n)+dexp(x(n))
       alph(1,n)=yknot(1,n)
       bet(1,n)=(yknot(1,n+1)-yknot(1,n))/(knot(1,n+1)-knot(1,n))
       yknot(2,n+1)=yknot(1,n+1)
       alph(2,n)=yknot(2,n)
       bet(2,n)=(yknot(2,n+1)-yknot(2,n))/(knot(2,n+1)-knot(2,n))


          n=n+1
          oldknot=knot
          knot(1,n+1-robo)=(knot(1,n-robo)+knot(1,n+1-robo))/2.0d0
          knot(2,n+1-robo)=knot(1,n+1-robo)
          do 202 l=n+2-robo,n+1
             knot(1,l)=oldknot(1,l-1)
             knot(2,l)=knot(1,l)
  202     continue
         robo=robo+2

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

         do 214 l=1,n
            aux=(yknot(2,l+1)-yknot(2,l))/(yknot(1,l+1)-yknot(2,l))
            x(n+l)=dlog(aux/(1.0d0-aux))
  214    continue
         xh=x

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


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

      subroutine funcv(x,fvec,n2)

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

      real*8 alph(2,20),bet(2,20),da(40),aux,bi,ti,erre
      real*8 dbdt1,dbt2,dbdp,cumu,dcumu,cumu1,dcumu1
      real*8 bbi,invbi,cumu2,dcumu2
      real*8 yknot(2,18)
      integer m,i,l,ll,k
      character rep
      logical check
      real*8 dbetdf,dbeta,const1,const2
      external dbetdf,dbeta

      real*8 knot(2,17),p(4),u(100000,4),slp
      integer mc,np
      logical cvrg
      common knot,p,u,slp,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

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

       const1=dbeta(p(1),p(2))
       const2=dbeta(p(3),p(4))

      do 898 m=1,mc

c*********************************************************
c             FOC for Bidder 1
c*********************************************************

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

                cumu1=dbetdf(ti,p(1),p(2))
                dcumu1=(ti**(p(1)-1))*((1-ti)**(p(2)-1))/const1
            endif
            goto 600
         endif

c      write(6,*)ti,bi,dbdt1,dbdp,cumu1,dcumu1
c      read *,rep

         check=.true.
         ll=0
         invbi=1.0d0
  697    if (check) then
            ll=ll+1
            if (bi.ge.yknot(2,ll).and.bi.le.yknot(2,ll+1).
     &          and.ll.le.n) then
                invbi=(bi-alph(2,ll))/bet(2,ll)+knot(2,ll)
                dbdt2=bet(2,ll)
                check=.false.

                cumu2=dbetdf(invbi,p(3),p(4))
                dcumu2=(invbi**(p(3)-1))*((1-invbi)**(p(4)-1))/const2
            endif
            if (ll.lt.n) then
               goto 697
            endif
         endif

c      write(6,*)bi,invbi,dbdt2,cumu2,dcumu2
c      read *,rep

         cumu=cumu1*cumu2**2
         dcumu=dcumu1*(cumu2**2)/dbdt1+2*dcumu2*cumu1*cumu2/dbdt2
         da(l)=da(l)+dbdp*((ti-bi)*dcumu-cumu)

c       write(6,*)cumu,dcumu,l,da(l)
c       read *,rep

c*********************************************************
c             FOC for Bidder 4
c*********************************************************

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

                cumu1=dbetdf(ti,p(3),p(4))
                dcumu1=(ti**(p(3)-1))*((1-ti)**(p(4)-1))/const2
            endif
            goto 699
         endif

         check=.true.
         ll=0
         invbi=1.0d0
  698    if (check) then
            ll=ll+1
            if (bi.ge.yknot(1,ll).and.bi.le.yknot(1,ll+1).
     &          and.ll.le.n) then
                invbi=(bi-alph(1,ll))/bet(1,ll)+knot(1,ll)
                dbdt2=bet(1,ll)
                check=.false.

                cumu2=dbetdf(invbi,p(1),p(2))
                dcumu2=(invbi**(p(1)-1))*((1-invbi)**(p(2)-1))/const1
            endif
            if (ll.lt.n) then
               goto 698
            endif
         endif

         cumu=cumu1*cumu2**2
         dcumu=dcumu1*(cumu2**2)/dbdt1+2*dcumu2*cumu1*cumu2/dbdt2
         da(n+l)=da(n+l)+dbdp*((ti-bi)*dcumu-cumu)

  898 continue

       do 763 l=1,n-1
          fvec(l)=da(l)*yknot(1,l+1)/mc
c          fvec(l)=da(l)/mc
          aux=dexp(x(n+l))
          aux=aux/(1.0d0+aux)**2
          fvec(n+l)=da(n+l)*aux/mc
c          fvec(n+l)=da(n+l)/mc
  763  continue
          fvec(n)=da(n)*yknot(1,n+1)/mc

c*********************************************************
      


      if (cvrg) then

         do 684 m=1,2*n-1
            write(6,*)'       x(',m,')=',x(m)
  684    continue

         do 685 m=1,n+1
            write(6,*)'       knot(1,',m,')=',knot(1,m)
  685    continue

         write(6,*)''
         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(1,l).and.ti.le.knot(1,l+1)) then
                   bi=alph(1,l)+bet(1,l)*(ti-knot(1,l))
                   check=.false.
               endif
               goto 700
            endif
            check=.true.
            l=0
  730       if (check) then
               l=l+1
               if (ti.ge.knot(2,l).and.ti.le.knot(2,l+1)) then
                   bbi=alph(2,l)+bet(2,l)*(ti-knot(2,l))
                   check=.false.
               endif
               goto 730
            endif
            write(6,151)ti,bi,bbi
  610    continue
      write(6,*)' '

      write(6,*)'PLAYERS 1 and 2'
      write(6,*)
     &'     l           knot(1,l-1)     knot(1,l)    CSE[knot(1,l)]'
      do 608 l=1,n-1
         write(6,151)l*1.0d0,knot(1,l),knot(1,l+1),alph(1,l+1)
  608 continue
      write(6,151)l*1.0d0,knot(1,n),knot(1,n+1),alph(1,n)
     &+bet(1,n)*(1-knot(1,n))
      write(6,*)' '
      write(6,*)'       l             alph(1,l)         bet(1,l)'
      do 609 l=1,n
         write(6,151)l*1.0d0,alph(1,l),bet(1,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

      write(6,*)''
      write(6,*)''
      write(6,*)'PLAYERS 3 and 4'
      write(6,*)
     &'     l           knot(2,l-1)     knot(2,l)    CSE[knot(2,l)]'
      do 908 l=1,n-1
         write(6,151)l*1.0d0,knot(2,l),knot(2,l+1),alph(2,l+1)
  908 continue
      write(6,151)l*1.0d0,knot(2,n),knot(2,n+1),alph(2,n)
     &+bet(2,n)*(1-knot(2,n))
      write(6,*)' '
      write(6,*)'       l             alph(2,l)         bet(2,l)'
      do 909 l=1,n
         write(6,151)l*1.0d0,alph(2,l),bet(2,l)
  909 continue
      write(6,*)' '
      erre=0.0d0
      write(6,*)'       l              x(l)            f(l)'
      do 907 l=1,n-1
         write(6,151)l*1.0d0,x(n+l),fvec(n+l)
         erre=erre+(fvec(n+l))**2
  907 continue
      write(6,*)' '
      write(6,*)'Norm of derivatives:',erre**0.5

      endif

c      if (n.gt.15) then
c         read *,rep
c      endif

  151 format (4d16.8)
      return
      end


