!This code simulates an SIR model with a quarantine strategy for tr=1 in networks where:
!1) cliques' sizes follow a Kronecker delta distribution and 
!2) each individual belongs to a number of cliques that also follows a Kronecker delta distribution.

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Input
!N_nodeC: total number of cliques
!kminC,kmaxC:   minimum and maximum size for the group of cliques
!kminI,kmaxI:   minimum and maximum "degree", i.e., the number of cliques to which an individual belongs
!nrea  : number of realizations

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!Output
!File: !fraction of recovered population in the plane beta-f

module globals
  implicit none
  save

  integer kminC,kmaxC,kminI,kmaxI
  integer N_node,N_nodeC,N_nodeI
  integer Rtot

  real(8) lambdaC,lambdaI
  real(8) kmedC,kmedI  
  real(8) r  !random number
  real(8) ff,beta

  integer,allocatable,dimension(:)::edge,edgep
  integer,allocatable,dimension(:)::node,nodep
  integer,allocatable,dimension(:)::kk,kkp
  integer,allocatable,dimension(:)::State
  
  real(8),allocatable,dimension(:)::PkC,PkI
    
end module globals

module random
   save
   integer::q1=1,q2=104
   integer ir(670)
   integer::sem(670)
   real(8)::nmax=2*(2**30-1)+1
   real(8) sig
end module random


