# raw data

t=c(5,4,6,3.5,4.5,7,0.4,6,8,6,7,4.8,5.2)
married=c(rep(1,11),rep(0,2))

library(survival)

km=survfit(Surv(t,married) ~ 1 );
str(km)

plot(km)

km$surv

# k-m by hand

s=c(12,11,10,9,7,3,1,0)
n=c(13,12,11,10,8,6,3,1)
y=c(0.4,3.5,4,4.5,5,6,7,8)
si.hat = s/n

s.hat=cumprod(si.hat)

cbind(y,s.hat)

plot(c(0,y),c(1,s.hat), xlim=c(0,10),ylim=c(0,1),type="s")


setwd("/Users/jameshanley/Documents/work/Espionage")

source("ExtractLinesFromLargeFileFunction.R")

result = ExtractLinesFromLargeFileFunction("oncologyFellowsMarriage.ps",90000,10000)

LENGTH=unlist(result[3]) ; LENGTH ; summary(LENGTH)

n.l = length(LENGTH) ; n.l

X = result[1][[1]] ; Y = result[2][[1]] ; length(X)

x.min=unlist(result[4]) ; y.min=unlist(result[5])
x.max=unlist(result[6]) ; y.max=unlist(result[7])
c(x.min,y.min,x.max,y.max)

x.range=c(40,255) ; y.range = c(85,215)

plot(x.range,y.range,col="white", xlim=x.range, ylim=y.range)
for(i in 1:length(X) )  lines(X[[i]],Y[[i]] )


# focus on tick marks 


plot(x.range,y.range,col="white", xlim=x.range, ylim=y.range)
for(i in 1:18 ) {
	L = LENGTH[i]	
	xa=0; if (X[[i]][1] > X[[i]][L]) xa=1
	ya=0; if (Y[[i]][1] > Y[[i]][L]) ya=1
	if(L==2) text(X[[i]][1],Y[[i]][1], toString(i), cex=0.5,adj=c(xa,ya) ) 
	if(L==2) lines(X[[i]],Y[[i]]) ;
	j = 1 + (i - 5*floor(i/5))
	lines(X[[i]],Y[[i]], lwd=0.5) 
   if(L==2) text(X[[i]][L],Y[[i]][L], toString(i), cex=0.5,adj=c(xa,ya) ) 
}

i.x.10 = 13; i.x.2 = 9
x1 = ( X[[i.x.10]][1] - X[[i.x.2]][1] ) / 8 ; x1

i.y.10=18; i.y.8 = 17; i.y.2=14 ; i.y.0=7

y1 =  (Y[[i.y.8]][1] - Y[[i.y.2]][1])/ 0.6  ; y1

y1 =  (Y[[i.y.10]][1] - Y[[i.y.0]][1])/ 1  ; y1

Y.s.0 = Y[[i.y.8]][1] - 0.8*y1 ; Y.s.0

Y.s.0 = Y[[i.y.0]][1]; Y.s.0

X.t.0 = X[[i.x.10]][1] - 10*x1 ; X.t.0

s = (Y[[19]] - Y.s.0)/y1 ; s = s[seq(1,17,2)] ; s ; length(s)

floor(cbind(s, c(1,s.hat)),2)
floor(100*cbind(s, c(1,s.hat)))

yr = (X[[19]] - X.t.0)/x1 ; yr=yr[seq(1,17,2)] ; yr

plot(yr,s,type="l",col="grey70")

jump = s[seq(1,15,2)]  - s[seq(3,17,2)] ;jump ; 1/jump

ratio = s[2:9] / s[1:8] ; 1/(1-ratio)

log.s = - log( s[seq(1,17,2)] ) ; log.s

j =  log.s[2:9] - log.s[1:8] ; j ; 1/j



1/13