Normalizing synthetic BT-Settl spectra

Introduction

The main of this work is to produce reliable set of spectra valid to train a ML approach making it possible to estimate physical parameter for class M stars, i.e. according to the Hardvard stellar classification it includes stars with surface temperature between 2000 and 3500 K and they represent about the 76% of all the main sequence's stars see ref..

Study of GRID based syntetic spectra

The PHOENIX Grid has been used (called as BT-SETTL) and some of the spectra fit into the M class type. First step will be to process all the filenames in order to identify physical parameters. Then it will be possible to identify the interesting subset.

Main changes in connection with the previous release of this file (M_prep_rep_02v_01.Rmd) was to consider not each window as independent but to consider part of the I-band as a continuous region and to consider all the available spectra to select the normallization factor. Each spectrum is scaled to cover a unit area below it.

setwd("~/M")
files <- list.files(path = "./BT-SETTL/SPECTRA_2011/", pattern = ".*BT-Settl.*7.[gz|bz2]", 
    ignore.case = TRUE)
nf <- length(files)
body <- apply(as.data.frame(files), 1, function(x, ct) {
    a <- strsplit(x, ct)
    return(a[[1]][1])
}, "\\.BT")
pest <- apply(as.data.frame(body), 1, function(x, ct) {
    a <- strsplit(x, ct)
    return(substr(a[[1]][1], 4, nchar(a[[1]][1])))
}, "a")
tas <- as.numeric(apply(as.data.frame(pest), 1, function(x, ct, pos = 1) {
    a <- strsplit(x, ct)
    return(a[[1]][pos])
}, "[-+]", 1)) * 100
sgn1 <- as.numeric(apply(as.data.frame(pest), 1, function(x, ct, pos = 1) {
    a <- paste(unlist(strsplit(x, ct)), collapse = "")
    return(ifelse(substr(a, pos, pos) == "-", -1, 1))
}, "[0-9\\.]", 1))
sgn2 <- as.numeric(apply(as.data.frame(pest), 1, function(x, ct, pos = 1) {
    a <- paste(unlist(strsplit(x, ct)), collapse = "")
    return(ifelse(substr(a, pos, pos) == "-", -1, 1))
}, "[0-9\\.]", 2))
ges <- as.numeric(apply(as.data.frame(pest), 1, function(x, ct, pos = 1) {
    a <- strsplit(x, ct)
    return(a[[1]][pos])
}, "[-+]", 2)) * sgn1
met <- as.numeric(apply(as.data.frame(pest), 1, function(x, ct, pos = 1) {
    a <- strsplit(x, ct)
    return(a[[1]][pos])
}, "[-+]", 3)) * sgn2

gg <- as.data.frame(cbind(T = tas, LgG = ges, M = met, nf = files))
gg[, 1] <- as.numeric(as.character(gg[, 1]))
gg[, 2] <- as.numeric(as.character(gg[, 2]))
gg[, 3] <- as.numeric(as.character(gg[, 3]))
grd <- gg[gg$T < 3600 & gg$T > 1999, ]
rm(gg)
d_bt <- as.data.frame(matrix(0, nrow(grd), 14))
d_bt[, 1] <- grd[, 4]
d_bt[, 2:4] <- grd[, 1:3]

colnames(d_bt) <- c("FileName", "T[K]", "Log(g)[dex]", "Met [dex]")
tab01 <- xtable(d_bt[, 1:4])
print(tab01, "html")

Now we start to process the interested syntetic files. They are 346 files.

The first planned step will be to load compressed files to reorder the data and to take them out of the global files:

