/*
 
Function: Compute  NPMLE
              interval censored data.

*/

#include <sys/time.h>
#include <sys/resource.h>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <string.h>


#define MAX_N_Obs  3000

static int comp(double *f1, double *f2)
{
 if (*f1 < *f2) return(-1);
    else if (*f1 > *f2) return(1);
            else return(0);
}

/* L_infinite distance */
double Linf_dist(double *a, double *b, int n)
{
 int i;
 double s1, s=0;

 for(i=1; i<=n; i++) {
    s1 = (s1=a[i]-b[i])>0? s1 : -s1;
    if (s < s1) s=s1;
   }
/*
    s += (a[i]-b[i])*(a[i]-b[i]);
 s = (float)sqrt((double)s);
*/
 return(s);
}

/* copied from mono mle1.2.c, Nov 8, 96
   make a Greatest Convex Minorant from
   interval hazards from output of N-C-P.

   March 14, 96
*/

gcm(q,nh,n)
double *q, *nh;
int n;
{
 int i,j,k,last;
 double  min_slope, slope;

 nh[0]=q[0]=0;
 for(i=1; i<=n; i++){
   nh[i] = nh[i]*q[i] + nh[i-1];
   q[i] += q[i-1];
   }
 i=0;

 while(i<n-1){
    min_slope=(nh[i+1]-nh[i])/(q[i+1]-q[i]);
    last=i+1;
    for(j=i+2; j<=n; j++){
       slope=(nh[j]-nh[i])/(q[j]-q[i]);
       if (slope<min_slope){
          /*
          printf("i=%d j=%d min_slope=%f    slope=%f\n",i,j,min_slope,slope);
          */
          min_slope=slope;
          last=j;
          }
       }

    for(k=i+1; k<last; k++)
       nh[k]=nh[i]+ (double)min_slope*(q[k]-q[i]);

    i=last;
    }

 for(i=n; i>=1; i--){
   q[i] -= q[i-1];
   nh[i] = (nh[i]-nh[i-1])/q[i];
   }
}

#define TB_EPSILON     1.0e-100
#define LM_EPSILON      1.0e-10
#define RIGHTCENSORING     9999
#define LARGENUM       1.0e+10
#define TB_MAX_ITER    500 

#define NdiscreteT 11

