First attent to upload plum

.Rbuildignore

@@ -0,0 +1,2 @@
+^.*\.Rproj$
+^\.Rproj\.user$

.gitignore

@@ -0,0 +1,4 @@
+.Rproj.user
+.Rhistory
+.RData
+.Ruserdata

DESCRIPTION

@@ -0,0 +1,14 @@
+Package: Plum
+Type: Package
+Title: 210Pb chronology
+Version: 0.1.0
+Author: Marco A. Aquino-Lopez
+Maintainer: Marco A. Aquino-Lopez <maquinolopez01@qub.ac.uk>
+Depends:
+    R (>= 3.0),
+    rPython
+Description: More about what it does (maybe more than one line)
+    Use four spaces when indenting paragraphs within the Description.
+License: What license is it under?
+Encoding: UTF-8
+LazyData: true

NAMESPACE

@@ -0,0 +1,3 @@
+exportPattern("^[[:alpha:]]+")
+importFrom(rPython, "python.load")
+

Plum.Rproj

@@ -0,0 +1,17 @@
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source

R/Data_creator.R

@@ -0,0 +1,124 @@
+
+
+timefunction<- function (x)(  (x^2)/3 + x/2)
+#'export
+Data_sim=function(tfuntion=timefunction, depths=1:30,supp=15,thick=1,Fi=100,rho=0,errors=TRUE){
+  dptS=depths[length(depths)]
+  lambda=0.03114
+  if (typeof(supp)=="character"){
+    supported=supp(depths)
+  }else{
+    supported=rep(supp,length(depths))}
+  ########Functions #########
+
+  t1 =tfuntion#<- function (x)(  (x^2)/3 + x/2) #  8*x-20*sin(x/(1.*pi)) )#
+  curve(t1,0,(dptS*.8))
+  abline(h=250)
+  abline(v=30)
+
+  dt <- function(x){}
+  body(dt) <- D(body(t1), 'x')
+  if(typeof(rho)=="double"){
+    rho <- function(x){
+      1.5-.5*cos(x/30*pi)
+    }
+  }else{rho=rho}
+  curve(rho,0,dptS)
+
+  r.sed <- function(x){ rho(x)/dt(x)}
+  curve(r.sed,0,dptS)
+
+  C0<- function(x){  Fi/r.sed(x)}
+  curve(C0,0,dptS)
+
+
+  f <- function(x){
+
+    return ((C0(x)*exp(-lambda*t1(x))))
+  }
+  curve(f,0,dptS)
+
+  f2 <- function(x){
+    return  ( (C0(x)*exp(-lambda*t1(x)))*rho(x))
+  }
+  curve(f2,0,dptS)
+
+  ##########Sampling #######
+
+  sample1=c()
+  density=c()
+  k=1
+  for (i in depths){
+    sample1=c(sample1,(integrate(f,(i-thick),i)$value+supported[k]))
+    density=c(density,integrate(rho,(i-thick),i)$value)
+
+    k=k+1
+  }
+
+  if(errors==TRUE){
+    uncer=seq(10,1,length.out = length(sample1))
+    sample1un=c()
+    unc=c()
+    suppo=abs(rnorm(length(supported),supported,rep(3,length(supported))) )
+    unc.cons=3
+    k=1
+    for (i in sample1){
+      sample1un=c(sample1un,i+rnorm(1,0,uncer[k]))
+      k=k+1
+    }
+  }else{
+    if(any(typeof(errors)=="double",typeof(errors)=="integer") ){
+      if (length(errors)==1){
+        sample1un=c()
+        suppo=abs(rnorm(length(supported),supported,rep(3,length(supported))) )
+        unc.cons=3
+        k=1
+        uncer=rep(errors,length(depths))
+        for (i in sample1){
+          sample1un=c(sample1un,i+rnorm(1,0,errors))
+          k=k+1
+        }
+      } else{if(length(errors)==length(depths)){
+        sample1un=c()
+        suppo=abs(rnorm(length(supported),supported,rep(3,length(supported))) )
+        unc.cons=3
+        uncer=errors
+        k=1
+        for (i in sample1){
+          sample1un=c(sample1un,i+rnorm(1,0,errors[k]))
+          k=k+1
+        }
+    }else{
+      print("Errors should either have the same lenth as the depths or be a constant")
+        }
+      }
+    }else{
+      sample1un=sample1
+      uncer=rep(0,length(depths))
+      suppo=abs(rnorm(length(supported),supported,rep(3,length(supported))) )
+    }
+  }
+
+
+  plot(depths,sample1un,pch=16,cex=.6,main="Simulated data",xlab="Depth (cm)",ylab="Bq/kg",ylim=c(0,(max(sample1un)+max(uncer))) )
+  segments(depths, sample1un+uncer, x1 = depths, y1 = sample1un-uncer)
+  points(depths,suppo,pch=16,col="red",cex=.6)
+  segments(depths, suppo+3, x1 = depths, y1 = suppo-3,col="red")
+
+
+  DATOS=matrix(c(depths,(density/10),sample1un,uncer,rep(thick,length(depths)),suppo,rep(3,length(suppo)) ),nrow=length(depths),byrow = FALSE)
+  return(DATOS)
+
+}
+
+
+#
+#
+# t <- function (x)(  (x^2)/3 + x/2) #  8*x-20*sin(x/(1.*pi)) )#
+# supp<- function(x)(10+abs(20*sin(x/pi))+.9*x/10)
+# curve(supp,0,30)
+# ro=function(x)(.9+cos(x/5.80*pi)+.085*x)
+# curve(ro,0,30)
+# Data_sim(t,seq(2,30,1),"supp",thick=2,Fi=50,rho="ro")
+#
+

