***  Renyi conditional mutual information I(x;y_tau|y) for 2 scalar variables
*** yy(.,.) input time series, yy(1,.) var 1, yy(2,.) var 2 in common block cirtq
*** input parameters:
**** nse series length, here max 131072
**** kmin,kmax,kstep range of forward lags tau, e.g 1,1,1 for just lag 1 sample
***** 1,5,1 for lags 1-5, 1,5,2 for lags 1,3,5. max 256 lags
**** nq number of equidistant bins in one dimension, here max 16
***** this is 3-dim estimator using nq**3 bins, be sure nq**3 << nse
*** qe(.) Renyi parameters q, or alpha, values filled below
*** error logical, true value indicates problems with binning, result invalid
*** r3(.,.,.) resulting RCMI in common block cirtq:
**** r3(direction,Renyi param index,forward lag index)
**** r3(1,.,.) causal direction 2->1, r3(2,.,.) 1->2



      subroutine r0cirt2u1lg(nse,kmin,kmax,kstep,nq,
     & qe,error)

      common /cirtq/ yy(2,131072),r3(2,64,256)

      real y0(2,131072)
      dimension ra(131072),rb(131072),ctp(2,16),
     &vvvr3(2,256),r2(2,256),r4(2,256)

      integer mmm(16,16,16),m12(16,16),m13(16,16),m23(16,16),
     & m1(16),m2(16),m3(16),m4(16),m14(16,16)
       logical error
      real qe(64) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      error=.false.
      

      qe(2)=0.1
      qe(3)=0.2
      qe(4)=0.3
      qe(5)=0.4
      qe(6)=0.5
      qe(7)=0.6
      qe(8)=0.7
      qe(9)=0.8
      qe(10)=0.9
      qe(1)=0.05
      qe(11)=0.95
      qe(12)=1.05
      qe(13)=1.1
      qe(14)=1.2
      qe(15)=1.3
      qe(16)=1.4
      qe(17)=1.5
      qe(18)=1.6
      qe(19)=1.7
      qe(20)=1.8
      qe(21)=1.9
      qe(22)=2.0
      qe(23)=2.1
      qe(24)=2.2
      qe(25)=2.3
      qe(26)=2.4
      qe(27)=2.5
      qe(28)=2.6
      qe(29)=2.7
      qe(30)=2.8
      qe(31)=2.9
      qe(32)=3.0

cc      

      do ivar=1,2

      xx=yy(ivar,1)
      xm=xx      
      do i=1,nse
        x1=yy(ivar,i)
        if (x1.gt.xx) xx=x1
        if (x1.lt.xm) xm=x1
        enddo
       xd=(xx-xm)/float(nq)
       do i=1,nse
        x1=yy(ivar,i)-xm
        x1=x1/xd
        ix=ifix(x1) +1
        if (ix.gt.nq) ix=nq
        y0(ivar,i)=ix
        enddo
        enddo  !ivar

**********************************
       ncut=0
       lagov=(kmax-kmin+kstep)/kstep
***** 
       nsea=nse-(kmax+iabs(kmin))-2*ncut
       pfa=float(nsea)


        i0=1
        if (kmin.lt.0) i0=1-kmin
        i9=0
        if (kmax.gt.0) i9=kmax
**********
********** lag cyc
**********
       ik=0
       do kk=kmin,kmax,kstep
        ik=ik+1

        do i1=1,nq
         do i2=1,nq
           do i3=1,nq
            mmm(i1,i2,i3)=0
            enddo
           m12(i1,i2)=0
           m23(i1,i2)=0
           m13(i1,i2)=0
           m14(i1,i2)=0
           enddo
          m1(i1)=0
          m2(i1)=0
          m3(i1)=0
          m4(i1)=0
          enddo


        do i=ncut+i0,nse-i9-ncut
         i1=y0(1,i)
         i2=y0(1,i+kk)
         i3=y0(2,i)
         i4=y0(2,i+kk)
         m1(i1)=m1(i1)+1
         m2(i2)=m2(i2)+1
         m3(i3)=m3(i3)+1
         m4(i4)=m4(i4)+1

         m12(i1,i2)=m12(i1,i2)+1
         m13(i1,i3)=m13(i1,i3)+1
         m23(i2,i3)=m23(i2,i3)+1
         m14(i1,i4)=m14(i1,i4)+1

         mmm(i1,i2,i3)=mmm(i1,i2,i3)+1
         enddo


          do iq=1,32
           q=qe(iq)