void IcmNPMLE(double * LC, double *RC, long *Npt)
{
char OUTfname[100],OUTfname2[100];
int i,j,k,N,n,n0, iter, iter2;
int INDlt[MAX_N_Obs], INDlc[MAX_N_Obs], INDrc[MAX_N_Obs];
    /* indices of the corrsponding LT,LC,RC in ordered y[] */
long perm[MAX_N_Obs];
double Tolerence=0.00001, M ;
double pTolerence=0.01 ;
double *tp1,
      y[MAX_N_Obs],
      wt[MAX_N_Obs],
      dW[MAX_N_Obs], dG[MAX_N_Obs], dpis[MAX_N_Obs];
double lr[2*MAX_N_Obs], d1[2*MAX_N_Obs], d2[2*MAX_N_Obs], 
       pis[2*MAX_N_Obs],npis[3*MAX_N_Obs];
double L0, Li, li, nu, dn, stepsize;
FILE *fp;
struct rusage ru; /* resource utilization, such as time */

N= *Npt;
j=1;
for(i=1; i<=N; i++){
   lr[i]=LC[i];
   lr[i+N]=RC[i];
  }

tp1=&(lr[1]);
qsort(tp1, N,sizeof(double), comp);
tp1=&(lr[1+N]);
qsort(tp1, N,sizeof(double), comp);

n=0; i=1; j=N+1;
y[0]=0;
while (i<=N && j<=2*N ){
  if (lr[i]>=lr[j]) j++;
     else if (i<N && lr[i+1]<lr[j]) i++; /* ==? */
	     else {
		   n=n+1;
                   y[n]=lr[j];
                   i++; j++;
                   /* added on Sun Nov 10 18:49:49 CST 1996 */
                   while (lr[j]== y[n] && j<=2*N) j++;
                  }
  }

for(i=1; i<=N; i++){
 /* INDlt[i]=-1; */
 INDlc[i]=INDrc[i]=-1;
 for(j=1; j<=n; j++) {
   /*
   following ``=" in LC[i]<=y[j] added on April 22, 97: [LC, RC)
   */
   if (INDlc[i]==-1 && LC[i]<=y[j] ) INDlc[i]=j-1;
   if (INDrc[i]==-1 && RC[i]<y[j] ) {
       INDrc[i]=j-1;
       break;
       }
   }
   if (INDrc[i]==-1) INDrc[i]=n;
if(INDlc[i]==-1 || INDrc[i]==-1)
   perror("Indices wrong!!!!"); 
 } 
/* initial cumulative hazards */
npis[0]=pis[0]=0;
for(j=0; j<=n; j++) {
  pis[j]=j*3.0/(n);
  wt[j]=j;
  }

L0=0;
for(i=1; i<=N; i++){
  dn = exp(-pis[INDlc[i]]);
  if (RC[i]<RIGHTCENSORING) 
          dn -= exp(-pis[INDrc[i]]);
  li=dn;
  if (li<TB_EPSILON) {
/*
     printf("li from obs %d = %lf\n", i, li);
*/
     li=TB_EPSILON;
     }
  L0 += log(li);
  }

iter=0;
while(iter<TB_MAX_ITER) {
iter=iter+1;

for(j=1; j<=n; j++){
 dW[j]=dG[j]=0;
 for(i=1; i<=N; i++){
   if (RC[i]>=RIGHTCENSORING){
      if (INDlc[i]==j) dW[j] -= 1;
      }
   else {
         if (INDlc[i]==j ) {
             dW[j] -= exp(-pis[j])/(exp(-pis[j])-exp(-pis[INDrc[i]]));
             dG[j] += exp(-pis[j]-pis[INDrc[i]])/((exp(-pis[j])-exp(-pis[INDrc[i]]))*(exp(-pis[j])-exp(-pis[INDrc[i]])));
             }
         if (INDrc[i]==j ) {
             dW[j] += exp(-pis[j])/(exp(-pis[INDlc[i]])-exp(-pis[j]));
             dG[j] += exp(-pis[j]-pis[INDlc[i]])/((exp(-pis[j])-exp(-pis[INDlc[i]]))*(exp(-pis[j])-exp(-pis[INDlc[i]])));
             }
        } 
   }
 }

for(j=1; j<=n; j++) {
  if (dG[j]>LM_EPSILON) dpis[j] = dW[j]/dG[j];
  else {
        dpis[j] =  dW[j]/LM_EPSILON ;
        dG[j]=LM_EPSILON;
        }
  }
	
stepsize=1;
while (stepsize>1e-10){
for(j=1; j<=n; j++) {
  npis[j] = pis[j]+ stepsize*dpis[j];
  }

gcm(dG, npis, n);
for(j=1; j<=n; j++) {
  if (npis[j]<0) npis[j] = 0;
  if (isinf(npis[j]) || isnan(npis[j])) npis[j] = 999;
/*
  if (isinf(npis[j])) npis[j] = 999;
*/
  }

npis[0]=0;

Li=0;
for(i=1; i<=N; i++){
  dn = exp(-npis[INDlc[i]]);
  if (RC[i]<RIGHTCENSORING) 
          dn -= exp(-npis[INDrc[i]]);
  li=dn;
  if (li<TB_EPSILON) li=TB_EPSILON;
  Li += log(li);
  }
/*
printf("Iteration=%d : DeltalogL=%lf Li=%lf  L0=%lf  \n", iter, Li-L0, Li,L0);
*/
if (Li>L0) break;
stepsize /= 2;
}
/*
for(j=1; j<=n; j++) printf(" %f ", npis[j]);
printf(" npis[n]= %lf \n", npis[n]);
*/

if (Li-L0<Tolerence && Linf_dist(pis,npis,n)<pTolerence) break;
   else 
       for(j=1; j<=n; j++)
        pis[j]=npis[j];
L0=Li;
}
if (Li>L0){
    for(j=1; j<=n; j++)
        pis[j]=npis[j];
    L0=Li;
    }

if (iter>=TB_MAX_ITER){
   printf("\n\nWARNING: No convergence after %d iterations!!\n\n",iter);
   for(j=1; j<=n; j++) printf("\n %lf ", pis[j]);
   printf("\n");
   }

/* returning values */
pis[0]=0;
for(j=1; j<=n; j++){
  LC[j]=y[j];
  RC[j] = exp(-pis[j-1])-exp(-pis[j]);
  if (RC[j]<1.e-10) RC[j]=0;
  }
*Npt = n;
}