IPAC (Met) based on class of features 1-Fs/Fc. Modeled with BT-Settl 2011-13
Introduction
After the analysis carried out for M stars within the IPAC analysis, we wnat to apply the band I features discovered there to estimate physical parameters from IRTF dataset. Target values came from Table III from Cesetti et al.
## 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 4.3
##
## 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.3-0
##
## 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-3
##
## 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] TRUETarget 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 BT-SETTL preprocesed data
#
# 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)
}
rm(list=c("model","model0","model2","model3","model4","model5","model6",
"model7","model8","ref","reftot"))
#Now, we will read the M stars from the dataset we have prepared when we forecasted the Temperature.
#
siz=10
if ( file.exists("~/git/M_prep_IPAC/Features_BT-settl_v3.RData")) {
load("~/git/M_prep_IPAC/Features_BT-settl_v3.RData")
} else {
stop(paste('File ~/git/M_prep_IPAC/Features_BT-settl_v3.RData does not exist', sep=''))
}Distance $ ^2 $ between spectra
Let’s find the closest spectra at different SNR ratio per IPAC spectum.
Closest spectrum for BT-Settl SNR=50
minchid=function(x,lib){
chid=function(x,y){
if ( nrow(x$data[[1]]) == nrow(y$data[[1]])) {
return(sum((x$data[[1]][,2]-y$data[[1]][,2])^2)/nrow(x$data[[1]]))
} else {
return (5000) # Enormous distance as it is not valid
}
}
ldist=ldply(lib,chid,x)
idx = which.min(ldist[,1])[1]
res = data.frame(T=lib[[idx]]$stellarp[1],logg=lib[[idx]]$stellarp[2],
met=lib[[idx]]$stellarp[3],name=lib[[idx]]$name,
dist=ldist[idx,1])
return(res)
}
#
chi2d_50=ldply(bf_clean,minchid,bp_50)
#And the same with SNR=10:
Closest spectrum for BT-Settl SNR=10
#
chi2d_10=ldply(bf_clean,minchid,bp_10)
#
#
# Reading the IPAC stellar type
load("~/git/M_prep/IPAC_stellar_stype.RData")
#
ref=data.frame(Name=as.vector(do.call(rbind,lapply(bf_clean,function(x){
return(sub('p','+',sub('.txt','',x$starname)))}))),
chi2_50=chi2d_50[,4],chi2_10=chi2d_10[,4])
ref$Name=as.character(ref$Name)
refm=ref
#ICA analysis for Met
Looking into the gravity properties ## SNR=10
# We prepare the coefs
if ( file.exists("~/git/M_prep_IPAC/BT-Settl-2013-PPR_bp10-M.RData")) {
load("~/git/M_prep_IPAC/BT-Settl-2013-PPR_bp10-M.RData")
} else {
dd2=as.data.frame(do.call(rbind,lapply(bp_10,function(x){return(x$data[[1]][,2])})))
jd2=JADE(dd2,10)
ss2=as.data.frame(t(jd2$W %*% t(dd2-jd2$Xmu)))
ss2[,(ncol(ss2)+1)]=unlist(lapply(bp_10,function(x){return(x$stellarp[3])}))
colnames(ss2)[ncol(ss2)]="Met"
#
folds=5
repeats=1
myControl = trainControl(method='cv', number=folds,
repeats=repeats, returnResamp='none',
returnData=FALSE, savePredictions=TRUE,
verboseIter=FALSE, allowParallel=TRUE,
index=createMultiFolds(ss2[,ncol(ss2)],
k=folds, times=repeats))
model2=train(ss2[,-ncol(ss2)], ss2[,ncol(ss2)], method='ppr',
trControl=myControl,tuneGrid=expand.grid(.nterms=3:(ncol(ss2)-1)))
d2=as.data.frame(do.call(rbind,lapply(bf_clean,
function(x){return(ldply(x$data,
function(x){return(x)})[,2])})))
si2=as.data.frame(t(jd2$W %*% t(d2-jd2$Xmu)))
M_coef2=predict(model2,si2)
save(M_coef2,si2,d2,model2,myControl,folds,repeats,ss2,jd2,dd2,
file="~/git/M_prep_IPAC/BT-Settl-2013-PPR_bp10-M.RData")
}
#SNR=50
# We prepare the coefs
if ( file.exists("~/git/M_prep_IPAC/BT-Settl-2013-PPR_bp50-M.RData")) {
load("~/git/M_prep_IPAC/BT-Settl-2013-PPR_bp50-M.RData")
} else {
dd=as.data.frame(do.call(rbind,lapply(bp_50,function(x){return(x$data[[1]][,2])})))
jd=JADE(dd,10)
ss=as.data.frame(t(jd$W %*% t(dd-jd$Xmu)))
ss[,(ncol(ss)+1)]=unlist(lapply(bp_50,function(x){return(x$stellarp[3])}))
colnames(ss)[ncol(ss)]="Met"
#
folds=5
repeats=1
myControl = trainControl(method='cv', number=folds,
repeats=repeats, returnResamp='none',
returnData=FALSE, savePredictions=TRUE,
verboseIter=FALSE, allowParallel=TRUE,
index=createMultiFolds(ss[,ncol(ss)],
k=folds, times=repeats))
model=train(ss[,-ncol(ss)], ss[,ncol(ss)], method='ppr',
trControl=myControl,tuneGrid=expand.grid(.nterms=3:(ncol(ss)-1)))
#
d1=as.data.frame(do.call(rbind,lapply(bf_clean,function(x){return(x$data[[1]][,2])})))
si=as.data.frame(t(jd$W %*% t(d1-jd$Xmu)))
M_coef=predict(model,si)
save(M_coef,si,dd,model,myControl,folds,repeats,ss,jd,dd,
file="~/git/M_prep_IPAC/BT-Settl-2013-PPR_bp50-M.RData")
}
#T foreseen by Luminosity models
# Loading IPAC prediction from Dr Sarro (objeto ipac_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]
}
}
#
#
# Leemos las estimaciones de Gravedad para estrellas IPAC de Rojas y Ayala con Incertidumbre
# leemos el bf_clean original para buscar nombres no truncados
if ( file.exists("~/git/M_prep_IPAC/IPAC_Names_TeorPar-M.RData")) {
load("~/git/M_prep_IPAC/IPAC_Names_TeorPar-M.RData")
} else {
load("~/git/M_prep_Noise/M_prep_clean_IPAC.RData")
rm(list=c("bq","d_ipac"))
name_IPAC=as.character(lapply(bf_clean, function(x){return(sub('.txt','',x$name))}))
rm("bq_clean")
ref_ipc=data.frame(Name=as.vector(do.call(rbind,lapply(bf_clean,function(x){
return(sub('p','+',sub('.txt','',x$starname)))}))))
ref_ipc$ICA_Met_10 = M_coef2
ref_ipc$ICA_Met_50 = M_coef
ref_ipc$SpT=NA
ref_ipc$M_teo=NA
ref_ipc$DM = NA
ref_ipc$T_teo = NA
ref_ipc$DT = NA
#
ddk=read.table(file="IPAC-NevesIII.tsv",header=FALSE,stringsAsFactors=FALSE)
ddj=read.table(file="U-IPAC-Rojas-uncert.tsv",sep="|",header=FALSE,stringsAsFactors=FALSE)
colnames(ddj)[c(1,11,12,13,14)]=c("Name","T","DT","Met","Dmet")
#
nombres=data.frame(orig=name_IPAC,LHS=rep(NA,length(name_IPAC)),GL=rep(NA,length(name_IPAC)),
OT=rep(NA,length(name_IPAC)),M2=rep(NA,length(name_IPAC)))
for ( i in (1:length(name_IPAC))) {
url=paste("http://ldwarf.ipac.caltech.edu/archive/mdwarf_spectra/",name_IPAC[i],'.txt',sep="")
ff =readLines(url,n=6)
lhs=trim(gsub("'","",strsplit(ff[grep("LHS_NAME",ff)],'=',fixed=TRUE)[[1]][2]))
gl =trim(gsub("'","",strsplit(ff[grep("GL_NAME",ff)],'=',fixed=TRUE)[[1]][2]))
ot =trim(gsub("'","",strsplit(ff[grep("OTHER_NM",ff)],'=',fixed=TRUE)[[1]][2]))
m2 =trim(gsub("'","",strsplit(ff[grep("2M_DESIG",ff)],'=',fixed=TRUE)[[1]][2]))
nombres[i,]=c(name_IPAC[i],lhs,gl,ot,m2)
}
#
vv = rep(NA,nrow(ref_ipc))
for (i in (1:nrow(ref_ipc))) {
pos=which(nombres[i,5] == stellar_stype[,1])
if ( length(pos) > 0 ) {
ref_ipc$SpT[i] = as.character(stellar_stype[pos[1],4])
}
pos=which(nombres[i,5] == trim(ddj[,1]))
if ( length(pos) > 0 ) {
vv[i]=pos[1]
ref_ipc$T_teo[i] = as.numeric(ddj[pos[1],"T"])
ref_ipc$DT[i] = as.numeric(ddj[pos[1],"DT"])
ref_ipc$M_teo[i] = as.numeric(ddj[pos[1],"Met"])
ref_ipc$DM[i] = as.numeric(ddj[pos[1],"Dmet"])
}
}
ref_ipc$Name=nombres[,1]
save(nombres,ref_ipc,file="~/git/M_prep_IPAC/IPAC_Names_TeorPar-M.RData")
}
refm$Name =nombres[,1]
#Preparation of potential features for Met
#
#
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: doParallelrequire(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)
}
fitnessM = function(string) {
par = decode(string)
if ( 0 %in% unlist(apply(par,1,sd)) ) {
return (-10000*NBLK*NBITS)
}
if ( length(unique(sort(par[,1]))) < nrow(par) |
length(unique(sort(par[,2]))) < nrow(par) ) {
return (-10000*NBLK*NBITS)
}
X = as.data.frame(matrix(0,nrow=nrow(bandas[[cnjts[1]]]$mat),ncol=(nrow(par)+1)))
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="-")
}
X[,ncol(X)]= YPT
colnames(X)[ncol(X)]="T"
mod = lm.fit(as.matrix(X), y)
class(mod) = "lm"
return( -BIC(mod))
}
fitnessM10 = function(string) {
par = decode(string)
if ( 0 %in% unlist(apply(par,1,sd)) ) {
return (-10000*NBLK*NBITS)
}
if ( length(unique(sort(par[,1]))) < nrow(par) |
length(unique(sort(par[,2]))) < nrow(par) ) {
return (-10000*NBLK*NBITS)
}
X = as.data.frame(matrix(0,nrow=nrow(bandas10[[cnjts[1]]]$mat),ncol=(nrow(par)+1)))
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="-")
}
X[,ncol(X)]= YPT
colnames(X)[ncol(X)]="T"
mod = lm.fit(as.matrix(X), y)
class(mod) = "lm"
return( -BIC(mod))
}
fitnessM50 = function(string) {
par = decode(string)
if ( 0 %in% unlist(apply(par,1,sd)) ) {
return (-10000*NBLK*NBITS)
}
if ( length(unique(sort(par[,1]))) < nrow(par) |
length(unique(sort(par[,2]))) < nrow(par) ) {
return (-10000*NBLK*NBITS)
}
X = as.data.frame(matrix(0,nrow=nrow(bandas50[[cnjts[1]]]$mat),ncol=(nrow(par)+1)))
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="-")
}
X[,ncol(X)]= YPT
colnames(X)[ncol(X)]="T"
mod = lm.fit(as.matrix(X), y)
class(mod) = "lm"
return( -BIC(mod))
}
#
SelFeatures=function(x,sets=cnjts,dat=bandas,fitness.=fitnessM,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.=fitnessM,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_NM11F2-2013_data.RData")) {
j=1
load("~/git/M_prep_IPAC/GA_NM11F2-2013_data.RData")
tt50=list()
GAs50 = list()
solsM50=list()
tt10=list()
GAs10 = list()
solsM10=list()
tt=list()
GAs=list()
solsM=list()
while(j <= NITER) {
cat(paste("Reading Iteration:",j,"<br>",sep=""))
flush.console()
load(paste("~/git/M_prep_IPAC/GA_2013_M_",sprintf("%02d",j),"_50_f2_data.RData",sep=""))
tt50[[j]]= t_tot
GAs50[[j]]=GA1
solsM50[[j]]=FeaturesExt(GA1,12)
load(paste("~/git/M_prep_IPAC/GA_2013_M_",sprintf("%02d",j),"_10_f2_data.RData",sep=""))
tt10[[j]]= t_tot
GAs10[[j]]=GA1
solsM10[[j]]=FeaturesExt(GA1,12)
load(paste("~/git/M_prep_IPAC/GA_2013_M_",sprintf("%02d",j),"_f2_data.RData",sep=""))
tt[[j]]= t_tot
GAs[[j]]=GA1
solsM[[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=6")
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=10
bandas=list()
bfndas=list()
bandas10=list()
bandas50=list()
for (i in seq(1,10,5)) { # 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]])
}
}
#
# BE CAREFUL !!! y should be in retionship with the variable being predicted
y=as.numeric(unlist(YPM))
cnjts = c("s=1","s=6")
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,fitnessM,fitnessM10,fitnessM50, encode,decode,
SelFeatures,file="~/git/M_prep_IPAC/GA_NM11F2-2013_data.RData")
#
tt=list()
GAs = list()
solsM=list()
j=1
while(j <= NITER) {
cat(paste("Starting Iteration:",j,"<br>",sep=""))
flush.console()
tt[[j]]= system.time({GA1 = ga(type="binary", fitness = fitnessM, nBits = NBITS*2*NG,
monitor = FALSE,parallel=18,maxiter=1000,popSize=8000,keepBest=TRUE)})
t_tot=tt[[j]]
GAs[[j]]=GA1
save(fitnessM,encode,decode,SelFeatures,GA1,t_tot,j,
file=paste("~/git/M_prep_IPAC/GA_2013_M_",sprintf("%02d",j),"_f2_data.RData",sep=""))
solsM[[j]]=FeaturesExt(GA1,12)
j=j+1
}
#
print(xtable(ldply(solsG,function(x){return(x$vpos[1])})),type="html")
#
j=1
tt10=list()
GAs10 = list()
solsM10=list()
while(j <= NITER) {
cat(paste("Starting Iteration:",j,"<br>",sep=""))
flush.console()
tt10[[j]]= system.time({GA1 = ga("binary", fitness = fitnessM10, nBits = NBITS*2*NG,
monitor = FALSE,parallel=18,maxiter=1000,popSize=8000,keepBest=TRUE)})
t_tot=tt10[[j]]
GAs10[[j]]=GA1
save(fitnessM10,encode,decode,SelFeatures,GA1,t_tot,j,
file=paste("~/git/M_prep_IPAC/GA_2013_M_",sprintf("%02d",j),"_10_f2_data.RData",sep=""))
solsM10[[j]]=FeaturesExt(GA1,12)
j=j+1
}
#
print(xtable(ldply(solsG10,function(x){return(x$vpos[1])})),type="html")
#
j=1
tt50=list()
GAs50 = list()
solsM50=list()
while(j <= NITER) {
cat(paste("Starting Iteration:",j,"<br>",sep=""))
flush.console()
tt50[[j]]= system.time({GA1 = ga("binary", fitness = fitnessM50, nBits = NBITS*2*NG,
monitor = FALSE,parallel=18,maxiter=1000,popSize=8000,keepBest=TRUE)})
t_tot=tt50[[j]]
GAs50[[j]]=GA1
save(fitnessM50,encode,decode,SelFeatures,GA1,t_tot,j,
file=paste("~/git/M_prep_IPAC/GA_2013_M_",sprintf("%02d",j),"_50_f2_data.RData",sep=""))
solsM50[[j]]=FeaturesExt(GA1,12)
j=j+1
}
#
print(xtable(ldply(solsM50,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(solsM,function(x){return(x$vpos[1])})),type="html")| V1 | |
|---|---|
| 1 | -233.58 |
| 2 | -210.08 |
| 3 | -314.10 |
| 4 | -290.51 |
| 5 | -160.79 |
| 6 | -91.53 |
| 7 | -280.71 |
| 8 | -301.00 |
| 9 | -274.60 |
# ISOL holds the best solution !
isol = which.max(as.numeric(ldply(solsM,function(x){return(x$vpos[1])})[,1]))
isol[1] 6
print(xtable(ldply(solsM10,function(x){return(x$vpos[1])})),type="html")| V1 | |
|---|---|
| 1 | -845.65 |
| 2 | -860.74 |
| 3 | -750.07 |
| 4 | -688.50 |
| 5 | -860.40 |
| 6 | -789.53 |
| 7 | -792.75 |
| 8 | -816.78 |
| 9 | -818.20 |
# ISOL holds the best solution !
isol10 = which.max(as.numeric(ldply(solsM10,function(x){return(x$vpos[1])})[,1]))
isol10 [1] 4
#
print(xtable(ldply(solsM50,function(x){return(x$vpos[1])})),type="html")| V1 | |
|---|---|
| 1 | -449.92 |
| 2 | -354.46 |
| 3 | -411.27 |
| 4 | -381.72 |
| 5 | -481.90 |
| 6 | -335.96 |
| 7 | -283.83 |
| 8 | -433.53 |
| 9 | -350.14 |
# ISOL holds the best solution !
isol50 = which.max(as.numeric(ldply(solsM50,function(x){return(x$vpos[1])})[,1]))
isol50 [1] 7
##
clase_luminosidad = function(x) {
# we look over Spectral SubClass
# for roman numbers V=>dwarft; III => giants I => Supergiants (II => I)
cl=c("V","IV","III","II","I")
i=1
res=NA
if ( ! is.na(x)) {
while(i <= length(cl) & is.na(res)) {
j=regexpr(cl[i],x)
if (length(j) > 0 ) {
if ( j[1] > 0 ) {
res=cl[i]
} else {
i=i+1
}
}
}
}
return(res)
}
clase_spt = function(x) {
# we look over Spectral SubType
# for roman numbers M0 ... M9
cl=c("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")
i=1
res=NA
if ( ! is.na(x)) {
while(i <= length(cl) & is.na(res)) {
j=regexpr(cl[i],x)
if (length(j) > 0 ) {
if ( j[1] > 0 ) {
res=cl[i]
} else {
i=i+1
}
}
}
}
return(res)
}
#
prep_datosM=function(res,T){
X=matrix(NA,nrow=length(res[[1]]$x),ncol=(length(res)+1))
for (j in 1:length(res)) {
wini=-eval(parse(text=res[[j]]$name))
winj=-eval(parse(text=res[[j]]$name0))
X[,j]=wini-res[[j]]$x/(res[[j]]$x0/winj)
}
X[,(length(res)+1)] = T
colnames(X)=c(paste("C",1:length(res),sep=""),"T")
return(X)
}
#
# 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)
}
#
#
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)
}
#
# Loading the predicted temperature
load(file="~/git/M_prep_IPAC/GA_predict_T_11F2.RData")
if ( file.exists("~/git/M_prep_IPAC/GA_IPAC_data_MM-2013_GF2_Models.RData")) {
load("~/git/M_prep_IPAC/GA_IPAC_data_MM-2013_GF2_Models.RData")
} else {
#
# Objects YTpknn00,YTpknn11,YTpknn55 are read from the file
# They were determined by the script prep_GA_case01_NT11F2_v1.Rmd
#
res00=SelFeatures(GAs[[isol]],cnjts,bandas)
res10=SelFeatures(GAs10[[isol10]],cnjts,bandas10)
res50=SelFeatures(GAs50[[isol50]],cnjts,bandas50)
#
X00 = prep_datosM(res00,as.numeric(YPT[,1]))
md5 = modelado(X00,as.numeric(YPM[,1]))
X10 = prep_datosM(res10,as.numeric(YPT[,1]))
md51= modelado(X10,as.numeric(YPM[,1]))
X50 = prep_datosM(res50,as.numeric(YPT[,1]))
md55= modelado(X50,as.numeric(YPM[,1]))
#
xf0 = prep_datosM(SelFeatures(GAs[[isol]],cnjts,bfndas),YTpknn00)
xf1 = prep_datosM(SelFeatures(GAs10[[isol10]],cnjts,bfndas),YTpknn11)
xf5 = prep_datosM(SelFeatures(GAs50[[isol50]],cnjts,bfndas),YTpknn55)
#
YMn00=predict(md5[[which(names(md5)=="rf")]],xf0)
YMn01=predict(md51[[which(names(md51)=="rf")]],xf0)
YMn05=predict(md55[[which(names(md55)=="rf")]],xf0)
YMn10=predict(md5[[which(names(md5)=="rf")]],xf1)
YMn11=predict(md51[[which(names(md51)=="rf")]],xf1)
YMn15=predict(md55[[which(names(md55)=="rf")]],xf1)
YMn50=predict(md5[[which(names(md5)=="rf")]],xf5)
YMn51=predict(md51[[which(names(md51)=="rf")]],xf5)
YMn55=predict(md55[[which(names(md55)=="rf")]],xf5)
# Correciones a las estimaciones de la tabla de Cesetti et al ya aplicadas
# ddj=ddj[! ddj$Name %in% toberemoved,]
# for ( i in 1:nrow(tobeupdated)) {
# ddj[ddj$Name %in% tobeupdated[i,"name"],"Teff"]=tobeupdated[i,"Teff"]
# }
ref=refm
ref_tot0i=data.frame(Name=as.character(ref$Name),
chi2d_10=ref$chi2_10, chi2d_50=ref$chi2_50,
Met_ICA_10=M_coef2,Met_ICA_50=M_coef,Met_rf_00=YMn00,Met_rf_11=YMn11,Met_rf_55=YMn55)
ref_tot1i=ref_ipc
ref_tot2i=merge(ref_tot0i,ref_tot1i,by="Name")
#
save(bpy,bfy,YPT,YPG,YPM,vf,siz,bandas,bfndas, bandas10,bandas50,cnjts,ref,
NG,NBLK,NBITS,MAXV,md5,md51,md55,res00,res10,res50,
X00,X10,X50,YMn00,YMn01,YMn05,YMn10,YMn11,YMn15,YMn50,YMn51,YMn55,xf0,xf1,xf5,
ref_tot0i,ref_tot1i,ref_tot2i,
file="~/git/M_prep_IPAC/GA_IPAC_data_MM-2013_GF2_Models.RData")
}
#
res00=SelFeatures(GAs[[isol]],cnjts,bandas)
res10=SelFeatures(GAs10[[isol10]],cnjts,bandas10)
res50=SelFeatures(GAs50[[isol50]],cnjts,bandas50)
ref=refm
ref_tot0=data.frame(Name=as.character(ref$Name),
chi2d_10=ref$chi2_10, chi2d_50=ref$chi2_50,
Met_ICA_10=M_coef2,Met_ICA_50=M_coef,Met_rf_00=YMn00,Met_rf_11=YMn11,Met_rf_55=YMn55)
ref_tot1=ref_ipc
ref_tot2=merge(ref_tot0,ref_tot1,by="Name")
colnames(ref_tot2)[which(colnames(ref_tot2)=="M_teo")]="Met"
ref_tot2$LC=apply(as.data.frame(ref_tot2[,"SpT"]),1,clase_luminosidad)
ref_tot2$Sp=apply(as.data.frame(ref_tot2[,"SpT"]),1,clase_spt)
#
print(xtable(ldply(res00,function(x){return(c(x$name,x$name0))})),type="html")| V1 | V2 | |
|---|---|---|
| 1 | 7188-7220.4 | 7854-7886.4 |
| 2 | 7080-7112.4 | 7926-7958.4 |
| 3 | 7116-7148.4 | 7098-7130.4 |
| 4 | 7422-7454.4 | 7836-7868.4 |
| 5 | 7350-7382.4 | 7998-8030.4 |
| 6 | 7224-7256.4 | 7818-7850.4 |
| 7 | 7710-7742.4 | 7062-7094.4 |
| 8 | 7476-7508.4 | 7944-7976.4 |
| 9 | 7134-7166.4 | 7584-7616.4 |
| 10 | 7836-7868.4 | 7278-7310.4 |
print(xtable(ldply(res10,function(x){return(c(x$name,x$name0))})),type="html")| V1 | V2 | |
|---|---|---|
| 1 | 7692-7724.4 | 7026-7058.4 |
| 2 | 6900-6932.4 | 7008-7040.4 |
| 3 | 7350-7382.4 | 7908-7940.4 |
| 4 | 6918-6950.4 | 6900-6932.4 |
| 5 | 7098-7130.4 | 7314-7346.4 |
| 6 | 7440-7472.4 | 7872-7904.4 |
| 7 | 7134-7166.4 | 7962-7994.4 |
| 8 | 7368-7400.4 | 7926-7958.4 |
| 9 | 7080-7112.4 | 7044-7076.4 |
| 10 | 7044-7076.4 | 7980-8012.4 |
print(xtable(ldply(res50,function(x){return(c(x$name,x$name0))})),type="html")| V1 | V2 | |
|---|---|---|
| 1 | 7098-7130.4 | 7926-7958.4 |
| 2 | 7188-7220.4 | 7962-7994.4 |
| 3 | 7368-7400.4 | 7980-8012.4 |
| 4 | 7116-7148.4 | 7872-7904.4 |
| 5 | 7062-7094.4 | 7206-7238.4 |
| 6 | 7584-7616.4 | 7170-7202.4 |
| 7 | 6936-6968.4 | 6918-6950.4 |
| 8 | 7692-7724.4 | 7890-7922.4 |
| 9 | 7134-7166.4 | 7548-7580.4 |
| 10 | 7494-7526.4 | 7998-8030.4 |
#
ggplot(data=ref_tot2) +
geom_point(aes(x=Met,y=chi2d_10,shape=LC),size=3) +
xlab("Theoretical Met [dex]") + ylab("Chi2 10 [dex]") +
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=Met,y=chi2d_50,shape=LC),size=3) +
xlab("Theoretical Met [dex]") + ylab("Chi2 50 [dex]") +
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=Met,y=Met_rf_00,shape=LC),size=3) +
xlab("Theoretical Met [dex]") + 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=Met,y=Met_rf_11,shape=LC),size=3) +
xlab("Theoretical Met [dex]") + 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=Met,y=Met_rf_55,shape=LC),size=3) +
xlab("Theoretical Met [dex]") + 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) +
#
#
lmod=c("rf","gbm","svmRadial","nnet","knn","bagEarth","kernelpls","cubist")
nmod=c("RF","GB","SVR","NNR","KNN","MARS","PLS","Rule-Regression")
for (i in 1:length(lmod)) {
YMm00=predict(md5[[which(names(md5)==lmod[i])[1]]],xf0)
YMm11=predict(md51[[which(names(md51)==lmod[i])[1]]],xf1)
YMm55=predict(md55[[which(names(md55)==lmod[i])[1]]],xf5)
#
ref=refm
ref_tot0m=data.frame(Name=as.character(ref$Name),
chi2d_10=ref$chi2_10, chi2d_50=ref$chi2_50,
Met_ICA_10=M_coef2,Met_ICA_50=M_coef,
YMm00,YMm11,YMm55,Tknn0=YTpknn00,Tknn1=YTpknn11,Tknn5=YTpknn55)
ref_tot1m=ref_ipc
ref_tot2m=merge(ref_tot0m,ref_tot1m,by="Name")
diffs=apply(ref_tot2m[,! colnames(ref_tot2m) %in% c("Name","SpT","Teff","LC","M_teo",
"DM","T_teo","DT")],2,FUN="-",ref_tot2m[,"M_teo"])
rownames(diffs)=ref_tot2m[,"Name"]
colnames(ref_tot2m)[which(colnames(ref_tot2m)=="M_teo")]="Met"
ref_tot2m$LC=apply(as.data.frame(ref_tot2m[,"SpT"]),1,clase_luminosidad)
ref_tot2m$Sp=apply(as.data.frame(ref_tot2m[,"SpT"]),1,clase_spt)
#
cat(paste("Metallicity 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("Met Modelling with: ",nmod[i],". SNR=oo",sep=""))
print(ggplot(data=ref_tot2m) +
geom_point(aes(x=Met,y=YMm00,shape=LC),size=3) +
xlab("Theoretical Met [dex]") +
ylab(paste(nmod[i]," predicted SNR=oo [dex]","")) +
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("Met Modelling with: ",nmod[i],". SNR=10",sep=""))
print(ggplot(data=ref_tot2m) +
geom_point(aes(x=Met,y=YMm11,shape=LC),size=3) +
xlab("Theoretical Met [dex]") +
ylab(paste(nmod[i]," predicted SNR=10 [dex]","")) +
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("Met Modelling with: ",nmod[i],". SNR=50",sep=""))
print(ggplot(data=ref_tot2m) +
geom_point(aes(x=Met,y=YMm55,shape=LC),size=3) +
xlab("Theoretical Met [dex]") +
ylab(paste(nmod[i]," predicted SNR=50 [dex]","")) +
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.| drmse | dmae | |
|---|---|---|
| chi2d_10 | 0.55 | 0.41 |
| chi2d_50 | 0.51 | 0.38 |
| Met_ICA_10 | 0.81 | 0.72 |
| Met_ICA_50 | 0.73 | 0.60 |
| YMm00 | 0.24 | 0.20 |
| YMm11 | 0.51 | 0.43 |
| YMm55 | 0.69 | 0.63 |
| Tknn0 | 3283.32 | 3262.90 |
| Tknn1 | 3295.15 | 3279.57 |
| Tknn5 | 3262.69 | 3241.80 |
| ICA_Met_10 | 0.81 | 0.72 |
| ICA_Met_50 | 0.73 | 0.60 |
Met Modelling with: RF. SNR=oo
Met Modelling with: RF. SNR=10
Met Modelling with: RF. 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
Metallicity Modelling with gbm. Error analysis follows.
| drmse | dmae | |
|---|---|---|
| chi2d_10 | 0.55 | 0.41 |
| chi2d_50 | 0.51 | 0.38 |
| Met_ICA_10 | 0.81 | 0.72 |
| Met_ICA_50 | 0.73 | 0.60 |
| YMm00 | 0.32 | 0.25 |
| YMm11 | 0.64 | 0.50 |
| YMm55 | 0.82 | 0.76 |
| Tknn0 | 3283.32 | 3262.90 |
| Tknn1 | 3295.15 | 3279.57 |
| Tknn5 | 3262.69 | 3241.80 |
| ICA_Met_10 | 0.81 | 0.72 |
| ICA_Met_50 | 0.73 | 0.60 |
Met Modelling with: GB. SNR=10
Met Modelling with: GB. SNR=50
Metallicity Modelling with svmRadial. Error analysis follows.
| drmse | dmae | |
|---|---|---|
| chi2d_10 | 0.55 | 0.41 |
| chi2d_50 | 0.51 | 0.38 |
| Met_ICA_10 | 0.81 | 0.72 |
| Met_ICA_50 | 0.73 | 0.60 |
| YMm00 | 3.51 | 2.82 |
| YMm11 | 0.45 | 0.36 |
| YMm55 | 0.60 | 0.50 |
| Tknn0 | 3283.32 | 3262.90 |
| Tknn1 | 3295.15 | 3279.57 |
| Tknn5 | 3262.69 | 3241.80 |
| ICA_Met_10 | 0.81 | 0.72 |
| ICA_Met_50 | 0.73 | 0.60 |
Met Modelling with: SVR. SNR=10
Met Modelling with: SVR. SNR=50
Metallicity Modelling with nnet. Error analysis follows.
| drmse | dmae | |
|---|---|---|
| chi2d_10 | 0.55 | 0.41 |
| chi2d_50 | 0.51 | 0.38 |
| Met_ICA_10 | 0.81 | 0.72 |
| Met_ICA_50 | 0.73 | 0.60 |
| YMm00 | 0.54 | 0.52 |
| YMm11 | 0.31 | 0.24 |
| YMm55 | 0.74 | 0.53 |
| Tknn0 | 3283.32 | 3262.90 |
| Tknn1 | 3295.15 | 3279.57 |
| Tknn5 | 3262.69 | 3241.80 |
| ICA_Met_10 | 0.81 | 0.72 |
| ICA_Met_50 | 0.73 | 0.60 |
Met Modelling with: NNR. SNR=10
Met Modelling with: NNR. SNR=50
Metallicity Modelling with knn. Error analysis follows.
| drmse | dmae | |
|---|---|---|
| chi2d_10 | 0.55 | 0.41 |
| chi2d_50 | 0.51 | 0.38 |
| Met_ICA_10 | 0.81 | 0.72 |
| Met_ICA_50 | 0.73 | 0.60 |
| YMm00 | 0.59 | 0.43 |
| YMm11 | 0.35 | 0.30 |
| YMm55 | 1.06 | 0.90 |
| Tknn0 | 3283.32 | 3262.90 |
| Tknn1 | 3295.15 | 3279.57 |
| Tknn5 | 3262.69 | 3241.80 |
| ICA_Met_10 | 0.81 | 0.72 |
| ICA_Met_50 | 0.73 | 0.60 |
Met Modelling with: KNN. SNR=oo
Met Modelling with: KNN. SNR=10
Met Modelling with: KNN. 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
Metallicity Modelling with bagEarth. Error analysis follows.
| drmse | dmae | |
|---|---|---|
| chi2d_10 | 0.55 | 0.41 |
| chi2d_50 | 0.51 | 0.38 |
| Met_ICA_10 | 0.81 | 0.72 |
| Met_ICA_50 | 0.73 | 0.60 |
| YMm00 | 0.86 | 0.70 |
| YMm11 | 0.68 | 0.55 |
| YMm55 | 0.75 | 0.64 |
| Tknn0 | 3283.32 | 3262.90 |
| Tknn1 | 3295.15 | 3279.57 |
| Tknn5 | 3262.69 | 3241.80 |
| ICA_Met_10 | 0.81 | 0.72 |
| ICA_Met_50 | 0.73 | 0.60 |
Met Modelling with: MARS. SNR=10
Met Modelling with: MARS. SNR=50
Metallicity Modelling with kernelpls. Error analysis follows.
| drmse | dmae | |
|---|---|---|
| chi2d_10 | 0.55 | 0.41 |
| chi2d_50 | 0.51 | 0.38 |
| Met_ICA_10 | 0.81 | 0.72 |
| Met_ICA_50 | 0.73 | 0.60 |
| YMm00 | 1.22 | 1.12 |
| YMm11 | 0.71 | 0.67 |
| YMm55 | 0.64 | 0.55 |
| Tknn0 | 3283.32 | 3262.90 |
| Tknn1 | 3295.15 | 3279.57 |
| Tknn5 | 3262.69 | 3241.80 |
| ICA_Met_10 | 0.81 | 0.72 |
| ICA_Met_50 | 0.73 | 0.60 |
Met Modelling with: PLS. SNR=10
Met Modelling with: PLS. SNR=50
Metallicity Modelling with cubist. Error analysis follows.
| drmse | dmae | |
|---|---|---|
| chi2d_10 | 0.55 | 0.41 |
| chi2d_50 | 0.51 | 0.38 |
| Met_ICA_10 | 0.81 | 0.72 |
| Met_ICA_50 | 0.73 | 0.60 |
| YMm00 | 1.20 | 1.11 |
| YMm11 | 0.45 | 0.37 |
| YMm55 | 0.52 | 0.42 |
| Tknn0 | 3283.32 | 3262.90 |
| Tknn1 | 3295.15 | 3279.57 |
| Tknn5 | 3262.69 | 3241.80 |
| ICA_Met_10 | 0.81 | 0.72 |
| ICA_Met_50 | 0.73 | 0.60 |
Met Modelling with: Rule-Regression. SNR=oo
Met Modelling with: Rule-Regression. SNR=10
Met Modelling with: Rule-Regression. SNR=50
#That’s all.
```