bt <- list()
if (file.exists("~/git/M_prep/bt-settl-2011-full.RData")) {
    load("~/git/M_prep/bt-settl-2011-full.RData")
} else {
    setwd("~/M/BT-SETTL/SPECTRA_2011/")
    for (i in 1:(nrow(grd))) {
        # for (i in 1:10) {
        ff <- as.character(grd[i, 4])
        if (grepl(".7.bz2", ff) > 0) {
            try(system(paste("bunzip2", as.character(grd[i, 4]), sep = " ")))
            j <- nchar(as.character(grd[i, 4])) - 4
        } else {
            try(system(paste("gunzip", as.character(grd[i, 4]), sep = " ")))
            j <- nchar(as.character(grd[i, 4])) - 3
        }
        f <- substr(as.character(grd[i, 4]), 1, j)
        dat <- read.fwf(file = f, c(13, 13))
        j <- cbind(str_trim(as.character(dat[[1]])), str_trim(as.character(dat[[2]])))
        dat <- str_replace(j, "D", "E")
        ss <- apply(dat, 2, as.numeric)
        ss[, 2] <- 10^ss[, 2]
        ss <- ss[order(ss[, 1]), ]
        try(system(paste("bzip2", f, sep = " ")))
        bt[[i]] <- list(data = ss, name = f)
        d_bt[i, 5:6] <- range(ss[, 2], na.rm = T)
        d_bt[i, 7:8] <- range(ss[, 1], na.rm = T)
        d_bt[i, 9:14] <- summary(diff(ss[, 1], 1))
    }
    rm(j, ss, dat, f, ff)
    save.image(file = "~/git/M_prep/bt-settl-2011-full.RData")
}
# d_bt[,12]<-tas d_bt[,13]<-ges d_bt[,14]<-met table04<-xtable(d_bt)
# print(table04,type='html') nfs[4]<-nf

As far as some diagrams are held several times in some files, a clean process is required, mostly to fix them and to have a real curve diagram.

Relevant literature paper adopted as guideline for some ideas in this work is: @article{Rayner:2009du, author = “Rayner, J.T. and Cushing, M.C. and Vacca, W.D.”, title = “{The Infrared Telescope Facility (IRTF) Spectral Library: Cool Stars}”, journal = “Astrophys.J.Suppl.”, volume = “185”, pages = “289-432”, doi = “10.1088/0067-0049/185/2/289”, year = “2009”, eprint = “0909.0818”, archivePrefix = “arXiv”, primaryClass = “astro-ph.SR”, SLACcitation = “%%CITATION = ARXIV:0909.0818;%%”, }

Some adjustments are going to be applied to the data like:

• Cutting into a window defined in the literature as most probable for getting answer from the spectra.
• Resampling the downgraded spectra to the resolution coming from observed spectra. In this particular case, as it is IPAC, the choice is 3.6A
• Consider relevant spectral windows for atoms exhibiting relevant activity as depicted in the literature