*0000000000000
        xxs=0.
        x1s=0.
        x4s=0.
        x12s=0.
        x13s=0.
        do i1=1,nq
           x1=float(m1(i1))/pfa
           x1s=x1s+x1**q
         do i2=1,nq

           x12=float(m12(i1,i2))/pfa
           x12s=x12s+x12**q
           x13=float(m13(i1,i2))/pfa
           x13s=x13s+x13**q
          do i3=1,nq

           xxx=float(mmm(i1,i2,i3))/pfa

           if (xxx.gt.0.) then

             xxs=xxs+xxx**q

             endif

            enddo


           enddo
          enddo
          
          xxs=(x12s*x13s)/(xxs*x1s)
                   if (xxs.gt.0.) then
          xxs=alog(xxs)/(1.-q)
          else
          xxs=0.
          endif



         r3(1,iq,ik)=xxs
         enddo !iq

         enddo !lag cyc


***********************
*** opposite direction
*******
      do i=1,nse
       xxx=y0(1,i)
       y0(1,i)=y0(2,i)
       y0(2,i)=xxx
       enddo
***********************

        i0=1
        if (kmin.lt.0) i0=1-kmin
        i9=0
        if (kmax.gt.0) i9=kmax
**********
********** lag cyc
**********
       ik=0
       do kk=kmin,kmax,kstep
        ik=ik+1

        do i1=1,nq
         do i2=1,nq
           do i3=1,nq
            mmm(i1,i2,i3)=0
            enddo
           m12(i1,i2)=0
           m23(i1,i2)=0
           m13(i1,i2)=0
           m14(i1,i2)=0
           enddo
          m1(i1)=0
          m2(i1)=0
          m3(i1)=0
          m4(i1)=0
          enddo


        do i=ncut+i0,nse-i9-ncut
         i1=y0(1,i)
         i2=y0(1,i+kk)
         i3=y0(2,i)
         i4=y0(2,i+kk)
         m1(i1)=m1(i1)+1
         m2(i2)=m2(i2)+1
         m3(i3)=m3(i3)+1
         m4(i4)=m4(i4)+1

         m12(i1,i2)=m12(i1,i2)+1
         m13(i1,i3)=m13(i1,i3)+1
         m23(i2,i3)=m23(i2,i3)+1
         m14(i1,i4)=m14(i1,i4)+1

         mmm(i1,i2,i3)=mmm(i1,i2,i3)+1
         enddo


         do iq=1,32
          q=qe(iq)
          

*0000000000000
        xxs=0.
        x1s=0.
        x4s=0.
        x12s=0.
        x13s=0.
        do i1=1,nq
           x1=float(m1(i1))/pfa
           x1s=x1s+x1**q
         do i2=1,nq

           x12=float(m12(i1,i2))/pfa
           x12s=x12s+x12**q
           x13=float(m13(i1,i2))/pfa
           x13s=x13s+x13**q
          do i3=1,nq

           xxx=float(mmm(i1,i2,i3))/pfa

           if (xxx.gt.0.) then

             xxs=xxs+xxx**q

             endif

            enddo

           enddo
          enddo
          
          xxs=(x12s*x13s)/(xxs*x1s)
                   if (xxs.gt.0.) then
          xxs=alog(xxs)/(1.-q)
          else
          xxs=0.
          endif

         r3(2,iq,ik)=xxs
         enddo !iq

         enddo !lag cyc


      return
      end