R/Plots.R

@@ -0,0 +1,209 @@
+#' @export
+
+fullchronology= function(folder,Data,resolution=200,supp_type=1,
+                         memory_shape=4., memory_mean=.7,
+                         acc_shape=1.5,acc_mean=20,
+                         fi_mean=50,fi_acc=2,
+                         As_mean=20,As_acc=2){
+
+
+  par(mfrow=c(1,1))
+  Ages=read.table(paste(folder,"Results/dates.csv",sep=""),sep=" ")
+  intervals=read.table(paste(folder,"Results/intervals.csv",sep=""),sep=",")
+  Depths=as.numeric(read.table(paste(folder,"Results/depths.csv",sep=""),sep=",") )
+  Output=read.table(paste(folder,"Results/Results_output.csv",sep=""),sep=",")
+  Lead=read.table(paste(folder,Data,sep=""),sep=",")
+  Plotval=read.table(paste(folder,"Results/Graphs.csv",sep=""),sep=",")
+  Slopes=read.table(paste(folder,"Results/Slopes.csv",sep=""),sep=",")
+  num_var=length(Output[0,])
+  if(length(Lead[1,])==5){
+    Ran=1  } else if(length(Lead[1,])==7){
+      if(supp_type==1){Ran=1}else{
+      Ran=length(Lead[,1])}}
+
+  maxA=max(Ages[,length(Ages[1,])])+.10
+  ageSeq=seq(from=0,to=maxA,maxA/resolution)
+  deptsSeq=seq(from=0,to=Depths[length(Depths)],Depths[length(Depths)]/resolution)
+  deptsSeq=deptsSeq
+  diffSep=(deptsSeq[2]-deptsSeq[1])/2
+  TotSeq=length(Ages[,1])
+  memory_shape2=(memory_shape*(1-memory_mean) )/memory_mean
+
+
+  layout(matrix(c(1,1,2,2,3,3,4,4,
+                  1,1,2,2,3,3,4,4,
+                  5,5,5,5,5,5,5,5,
+                  5,5,5,5,5,5,5,5,
+                  5,5,5,5,5,5,5,5,
+                  5,5,5,5,5,5,5,5,
+                  5,5,5,5,5,5,5,5), 7, 8, byrow = TRUE))
+
+  d <- density(as.numeric(Output[-1,(Ran+2)]))
+  plot(d, xlab="",main="Memory",ylab = "",xlim=c(0,1),xaxs="i",yaxs="i")
+  polygon(d, col=gray(.6))
+  lines(seq(0,1,.01),dbeta(x = seq(0,1,.01),shape1 = memory_shape,shape2 = memory_shape2),col="green")
+
+
+  d <- density(as.numeric(unlist(Output[-1,-c(c(1:(Ran+2)),num_var)])))
+  plot(d,xlab="",main="Acc",ylab = "",xlim=c(0,80),xaxs="i",yaxs="i")
+  polygon(d, col=gray(.6))
+  lines(seq(0,100,.5),dgamma(seq(0,100,.5),shape=acc_shape,scale=acc_mean/acc_shape),col="green")
+
+  if(Ran==1){
+  d <- density(as.numeric(Output[-1,2]))
+  plot(d,xlab="",main="Supported Act",ylab = "",xaxs="i",yaxs="i")
+  polygon(d, col=gray(.6))
+  lines(seq(0,100,.05),dgamma(seq(0,100,.05),shape=As_acc,scale=As_mean/As_acc),col="green")
+  }else {
+    min1=min(as.numeric(unlist(Output[-1,2:(Ran+1)])))+.1
+    max1=max(as.numeric(unlist(Output[-1,2:(Ran+1)])))+.1
+    plot(0,0,xlim=c(Lead[1,1],Lead[Ran,1]),ylim=c(min1,max1),xlab="Depth (cm)",ylab="",main="Supported 210Pb",xaxs="i",yaxs="i")
+    for (k in 1:Ran) {
+      points(rep(Lead[k,1],length(Output[-1,1+k])), Output[-1,1+k],pch=19,col=rgb(0,0,0,.03) )
+      points(Lead[k,1],Lead[k,6],col="red",pch=18)
+
+    }
+  }
+
+  d <- density(as.numeric(Output[-1,1]))
+  plot(d,xlab="",main="Supply of 210Pb",ylab = "",xaxs="i",yaxs="i")
+  polygon(d, col=gray(.6))
+  lines(seq(0,350,.05),dgamma(seq(0,350,.05),shape=fi_acc,scale=fi_mean/fi_acc),col="green")
+
+  plot(-1,-1, xlim=c(0,Depths[length(Depths)]),ylim = c(0, maxA),xlab = "Depth (cm)",ylab="Age (years)" ,
+       xaxs="i",yaxs="i",main= gsub("\\..*","",Data))
+  for (i2 in 1:(resolution-1)){
+    for (i in 1:(resolution-1) ){
+      rect(deptsSeq[i2], ageSeq[i], deptsSeq[i2+1], ageSeq[i+1], density = NA, border = NA,
+           col = gray((1-Plotval[i2,i])))
+    }
+  }
+  lines(Depths,c(0,intervals[,2]),type="l", lty=2, lwd=1,col="red")
+  lines(Depths,(c(0,intervals[,4])),type="l", lty=2, lwd=1,col="red")
+  lines(Depths,(c(0,intervals[,3])),type="l", lty=2, lwd=1,col="red")
+  for (i in 1:length(Lead[,1])){
+    rug(Lead[i,1], lwd = Lead[i,5],col = "blue")
+  }
+
+  par(mfrow=c(1,1))
+}
+
+
+#' @export
+chronologylinesP= function(folder,Data,...){
+  Ages=read.table(paste(folder,"Results/dates.csv",sep=""),sep=" ")
+  intervals=read.table(paste(folder,"Results/intervals.csv",sep=""),sep=",")
+  Depths=as.numeric(read.table(paste(folder,"Results/depths.csv",sep=""),sep=",") )
+  Output=read.table(paste(folder,"Results/Results_output.csv",sep=""),sep=",")
+  Lead=read.table(paste(folder,Data,sep=""),sep=",")
+  Plotval=read.table(paste(folder,"Results/Graphs.csv",sep=""),sep=",")
+  Slopes=read.table(paste(folder,"Results/Slopes.csv",sep=""),sep=",")
+  num_var=length(Output[0,])
+  iterations=length(Ages[,1])
+
+
+  plot(Depths,c(0,Ages[2,]),type="l",col=rgb(0,0,0,.01), lwd=2,ylim = c(0,max(Ages[,length(Ages[1,])])),
+       xlab = "Depth (cm)",ylab="Age (years)",...)
+  for (i in 1:(iterations-1)){
+    lines(Depths,c(0,Ages[i,]),type="l",col=rgb(0,0,0,.01), lwd=2)
+  }
+  lines(Depths,c(0,intervals[,2]),type="l", lty=2, lwd=1)
+  lines(Depths,(c(0,intervals[,4])),type="l", lty=2, lwd=1)
+  lines(Depths,(c(0,intervals[,3])),type="l", lty=2, lwd=1)
+  for (i in 1:length(Lead[,1])){
+    rug(Lead[i,1], lwd = Lead[i,5],col = "blue")
+  }
+
+
+}
+
+
+#' @export
+chronologylines= function(folder,Data,...){
+  Ages=read.table(paste(folder,"Results/dates.csv",sep=""),sep=" ")
+  intervals=read.table(paste(folder,"Results/intervals.csv",sep=""),sep=",")
+  Depths=as.numeric(read.table(paste(folder,"Results/depths.csv",sep=""),sep=",") )
+  Output=read.table(paste(folder,"Results/Results_output.csv",sep=""),sep=",")
+  Lead=read.table(paste(folder,Data,sep=""),sep=",")
+  Plotval=read.table(paste(folder,"Results/Graphs.csv",sep=""),sep=",")
+  Slopes=read.table(paste(folder,"Results/Slopes.csv",sep=""),sep=",")
+  num_var=length(Output[0,])
+  iterations=length(Ages[,1])
+
+
+plot(Depths,c(0,Ages[2,]),type="l",col=rgb(0,0,0,.01), lwd=2,ylim = c(0,max(Ages[,length(Ages[1,])])),
+     xlab = "Depth (cm)",ylab="Age (years)",main= gsub("\\..*","",Data),...)
+for (i in 1:(iterations-1)){
+  lines(Depths,c(0,Ages[i,]),type="l",col=rgb(0,0,0,.01), lwd=2)
+}
+lines(Depths,c(0,intervals[,2]),type="l", lty=2, lwd=1)
+lines(Depths,(c(0,intervals[,4])),type="l", lty=2, lwd=1)
+lines(Depths,(c(0,intervals[,3])),type="l", lty=2, lwd=1)
+for (i in 1:length(Lead[,1])){
+  rug(Lead[i,1], lwd = Lead[i,5],col = "blue")
+}
+
+
+}
+
+
+
+#' @export
+slopes= function(folder,Data){
+    Depths=as.numeric(read.table(paste(folder,"Results/depths.csv",sep=""),sep=",") )
+    Slopes=read.table(paste(folder,"Results/Slopes.csv",sep=""),sep=",")
+    iterations=length(Slopes[,1])
+  maxS=max(Slopes)+.10
+  plot(-10,-10, xlim=c(Depths[length(Depths)]),ylim = c(0, maxS),xlab = "Depth (cm)",ylab="Slopes (Accumulations)" )
+  for (i in 1:(iterations-1)){
+    lines(Depths,as.numeric(c(Slopes[i,])),type="l",col=rgb(0,0,0,.01), lwd=2)
+  }
+  for (i in 1:length(Lead[,1])){
+    rug(Lead[i,1], lwd = Lead[i,5],col = "blue")
+  }
+}
+
+#' @export
+ageof=function(x,interval=.95,folder){
+  if(x!=0){
+    Ages=read.table(paste(folder,"Results/dates.csv",sep=""),sep=" ")
+    Depths=as.numeric(read.table(paste(folder,"Results/depths.csv",sep=""),sep=",") )
+    Slopes=read.table(paste(folder,"Results/Slopes.csv",sep=""),sep=",")
+    depfix=which(Depths<x)
+    depfix=depfix[length(depfix)]
+    m2=Slopes[,depfix]
+    sumages=c()
+    if(depfix!=1){
+      for (i in 1:length(Ages[,1])){
+        sumages=c(sumages, Ages[i,(depfix-1)] + Slopes[i,depfix]* (x-Depths[depfix]) )
+
+      }
+    }else{sumages= Slopes[,depfix]* (x)}
+    n=length(Ages[,1])
+    mean1=mean(sumages)
+    inter=(1-interval)/2
+    lim1=sort(sumages)[as.integer(n*inter)]
+    lim2=sort(sumages)[as.integer(n*(inter+interval))]
+
+    d<- density(sumages)
+    plot(d, main=paste("Age at ",x, "cm"),ylab = "",xlab = "Age",xaxs="i",yaxs="i")
+    polygon(d, col=gray(.6))
+    points(x=lim1,y=(max(d$y)*.015),col="red",pch="(")
+    points(x=lim2,y=(max(d$y)*.015),col="red",pch=")")
+    points(x=mean1,y=(max(d$y)*.015),col="red",pch=16)
+    print(paste("the age of depth",x,"is between") )
+    print(c(lim1,lim2))
+    print(paste("with a ", interval, "%"," condifence interval and a mean of:",sep = ""))
+    print(mean1)
+
+
+  }else{
+    print("For depth 0, the age is equal to the collection date.")
+  }
+
+
+}
+
+
+
+