-
Notifications
You must be signed in to change notification settings - Fork 5
/
fit.thinplate_2d.R
60 lines (48 loc) · 2.08 KB
/
fit.thinplate_2d.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#' Function to fit a local radial basis function (CReSS) as a two dimensional smooth
#'
#'
#' @export
#'
"fit.thinPlate_2d" <- function(fitnessMeasure, dists,aR,radii,baseModel,radiusIndices,models, currentFit, interactionTerm, data, initDisp, cv.opts, basis='gaussian', printout) {
if (isS4(baseModel)){
attributes(baseModel@misc$formula)$.Environment<-environment()
baseModel@splineParams[[1]]$knotPos<-aR
baseModel@splineParams[[1]]$radiusIndices<-radiusIndices
} else {
attributes(baseModel$formula)$.Environment<-environment()
baseModel$splineParams[[1]]$knotPos<-aR
#baseModel$splineParams[[1]]$knotPos<-baseModel$splineParams[[1]]$mapInd[aR]
baseModel$splineParams[[1]]$radiusIndices<-radiusIndices
}
#print("ooooooooooooooooooooooooooooooooooooooo")
#print("Fitting Model...")
#print("ooooooooooooooooooooooooooooooooooooooo")
#aR<- aR
#radiusIndices<- radiusIndices
#cat("R Indices: ", radiusIndices, "\n")
# print(aR)
if(basis=='gaussian'){
bspl<- "LRF.g(radiusIndices, dists, radii, aR)"
}
if(basis=='exponential'){
bspl<- "LRF.e(radiusIndices, dists, radii, aR)"
}
if(is.null(interactionTerm)){
test<-paste("update(baseModel, . ~ . + ",bspl, ")",sep="")
currentModel<-eval(parse(text=test))
}else{
test<-paste("update(baseModel, . ~ . + ",bspl, "*",interactionTerm, ")", sep="")
currentModel<-eval(parse(text=test))
}
tempFit <- get.measure_2d(fitnessMeasure, currentFit, currentModel,data, dists,aR,radii,radiusIndices, initDisp, cv.opts, printout)$fitStat
# if(tempFit <= (currentFit+10)){
# models[[length(models)+1]] = list(aR,radiusIndices, radii, tempFit)
# }
models<-NULL
#modelinprogress<<-currentModel
#print("ooooooooooooooooooooooooooooooooooooooo")
#print("Model fitted...")
#print(paste('disp= ', summary(currentModel)$dispersion, ', num knots: ', length(aR), ', fitstat: ',tempFit,sep=''))
#print("ooooooooooooooooooooooooooooooooooooooo")
return(list(currentModel=currentModel,models=models, fitStat=tempFit))
}