#
ordenar <- function(x) {
    x$data <- x$data[order(x$data[, 1]), ]
    return(x)
}
#
corte_bt <- function(x, lim = 1) {
    xx <- c(which(diff(x$data[, 1], 1) < 0), length(x$data[, 1]))
    if (lim == 0) {
        return(xx)
    } else {
        return(xx[1:lim])
    }
}
#
interp_seg <- function(x, filtro, paso = 1, cortes = list(), fnorm = 1, convolucion = TRUE) {
    y <- list()
    y$name <- x$name
    y$nch <- length(cortes)
    y$cuts <- cortes
    y$factors <- rep(0, y$nch)
    y$data <- list()
    y$err <- 0
    lag <- ifelse(convolucion, 4 * paso, 0)
    xmin <- which(x$data[, 1] >= (cortes[[1]][1] - lag))[1]
    xmax <- which(x$data[, 1] >= (cortes[[y$nch]][2] + lag))[1]
    if (is.na(xmin) || is.na(xmax)) {
        for (i in 1:length(cortes)) {
            nx <- seq(cortes[[i]][1], cortes[[i]][2], paso)
            ndat <- matrix(0, nrow = length(nx), ncol = 2)
            ndat[, 1] <- nx
            ndat[, 2] <- NA
            y$data[[i]] <- ndat
            y$factors[i] <- 0
        }
        y$err <- y$nch
        return(y)
    }
    if (convolucion && xmin > length(filtro)) {
        xmin <- xmin - length(filtro)
    }
    if (convolucion && xmax < length(x$data[, 1]) - length(filtro)) {
        xmax <- xmax + length(filtro)
    }
    x$data <- x$data[xmin:xmax, ]
    if (convolucion) {
        xtmp <- x$data[, 2] * fnorm
        ytmp <- conv(xtmp, filtro)
        ytmp <- ytmp[(floor(length(filtro)/2) + 1):(length(xtmp) + floor(length(filtro)/2))]
        x$data[, 2] <- ytmp
    }
    vec <- 3
    for (i in 1:length(cortes)) {
        xl <- cortes[[i]][1] - vec * paso
        xu <- cortes[[i]][2] + vec * paso
        dat <- as.matrix(x$data[(x$data[, 1] >= xl) & (x$data[, 1] <= xu), ])
        nx <- seq(cortes[[i]][1], cortes[[i]][2], paso)
        ndat <- matrix(0, nrow = length(nx), ncol = 2)
        ndat[, 1] <- nx
        for (k in 1:length(nx)) {
            l <- which(abs(dat[, 1] - nx[k]) == min(abs(dat[, 1] - nx[k])))
            xtmp <- dat[max(1, l - vec):(min(l + vec, dim(dat)[1])), ]
            if (sum(is.na(xtmp[, 2])) == 0) {
                ytmp <- spline(xtmp[, 1], xtmp[, 2], n = dim(xtmp)[1], xout = nx[k])
                ndat[k, 2] <- ytmp$y
            } else {
                ndat[k, 2] <- NA
            }
        }
        y$data[[i]] <- ndat
    }
    return(y)
}
#
`?` <- function(x, y) eval(sapply(strsplit(deparse(substitute(y)), ":"), function(e) parse(text = e))[[2 - 
    as.logical(x)]])
#
fnspec <- function(x, fn = "max") {
    FUN <- match.fun(fn)
    return(unlist(lapply(x$data, function(x) FUN(x[, 2]))))
}
# cortes<-list(c(8450,8680),c(8700,8900),c(8980,9270))
# =================================================================================
# Processing starts
# =================================================================================
if (file.exists("~/git/M_prep/M_prep_res_BT-Settl.RData")) {
    load("~/git/M_prep/M_prep_res_BT-Settl.RData")
} else {
    bi <- read.csv(file = "~/git/M_sel/banda_I.csv", header = TRUE, sep = ";", 
        skip = 2)
    bi <- bi[, -c(2, 5, 8)]
    colnames(bi) <- c("Element", "From", "To", "ContFrom_1", "ContTo_1", "ContFrom_2", 
        "ContTo_2")
    bi[, 8] <- 0
    bi[, 9] <- 1
    # WARNING Huge number of IPAC spectra end by 9152A, so there is no matter
    # to consider the Pa6 subband proposed in the paper ...
    bi <- bi[1:9, ]
    border <- list(range(bi[, c(2:3)]))
    cortes <- border
    #
    iir <- dnorm(seq(-50, 50, 1), 0, 3.6)
    bp <- lapply(bt, interp_seg, iir, paso = 3.6, cortes = cortes, convolucion = TRUE)
    #
    for (i in 1:length(bp)) {
        bp[[i]]$stellarp <- d_bt[i, 2:4]
    }
    save(d_bt, bp, border, cortes, iir, bi, file = "~/git/M_prep/M_prep_res_BT-Settl.RData")
}
#
tab03 <- xtable(bi[, 1:3])
print(tab03, "html")

Now it is time to noralize to unit area below the absortion curve in the window 8461, 9030 A. Then we will write out the spectra as they were prepaired for being used.

