неверный тип индекса «ошибка» в коде R

Question

Я пытаюсь запустить мой код R, который должен прочитать несколько файлов текстовых документов и вывести график, аналогичный приведенному в ссылке ниже.

https://imgur.com/IOhX42L

Тем не менее, я продолжаю получать следующую ошибку:

Error in log10(1 - DF$q)[sapply(xticks, function(x) which.min(abs(x -  : 
  invalid subscript type 'list'

Я понимаю, что пользователи здесь ранее имели ту же проблему с этой проблемой, но мне было интересно, если кто-то мог взглянуть на мой укажите c код и предложите, в чем проблема может быть здесь

### LIBRARIES 
source('utils.R') 
options(scipen=10000)
library(plotrix)

logs_dir <- 'LOGS_3MS' 

prof.idx <- read.table(file.path(logs_dir, 'profiles.index'), skip=1, 
    col.names=c('mdl_num', 'priority', 'prof_num'))

hstry <- read.table(file.path(logs_dir, 'history.data'), header=1, skip=5)

decreasing_L <- with(hstry, 
    which((diff(log_L) < 0)
        & center_h1[-1] > 0.65))
if (any(decreasing_L)) {
    goes_back_up <- diff(decreasing_L) > 1
    pms <- max(decreasing_L)
    hstry <- hstry[-1:-pms,]
}
decreasing_Teff <- with(hstry, 
    which((diff(log_Teff) < 0)
        & center_h1[-1] > 0.65))
if (any(decreasing_Teff)) {
    goes_back_up <- diff(decreasing_Teff) > 1
    pms <- max(decreasing_Teff)
    hstry <- hstry[-1:-pms,]
}

prof_num <- 1360

for (prof_num in c(1360, 3000)) {#prof.idx$prof_num) { 

DF <- read.table(file.path(logs_dir, paste0('profile', prof_num, '.data')), 
    header=1, skip=5)
hstry. <- hstry[hstry$model_number == 
    prof.idx[prof.idx$prof_num == prof_num,]$mdl_num,]

plot_hif <- function(..., make.x=T, make.y=T,
        text.cex=1, mgp=utils.mgp, font=utils.font, mar=utils.mar, short=F) {

    par(mar=mar+c(0.3, -0.5, 2, -0.1), lwd=1.66, las=1, cex.axis=text.cex,
        mfrow=c(1,2))

    ### HRD 

    xlim <- log10(c(8700, 3300)) 
    ylim <- c(-0.5, 4.6)

    plot(NA, 
        axes=F, xaxs='i', yaxs='i', 
        type='l', lwd=3, col=1, 
        xlim=xlim, 
        ylim=ylim, 
        xlab="", 
        ylab="")

    spectral.divs <- log10(c(30000, 10000, 7500, 6000, 5200, 3700, 2400))
    rs <- c(175/255, 199/255, 1, 1, 1, 1, 1, 1)
    gs <- c(201, 216, 244, 229, 217, 199, 166)/255
    bs <- c(1, 1, 243/255, 207/255, 178/255, 142/255, 81/255)
    cols <- c(
        rgb(175/255, 201/255, 1),       # O
        rgb(199/255, 216/255, 1),       # B
        rgb(1,       244/255, 243/255), # A 
        rgb(1,       229/255, 207/255), # F 
        rgb(1,       217/255, 178/255), # G 
        rgb(1,       199/255, 142/255), # K 
        rgb(1,       166/255, 81/255))  # M
    for (ii in 1:length(spectral.divs)) {
        div <- spectral.divs[ii]
        if (div > xlim[1]) next  
        if (div < xlim[2]) div <- xlim[2]
        if (ii == 1) {
            rect(xlim[1], ylim[1], div, ylim[2], col=cols[ii], border=NA)
        } else {
            prev <- spectral.divs[ii-1]
            if (prev > xlim[1]) prev <- xlim[1]
            rect(prev, ylim[1], div, ylim[2], col=cols[ii], border=NA)
        }
    }
    for (ii in 2:(length(spectral.divs)-1)) {
        div <- spectral.divs[ii]
        if (div > xlim[1]) next  
        if (div < xlim[2]) div <- xlim[2]
        gradient.rect(div+0.0025, ylim[1], div-0.0025, ylim[2],
            nslices=10, border=NA, 
            reds=c(rs[ii], rs[ii+1]), 
            greens=c(gs[ii], gs[ii+1]),
            blues=c(bs[ii], bs[ii+1]))
    }

    lines(hstry$log_Teff, hstry$log_L, type='l', lwd=2, col=1)
    points(hstry$log_Teff[1], hstry$log_L[1], pch=20, cex=0.5, col=1)

    mdl.col <- cols[which.min(hstry.$log_Teff < spectral.divs)]
    mdl.cex <- 1+hstry.$log_R
    points(hstry.$log_Teff, hstry.$log_L, pch=21, cex=mdl.cex, 
        col="#FFFFFF", lwd=par()$lwd,
        bg=with(DF[nrow(DF),], 
            rgb(red=y*.8, green=x*.8, blue=z*.8)))

    age <- round(signif(hstry.$star_age/10**9, 4), 5)
    mass <- signif(hstry.$star_mass, 3)
    par(family="Helvetica LT Std Light")
    legend('topleft', 
        cex=0.8*text.cex,
        bty='n', inset=c(-0.05, 0),
        legend=c(as.expression(bquote(tau/Gyr==.(age))),
                 as.expression(bquote(M/M["sun"]==.(mass)))))
    par(family=font)

    nxticks <- 3
    nyticks <- 4
    nxminor <- 5
    nyminor <- 4
    xticks <- pretty(xlim, n=nxticks)
    yticks <- pretty(ylim, n=nyticks)
    xticks.minor <- pretty(xlim, n=nxticks*nxminor)
    yticks.minor <- pretty(ylim, n=nyticks*nyminor)
    xticks.minor <- xticks.minor[!xticks.minor %in% xticks]
    yticks.minor <- yticks.minor[!yticks.minor %in% yticks]
    par(mgp=mgp+c(0, 0.3, 0))
    xpos <- seq(10**xlim[2], 10**xlim[1], 2000)
    xpos2 <- seq(10**xlim[2], 10**xlim[1], 400)
    axis(side=1, tcl=-0.346/2, at=log10(xpos2), 
        labels=F, lwd.ticks=par()$lwd)
    axis(side=1, tcl=-0.346, at=log10(xpos), 
        labels=xpos, cex.axis=text.cex, lwd.ticks=par()$lwd)
    par(mgp=mgp+c(0, 0.43, 0))

    magaxis(2, tcl=-0.346, lwd=0, lwd.ticks=par()$lwd, unlog=T, family=font,
        mgp=mgp+c(0, 0.43, 0))
    #axis(2, tcl=-0.346, lwd=0, lwd.ticks=par()$lwd, tick=T, at=yticks,
    #    labels=as.logical(make.y))
    #axis(2, tcl=-0.346/2, lwd=0, lwd.ticks=par()$lwd, tick=T, 
    #    at=yticks.minor, labels=F)

    box(lwd=par()$lwd)

    par(mgp=mgp+c(0.55, 0, 0))
    if (make.x) title(xlab=expression(T["eff"]/K))
    par(mgp=mgp+c(0.55, 0, 0))
    if (make.y) title(ylab=expression(L/L["sun"]))
    if (make.x) mtext(expression("Spectral Type"), 
        side=3, line=1.5, cex=text.cex)

    spectral.labs <- c("O", "B", "A", "F", "G", "K", "M")
    selector <- 1:(length(spectral.divs))
    spectral.Teffs <- sapply(selector, 
        function(ii) {
            div <- spectral.divs[ii]
            if (div > xlim[1]) return(Inf)
            if (div < xlim[2]) div <- xlim[2]
            if (ii == 1) return((xlim[1]+div)/2)
            prev <- spectral.divs[ii-1]
            if (prev > xlim[1]) prev <- xlim[1]
            (div+prev)/2
        })
    axis(3, at=spectral.divs, tcl=-0.346, labels=F, cex.axis=text.cex,
        lwd.ticks=par()$lwd)
    par(mgp=mgp+c(0, -0.1, 0))
    axis(3, at=spectral.Teffs, labels=spectral.labs[selector], 
        cex.axis=text.cex, tcl=0)


    ### ENVELOPE STRUCTURE 

    xlim <- c(-10, 0)
    ylim <- c(0, 12)

    plot(NA, axes=F, 
        xaxs='i',  yaxs='i', 
        xlim=xlim, ylim=ylim,
        xlab="",   ylab="")

    lines(log10(1-DF$q), DF$neutral_fraction_H, lwd=3)
    lines(log10(1-DF$q), DF$temperature/10**4, lwd=3, col=blue)
    lines(log10(1-DF$q), DF$opacity/10**1, lwd=3, col=orange)

    nxticks <- 6
    nyticks <- 5
    nxminor <- 4
    nyminor <- 4
    xticks <- pretty(xlim, n=nxticks)
    yticks <- pretty(ylim, n=nyticks)
    xticks.minor <- pretty(xlim, n=nxticks*nxminor)
    yticks.minor <- pretty(ylim, n=nyticks*nyminor)
    xticks.minor <- xticks.minor[!xticks.minor %in% xticks]
    yticks.minor <- yticks.minor[!yticks.minor %in% yticks]

    par(mgp=mgp+c(0, 0.3, 0))
    axis(1, tcl=-0.346, lwd=0, lwd.ticks=par()$lwd, tick=T, at=xticks,
        labels=as.logical(make.x))
    axis(1, tcl=-0.346/2, lwd=0, lwd.ticks=par()$lwd, tick=T, 
        at=xticks.minor, labels=F)

    par(mgp=mgp+c(0, 0.43, 0))
    axis(2, tcl=-0.346, lwd=0, lwd.ticks=par()$lwd, tick=T, at=yticks,
        labels=as.logical(make.y))
    axis(2, tcl=-0.346/2, lwd=0, lwd.ticks=par()$lwd, tick=T, 
        at=yticks.minor, labels=F)

    par(mgp=mgp+c(0, 0.3, 0))
    axis(3, tcl=-0.346, lwd=0, lwd.ticks=par()$lwd, tick=T, 
        at=log10(1-DF$q)[sapply(xticks, function(x) 
            which.min(abs(x-log10(1-DF$radius/10**hstry.$log_R))))],
        labels=xticks)
    axis(3, tcl=-0.346/2, lwd=0, lwd.ticks=par()$lwd, tick=T, 
        at=log10(1-DF$q)[sapply(xticks.minor, function(x) 
            which.min(abs(x-log10(1-DF$radius/10**hstry.$log_R))))], 
        labels=F)

    par(family="Helvetica LT Std Light")
    leg.spot <- if (prof_num == 3000) "topleft" else "bottomright"
    inset <- if (prof_num == 3000) c(0, 0) else c(0, 0.08)
    legend(leg.spot, inset=inset, 
        col=c(orange, blue, 1), lwd=3, cex=0.8*text.cex,
        bty='n', lty=c(1,1,1),
        legend=c(expression(Opacity~kappa/(10~cm^2/g)),
                 expression(Temp.~T/(10^4~K)),
                 expression(Neutral~fraction~H)
            ))
    par(family=font)

    box(lwd=par()$lwd)

    par(mgp=mgp+c(0.55, 0, 0))
    if (make.x) title(xlab=expression(Fractional~mass~log[10](1-m/M)))

    par(mgp=mgp+c(0.55, 0, 0))
    if (make.y) title(ylab=expression(Envelope~structure)) 
    if (make.x) mtext(expression(Fractional~radius~log[10](1-r/R)), 
        side=3, line=1.5, cex=text.cex)
}

#plot_hif()

make_plots(plot_hif, paste0('hif', sprintf("%06d", prof_num)),
        filepath=file.path('plots'),
        cex.paper=0.93, 
        wide=T, thin=F, short=F, slides=F, make_png=T, #make_pdf=F,
        make.x=T,
        make.y=T,
        font="Palatino Linotype", 
        use.cairo=T)

}




### LIBRARIES 
#source('../scripts/utils.R') 
#options(scipen=10000)

logs_dir <- 'LOGS_3MS'
prof.idx <- read.table(file.path(logs_dir, 'profiles.index'), skip=1, 
    col.names=c('mdl_num', 'priority', 'prof_num'))
hstry <- read.table(file.path(logs_dir, 'history.data'), header=1, skip=5)



taus <- c()
ages <- c()
logqs <- c()
for (prof_num in prof.idx$prof_num) { 

    prof.path <- file.path(logs_dir, paste0('profile', prof_num, '.data'))
    if (!file.exists(prof.path)) next 
    print(prof.path)
    DF <- read.table(prof.path, header=1, skip=5)
    hstry. <- hstry[hstry$model_number == 
        prof.idx[prof.idx$prof_num == prof_num,]$mdl_num,]

    hif.idx <- min(which(DF$neutral_fraction_H <= 0.5))

    logq <- log10(1-DF$mass[hif.idx] / hstry.$star_mass)
    tau <- DF$opacity[hif.idx]

    ages <- c(ages, hstry.$star_age / 10**9)
    logqs <- c(logqs, logq)
    taus <- c(taus, tau)

}


tau.DF <- data.frame(ages, taus, logqs)
write.table(tau.DF, file="taus.dat")
tau.DF <- read.table('taus.dat', header=1)

options(scipen=1)
plot_hif_tau <- function(..., make.x=T, make.y=T, 
        text.cex=1, mgp=utils.mgp, font=utils.font, mar=utils.mar, short=F) {

    par(mar=mar+c(2, 1, 2, 0.5), lwd=3.66, las=1, cex.axis=text.cex)

    xlim <- c(10, 12.5)
    #xlim <- c(0, max(hstry$star_age)/10**9)
    ylim <- c(0.1, 1e6) #log10(c(0.1, 1e7))

    plot(NA, 
        axes=F, xaxs='i', yaxs='i', log='y', 
        type='l', lwd=3, col=1, 
        xlim=xlim, 
        ylim=ylim, 
        xlab="", 
        ylab="")

    ## EVOLUTIONARY STAGE 
    TAMS <- with(hstry, star_age[min(which(center_h1 <= 1e-5))]/10**9)
    #abline(v=TAMS, lty=2, lwd=par()$lwd)

    rgb <- with(hstry, center_h1 <= 1e-5 & center_he4 > 1e-5)
    RGB.TIP <- with(hstry[rgb,], star_age[which.max(log_L)]/10**9)
    #abline(v=RGB.TIP, lty=2, lwd=par()$lwd)

    AGB <- with(hstry, 
        star_age[min(which(center_h1 <= 1e-5 & center_he4 <= 1e-5))]/10**9)
    #abline(v=AGB, lty=2, lwd=par()$lwd)


    spectral.divs <- c(TAMS, RGB.TIP, AGB, xlim[2])
    #rs <- c(175/255, 199/255, 1, 1, 1, 1, 1, 1)
    #gs <- c(201, 216, 244, 229, 217, 199, 166)/255
    #bs <- c(1, 1, 243/255, 207/255, 178/255, 142/255, 81/255)
    rs <- c(175/255, 199/255, 1, 175/255, 1, 1, 1, 1)
    gs <- c(201, 216, 244, 201, 217, 199, 166)/255
    bs <- c(1, 1, 243/255, 1, 178/255, 142/255, 81/255)
    cols <- c(
        rgb(175/255, 201/255, 1),       # O 1
        rgb(199/255, 216/255, 1),       # B 2
        rgb(1,       244/255, 243/255), # A 3
        rgb(175/255, 201/255, 1),       # O 1#rgb(1,       229/255, 207/255), # F 4
        rgb(1,       217/255, 178/255), # G 5
        rgb(1,       199/255, 142/255), # K 6
        rgb(1,       166/255, 81/255))  # M 7
    #idxs <- c(4, 6, 2, 5)
    #rs <- rs[idxs]
    #gs <- gs[idxs]
    #bs <- bs[idxs]
    #cols <- cols[idxs]
    #cols <- c("#E8E9EB", "#E4B363", "#EF6461", "#E0DFD5")
    for (ii in 1:length(spectral.divs)) {
        div <- spectral.divs[ii]
        if (div < xlim[1]) next 
        if (div > xlim[2]) div <- xlim[2]
        if (ii == 1) {
            rect(xlim[1], ylim[1], div, ylim[2], col=cols[ii], border=NA)
        } else {
            prev <- spectral.divs[ii-1]