BT-Settl preparation for IPAC (T) based on class of features 1-Fs/Fc.
Introduction
After the pre-processing analysis carried out for M stars within the IPAC dataset, we will proceed with a similar methodology to the one ran for IRTF.
## Loading required package: MASS
## Loading required package: Cubist
## Loading required package: lattice
## Loading required package: caret
## Loading required package: ggplot2
## Loading required package: xtable
## Loading required package: Peaks
## Loading required package: magic
## Loading required package: abind
## Loading required package: JADE
## Loading required package: doMC
## Loading required package: foreach
## Loading required package: iterators
## Loading required package: parallel
## Loading required package: segmented
## Loading required package: fftw
## Loading required package: FITSio
## Loading required package: stringr
## Loading required package: ztable
## Welcome to package ztable ver 0.1.5
## Loading required package: signal
##
## Attaching package: 'signal'
##
## The following objects are masked from 'package:stats':
##
## filter, poly
##
## Loading required package: e1071
## Loading required package: class
## Loading required package: zoo
##
## Attaching package: 'zoo'
##
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
##
## Loading required package: plyr
## Loading required package: compare
##
## Attaching package: 'compare'
##
## The following object is masked from 'package:base':
##
## isTRUE
##
## Loading required package: mgcv
## Loading required package: nlme
## This is mgcv 1.8-6. For overview type 'help("mgcv-package")'.
##
## Attaching package: 'mgcv'
##
## The following object is masked from 'package:magic':
##
## magic
##
## Loading required package: elasticnet
## Loading required package: lars
## Loaded lars 1.2
##
## Loading required package: pbapply
## Loading required package: nnet
## Loading required package: kernlab
## Loading required package: pls
##
## Attaching package: 'pls'
##
## The following object is masked from 'package:caret':
##
## R2
##
## The following object is masked from 'package:stats':
##
## loadings
##
## Loading required package: multicore
##
## Attaching package: 'multicore'
##
## The following objects are masked from 'package:parallel':
##
## mclapply, mcparallel, pvec
##
## The following object is masked from 'package:lattice':
##
## parallel
##
## Loading required package: devtools
## Loading required package: caretEnsemble
## Loading required package: mlbench
## Loading required package: mclust
## Package 'mclust' version 5.0.1
## Type 'citation("mclust")' for citing this R package in publications.
##
## Attaching package: 'mclust'
##
## The following object is masked from 'package:mgcv':
##
## mvn
##
## Loading required package: analogue
## Loading required package: vegan
## Loading required package: permute
##
## Attaching package: 'permute'
##
## The following object is masked from 'package:devtools':
##
## check
##
## The following object is masked from 'package:kernlab':
##
## how
##
## This is vegan 2.2-1
##
## Attaching package: 'vegan'
##
## The following object is masked from 'package:pls':
##
## scores
##
## The following object is masked from 'package:caret':
##
## tolerance
##
## This is analogue 0.16-0
##
## Attaching package: 'analogue'
##
## The following objects are masked from 'package:pls':
##
## RMSEP, crossval, pcr
##
## The following object is masked from 'package:compare':
##
## compare
##
## The following object is masked from 'package:plyr':
##
## join
##
## Loading required package: cluster
## [[1]]
## [1] TRUE
##
## [[2]]
## [1] TRUE
##
## [[3]]
## [1] TRUE
##
## [[4]]
## [1] TRUE
##
## [[5]]
## [1] TRUE
##
## [[6]]
## [1] TRUE
##
## [[7]]
## [1] TRUE
##
## [[8]]
## [1] TRUE
##
## [[9]]
## [1] TRUE
##
## [[10]]
## [1] TRUE
##
## [[11]]
## [1] TRUE
##
## [[12]]
## [1] TRUE
##
## [[13]]
## [1] TRUE
##
## [[14]]
## [1] TRUE
##
## [[15]]
## [1] TRUE
##
## [[16]]
## [1] TRUE
##
## [[17]]
## [1] TRUE
##
## [[18]]
## [1] TRUE
##
## [[19]]
## [1] TRUE
##
## [[20]]
## [1] TRUE
##
## [[21]]
## [1] TRUE
##
## [[22]]
## [1] TRUE
##
## [[23]]
## [1] TRUE
##
## [[24]]
## [1] TRUE
##
## [[25]]
## [1] TRUE
##
## [[26]]
## [1] TRUE
##
## [[27]]
## [1] TRUE
##
## [[28]]
## [1] TRUE
##
## [[29]]
## [1] TRUE
##
## [[30]]
## [1] TRUE
##
## [[31]]
## [1] TRUE
##
## [[32]]
## [1] TRUE
##
## [[33]]
## [1] TRUE
##
## [[34]]
## [1] TRUE
##
## [[35]]
## [1] TRUE
##
## [[36]]
## [1] TRUE
##
## [[37]]
## [1] TRUE
Target T will be between 2000 K and 4200 K
Prediction for Temperature
We will use the I-Band according to GA models noised by SNR=50 and SNR=10 as well as the features proposed by Cestti et al. also noised by SNR=50 and SNR=10.
Recovering already built models
# Procesado del genético de Tª
range_elem<-function(z,l,bi){
look=function(x,l,y){
j<-which(y>= as.numeric(x[l]))[1];
k<-(which(y>as.numeric(x[(l+1)]))[1])-1;
return(c(j,k))
}
y=z$data[[1]][,1]
lm<-t(apply(bi,1,look,l,y))
rownames(lm)<-bi[,1]
colnames(lm)<-c("From","To")
return(lm)
}
int_spec<-function(x,idx,norm=0){
y<-x$data[[1]][eval(parse(text=idx)),]
xz<-diff(as.numeric(y[,1]),1)
yz<-as.numeric(y[,2])
if ( norm > 0 ){
z<-sum(xz)
} else {
z<-sum(xz*rollmean(yz,2))
}
return(z)
}
#
feature_extr<-function(sn,bp){
sig<-sn[1]
noi<-sn[2]
nof<-sn[3] # Allowing to consider 2 bands as the paper
Fcont<-unlist(lapply(bp,int_spec,noi,0))/
unlist(lapply(bp,int_spec,noi,1))
if (! is.na(nof)){ # In case of two bands => average
Fcont2=unlist(lapply(bp,int_spec,nof,0))/
unlist(lapply(bp,int_spec,nof,1))
Fcont = (Fcont + Fcont2) /2.
}
fea<-unlist(lapply(bp,int_spec,sig,1))-
unlist(lapply(bp,int_spec,sig,0))/Fcont
return(fea)
}
#
sacar=function(x,pos=2) {
if ( length(x) < pos ) {
return(NA);
}
return(x[pos]);
}
#
######################
#
# Loading the object BT-Settl
if (file.exists("~/git/M_prep_IPAC/Models_BT_settl_IPAC33_v2.RData")) {
load( file="~/git/M_prep_IPAC/Models_BT_settl_IPAC33_v2.RData")
} else {
cat("ERROR: ~/git/M_prep_IPAC/Models_BT_settl_IPAC33_v2.RData NOT FOUND !!!")
exit(2)
}
#
Now, we will read the M stars from the IPAC dataset in fits format.
Reading the IPAC dataset
#
scaling=function(x,xlim=NULL,xp=NULL) {
y=x
if (! is.null(y$factors)) {
y$data[,2]=y$data[,2]*y$factors
}
if ( is.null(xlim)){
ixp=rep(TRUE,nrow(y$data))
xlim=range(y$data[,1])
}
if ( ! is.null(xp)) {
xp = xp[(xp > xlim[1]) & (xp < xlim[2])]
xli0 = y$data[rev(which(y$data[,1] < xlim[1]))[1],1] - 0.001 # Rounding ...
xli1 = y$data[which(y$data[,1]>xlim[2])[1],1] + 0.001 # Rounding ...
ixp = (y$data[,1] >= xli0 & y$data[,1] <= xli1)
} else {
ixp = (y$data[,1]>= xlim[1] & y$data[,1] <= xlim[2])
}
z=y$data[ixp,]
jp=rep(0,nrow(z)+1)
# jp[is.finite(z[,2])]=1 # LSB has requested to interpolate the holes 30/7/2014
jp[nrow(z)+1]=1
zp=diff(jp)
ip=which(zp != 0)
ini=1
j=0
y$data=list()
for (i in ip) {
if (zp[i] != -1 ) {
j=j+1
y$data[[j]]=z[ini:i,]
} else {
ini=i+1
}
}
y$nch=j
y$factors=rep(NA,y$nch)
for (i in 1:y$nch) {
if ( is.na(y$data[[i]][nrow(y$data[[i]]),2])) {
kk = rev(which(! is.na(y$data[[i]][,2])))[1]
y$data[[i]][nrow(y$data[[i]]),2]=y$data[[i]][kk,2]
}
if ( is.na(y$data[[i]][1,2])) {
kk = which(! is.na(y$data[[i]][,2]))[1]
y$data[[i]][1,2]=y$data[[i]][kk,2]
}
if ( ! is.null(xp)) {
yp=(approx(y$data[[i]][,1],y$data[[i]][,2],xout=xp))$y
y$data[[i]]=data.frame(x=xp,y=yp)
} else {
yp=(approx(y$data[[i]][,1],y$data[[i]][,2],xout=y$data[[i]][,1]))$y
y$data[[i]]=data.frame(x=y$data[[i]][,1],y=yp)
}
y$factors[i]=sum(rollmean(y$data[[i]][,2],2)*diff(y$data[[i]][,1]))
}
ss=sum(y$factors)
for (i in 1:y$nch) {
y$data[[i]][,2]=y$data[[i]][,2]/ss
}
return(y)
}
#
pinta=function(x,cortes=NULL,borders=FALSE, boundary=FALSE, close=FALSE,
xlim=NULL, ylim=NULL,plt=NULL,width=12,height=8,...){
if (is.null(x$nch)){
nch=1
} else {
nch=x$nch
}
if ( length(xlim) < 2 ) {
xmin=min(x$data[[1]][,1])
xmax=max(x$data[[nch]][,1])
xl=c(xmin,xmax)
} else {
xl=xlim
}
if ( length(ylim) < 2 ) {
ymin=min(x$data[[1]][,2])
ymax=max(x$data[[1]][,2])
if ( nch > 1 ) {
for (i in 2:nch){
ymin=min(ymin,x$data[[i]][,2])
ymax=max(ymax,x$data[[i]][,2])
}
}
yl=c(ymin,ymax)
} else {
yl = ylim
}
if (! is.null(plt) ) {
pdf(file=plt,width,height)
}
if ( boundary) {
plot(0,xlim=xl,ylim=yl,...)
}
for (i in 1:nch) {
lines(x$data[[i]],...)
}
for (i in cortes) {
abline(v=i[1],col=2)
abline(v=i[2],col=2)
}
if ( borders) {
for (i in nch) {
abline(v=min(x$data[[i]][,1]),col=3)
abline(v=max(x$data[[i]][,1]),col=3)
}
}
if ( close ) {
dev.off()
}
}
setwd('~/git/M_prep_IPAC')
#
A number of 595 have been readed from the individual spectra dataset.
Preparation of potential features for T
#
if ( file.exists("~/git/M_prep_IPAC/Features_BT-settl_v30.RData")) {
load("~/git/M_prep_IPAC/Features_BT-settl_v30.RData")
} else {
bpy=do.call(rbind,lapply(bp_clean,function(x){return(x$data[[1]][,2])}))
bpy10=do.call(rbind,lapply(bp_10,function(x){return(x$data[[1]][,2])}))
bpy50=do.call(rbind,lapply(bp_50,function(x){return(x$data[[1]][,2])}))
YPT=do.call(rbind,lapply(bp_clean,function(x){return(x$stellarp[1])}))
YPG=do.call(rbind,lapply(bp_clean,function(x){return(x$stellarp[2])}))
YPM=do.call(rbind,lapply(bp_clean,function(x){return(x$stellarp[3])}))
vf=as.character(bp_clean[[1]]$data[[1]][,1])
colnames(bpy)=vf
colnames(bpy10)=vf
colnames(bpy50)=vf
bfy=do.call(rbind,lapply(bf_clean,function(x){return(x$data[[1]][,2])}))
colnames(bfy)=vf
#
area=function(v,nv) {
y=data.frame(x=as.numeric(nv),y=v)
xz = diff(as.numeric(y[,1]),1)
yz = as.numeric(y[,2])
z = sum(xz*rollmean(yz,2))
return(z)
}
#
######################################################################
#
siz=30
bandas=list()
bfndas=list()
bandas10=list()
bandas50=list()
for (i in seq(1,siz,15)) { # LSB asked for steps of 5 pixels
idx=paste("s=",i,sep="")
lvf=vf
if ( i > 1) {
lvf=vf[-c(1:(i-1))]
}
ff=as.factor(sort(rep(1:ceiling(length(lvf)/siz),siz))[1:length(lvf)])
st=split(lvf,ff)
for (j in 1:length(st)) {
if (length(st[[j]]) < siz)
st[[j]] = NULL
}
lst = length(st)
bandas[[idx]]=list(names= as.vector(unlist(lapply(st,
function(x){return(paste(range(x),collapse="-"))}))),
mat=matrix(NA,ncol=lst,nrow=nrow(bpy)))
bfndas[[idx]]=list(names= as.vector(unlist(lapply(st,
function(x){return(paste(range(x),collapse="-"))}))),
mat=matrix(NA,ncol=lst,nrow=nrow(bfy)))
bandas10[[idx]]=list(names= as.vector(unlist(lapply(st,
function(x){return(paste(range(x),collapse="-"))}))),
mat=matrix(NA,ncol=lst,nrow=nrow(bpy10)))
bandas50[[idx]]=list(names= as.vector(unlist(lapply(st,
function(x){return(paste(range(x),collapse="-"))}))),
mat=matrix(NA,ncol=lst,nrow=nrow(bpy50)))
for (j in 1:lst) {
bandas[[idx]]$mat[,j]=apply(bpy[,st[[j]]],1,area,st[[j]])
bandas10[[idx]]$mat[,j]=apply(bpy10[,st[[j]]],1,area,st[[j]])
bandas50[[idx]]$mat[,j]=apply(bpy50[,st[[j]]],1,area,st[[j]])
}
for (j in 1:lst) {
bfndas[[idx]]$mat[,j]=apply(bfy[,st[[j]]],1,area,st[[j]])
}
}
#
save(bpy,bfy,bandas,bfndas,vf, bpy10,bpy50,st,
bandas10,bandas50,YPT,YPG,YPM,
file="~/git/M_prep_IPAC/Features_BT-settl_v30.RData")
}
#
Otras features de genéticos que evolucionan construyendo las features (Población 8000 individuos, 1000 evoluciones)
#
#
require(GA)
GA_r Loading required package: GA
GA_r Package 'GA' version 2.2
GA_r Type 'citation("GA")' for citing this R package in publications.
require(doParallel)
GA_r Loading required package: doParallel
require(plyr)
#
#
area=function(v,nv) {
y=data.frame(x=as.numeric(nv),y=v)
xz = diff(as.numeric(y[,1]),1)
yz = as.numeric(y[,2])
z = sum(xz*rollmean(yz,2))
return(z)
}
fitness = function(string) {
par = decode(string)
if ( 0 %in% unlist(apply(par,1,sd)) ) {
return (-1000*NBLK*NBITS)
}
if ( length(unique(sort(par[,1]))) < nrow(par) |
length(unique(sort(par[,2]))) < nrow(par) ) {
return (-1000*NBLK*NBITS)
}
X = as.data.frame(matrix(0,nrow=nrow(bandas[[cnjts[1]]]$mat),ncol=(nrow(par))))
colnames(X)[1] ="Intercept"
for ( i in 1:nrow(par)) {
bi=par[i,1] %% NBLK
bj=par[i,2] %% NBLK
if ( bi==0 ) bi=NBLK
if ( bj==0 ) bj=NBLK
idxi=(par[i,1] %/% NBLK) +1
idxj=(par[i,2] %/% NBLK) +1
wini=-eval(parse(text=bandas[[cnjts[bi]]]$names[idxi]))
winj=-eval(parse(text=bandas[[cnjts[bj]]]$names[idxj]))
X[,(i)]=wini-bandas[[cnjts[bi]]]$mat[,idxi] / (bandas[[cnjts[bj]]]$mat[,idxj]/winj)
colnames(X)[(i)] = paste(par[i,],collapse="-")
}
mod = lm.fit(as.matrix(X), y)
class(mod) = "lm"
return( -BIC(mod))
}
fitness10 = function(string) {
par = decode(string)
if ( 0 %in% unlist(apply(par,1,sd)) ) {
return (-1000*NBLK*NBITS)
}
if ( length(unique(sort(par[,1]))) < nrow(par) |
length(unique(sort(par[,2]))) < nrow(par) ) {
return (-1000*NBLK*NBITS)
}
X = as.data.frame(matrix(0,nrow=nrow(bandas10[[cnjts[1]]]$mat),ncol=(nrow(par))))
colnames(X)[1] ="Intercept"
for ( i in 1:nrow(par)) {
bi=par[i,1] %% NBLK
bj=par[i,2] %% NBLK
if ( bi==0 ) bi=NBLK
if ( bj==0 ) bj=NBLK
idxi=(par[i,1] %/% NBLK) +1
idxj=(par[i,2] %/% NBLK) +1
wini=-eval(parse(text=bandas10[[cnjts[bi]]]$names[idxi]))
winj=-eval(parse(text=bandas10[[cnjts[bj]]]$names[idxj]))
X[,(i)]=wini-bandas10[[cnjts[bi]]]$mat[,idxi] / (bandas10[[cnjts[bj]]]$mat[,idxj]/winj)
colnames(X)[(i)] = paste(par[i,],collapse="-")
}
mod = lm.fit(as.matrix(X), y)
class(mod) = "lm"
return( -BIC(mod))
}
fitness50 = function(string) {
par = decode(string)
if ( 0 %in% unlist(apply(par,1,sd)) ) {
return (-1000*NBLK*NBITS)
}
if ( length(unique(sort(par[,1]))) < nrow(par) |
length(unique(sort(par[,2]))) < nrow(par) ) {
return (-1000*NBLK*NBITS)
}
X = as.data.frame(matrix(0,nrow=nrow(bandas50[[cnjts[1]]]$mat),ncol=(nrow(par))))
colnames(X)[1] ="Intercept"
for ( i in 1:nrow(par)) {
bi=par[i,1] %% NBLK
bj=par[i,2] %% NBLK
if ( bi==0 ) bi=NBLK
if ( bj==0 ) bj=NBLK
idxi=(par[i,1] %/% NBLK) +1
idxj=(par[i,2] %/% NBLK) +1
wini=-eval(parse(text=bandas50[[cnjts[bi]]]$names[idxi]))
winj=-eval(parse(text=bandas50[[cnjts[bj]]]$names[idxj]))
X[,(i)]=wini-bandas50[[cnjts[bi]]]$mat[,idxi] / (bandas50[[cnjts[bj]]]$mat[,idxj]/winj)
colnames(X)[(i)] = paste(par[i,],collapse="-")
}
mod = lm.fit(as.matrix(X), y)
class(mod) = "lm"
return( -BIC(mod))
}
#
SelFeatures=function(x,sets=cnjts,dat=bandas,fitness.=fitness,pp=FALSE){
pos=which.max(as.numeric(unlist(ldply(x@bestSol,fitness.,.parallel=pp))))
sol=x@bestSol[[pos]]
par = decode(sol)
res = list()
for ( i in 1:nrow(par)) {
bi=par[i,1] %% NBLK
bj=par[i,2] %% NBLK
if ( bi==0 ) bi=NBLK
if ( bj==0 ) bj=NBLK
idxi=(par[i,1] %/% NBLK) +1
idxj=(par[i,2] %/% NBLK) +1
res[[i]]=list(set1=sets[bi],pos1=idxi,set2=sets[bj],pos2=idxj,par=par[i,],
x=dat[[sets[bi]]]$mat[,idxi],
x0=dat[[sets[bj]]]$mat[,idxj],
name=dat[[sets[bi]]]$nam[idxi],
name0=dat[[sets[bj]]]$nam[idxj])
}
return(res)
}
#
FeaturesExt=function(x,num,sets=cnjts,dat=bandas,fitness.=fitness,pp=FALSE){
vals=as.numeric(unlist(ldply(x@bestSol,fitness.,.parallel=pp)))
fiveval=summary(vals)
uv = unique(sort(vals,decreasing=TRUE))
uvals = uv[1:min(num,length(uv))]
isols=rep(0,length(uvals))
for (i in 1:length(uvals)) {
isols[i] = which(vals==uvals[i])[1]
}
res=NA
if ( length(isols) > 0) {
par=decode(x@bestSol[[1]])
res = as.data.frame(matrix(NA,nrow=nrow(par),ncol=(2*length(isols))))
nam = as.data.frame(matrix(NA,nrow=nrow(par),ncol=(2*length(isols))))
for (j in 1:length(isols)) {
sol=x@bestSol[[isols[j]]]
par = decode(sol)
res[,(2*(j-1)+1):(2*j)]=par
names(res)[(2*(j-1)+1):(2*j)] = c(paste(j,":s",sep=""),paste(j,":c",sep=""))
for ( i in 1:nrow(par)) {
bi=par[i,1] %% NBLK
bj=par[i,2] %% NBLK
if ( bi==0 ) bi=NBLK
if ( bj==0 ) bj=NBLK
idxi=(par[i,1] %/% NBLK) +1
idxj=(par[i,2] %/% NBLK) +1
nam[i,(2*(j-1)+1)]=dat[[sets[bi]]]$nam[idxi]
nam[i,(2*j)]=dat[[sets[bj]]]$nam[idxj]
}
names(nam)[(2*(j-1)+1):(2*j)] = c(paste(j,":s",sep=""),paste(j,":c",sep=""))
}
}
return(list(par=res,name=nam,pos=isols,vpos=vals[isols],fivevals=fiveval))
}
#
decode = function(string) {
string = gray2binary(string)
n = binary2decimal(string[1:NBITS])
c = binary2decimal(string[(NBITS + 1):(2*NBITS)])
res = data.frame(p1=n%%MAXV,p2=c%%MAXV)
ini = 2 * NBITS
if ( NG > 1) {
for (i in 2:NG) {
n = binary2decimal(string[(ini+1):(ini+NBITS)])
c = binary2decimal(string[(ini + NBITS + 1):(ini + 2*NBITS)])
res = rbind(res, c(n%%MAXV,c%%MAXV))
ini = ini + 2*NBITS
}
}
return(res)
}
#
encode = function(res,nbits=NBITS) {
if (! is.data.frame(res)) {
warning(paste("ENCODE: First parameter is not data.frame:",str(res),sep=" "))
return(NULL)
}
cad=rep(0,ncol(res)*nrow(res)*nbits)
for (i in 1:nrow(res)) {
for (j in 1:ncol(res)) {
k=((i-1)*ncol(res) + (j - 1))*nbits+1
cad[k:(k+nbits-1)]=decimal2binary(res[i,j],nbits)
}
}
return(binary2gray(cad))
}
#
#
NITER=9
if ( file.exists("~/git/M_prep_IPAC/GA_NT11F2_IPAC_30.RData")) {
j=1
load("~/git/M_prep_IPAC/GA_NT11F2_IPAC_30.RData")
tt50=list()
GAs50 = list()
solsT50=list()
tt10=list()
GAs10 = list()
solsT10=list()
tt=list()
GAs=list()
solsT=list()
while(j <= NITER) {
cat(paste("Reading Iteration:",j,"<br>",sep=""))
flush.console()
load(paste("~/git/M_prep_IPAC/GA_T_",sprintf("%02d",j),"_50_f2_IPAC_30.RData",sep=""))
tt50[[j]]= t_tot
GAs50[[j]]=GA1
solsT50[[j]]=FeaturesExt(GA1,12)
load(paste("~/git/M_prep_IPAC/GA_T_",sprintf("%02d",j),"_10_f2_IPAC_30.RData",sep=""))
tt10[[j]]= t_tot
GAs10[[j]]=GA1
solsT10[[j]]=FeaturesExt(GA1,12)
load(paste("~/git/M_prep_IPAC/GA_T_",sprintf("%02d",j),"_f2_IPAC_30.RData",sep=""))
tt[[j]]= t_tot
GAs[[j]]=GA1
solsT[[j]]=FeaturesExt(GA1,12)
j = j+1
}
} else {
# cnjts = c("s=1","s=6","s=11","s=16","s=21","s=26")
cnjts = c("s=1","s=16")
bpy=do.call(rbind,lapply(bp_clean,function(x){return(x$data[[1]][,2])}))
YPT=do.call(rbind,lapply(bp_clean,function(x){return(x$stellarp[1])}))
YPG=do.call(rbind,lapply(bp_clean,function(x){return(x$stellarp[2])}))
YPM=do.call(rbind,lapply(bp_clean,function(x){return(x$stellarp[3])}))
vf=as.character(bp_clean[[1]]$data[[1]][,1])
colnames(bpy)=vf
bfy=do.call(rbind,lapply(bf_clean,function(x){return(x$data[[1]][,2])}))
colnames(bfy)=vf
#
siz=30
bandas=list()
bfndas=list()
bandas10=list()
bandas50=list()
for (i in seq(1,30,15)) { # LSB asked for steps of 5 pixels
idx=paste("s=",i,sep="")
lvf=vf
if ( i > 1) {
lvf=vf[-c(1:(i-1))]
}
ff=as.factor(sort(rep(1:ceiling(length(lvf)/siz),siz))[1:length(lvf)])
st=split(lvf,ff)
for (j in 1:length(st)) {
if (length(st[[j]]) < siz)
st[[j]] = NULL
}
lst = length(st)
bandas[[idx]]=list(names= as.vector(unlist(lapply(st,
function(x){return(paste(range(x),collapse="-"))}))),
mat=matrix(NA,ncol=lst,nrow=nrow(bpy)))
bfndas[[idx]]=list(names= as.vector(unlist(lapply(st,
function(x){return(paste(range(x),collapse="-"))}))),
mat=matrix(NA,ncol=lst,nrow=nrow(bfy)))
bandas10[[idx]]=list(names= as.vector(unlist(lapply(st,
function(x){return(paste(range(x),collapse="-"))}))),
mat=matrix(NA,ncol=lst,nrow=nrow(bpy10)))
bandas50[[idx]]=list(names= as.vector(unlist(lapply(st,
function(x){return(paste(range(x),collapse="-"))}))),
mat=matrix(NA,ncol=lst,nrow=nrow(bpy50)))
for (j in 1:lst) {
bandas[[idx]]$mat[,j]=apply(bpy[,st[[j]]],1,area,st[[j]])
bandas10[[idx]]$mat[,j]=apply(bpy10[,st[[j]]],1,area,st[[j]])
bandas50[[idx]]$mat[,j]=apply(bpy50[,st[[j]]],1,area,st[[j]])
}
for (j in 1:lst) {
bfndas[[idx]]$mat[,j]=apply(bfy[,st[[j]]],1,area,st[[j]])
}
}
#
#
y=as.numeric(unlist(YPT))
cnjts = c("s=1","s=16")
NG=10 # Número de Features
NBLK = length(cnjts)
NBITS= ceiling(log(NBLK*ncol(bandas[[cnjts[1]]]$mat),2))
MAXV = ncol(bandas[[cnjts[1]]]$mat)
#
save(bpy,bpy10,bpy50,bfy,YPT,YPG,YPM,vf,siz,bandas,bfndas, bandas10,bandas50,
cnjts,st,NG,NBLK,NBITS,MAXV,y,fitness,encode,decode,
SelFeatures,file="~/git/M_prep_IPAC/GA_NT11F2_IPAC_30.RData")
#
tt=list()
GAs = list()
solsT=list()
j=1
while(j <= NITER) {
cat(paste("Starting Iteration:",j,"<br>",sep=""))
flush.console()
tt[[j]]= system.time({GA1 = ga("binary", fitness = fitness, nBits = NBITS*2*NG,
monitor = FALSE,parallel=18,maxiter=1000,popSize=8000,keepBest=TRUE)})
t_tot=tt[[j]]
GAs[[j]]=GA1
save(fitness,encode,decode,SelFeatures,GA1,t_tot,j,
file=paste("~/git/M_prep_IPAC/GA_T_",sprintf("%02d",j),"_f2_IPAC_30.RData",sep=""))
solsT[[j]]=FeaturesExt(GA1,12)
j=j+1
}
#
print(xtable(ldply(solsT,function(x){return(x$vpos[1])})),type="html")
#
j=1
tt10=list()
GAs10 = list()
solsT10=list()
while(j <= NITER) {
cat(paste("Starting Iteration:",j,"<br>",sep=""))
flush.console()
tt10[[j]]= system.time({GA1 = ga("binary", fitness = fitness10, nBits = NBITS*2*NG,
monitor = FALSE,parallel=18,maxiter=1000,popSize=8000,keepBest=TRUE)})
t_tot=tt10[[j]]
GAs10[[j]]=GA1
save(fitness,encode,decode,SelFeatures,GA1,t_tot,j,
file=paste("~/git/M_prep_IPAC/GA_T_",sprintf("%02d",j),"_10_f2_IPAC_30.RData",sep=""))
solsT10[[j]]=FeaturesExt(GA1,12)
j=j+1
}
#
print(xtable(ldply(solsT10,function(x){return(x$vpos[1])})),type="html")
#
j=1
tt50=list()
GAs50 = list()
solsT50=list()
while(j <= NITER) {
cat(paste("Starting Iteration:",j,"<br>",sep=""))
flush.console()
tt50[[j]]= system.time({GA1 = ga("binary", fitness = fitness50, nBits = NBITS*2*NG,
monitor = FALSE,parallel=18,maxiter=1000,popSize=8000,keepBest=TRUE)})
t_tot=tt50[[j]]
GAs50[[j]]=GA1
save(fitness,encode,decode,SelFeatures,GA1,t_tot,j,
file=paste("~/git/M_prep_IPAC/GA_T_",sprintf("%02d",j),"_50_f2_IPAC_30.RData",sep=""))
solsT50[[j]]=FeaturesExt(GA1,12)
j=j+1
}
#
print(xtable(ldply(solsT50,function(x){return(x$vpos[1])})),type="html")
}
Reading Iteration:1
Reading Iteration:2
Reading Iteration:3
Reading Iteration:4
Reading Iteration:5
Reading Iteration:6
Reading Iteration:7
Reading Iteration:8
Reading Iteration:9
#
print(xtable(ldply(solsT,function(x){return(x$vpos[1])})),type="html")
|
|
V1
|
|
1
|
-11161.73
|
|
2
|
-11284.13
|
|
3
|
-11281.02
|
|
4
|
-11183.67
|
|
5
|
-11181.73
|
|
6
|
-11176.48
|
|
7
|
-11157.40
|
|
8
|
-11249.65
|
|
9
|
-11102.38
|
# ISOL holds the best solution !
isol = which.max(as.numeric(ldply(solsT,function(x){return(x$vpos[1])})[,1]))
isol
[1] 9
print(xtable(ldply(solsT10,function(x){return(x$vpos[1])})),type="html")
|
|
V1
|
|
1
|
-11310.37
|
|
2
|
-11353.04
|
|
3
|
-11361.24
|
|
4
|
-11394.16
|
|
5
|
-11346.44
|
|
6
|
-11327.27
|
|
7
|
-11384.69
|
|
8
|
-11262.78
|
|
9
|
-11285.79
|
# ISOL holds the best solution !
isol10 = which.max(as.numeric(ldply(solsT10,function(x){return(x$vpos[1])})[,1]))
isol10
[1] 8
#
print(xtable(ldply(solsT50,function(x){return(x$vpos[1])})),type="html")
|
|
V1
|
|
1
|
-11321.64
|
|
2
|
-11247.15
|
|
3
|
-11104.20
|
|
4
|
-11169.37
|
|
5
|
-11202.69
|
|
6
|
-11266.13
|
|
7
|
-11213.65
|
|
8
|
-11195.41
|
|
9
|
-11165.31
|
# ISOL holds the best solution !
isol50 = which.max(as.numeric(ldply(solsT50,function(x){return(x$vpos[1])})[,1]))
isol50
[1] 3
#
T foreseen by Luminosity models
# Loading IPAC prediction from Dr Sarro (objeto piac_temp)
load("~/git/M_prep_Noise/Teffs.RData")
ipac_temp[,1]=as.character(ipac_temp[,1])
ipac_temp[,2]=as.character(ipac_temp[,2])
ipac_temp[,3]=as.character(ipac_temp[,3])
ipac_temp[,4]=as.character(ipac_temp[,4])
ipac_temp[,5]=as.character(ipac_temp[,5])
ipac_temp[,6]=as.character(ipac_temp[,6])
ipac_temp[ipac_temp=='\xa0']=""
#
lnm<-unlist(lapply(bf_clean,function(x){return(x$name)}))
buscar=function(x,y,l=4){
for ( i in x[1:l]) {
if (nchar(i) > 1) {
stg = gsub('+','p',gsub(' ','_',i),fixed=TRUE)
stg = gsub('\\(','\\\\(',gsub('\\)','\\\\)',stg))
if (nchar(stg) > 0) {
pos=grep(stg,y)
if ( length(pos) == 1 ) {
return (pos)
}
}
}
}
return(-1)
}
#
idx=apply(ipac_temp[,-c(5,7)],1,buscar,lnm,6)
Teff_scs = rep(NA,max(idx))
Tipo_scs = rep(NA,max(idx))
idpos=data.frame(LSB=1:length(idx),BPG=idx)
for ( i in 1:nrow(idpos)) {
if (idpos[i,2]>0) {
Teff_scs[idpos[i,2]]= ipac_temp[idpos[i,1],7]
Tipo_scs[idpos[i,2]]= ipac_temp[idpos[i,1],5]
}
}
#
Modelado de T con esas Features
#
clase_luminosidad = function(x) {
# we look over Spectral SubClass
# for roman numbers V=>dwarft; III => giants I => Supergiants (II => I)
cl=c("V","III","I")
i=1
res=NA
if ( ! is.na(x)) {
while(i <= length(cl) & is.na(res)) {
j=regexpr(cl[i],x)
if (j[1] > 0) {
res=cl[i]
} else {
i=i+1
}
}
}
return(res)
}
clase_m = function(x) {
# we look over Spectral SubClass
# for roman numbers V=>dwarft; III => giants I => Supergiants (II => I)
# cl=c("M0.5","M0","M1.5","M1","M2.5","M2","M3.5","M3","M4.5","M4","M5.5","M5","M6.5","M6",
# "M7.5","M7","M8.5","M8","M9.5","M9")
cl=c("M0","M1","M2","M3","M4","M5","M6","M7","M8","M9")
i=1
res=NA
if ( ! is.na(x) ) {
while(i <= length(cl) & is.na(res)) {
j=regexpr(cl[i],x)
if (j[1] > 0) {
res=cl[i]
} else {
i=i+1
}
}
}
return(res)
}
#
#
prep_datos=function(res){
X=matrix(NA,nrow=length(res[[1]]$x),ncol=length(res))
for (j in 1:length(res)) {
X[,j]=res[[j]]$x/res[[j]]$x0
}
colnames(X)=paste("C",1:length(res),sep="")
return(X)
}
#
modelado=function(X,Y) {
#
trn=1:nrow(X)
folds=5
repeats=1
#
mseeds <- vector(mode = "list", length = (folds+1))
for(i in 1:folds) mseeds[[i]] <- sample.int(nrow(X), 100)
mseeds[[(folds+1)]] <- sample.int(1000, 1)
myControl = trainControl(method='cv', number=folds,
repeats=repeats, returnResamp='none',
returnData=FALSE, savePredictions=TRUE,
verboseIter=FALSE, allowParallel=TRUE,seeds=mseeds)
#Train some models
model_1.1 = "train(X[trn,], Y[trn], method='gbm',
trControl=myControl,
tuneGrid=expand.grid(n.trees=seq(200,500,100), n.minobsinnode=c(2,5),
interaction.depth=seq(3,15,2),
shrinkage = seq(0.01,0.1,0.02)),verbose=FALSE)"
model_1.2 = "train(X[trn,], Y[trn], method='blackboost',
trControl=myControl, tuneGrid=expand.grid(.mstop=10^(2:4),
.maxdepth=seq(5,15,3)))"
model_1.3 = "train(X[trn,], Y[trn], method='rf',
trControl=myControl,
tuneGrid=expand.grid(.mtry=(2:ncol(X))),
proximity=TRUE)"
model_1.4 = "train(X[trn,], Y[trn], method='knn',
trControl=myControl, trace=FALSE,
tuneGrid=expand.grid(k=(2:min(9,ncol(X)))))"
model_1.5 = "train(X[trn,], Y[trn], method='ppr',
trControl=myControl,tuneGrid=expand.grid(.nterms=10))"
model_1.6 = "train(X[trn,], Y[trn], method='earth',
trControl=myControl,tuneGrid=expand.grid(.degree = 3,
.nprune = (1:10) * 2), metric = 'RMSE',
maximize = FALSE)"
model_1.7 = "train(X[trn,], Y[trn], method='glm',
trControl=myControl)"
model_1.8 = "train(X[trn,], Y[trn], method='svmRadial',
trControl=myControl,tuneGrid=expand.grid(.C=10^(-3:3),
.sigma=c(10^(-4:3))))"
model_1.9 = "train(X[trn,], Y[trn], method='cubist',
trControl=myControl,commiittees=10)"
model_1.10 = "train(X[trn,], Y[trn], method='kernelpls',
trControl=myControl,tuneGrid=expand.grid(.ncomp=seq(2,ncol(X))))"
model_1.11 = "train(X[trn,], Y[trn], method='nnet',
trControl=myControl,tuneGrid=expand.grid(.size=seq(2,round(2*ncol(X)),3),
.decay=c(1.e-5,1.e-2)),
linout=TRUE, rang = 300, trace=FALSE,
maxit=10000, reltol=1.0e-11, abstol=1.0e-6)"
model_1.12 = "train(X[trn,], Y[trn], method='bagEarth',
trControl=myControl,tuneGrid=expand.grid(.nprune=(1:10)*3,.degree=1:5))"
model_1.13 = "train(X[trn,], Y[trn], method='cubist', trControl=myControl,
tuneGrid=expand.grid(.committees=3:50, .neighbors=1:7))"
#
#Make a list of all the models
lmodels=c("model_1.1", "model_1.2", "model_1.3", "model_1.4",
"model_1.5", "model_1.6", "model_1.7", "model_1.8",
"model_1.9", "model_1.10","model_1.11","model_1.12",
"model_1.13")
#
all.models_5 =list()
for (i in lmodels) {
all.models_5[[length(all.models_5)+1]] = eval(parse(text=gsub('\n','',get(i))))
}
names(all.models_5) = sapply(all.models_5, function(x) x$method)
#
ensam_5 <- caretList(X[trn,],Y[trn],methodList=c('svmRadial',
'gbm','blackboost','rf','earth','bagEarth'))
ens <- caretEnsemble(ensam_5)
#Make a linear regression ensemble
if (exists("ens")){
all.models_5[[length(all.models_5)+1]] = get("ens")
names(all.models_5)[length(all.models_5)]="ENS"
}
return(all.models_5)
}
#
# Function that returns Root Mean Squared Error
rmse <- function(error) {
sqrt(mean(error^2,na.rm=TRUE))
}
# Function that returns Mean Absolute Error
mae <- function(error) {
mean(abs(error),na.rm=TRUE)
}
#
#
if ( file.exists("~/git/M_prep_IPAC/GA_IPAC_MT_TF2_30_Models.RData")) {
load("~/git/M_prep_IPAC/GA_IPAC_MT_TF2_30_Models.RData")
} else {
res00=SelFeatures(GAs[[isol]],cnjts,bandas)
res10=SelFeatures(GAs10[[isol10]],cnjts,bandas10)
res50=SelFeatures(GAs50[[isol50]],cnjts,bandas50)
#
X00 = prep_datos(res00)
md5 = modelado(X00,as.numeric(YPT[,1]))
X10 = prep_datos(res10)
md51= modelado(X10,as.numeric(YPT[,1]))
X50 = prep_datos(res50)
md55= modelado(X50,as.numeric(YPT[,1]))
#
xf0 = prep_datos(SelFeatures(GAs[[isol]],cnjts,bfndas))
xf1 = prep_datos(SelFeatures(GAs10[[isol10]],cnjts,bfndas))
xf5 = prep_datos(SelFeatures(GAs50[[isol50]],cnjts,bfndas))
#
YTn00=predict(md5[[which(names(md5)=="rf")]],xf0)
YTn01=predict(md51[[which(names(md51)=="rf")]],xf0)
YTn05=predict(md55[[which(names(md55)=="rf")]],xf0)
YTn10=predict(md5[[which(names(md5)=="rf")]],xf1)
YTn11=predict(md51[[which(names(md51)=="rf")]],xf1)
YTn15=predict(md55[[which(names(md55)=="rf")]],xf1)
YTn50=predict(md5[[which(names(md5)=="rf")]],xf5)
YTn51=predict(md51[[which(names(md51)=="rf")]],xf5)
YTn55=predict(md55[[which(names(md55)=="rf")]],xf5)
#
Clase_lum=apply(as.data.frame(as.character(Tipo_scs),
stringsAsFactors=FALSE),1,clase_luminosidad)
Clase_m =apply(as.data.frame(as.character(Tipo_scs),
stringsAsFactors=FALSE),1,clase_m)
SpT=Clase_m
LC =Clase_lum
ref_tot0=data.frame(Name=as.character(ref$Name),
chi2d_10=ref$Chi2_T_10, chi2d_50=ref$Chi2_T_50,
YTn00,YTn01,YTn05,YTn10,YTn11,YTn15,YTn50,
YTn51,YTn55,Teff_LSB=Teff_scs,SpT=SpT,LC=LC)
#
save(bpy,bfy,YPT,YPG,YPM,vf,siz,bandas,bfndas, bandas10,
bandas50,cnjts, NG,NBLK,NBITS,MAXV,md5,md51,
md55,res00,res10,res50, X00,X10,X50,YTn00,YTn01,
YTn05,YTn10,YTn11,YTn15,YTn50,YTn51,YTn55,xf0,xf1,xf5,
ref_tot0,LC,SpT,
file="~/git/M_prep_IPAC/GA_IPAC_MT_TF2_30_Models.RData")
}
#
YTpknn00=predict(md5[[which(names(md5)=="knn")]],xf0)
YTpknn11=predict(md51[[which(names(md51)=="knn")]],xf1)
YTpknn55=predict(md55[[which(names(md55)=="knn")]],xf5)
save(YTpknn00,YTpknn11,YTpknn55,xf0,xf1,xf5,
file="~/git/M_prep_IPAC/GA_IPAC_predict_T_11F2_30.RData")
res00=SelFeatures(GAs[[isol]],cnjts,bandas)
res10=SelFeatures(GAs10[[isol10]],cnjts,bandas10)
res50=SelFeatures(GAs50[[isol50]],cnjts,bandas50)
ref_tot0=data.frame(Name=as.character(ref$Name),
chi2d_10=ref$Chi2_T_10, chi2d_50=ref$Chi2_T_50,
YTn00,YTn01,YTn05,YTn10,YTn11,YTn15,YTn50,
YTn51,YTn55,Teff_LSB=Teff_scs,SpT=SpT,LC=LC)
#
ref_tot2=ref_tot0
#
print(xtable(ldply(res00,function(x){return(c(x$name,x$name0))})),type="html")
|
|
V1
|
V2
|
|
1
|
7602-7706.4
|
7764-7868.4
|
|
2
|
7224-7328.4
|
7872-7976.4
|
|
3
|
7926-8030.4
|
7170-7274.4
|
|
4
|
7656-7760.4
|
7548-7652.4
|
|
5
|
6954-7058.4
|
7926-8030.4
|
|
6
|
6900-7004.4
|
7818-7922.4
|
|
7
|
7818-7922.4
|
7386-7490.4
|
|
8
|
7116-7220.4
|
7494-7598.4
|
|
9
|
7764-7868.4
|
7710-7814.4
|
|
10
|
7872-7976.4
|
7224-7328.4
|
print(xtable(ldply(res10,function(x){return(c(x$name,x$name0))})),type="html")
|
|
V1
|
V2
|
|
1
|
7602-7706.4
|
7008-7112.4
|
|
2
|
7656-7760.4
|
6954-7058.4
|
|
3
|
7872-7976.4
|
7224-7328.4
|
|
4
|
7062-7166.4
|
7494-7598.4
|
|
5
|
6900-7004.4
|
7818-7922.4
|
|
6
|
7278-7382.4
|
7548-7652.4
|
|
7
|
7224-7328.4
|
7872-7976.4
|
|
8
|
7332-7436.4
|
7926-8030.4
|
|
9
|
7926-8030.4
|
7170-7274.4
|
|
10
|
6954-7058.4
|
7278-7382.4
|
print(xtable(ldply(res50,function(x){return(c(x$name,x$name0))})),type="html")
|
|
V1
|
V2
|
|
1
|
7764-7868.4
|
7332-7436.4
|
|
2
|
6954-7058.4
|
7872-7976.4
|
|
3
|
7548-7652.4
|
7926-8030.4
|
|
4
|
7116-7220.4
|
7062-7166.4
|
|
5
|
7818-7922.4
|
7764-7868.4
|
|
6
|
7170-7274.4
|
7116-7220.4
|
|
7
|
6900-7004.4
|
7818-7922.4
|
|
8
|
7008-7112.4
|
7224-7328.4
|
|
9
|
7872-7976.4
|
7278-7382.4
|
|
10
|
7494-7598.4
|
7710-7814.4
|
#
ggplot(data=ref_tot2) +
geom_point(aes(x=Teff_LSB,y=chi2d_10,shape=LC),size=3) +
xlab("Theoretical T [K]") + ylab("Chi2 10 [K]") +
theme_bw() + # scale_colour_brewer(palette="Set1") +
geom_abline(position="identity", colour="gray") +
xlim(2000,4200) + ylim(2000,4200) +
guides(col=guide_legend(ncol=2))
#
ggplot(data=ref_tot2) +
geom_point(aes(x=Teff_LSB,y=chi2d_50,shape=LC),size=3) +
xlab("Theoretical T [K]") + ylab("Chi2 50 [K]") +
theme_bw() + # scale_colour_brewer(palette="Set1") +
geom_abline(position="identity", colour="gray") +
xlim(2000,4200) + ylim(2000,4200) +
guides(col=guide_legend(ncol=2))
#
ggplot(data=ref_tot2) +
geom_point(aes(x=Teff_LSB,y=YTn00,shape=LC),size=3) +
xlab("Theoretical T [K]") +
ylab("ML predicted Msnr=oo Fsnr=oo [K]") +
theme_bw() + # scale_colour_brewer(palette="Set1") +
geom_abline(position="identity", colour="gray") +
xlim(2000,4200) + ylim(2000,4200) +
guides(col=guide_legend(ncol=2))
#
ggplot(data=ref_tot2) +
geom_point(aes(x=Teff_LSB,y=YTn11,shape=LC),size=3) +
xlab("Theoretical T [K]") +
ylab("ML predicted Msnr=10 Fsnr=10 [K]") +
theme_bw() + # scale_colour_brewer(palette="Set1") +
geom_abline(position="identity", colour="gray") +
xlim(2000,4200) + ylim(2000,4200) +
guides(col=guide_legend(ncol=2))
#
ggplot(data=ref_tot2) +
geom_point(aes(x=Teff_LSB,y=YTn55,shape=LC),size=3) +
xlab("Theoretical T [K]") +
ylab("ML predicted Msnr=50 Fsnr=50 [K]") +
theme_bw() + # scale_colour_brewer(palette="Set1") +
geom_abline(position="identity", colour="gray") +
xlim(2000,4200) + ylim(2000,4200) +
guides(col=guide_legend(ncol=2))
#
lmod=c("rf","gbm","svmRadial","nnet","knn","bagEarth","kernelpls","cubist")
nmod=c("RF","GB","SVR","NNR","KNN","MARS","PLS","Rule-Regression")
for(ll in c(1:2)) {
if (ll == 2) {
pdf(file="~/git/M_prep_IPAC/GA_model_30_T.pdf",width=12,heigh=10)
}
for (i in 1:length(lmod)) {
YTm00=predict(md5[[which(names(md5)==lmod[i])[1]]],xf0)
YTm11=predict(md51[[which(names(md51)==lmod[i])[1]]],xf1)
YTm55=predict(md55[[which(names(md55)==lmod[i])[1]]],xf5)
#
ref_tot0m=data.frame(Name=as.character(ref$Name),
chi2d_10=ref$Chi2_T_10, chi2d_50=ref$Chi2_T_50,
YTm00,YTm11,YTm55,Teff_LSB=Teff_scs,SpT=SpT,LC=LC)
ref_tot2m=ref_tot0m
diffs=apply(ref_tot2m[,! colnames(ref_tot2m) %in% c("Name","SpT","Teff_LSB","LC")],
2,FUN="-",ref_tot2m[,"Teff_LSB"])
rownames(diffs)=ref_tot2m[,"Name"]
#
cat(paste("T Modelling with ",lmod[i],". Error analysis follows.",sep=""))
errors=data.frame(drmse = apply(diffs,2,rmse),dmae=apply(diffs,2,mae))
print(xtable(errors),type="html")
#
cat(paste("T Modelling with",nmod[i],". SNR=oo",sep=""))
print(ggplot(data=ref_tot2m) +
geom_point(aes(x=Teff_LSB,y=YTm00,shape=LC),size=3) +
xlab("Theoretical T [K]") +
ylab(paste(nmod[i]," predicted SNR=oo [K]","")) +
theme_bw() + # scale_colour_brewer(palette="Set1") +
geom_abline(position="identity", colour="red") +
xlim(2000,4200) + ylim(2000,4200) +
guides(col=guide_legend(ncol=2)))
#
cat(paste("T Modelling with",nmod[i],". SNR=10",sep=""))
print(ggplot(data=ref_tot2m) +
geom_point(aes(x=Teff_LSB,y=YTm11,shape=LC),size=3) +
xlab("Theoretical T [K]") +
ylab(paste(nmod[i]," predicted SNR=10 [K]","")) +
theme_bw() + # scale_colour_brewer(palette="Set1") +
geom_abline(position="identity", colour="red") +
xlim(2000,4200) + ylim(2000,4200) +
guides(col=guide_legend(ncol=2)))
#
cat(paste("T Modelling with",nmod[i],". SNR=50",sep=""))
print(ggplot(data=ref_tot2m) +
geom_point(aes(x=Teff_LSB,y=YTm55,shape=LC),size=3) +
xlab("Theoretical T [K]") +
ylab(paste(nmod[i]," predicted SNR=50 [K]","")) +
theme_bw() + # scale_colour_brewer(palette="Set1") +
geom_abline(position="identity", colour="red") +
xlim(2000,4200) + ylim(2000,4200) +
guides(col=guide_legend(ncol=2)))
}
}
GA_r3 Loading required package: randomForest
GA_r3 randomForest 4.6-7
GA_r3 Type rfNews() to see new features/changes/bug fixes.
T Modelling with rf. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
213.31
|
142.29
|
|
YTm11
|
252.26
|
181.39
|
|
YTm55
|
375.31
|
285.51
|
T Modelling withRF. SNR=oo
T Modelling withRF. SNR=10
T Modelling withRF. SNR=50
GA_r3 Loading required package: gbm
GA_r3 Loading required package: survival
GA_r3
GA_r3 Attaching package: 'survival'
GA_r3
GA_r3 The following object is masked from 'package:caret':
GA_r3
GA_r3 cluster
GA_r3
GA_r3 Loading required package: splines
GA_r3 Loaded gbm 2.1.1
T Modelling with gbm. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
252.36
|
195.10
|
|
YTm11
|
308.49
|
241.92
|
|
YTm55
|
439.70
|
354.19
|
T Modelling withGB. SNR=oo
T Modelling withGB. SNR=10
T Modelling withGB. SNR=50
T Modelling with svmRadial. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
1263.58
|
1168.71
|
|
YTm11
|
313.69
|
228.18
|
|
YTm55
|
504.40
|
433.35
|
T Modelling withSVR. SNR=oo
T Modelling withSVR. SNR=10
T Modelling withSVR. SNR=50
T Modelling with nnet. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
1741.56
|
1646.00
|
|
YTm11
|
386.57
|
313.57
|
|
YTm55
|
605.90
|
449.79
|
T Modelling withNNR. SNR=oo
T Modelling withNNR. SNR=10
T Modelling withNNR. SNR=50
T Modelling with knn. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
199.70
|
119.72
|
|
YTm11
|
226.99
|
169.73
|
|
YTm55
|
276.19
|
175.42
|
T Modelling withKNN. SNR=oo
T Modelling withKNN. SNR=10
T Modelling withKNN. SNR=50
GA_r3 Loading required package: earth
GA_r3 Loading required package: plotmo
GA_r3 Loading required package: plotrix
GA_r3 Loading required package: TeachingDemos
T Modelling with bagEarth. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
434.32
|
363.85
|
|
YTm11
|
301.41
|
207.20
|
|
YTm55
|
470.88
|
367.38
|
T Modelling withMARS. SNR=oo
T Modelling withMARS. SNR=10
T Modelling withMARS. SNR=50
T Modelling with kernelpls. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
618.61
|
537.56
|
|
YTm11
|
315.16
|
210.06
|
|
YTm55
|
436.10
|
314.34
|
T Modelling withPLS. SNR=oo
T Modelling withPLS. SNR=10
T Modelling withPLS. SNR=50
T Modelling with cubist. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
1149.43
|
1016.61
|
|
YTm11
|
275.52
|
183.72
|
|
YTm55
|
435.17
|
361.92
|
T Modelling withRule-Regression. SNR=oo
T Modelling withRule-Regression. SNR=10
T Modelling withRule-Regression. SNR=50T Modelling with rf. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
213.31
|
142.29
|
|
YTm11
|
252.26
|
181.39
|
|
YTm55
|
375.31
|
285.51
|
T Modelling withRF. SNR=ooT Modelling withRF. SNR=10T Modelling withRF. SNR=50T Modelling with gbm. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
252.36
|
195.10
|
|
YTm11
|
308.49
|
241.92
|
|
YTm55
|
439.70
|
354.19
|
T Modelling withGB. SNR=ooT Modelling withGB. SNR=10T Modelling withGB. SNR=50T Modelling with svmRadial. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
1263.58
|
1168.71
|
|
YTm11
|
313.69
|
228.18
|
|
YTm55
|
504.40
|
433.35
|
T Modelling withSVR. SNR=ooT Modelling withSVR. SNR=10T Modelling withSVR. SNR=50T Modelling with nnet. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
1741.56
|
1646.00
|
|
YTm11
|
386.57
|
313.57
|
|
YTm55
|
605.90
|
449.79
|
T Modelling withNNR. SNR=ooT Modelling withNNR. SNR=10T Modelling withNNR. SNR=50T Modelling with knn. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
199.70
|
119.72
|
|
YTm11
|
226.99
|
169.73
|
|
YTm55
|
276.19
|
175.42
|
T Modelling withKNN. SNR=ooT Modelling withKNN. SNR=10T Modelling withKNN. SNR=50T Modelling with bagEarth. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
434.32
|
363.85
|
|
YTm11
|
301.41
|
207.20
|
|
YTm55
|
470.88
|
367.38
|
T Modelling withMARS. SNR=ooT Modelling withMARS. SNR=10T Modelling withMARS. SNR=50T Modelling with kernelpls. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
618.61
|
537.56
|
|
YTm11
|
315.16
|
210.06
|
|
YTm55
|
436.10
|
314.34
|
T Modelling withPLS. SNR=ooT Modelling withPLS. SNR=10T Modelling withPLS. SNR=50T Modelling with cubist. Error analysis follows.
|
|
drmse
|
dmae
|
|
chi2d_10
|
168.61
|
122.99
|
|
chi2d_50
|
155.71
|
99.24
|
|
YTm00
|
1149.43
|
1016.61
|
|
YTm11
|
275.52
|
183.72
|
|
YTm55
|
435.17
|
361.92
|
T Modelling withRule-Regression. SNR=ooT Modelling withRule-Regression. SNR=10T Modelling withRule-Regression. SNR=50
dev.off()
png 2
#
That’s all folks.