# Looking for spectra without NA as y
mettering <- function(x, blq = 1) {
    y <- x$data[[blq]]
    yv <- sum(is.na(y[, 2])) == 0
    st <- sum(diff(y[, 1]) * rollmean(y[, 2], 2))
    md <- mean(y[, 2])
    sv <- sd(y[, 2])
    lg <- nrow(y)
    return(list(ok = yv, area = st, mean = md, sdev = sv, lngth = lg))
}
#
scspec <- function(x, area = 1) {
    lrt <- mettering(x, 1)
    for (i in x$nch) {
        x$factors[i] <- ifelse(x$factors[i] == 0, area/lrt$area, x$factors[i] * 
            area/lrt$area)
        x$data[[i]][, 2] <- x$data[[i]][, 2] * area/lrt$area
    }
    return(x)
}
#
prepare <- function(x, chunk = 1) {
    tmp <- t(x$data[[chunk]])
    colnames(tmp) <- paste("wl_", tmp[1, ], sep = "")
    dtmp <- tmp[2, ]
    if ("stellarp" %in% (attributes(x)$names)) {
        dtmp <- cbind(t(dtmp), x$stellarp)
    }
    return(dtmp)
}
#
# ===========================================================================
# Processing
# ===========================================================================
if (file.exists("~/git/M_prep/M_prep_clean_BT-Settl.RData")) {
    load("~/git/M_prep/M_prep_clean_BT-Settl.RData")
} else {
    ltp <- unlist(lapply(bp, function(x) sum(is.na(x$data[[1]][, 2])) == 0))
    bp_clean <- bp[ltp]
    m_bp <- lapply(bp_clean, mettering)
    bp_clean <- lapply(bp_clean, scspec, area = 1)
    save.image(file = "~/git/M_prep/M_prep_clean_BT-Settl.RData")
}
# Export data for training classical ICA models.
require(plyr)
df_data <- ldply(bp_clean, prepare, 1)
rownames(df_data) <- d_bt[, 1]
setwd("~/git/M_prep/")
write.table(df_data, "train_x_01.txt")
#

Processing IPAC dataset

Now it's time to apply the same procedeure to the test data from IPAC library

if (file.exists("~/git/M_prep/M_prep_cleanip_BT-Settl.RData")) {
    load("~/git/M_prep/M_prep_cleanip_BT-Settl.RData")
} else {
    setwd("~/M")
    files <- list.files(path = "./ipac", pattern = ".*txt$", ignore.case = TRUE)
    nf <- length(files)
    d_ipac <- as.data.frame(matrix(0, nf, 10))
    d_ipac[, 1] <- files
    colnames(d_ipac) <- c("FileName", "MinFL", "MaxFL", "MinWL", "MaxWL", "MinST", 
        "1Q_ST", "Mean", "3Q_ST", "Max_ST")
    #
    ipc <- list()
    for (i in files) {
        ss <- read.table(file = paste("ipac/", i, sep = ""), skip = 16, header = FALSE)
        ipc[[i]] <- list(data = ss, name = i)
    }
    d_ipac[, 6:10] <- ldply(ipc, function(x) {
        return(fivenum(x$data[, 1]))
    })[, -1]
    # Selection for valid spectra let's remove the spectra with NA
    bq <- lapply(ipc, interp_seg, iir, paso = 3.6, cortes = border, convolucion = FALSE)
    ltq <- unlist(lapply(bq, function(x) sum(is.na(x$data[[1]][, 2])) == 0))
    bq_clean <- bq[ltq]
    m_bq <- lapply(bq_clean, mettering)
    bq_clean <- lapply(bq_clean, scspec, area = 1)
    save.image(file = "~/git/M_prep/M_prep_cleanip_BT-Settl.RData")
}
tab02 <- xtable(d_ipac[, c(1, 6, 10)])
print(tab02, "html")

#
df_data <- t(unlist(sapply(bq_clean, prepare, 1)))
rownames(df_data) <- d_ipac[ltq, 1]
write.table(df_data, "test_x_01.txt")
#
save.image(file = "M_prep_BT-Settl.RData")

This finishes the data preparation step. It considers 346 training patterns from syntetic spectra and 595 to be interpreted inside the window 8461, 9030.