# -----------------------------------------------------------------------
# argomenti per la funzione:

# tau   vettore con il tempo di osservazione di ciascun soggetto
# tempi matrice (righe = soggetti) con i tempi di accadimento
# nr    vettore con il numero di ricorrenze per ciascun soggetto
# xcov      matrice delle covariate
# betastart (opzionale) TRUE se parte con la Poisson


# prepara le due matrici di lavoro
# per entrambe le RIGHE sono uguali al numero di SOGGETTI
# e le COLONNE al TEMPO di osservazione PIU' GRANDE
# ATTENZIONE: si considerano soltanto tempi interi e maggiori di zero

# nella PRIMA c'è UNO se il soggetto (riga) è sotto
# osservazione al tempo t (colonna)
# nella SECONDA c'è il numero di recidive osservate
# sul soggetto (riga) al tempo t (colonna)


mfreg <- function(tau,tempi,nr,xcov,betastart=FALSE)
{

    k      <- length(tau)
    taumax <- max(tau)



    delta <- matrix(0,nrow=k,ncol=taumax)
    enne  <- delta
    for (i in 1:k)
    {
        delta[i,(1:tau[i])] <- 1
        if (nr[i]>0) {
           ttt <- tempi[i,]
           abc <- ttt>0             # tempi di ricorrenza (ev. duplicati)
           def <- ttt[abc]          # posizione dei tempi
           ghi <- table(def)        # numero di ricorrenze ai vari tempi
           jkl <- sort(unique(def)) # tempi di ricorrenza distinti
           enne[i,jkl] <- ghi
        }
    }


    Delta <- apply(delta,2,sum) # totale soggetti sotto osservazione al tempo t
    Enne  <- apply(enne, 2,sum) # totale ricorrenze osservate al tempo t


    # -------------------------------------------------------------------
    idtempi <- c(1:taumax)
    ok <- which(Enne>0)
    delta <- delta[,ok]
    enne  <- enne[,ok]
    Delta <- Delta[ok]
    Enne  <- Enne[ok]
    Tempi <- ncol(delta)
    idtempi <- idtempi[ok]
    # -------------------------------------------------------------------
    p <- ncol(xcov)
    b <- rep(0,p)
    if (betastart) {
       fit <- glm(nr~xcov,family="poisson")
       b   <- fit$coef[-1]
    }
    lista <- list(xcov,delta,enne,Enne)
    out <- nsolve(fun,jac,b,lista)
    # -------------------------------------------------------------------
    # varianze e covarianze
    b <- out
    expxtb <- exp(xcov %*% b)
    Erre   <- apply(sweep(delta,1,expxtb,"*"),2,sum) # R(t)
    emmezero <- Enne/Erre
    p <- length(b)
    W <- array(dim=c(k,Tempi,p))
    B <- matrix(nrow=k,ncol=p)
    for (j in 1:p)
    {
        tmp <- xcov[,j] * expxtb
        tmp <- sweep(delta,1,tmp,"*")
        tmp <- apply(tmp,2,sum)
        tmp <- tmp/Erre
        aaa <- matrix(rep(xcov[,j],Tempi),ncol=Tempi,byrow=FALSE)
        bbb <- matrix(rep(tmp     ,k)     ,nrow=k     ,byrow=TRUE)
        W[,,j] <- aaa - bbb

        tmp <- enne - expxtb %*% t(emmezero)
        tmp <- delta * W[,,j] * tmp
        B[,j] <- apply(tmp,1,sum)
    }
    B1 <- (t(B) %*% B) / k

    tmp <- sweep(delta,2,emmezero,"*")
    A1 <- matrix(nrow=p,ncol=p)
    for (j in 1:p)
    {
        for (l in 1:p)
        {
            aaa <- matrix(rep((xcov[,j]*expxtb),Tempi),ncol=Tempi,byrow=FALSE)
            bbb <- W[,,l]
            ccc <- tmp * aaa * bbb
            A1[j,l] <- sum(ccc)/k
        }
    }

    invA1 <- solve(A1)
    asvar <- (invA1 %*% B1 %*% t(invA1)) / k
    se    <- sqrt(diag(asvar))
    z     <- abs(out/se)
    pval  <- pnorm(z,lower.tail=FALSE)*2

    ris <- list(estimates = cbind(out,se,z,pval),
                 asvar = asvar,
                 basemf = emmezero,
                 times = idtempi)

    ris
}
# -----------------------------------------------------------------------




# -----------------------------------------------------------------------
# calcola il valore numerico delle equazioni 3.5
# avendo come input un vettore di coefficienti di regressione
# e le matrici Xcov, delta, enne e il vettore Enne
fun <- function(b,lista)
{
    ncov  <- length(b)
    f     <- rep(NA,ncov)
    X     <- lista[[1]] # covariate
    delta <- lista[[2]] # a rischio sì/no
    enne  <- lista[[3]] # matrice soggetti x tempi
    Enne  <- lista[[4]] # totale eventi ai singoli tempi
    
    expxtb  <- exp(X %*% b) # risk score per i singoli soggetti
    Erre    <- apply(sweep(delta,1,expxtb,"*"),2,sum) # R(t)
    
    for (j in 1:ncov)
    {
        tmp <- enne - expxtb %*% t(Enne/Erre)
        tmp <- delta * tmp
        tmp <- sweep(tmp,1,X[,j],"*")
        f[j] <- sum(tmp)
    }
    f
}
# -----------------------------------------------------------------------