Program MainProgram
  use globals
  use random

  implicit none

  integer i,j,rea,nrea,nptosf,nptosb,ptosf,ptosb
  real(8) product,sc
  
  real(8),allocatable,dimension(:,:)::MatMax
  integer,allocatable,dimension(:,:)::Countsc
  
  character(8) Var01

  !Parameters
  !N_nodeC=2000000
  N_nodeI=1000000

  nrea     =1000

  sc       =0.01d0

  kminC    =7
  kmaxC    =kminC
  kminI    =3
  kmaxI    =kminI
  
  nrea     =1000
  nptosb   =100
  nptosf   =100

  
  write(Var01,'(I8)') N_nodeI
  call zeros(Var01)  
  
  allocate(PkC(kminC:kmaxC),PkI(kminI:kmaxI))
  allocate(MatMax(nptosb,nptosf),Countsc(nptosb,nptosf))
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !Degree distribution Cliques
  PkC           =0d0
  PkC(kminC)    =1d0
  
  kmedC=0d0
  do i=kminC,kmaxC
     kmedC=kmedC+i*PkC(i)
  enddo

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
  !Degree distribution Individuals
  PkI           =0d0
  PkI(kminI)    =1d0
   
  kmedI=0d0
  do i=kminI,kmaxI
     kmedI=kmedI+i*PkI(i)
  enddo
  
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    
  !N_nodeI=kmedC*(N_nodeC/kmedI)
  N_nodeC=kmedI*(N_nodeI/kmedC)
  N_node=N_nodeC+N_nodeI
    
  print*,"Number of cliques", N_nodeC,"Number of individuals",N_nodeI, "N_nodeC+N_nodeI",N_node


  allocate(kk(N_node),node(N_node),State(N_nodeI),kkp(N_nodeI),nodep(N_nodeI))  

  call initialize_random

  MatMax   =0d0
  Countsc  =0

  do rea=1,nrea
     print*,rea
     call ConfModel     
     do ptosb=1,nptosb
        beta=0.01d0*ptosb
        do ptosf=1,nptosf
           ff=0.01d0*ptosf
           State =1
           call SIR     
           if(Rtot/dble(N_nodeI)>sc)then
              MatMax(ptosb,ptosf)  =MatMax(ptosb,ptosf)+Rtot/dble(N_nodeI)
              Countsc(ptosb,ptosf) =Countsc(ptosb,ptosf)+1
           endif
        enddo
     enddo
     deallocate(edge,edgep)
     open(2,file='PhaseDiagram_beta_and_ff_'//Var01//'.dat')
     write(2,*)'#',rea
     write(2,*)'#','kminC,kmaxC,kminI,kmaxI,lambdaC,lambdaI'
     write(2,*)'#',kminC,kmaxC,kminI,kmaxI,lambdaC,lambdaI 
     do ptosb=1,nptosb
        beta=0.01d0*ptosb
        do ptosf=1,nptosf
           ff=0.01d0*ptosf
           if(Countsc(ptosb,ptosf)>0)then
              write(2,*)beta,ff,MatMax(ptosb,ptosf)/Countsc(ptosb,ptosf)
           else
              write(2,*)beta,ff,0d0
           endif
        enddo
     enddo       
     close(2)  
  enddo
  

end Program MainProgram

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!!   SIR
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine SIR
  use globals
  use random
  implicit none
  integer i,j
  integer sel,nod,vec
  integer Ninf,Niso
  integer Ninfnew
  integer tt
  integer, dimension(N_nodeI):: ListIndInf,ListIndIsol
  integer, dimension(N_nodeI):: ListIndInfNew
  
  !State=0  :susceptible
  !State=1  :infected
  !State=2  :recovered
  !State=3  :isolated
  
  
  State         =0
  Ninf          =0   !number of infected individuals
  Niso          =0   !number of isolated individuals

  ListIndInf    =0   !List of infected individuals
  ListIndIsol   =0   !List of isolated individuals


  Rtot          =0   !number of recovered individuals


  !Selecting the index case
  call rand
  sel  =N_nodeI*r+1
  Rtot  =Rtot+1
  call rand
  if(r<ff)then
     State(sel)=3
     Niso      =Niso+1
     ListIndIsol(Niso)=sel
  else
     State(sel)=1
     Ninf =Ninf+1
     ListIndInf(Ninf) =sel
  endif  
  
  tt=0
  
  do while(Ninf>0)


     !INFECTION STEP:
     !Those infected individuals who have not been isolated, will now transmit the disease to their susceptible neighbors.
     
     Ninfnew       =0
     ListIndInfNew =0
     do i=1,Ninf
        nod =ListIndInf(i)
        if(State(nod)==3)cycle
        do j=1,kkp(nod)
           vec =edgep(nodep(nod)+j-1)
           if(State(vec)==0)then
              call rand
              if(r<beta)then
                 State(vec)  =1
                 Ninfnew =Ninfnew+1
                 Rtot    =Rtot+1
                 ListIndInfNew(Ninfnew)=vec
              endif
           endif
        enddo               
        State(nod)=2
     enddo
     
     Ninf        =Ninfnew
     ListIndInf  =ListIndInfNew
     
     tt    =tt+1
     
     !ISOLATION STEP
     !The new infected individuals are detected and isolated with probability ff
     Niso         =0
     ListIndIsol  =0
     do i=1,Ninf
        nod  =ListIndInf(i)
        call rand
        if(r<ff)then
           Niso =Niso+1
           ListIndIsol(Niso)=nod
           State(nod)=3
        endif
     enddo

     !Isolating neighbors
     do i=1,Niso
        nod =ListIndIsol(i)
        do j=1,kkp(nod)
           vec =edgep(nodep(nod)+j-1)
           if(State(vec)==0.or.State(vec)==1)then
              State(vec) =3    
           endif
        enddo        
     enddo 


  enddo

end subroutine SIR




!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!!   Configuration Model
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!This subroutine creates three arrays: kk, edge, and node, which contain all the network topology information.
!kk(i) is the degree of node "i".
!The subset edge(node(i):node(i)+kk(i)-1) contains the list of neighbors of node "i".

subroutine ConfModel
  use globals
  use random
  implicit none
  integer i,j,k,ll
  integer nod1,lug,nod2,vec
  integer Aux,flag
  integer m,n,stubmaxpos,stubminpos
  integer counting, nstubs,nstubsC,nstubsI
  integer nstubAux,stubposC,stubposI

  real(8) ws

  integer, allocatable, dimension(:)::listStubC,listStubI,kkAux

  real(8),dimension(kminC-1:kmaxC)::CumulativePkC
  real(8),dimension(kminI-1:kmaxI)::CumulativePkI
  
  
  allocate(kkAux(n_node))

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!Cumulative distribution group 1 (cliques)
  
  CumulativePkC   =0d0
  ws             =0d0
  do i=kminC,kmaxC
     ws                = ws+PkC(i)
     CumulativePkC(i)   = ws
  enddo
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!Cumulative distribution group 2 (Individuals)
  
  CumulativePkI   =0d0
  ws             =0d0
  do i=kminI,kmaxI
     ws                = ws+PkI(i)
     CumulativePkI(i)   = ws
  enddo  
  
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!Assigning Connectivity to all nodes
2 kk=0 
  do i=1,N_nodeC
     call rand
     do j=kminC,kmaxC
        if(CumulativePkC(j-1)<r.and.r<=CumulativePkC(j)) then
           kk(i)    = j
           exit
        endif
     enddo
  enddo
  do i=N_nodeC+1,N_node
     call rand
     do j=kminI,kmaxI
        if(CumulativePkI(j-1)<r.and.r<=CumulativePkI(j)) then
           kk(i)    = j
           exit
        endif
     enddo
  enddo  

  nstubsC =sum(kk(1:N_nodeC))
  nstubsI =sum(kk(N_nodeC+1:N_node))
  nstubs  =nstubsC+nstubsI
  
  
  if(abs(nstubsC-nstubsI)>0.01d0*(kmedI*N_nodeI)) goto 2
  
  allocate(listStubC(nstubsC),listStubI(nstubsI))  
  allocate(edge(nstubs))

  listStubC     =0
  listStubI     =0
  nstubAux      =0

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!List of stubs

  do i=1,N_nodeC
     if(kk(i)==0) cycle
     do k=1,kk(i)
        nstubAux           =nstubAux+1
        listStubC(nstubAux)=i
     enddo
  enddo
  
  nstubAux      =0
  do i=N_nodeC+1,N_node
     if(kk(i)==0) cycle
     do k=1,kk(i)
        nstubAux           =nstubAux+1
        listStubI(nstubAux)=i
     enddo
  enddo  

  node(1)       =1
  do i=1,n_node-1
     node(i+1)  = node(i)+kk(i)
  enddo

  edge        =0
  kkAux       =0
  counting    =0
  do while(nstubsC>0.and.nstubsI>0)
3    if(counting==100) then
	deallocate(listStubC,listStubI,edge)
	go to 2
     end if
     
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !!Randomly choose two stubs
     
     call rand
     stubposC    =r*nstubsC+1
     call rand
     stubposI    =r*nstubsI+1

     m        =listStubC(stubposC)   !the first stub corresponds to node "m"
     n        =listStubI(stubposI)   !the second stub corresponds to node "n"

     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     !!Checking whether these two stubs can be connected or not

     if(m.ne.n) then
	do k=1,kkAux(m)
           if(edge(node(m)+k-1)==n)then
              counting     =counting+1
              go to 3   !if nodes "m" and "n" are already connected, we have to choose another pair of stubs
           endif
	enddo
	edge(node(m)+kkAux(m))   =n
	edge(node(n)+kkAux(n))   =m

	kkAux(m)                 =kkAux(m)+1
	kkAux(n)                 =kkAux(n)+1
	listStubC(stubposC)      =listStubC(nstubsC)
	listStubI(stubposI)      =listStubI(nstubsI)
	
	nstubsC             =nstubsC-1
	nstubsI             =nstubsI-1
	counting            =0
     else
        counting           =counting+1
        go to 3 !self-loops are forbidden
     endif
  enddo
  
  kk   =kkAux
  
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!!!Projection
  
  do i=N_nodeC+1,N_node
     nod1 =i-N_nodeC
     Aux  =0
     do j=1,kk(i)
        lug      =edge(node(i)+j-1)
        Aux      =Aux+kk(lug)
     enddo
     kkp(nod1)  =Aux-kk(i)
  enddo    
  
  nodep    =0
  nodep(1) =1

  do i=2,N_nodeI
     nodep(i)=nodep(i-1)+kkp(i-1)
  enddo
  
    
  allocate(edgep(sum(kkp)))
  
  kkp    =0
  edgep  =0

  do i=N_nodeC+1,N_node
     nod1 =i-N_nodeC
     
     do j=1,kk(i)
        lug      =edge(node(i)+j-1)
        do k=1,kk(lug)
           nod2  =edge(node(lug)+k-1)-N_nodeC
           if(nod2==nod1)cycle
           flag  =0
           do ll=1,kkp(nod1)
              vec   =edgep(nodep(nod1)+ll-1)
              if(vec==nod2)then
                 flag =1
                 exit
              endif
           enddo
           if(flag==0)then
              edgep(nodep(nod1)+kkp(nod1))   =nod2
              edgep(nodep(nod2)+kkp(nod2))   =nod1
              kkp(nod1)  =kkp(nod1)+1
              kkp(nod2)  =kkp(nod2)+1
           endif
        enddo
     enddo  
  enddo    
  
  
  deallocate(kkAux,listStubC,listStubI)
end subroutine ConfModel



!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!!   Seed for the Pseudorandom number generator
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine initialize_random
  use globals
  use random
  implicit none

  integer i,see(33)
  integer hour

  CALL SYSTEM_CLOCK(hour)
  !hour=26896048
  see=hour
  print*,hour
  !read(*,*)
  CALL RANDOM_SEED(put=see)		
  CALL RANDOM_SEED(get=see)

  do i=1,670
     call random_number(r)
     sem(i)   =r*nmax
     ir(i)    =i+1
  enddo
  ir(670)=1
  return
end subroutine initialize_random
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!!   Pseudo-random number generator
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine rand
  use globals
  use random
  implicit none

 1 q1=ir(q1)
  q2=ir(q2)
  sem(q1)= IEOR(sem(q1),sem(q2))
  r=dfloat(sem(q1))/nmax
  if(r==1d0) go to 1
  return
end subroutine rand


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!
!!!   Zero Character
!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

subroutine zeros(b)
 implicit none
 
 character(*):: b
 integer i
 
 do i=1,len(b)
   if (b(i:i)==' ') then  
      b(i:i)='0'
   endif
 end do
end subroutine zeros


