# aug 18

#Fig1A 

setwd("/Users/jameshanley/Documents/work/Espionage")
source("ExtractLinesFunction.R")

result = ExtractLinesFunction("AtrialFibrillation.ps")
X=result[[1]]; Y=result[[2]] ; LENGTH=result[[3]]

# flip vertical 

x.min=min(unlist(X));  x.max=max(unlist(X))
y.min=min(unlist(Y));  y.max=max(unlist(Y))

NP=length(X) ; for (i in 1:NP) Y[[i]] = y.max-Y[[i]]

y.min=min(unlist(Y));  y.max=max(unlist(Y))

c(x.min,x.max,y.min,y.max)


# plot selected sections of page
par(mfrow=c(1,1),mar = c(.01,.01,.1,.1))
 # x.range=c(0.01,0.99)*x.max ; y.range=c(0.01,0.99)*y.max
 # x.range=c(0.08,0.74)*x.max ; y.range=c(0.11,0.90)*y.max
 # x.range=c(0.05,0.60)*x.max ; y.range=c(0.71,0.90)*y.max
 # x.range=c(0.10,0.55)*x.max ; y.range=c(0.58,0.87)*y.max
x.range=c(0.08,0.92)*x.max ; y.range=c(0.40,0.90)*y.max

plot(x.range,y.range,col="white")

range=1:NP #range=1:400 range=400:800 range=c(1:50, 1402:1707)

for (i in range) {
	xa=0; if (X[[i]][1] > X[[i]][length(X[[i]])]) xa=1
	ya=0; if (Y[[i]][1] > Y[[i]][length(Y[[i]])]) ya=1
	text(X[[i]][1],Y[[i]][1], toString(i), cex=0.5,adj=c(xa,ya) ) 
	LL=LENGTH[i]
	if (LL > 4) {  
		lines(X[[i]],Y[[i]]) ; 
		text(X[[i]][LL],Y[[i]][LL], toString(i), cex=0.5,adj=c(xa,ya) ) }
}

LENGTH

T=seq(0,840,120)

res="hi";

# hi res

if(res=="hi") {
index=89; if(index==89) N=c(6958,6211,5786,5468,4406,3407,2472,1496)
index=90; if(index==90) N=c(7004,6327,5911,5542,4461,3478,2519,1538) 
i.x.zero=43 ; i.x.max=50 ; i.y.zero=42 ; i.y.max=40 ; yd=20
}

# lo res
if(res=="lo") {
index=87; if(index==87) N=c(6958,6211,5786,5468,4406,3407,2472,1496)
index=88; if(index==88) N=c(7004,6327,5911,5542,4461,3478,2519,1538)
i.x.zero=28 ; i.x.max=35 ; i.y.zero=27 ; i.y.max=17 ; yd=1
}

dx=(X[[i.x.max]][1]-X[[i.x.zero]][1])/840 ; dx; 
x.1all=(X[[index]] - X[[i.x.zero]][1])/dx ; x.1all

dy=yd*(Y[[i.y.max]][1]-Y[[i.y.zero]][1]) ; dy
c.1all= (Y[[index]] - Y[[i.y.zero]][1])/dy ; c.1all

yy = 1-c.1all ; nn=length(yy) ; nn

y=c(1,yy) ; length(y)
t=c(0,x.1all) ;  length(t)

I.hi=c()
I.lo=c()
I.t=c()
Y.hi=c()
Y.lo=c()
for (i in 2:length(y)) {
	if( ( y[i] < y[i-1]) & (t[i]==t[i-1]) ) {
		I.hi=c(I.hi,i-1) ; Y.hi=c(Y.hi,y[i-1])
		I.lo=c(I.lo,i) ;   Y.lo=c(Y.lo,y[i])
		I.t=c(I.t,t[i])
		}
	}

I.hi ; I.lo ; I.t ; Y.hi ; Y.lo
I= c(I.hi , I.lo ) ; I ; length(I)

ratio=Y.lo / Y.hi ;  n.ratios=length(ratio)
n.risk=round( 1/(1-ratio),0) ; n.risk

plot(I.t,n.risk,ylim=c(0,7500),pch=19,cex=0.5)
points(T,N,col="red",cex=1,pch=19)

jump = Y.hi - Y.lo ; multiple=c() 
for ( i in 2:length(jump) ) {
	hit=0;
	if( jump[i] > 1.5*jump[i-1] )  hit =1  
	if(i > 2) if( jump[i] > 1.5*jump[i-2] ) hit=1
    if(i == length(jump)) if (jump[i] > 1.5*jump[i-2] ) hit=1
     if(hit==1) multiple = c(multiple,i) 
}
n.mult=length(multiple)
multiple
n.mult

fitted=approx(I.t[-multiple], n.risk[-multiple], I.t[multiple] )$y

n.events = rep(1,length(jump)) ; n.events

n.events[multiple] = round(fitted/ n.risk[multiple],0); 
n.events

sum(n.events)

n.risk[multiple] = n.risk[multiple] * n.events[multiple]

plot(I.t,n.risk,ylim=c(0,7500),pch=19,cex=0.5)
# plot(I.t,n.risk,ylim=c(0,7500),pch=19,cex=0.5,type="l")
points(T,N,col="red",cex=2,pch=19)


d.t = I.t - c(0,I.t[1:(length(I.t)-1)]) ; d.t
pt=sum(n.risk*d.t) / 365.25

sum(n.events)/pt

####### to see some of the censored obsns (not all are given .. but even then too dense unless thin them out )

y.rest=y[-I] ; length(y.rest)
x.rest=t[-I] ; length(x.rest)
N=length(x.rest) ; selected=seq(1,N,20)
plot(x.rest,y.rest,col="white",pch=19,cex=0.25)
for (i in 1:length(Y.hi))  segments(I.t[i],Y.hi[i],I.t[i],Y.lo[i],lwd=1,col="red")
points(x.rest[selected],y.rest[selected],col="blue",pch=19,cex=0.25)