# -----------------------------------------------------------------------
# calcola il valore numerico degli elementi dello jacobiano
# avendo come input un vettore di coefficienti di regressione
# e le matrici Xcov e delta e il vettore Enne
# (in realtà la matrice enne non viene usata)
jac <- function(b,lista)
{
    ncov  <- length(b)
    jcb   <- matrix(nrow=ncov,ncol=ncov)
    X     <- lista[[1]] # covariate
    delta <- lista[[2]] # a rischio sì/no
    Enne  <- lista[[4]] # totale eventi ai singoli tempi
    
    expxtb  <- exp(X %*% b) # risk score per i singoli soggetti
    Erre    <- apply(sweep(delta,1,expxtb,"*"),2,sum) # R(t)
    
    for (j in 1:ncov)
    {
        for (l in 1:j)
        {
            tmp <- expxtb %*% t(Enne/Erre)
            tmp <- delta * tmp
            tmp <- sweep(tmp,1,(X[,j]*X[,l]),"*")
            add.1 <- sum(tmp)
            
            tmp <- sweep(delta,1,(X[,l]*expxtb),"*")
            Somma <- apply(tmp,2,sum)
            tmp <- Enne * Somma / (Erre^2)
            tmp <- expxtb %*% t(tmp)
            tmp <- delta * tmp
            tmp <- sweep(tmp,1,X[,j],"*")
            add.2 <- sum(tmp)
            jcb[j,l] <- add.2 - add.1
            jcb[l,j] <- jcb[j,l]
        }
    }
    jcb
}
# -----------------------------------------------------------------------




# -----------------------------------------------------------------------
# argomenti per la funzione:

# tau   vettore con il tempo di osservazione di ciascun soggetto
# tempi matrice (righe = soggetti) con i tempi di accadimento
# nr    vettore con il numero di ricorrenze per ciascun soggetto

# prepara le due matrici di lavoro
# per entrambe le RIGHE sono uguali al numero di SOGGETTI
# e le COLONNE al TEMPO di osservazione PIU' GRANDE
# ATTENZIONE: si considerano soltanto tempi interi e maggiori di zero

# nella PRIMA c'è UNO se il soggetto (riga) è sotto
# osservazione al tempo t (colonna)
# nella SECONDA c'è il numero di recidive osservate
# sul soggetto (riga) al tempo t (colonna)

cmfplot <- function(tau,tempi,nr,conf=0.95,out=FALSE,plot=TRUE,plot.conf=TRUE)
{
    k      <- length(tau)
    taumax <- max(tau)

    delta <- matrix(0,nrow=k,ncol=taumax)
    enne  <- delta
    for (i in 1:k)
    {
        delta[i,(1:tau[i])] <- 1
        if (nr[i]>0) {
           ttt <- tempi[i,]
           abc <- ttt>0             # tempi di ricorrenza (ev. duplicati)
           def <- ttt[abc]          # posizione dei tempi
           ghi <- table(def)        # numero di ricorrenze ai vari tempi
           jkl <- sort(unique(def)) # tempi di ricorrenza distinti
           enne[i,jkl] <- ghi
        }
    }

    Delta <- apply(delta,2,sum) # totale soggetti sotto osservazione al tempo t
    Enne  <- apply(enne, 2,sum) # totale ricorrenze osservate al tempo t
    m     <- Enne/Delta
    M     <- cumsum(m) # questa è la Cumulative Mean Function

# stima della varianza
   abc  <- sweep(delta,2,Delta,"/")
   def  <- sweep(enne,2,m    ,"-")
   ghi  <- abc*def
   ghi  <- t(apply(ghi,1,cumsum))
   ghi  <- ghi^2
   vart <- apply(ghi,2,sum)
   z    <- qnorm(conf+(1-conf)/2)
   ermi <- z*sqrt(vart)
   
   if (!plot) plot.conf=FALSE

   if (plot) {
      ymax <- max(M+ermi)
      plot(c(0,taumax),c(0,ymax),type="n",xlab="Times",ylab="Cumulative Mean Function")
      lines(c(0:taumax),c(0,M),)
   }

   if (plot.conf) {
      lines(c(0:taumax),c(0,M-ermi),lty=2)
      lines(c(0:taumax),c(0,M+ermi),lty=2)
   }
   
   if (out) cbind(c(1:taumax),M,sqrt(vart),M-ermi,M+ermi,Delta,Enne)
}




# -----------------------------------------------------------------------
nsolve <- function(f,jac,x,lista=NA,eps=10^(-6),maxit=20)
{
    error <- 2*eps
    niter <- 0

    repeat
    {
       F  <- f(x,lista)
       error <- max(abs(F))
       if (error<eps) break
       J  <- jac(x,lista)
       dx <- solve(J,-F)
       i <- 0
       repeat
       {
          xx <- x + 0.5^i*dx
          FF <- f(xx,lista)
          if (sum(FF^2)<sum(F^2)) break
          cat("Ciclo interno",i,xx,sum(FF^2),sum(F^2),"\n")
          i <- i+1
       }
       x <- x + 0.5^i*dx
       niter <- niter+1
       cat(niter,error,"\n")
       if (niter>=maxit) break
    }
    cat("--------------------\n",F,"\n")
    x
}
# -----------------------------------------------------------------------