# aug 26 setwd("/Users/jameshanley/Documents/work/Espionage") source("ExtractLinesFromLargeFileFunction.R") result = ExtractLinesFromLargeFileFunction("BisphosphonatesFig2.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(310,530) ; y.range = c(525,685) 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:40 ) { 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 text(X[[i]][1],Y[[i]][1], toString(i), cex=0.5,adj=c(xa,ya) ) lines(X[[i]],Y[[i]]) ; text(X[[i]][L],Y[[i]][L], toString(i), cex=0.5,adj=c(xa,ya) ) } # zoom on rightmost portion x.zoom=c(490,515) ; y.zoom = c(632,638) plot(x.range,y.range,col="white", xlim=x.zoom, ylim=y.zoom) for(i in 1:n.l ) lines(X[[i]],Y[[i]],lwd=0.25) x.zoom=c(490,515) ; y.zoom = c(632,638) for(i in 1:n.l ) if (X[[i]][1] > x.zoom[1] & X[[i]][1] < x.zoom[2] & Y[[i]][1] > y.zoom[1] & Y[[i]][1] < y.zoom[2]) { print( c( X[[i]], Y[[i]]) ) } i.x.84 = 20; i.x.0 = 14 x1 = ( X[[i.x.84]][1] - X[[i.x.0]][1] ) / 84 ; x1 i.y.1 = 8; i.y.0=13 y1 = Y[[i.y.1]][1] - Y[[i.y.0]][1] ; y1 X.sd = as.vector(sapply(X,function(x) sd(x))) Y.sd = as.vector(sapply(Y,function(x) sd(x))) dot = (Y.sd==0 & X.sd==0) ; sum(dot) vertical.line = (Y.sd > 0 & X.sd == 0 ) ; n.jumps=sum(vertical.line) ; n.jumps i.jumps = (1:n.l)[vertical.line] ; for(i in i.jumps ) { lines(X[[i]],Y[[i]], col="red") } h.line = (Y.sd == 0 & X.sd > 0 ) ; h.jumps=(1:n.l)[h.line] for(i in h.jumps ) { lines(X[[i]],Y[[i]], col="blue") } i.dot = i.jumps = (1:n.l)[dot] t.dot = ( unlist(X[i.dot]) - X[[i.y.0]][1] ) / x1 n.t = table(1+round(t.dot,1)) mean(t.dot > 10) p=rep(NA,168); for(i in 1:168) p[i]=mean(t.dot > i/12) plot(1:168,p) plot(1:length(n.t),n.t) y.jumps = Y[i.jumps] ; x.jumps = X[i.jumps] # plot jumps i.x.14 = 20; i.x.0 = 13 x1 = ( X[[i.x.14]][1] - X[[i.x.0]][1] ) / 14 ; x1 i.y.02=6; i.y.0=12 y1 = 50*( Y[[i.y.02]][1] - Y[[i.y.0]][1] ) ; y1 ci.jumps = ( y.jumps.1 - Y[[i.y.0]][1]) at.risk=c(); jump=c() ; t=c() ; ci=c() plot(x.range,y.range,col="white", xlim=c(330,530), ylim=y.range) for(I in 1:n.jumps ) { i = i.jumps[I] 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 lines(X[[i]],Y[[i]]) ; ci.0 = ( Y[[i]][1] - Y[[i.y.0]][1] ) / y1 ci.1 = ( Y[[i]][2] - Y[[i.y.0]][1] ) / y1 #print( c(ci.0, ci.1) ) ratio = (1-ci.1) / (1-ci.0) jmp = ci.1 - ci.0 year = round( ( X[[i]][1] - X[[i.y.0]][1] ) / x1, 3) at.r = 1/(1-ratio) #print(c(t,at.r)) if( year > 0 & year < 15 & jmp > 0 ) { at.risk=c(at.risk,at.r) jump=c(jump,jmp) t=c(t,year) ci=c(ci,ci.1) } } ci jump at.risk t length(t) i.1= 1:112 ; i.0=113:262 t.0=t[i.0] ; t.0 ; at.risk.0 = at.risk[i.0] ; ci.0=ci[i.0] t.1=t[i.1] ; t.1 ; at.risk.1 = at.risk[i.1] ; ci.1=ci[i.1] plot(t.0,ci.0,col="red",type="s",xlim=c(0,12),ylim=c(0,0.01)) points(t.1,ci.1,col="blue",type="s") plot(t.0,at.risk.0,col="blue") points(t.1,at.risk.1,col="red") N.1=c(65078,58902,20288) N.0=c(80101,73534,23758); T.0=c(5,7,10) plot(T.0,N.1,col="red",cex=2,pch=19) points(T.0,N.0,col="blue",cex=2,pch=19)