2

I have a pretty ugly looking script that I use to process evaluation files from simulations, it looks horrible and no I am not a coder, though it usually works, however it is currently not.

Just to clarify, the script normally iterates over multiple input files and it does works on my mac and on the cluster where I run simulations. I am now trying to run it on a VPS running Ubuntu server and it produces some strange output. I cannot figure out how to fix this.

This is the script in full:

#!/usr/bin/awk -f
FNR==1 && NR!=1 { endfile(); avgLT=totFrames=avgLTsq=avgFramessq=denom=0 }
FNR==1 { out1="analLT_"FILENAME; out2="sumLT_"FILENAME; out3="reportLT.txt"; print "-> Input file is: "FILENAME >> out3; next
       }
FNR==1 { next }

{
   avgLT+=$4; totFrames+=$5; ++denom;
   printf "%10.4f %10.1f\n",$4,$5 > out1
  }

END { endfile() }
function endfile()
{
  x="\nNO DATA POINTS IN INPUT => NO HYDROGEN BONDS DETECTED!"
  if (avgLT==0 && denom==0) {
    print x > out1; print x > out2; print x"\n\n----------------------------------------\n" >> out3;
    close(out1); close(out2); close(out3); return
  }
  if (avgLT>0) {
    avgAvgLT=avgLT/denom
    avgFrames=totFrames/denom
    while ((getline<out1)>0) {
      avgLTsq+=(($1-avgAvgLT)^2)
      avgFramessq+=(($2-avgFrames)^2)
    }
  close(out1)
    printf "\n   Summary data for hbond lifetime analysis:\n\n" > out2
    printf "   Summed Avg Lifetime:    %10.4f\n",avgLT > out2
    printf "   Average Lifetime:       %10.4f\n",avgAvgLT > out2
    printf "      Summed Frames:  %10.0f\n",totFrames > out2
    printf "      Average Frames:      %10.4f\n",avgFrames > out2
    printf "\n   Summary data for hbond lifetime analysis:\n\n" >> out3
    printf "   Summed Avg Lifetime:    %10.4f\n",avgLT >> out3
    printf "   Average Lifetime:       %10.4f\n",avgAvgLT >> out3
    printf "      Summed Frames:  %10.0f\n",totFrames >> out3
    printf "      Average Frames:      %10.4f\n",avgFrames >> out3

    if (denom>1) {
      sd_avgLT=sqrt(avgLTsq/(denom-1)); semAvgLT=(sd_avgLT/(sqrt(denom))); sd_totFrames=sqrt(avgFramessq/(denom-1)); semTotFrames=(sd_totFrames/(sqrt(denom)))
      printf "\n   SD lifetime:            %10.4f\n",sd_avgLT > out2
      printf "   SEM lifetime:           %10.4f\n",semAvgLT > out2
      printf "      SD Frames:           %10.4f\n",sd_totFrames > out2
      printf "      SEM Frames:          %10.4f\n\n",semTotFrames > out2
      printf "\n   SD lifetime:            %10.4f\n",sd_avgLT >> out3
      printf "   SEM lifetime:           %10.4f\n",semAvgLT >> out3
      printf "      SD Frames:           %10.4f\n",sd_totFrames > out3
      printf "      SEM Frames:          %10.4f\n\n",semTotFrames > out3
    } if (denom>1 && denom!=2) {print "----------------------------------------\n" >> out3 }
      if (denom==1) { print "   Single HBOND event, no SD or SEM calculation possible!" > out2;
             print "\n   Single HBOND event, no SD or SEM calculation possible!\n\n----------------------------------------\n" >> out3
           }
      if (denom==2) { print "\n   2 Hydrogen bond events found! No proper SD or SEM!" > out2;
             print "   2 Hydrogen bond events found! No proper SD or SEM!\n\n----------------------------------------\n" >> out3
           }
}
  close(out3)
  close(out2)
}

It takes a 5 column input file, works on 2 columns and puts the same columns in a separate file for later processing (out1). This file should then be processed to calculate some statistics though this does not happen on the VPS, all I get is 0.0000 values.

The problem seems to be in the while function somehow:

while ((getline<out1)>0) {
      avgLTsq+=(($1-avgAvgLT)^2)
      avgFramessq+=(($2-avgFrames)^2)
    }

As by the end of the script, when things get printed to files, it seems like I am getting reasonable values for the Sums and Averages calculated (avgLT, avgAvgLT, totFrames and avgFrames). When it gets to the Statistics part (sd_avgLT, semAvgLT, sd_totFrames and semTotFrames), all of these are printed to both out2 and out3 though all values are 0.0000 and not what they should be.

The "math" seems to be OK as running the commands separately on the out1 file:

awk ' BEGIN { avgAvgLT=1.4264 } { avgLTsq+=(($1-avgAvgLT)^2) } END { print avgLTsq }' analLT_multiple.out
awk ' BEGIN { avgFrames=4.4831 } { avgFramessq+=(($2-avgFrames)^2) } END { print avgFramessq }' analLT_multiple.out
awk ' BEGIN { avgLTsq=30.3478; denom=89 } { sd_avgLT=sqrt(avgLTsq/(denom-1)) } END { print sd_avgLT }' analLT_multiple.out
awk ' BEGIN { sd_avgLT=0.587249; denom=89 } {semAvgLT=(sd_avgLT/(sqrt(denom))) }  END { print semAvgLT }' analLT_multiple.out
awk ' BEGIN { avgFramessq=2040.22; denom=89 } { sd_totFrames=sqrt(avgFramessq/(denom-1)) } END { print sd_totFrames }' analLT_multiple.out
awk ' BEGIN { sd_totFrames=4.81501; denom=89 } { semTotFrames=(sd_totFrames/(sqrt(denom))) } END { print semTotFrames }' analLT_multiple.out

Gives me non-zero values and these values seem reasonable, however the script gives me all 0.0000 values for all. I also tried printing out the values of the variables in the script while running over multiple files and the denom variable works though sd_avgLT, semAvgLT, sd_totFrames and semTotFrames are all returned zero or empty.

My "conclusion" (I would say guess here) is, as stated that there is something wrong with the while function though I don't see what.

I put a sample input file up on pastebin https://pastebin.com/JsuTz0mD if you want to try to run the script yourself.

Any input/feedback or solutions that would make this script run on my VPS system would be enormously appreciated

2 Answers 2

4

It's either that the awk you are using (GNU awk or mawk) does not flush the data written to to the out1 file as you write it, or that the awk reads nothing from a file handle that you are keeping open for writing. This means that when you then read from that file in the END block, there is no data read. The BSD awk implementations don't seem to have this issue, and your code works as expected on e.g. OpenBSD and macOS.

The solution is simple, use close(out1) unconditionally in the END block before reading from it with getline. Currently, you close it after reading from it.

Also, consider being a bit more consistent with your > and >>. In this code, I believe that you can use > throughout.

2
  • Thank you so much, this would have taken me days to figure out by pure trial and error. Yes, moving the close(out1) before getline solved the issue on the VPS installation. Should work equally well on macOS. I'll test it later. My thought regarding the directs of output was that while iterating over files, using > would not append to the report file so I get one document for all files as intended. However, I could obviously make all the "prints" cleaner and I will look into that. Again, a super big thank you!
    – gusgrave
    Jun 19, 2019 at 6:59
  • 1
    @gusgrave Using > would be fine. It would truncate the output file to zero size if you close it and then write to it again. I don't think you are doing that. awk keeps filehandles open until you close them or the program exits. Doing multiple prints with > is therefore not an issue unless, as I mentioned, explicitly close the file and then open it for writing again.
    – Kusalananda
    Jun 19, 2019 at 7:10
1

Not an answer as @Kusalananda already told you what's wrong but let's tidy your script up a bit to make it readable and have a bit less duplication of code:

FNR == 1 {
    if ( NR != 1 ) {
        endfile()
    }
    avgLT = totFrames = denom = 0
    out1 = "analLT_" FILENAME
    out2 = "sumLT_" FILENAME
    out3 = "reportLT.txt"
    print "-> Input file is: " FILENAME > out3
    next
}

{
    avgLT += $4
    totFrames += $5
    ++denom
    printf "%10.4f %10.1f\n", $4, $5 > out1
}

END {
    endfile()
}

function endfile(       x, avgAvgLT, avgFrames, sd_avgLT,
                        semAvgLT, sd_totFrames, semTotFrames )
{
    if (avgLT == 0 && denom == 0 ) {
        x = "\nNO DATA POINTS IN INPUT => NO HYDROGEN BONDS DETECTED!"
        print x         > out1
        print x         > out2
        print x         > out3
    }
    else if (avgLT > 0) {
        avgAvgLT = avgLT / denom
        avgFrames = totFrames / denom

        close(out1)
        while ((getline < out1) > 0) {
            avgLTsq     += (($1 - avgAvgLT) ^ 2)
            avgFramessq += (($2 - avgFrames) ^ 2)
        }
        close(out1)

        printf "\n   Summary data for hbond lifetime analysis:\n\n"             > out2
        printf "   Summed Avg Lifetime:    %10.4f\n", avgLT                     > out2
        printf "   Average Lifetime:       %10.4f\n", avgAvgLT                  > out2
        printf "      Summed Frames:  %10.0f\n", totFrames                      > out2
        printf "      Average Frames:      %10.4f\n", avgFrames                 > out2

        printf "\n   Summary data for hbond lifetime analysis:\n\n"             > out3
        printf "   Summed Avg Lifetime:    %10.4f\n", avgLT                     > out3
        printf "   Average Lifetime:       %10.4f\n", avgAvgLT                  > out3
        printf "      Summed Frames:  %10.0f\n", totFrames                      > out3
        printf "      Average Frames:      %10.4f\n", avgFrames                 > out3

        if (denom == 1) {
            x = "   Single HBOND event, no SD or SEM calculation possible!"
            print x     > out2
            print ""    > out3
            print x     > out3
        }
        else if (denom > 1) {
            sd_avgLT = sqrt(avgLTsq / (denom - 1))
            semAvgLT = (sd_avgLT / (sqrt(denom)))
            sd_totFrames = sqrt(avgFramessq / (denom - 1))
            semTotFrames = (sd_totFrames / (sqrt(denom)))

            printf "\n   SD lifetime:            %10.4f\n", sd_avgLT            > out2
            printf "   SEM lifetime:           %10.4f\n", semAvgLT              > out2
            printf "      SD Frames:           %10.4f\n", sd_totFrames          > out2
            printf "      SEM Frames:          %10.4f\n\n", semTotFrames        > out2

            printf "\n   SD lifetime:            %10.4f\n", sd_avgLT            > out3
            printf "   SEM lifetime:           %10.4f\n", semAvgLT              > out3
            printf "      SD Frames:           %10.4f\n", sd_totFrames          > out3
            printf "      SEM Frames:          %10.4f\n\n", semTotFrames        > out3

            if (denom == 2) {
                x = "   2 Hydrogen bond events found! No proper SD or SEM!"
                print ""        > out2
                print x         > out2
                print x         > out3
        }
    }

    print "\n\n----------------------------------------\n"                      > out3

    close(out1)
    close(out2)
    close(out3)
}

That while getline loop on out1 obviously isn't really necessary since you could just store the data in an array instead of writing to out1 in the main body of the script, e.g.:

FNR == 1 {
    if ( NR != 1 ) {
        endfile()
    }
    avgLT = totFrames =  denom = 0
    out1 = "analLT_" FILENAME
    out2 = "sumLT_" FILENAME
    out3 = "reportLT.txt"
    print "-> Input file is: " FILENAME > out3
    next
}

{
    avgLT += $4
    totFrames += $5
    ++denom
    fnr2avgLT[FNR] = avgLT
    fnr2totFrames[FNR] = totFrames
}

END {
    endfile()
}

function endfile(       i, x, avgAvgLT, avgFrames, sd_avgLT,
                        semAvgLT, sd_totFrames, semTotFrames )
{
    if (avgLT == 0 && denom == 0 ) {
        x = "\nNO DATA POINTS IN INPUT => NO HYDROGEN BONDS DETECTED!"
        print x         > out1
        print x         > out2
        print x         > out3
    }
    else if (avgLT > 0) {
        avgAvgLT = avgLT / denom
        avgFrames = totFrames / denom

        for (i=1; i<=FNR; i++) {
            avgLT = fnr2avgLT[i]
            totFrames = fnr2totFrames[i]
            printf "%10.4f %10.1f\n", avgLT, totFrames > out1

            avgLTsq     += ((avgLT - avgAvgLT) ^ 2)
            avgFramessq += ((totFrames - avgFrames) ^ 2)
        }

        printf "\n   Summary data for hbond lifetime analysis:\n\n"             > out2
        printf "   Summed Avg Lifetime:    %10.4f\n", avgLT                     > out2
        printf "   Average Lifetime:       %10.4f\n", avgAvgLT                  > out2
        printf "      Summed Frames:  %10.0f\n", totFrames                      > out2
        printf "      Average Frames:      %10.4f\n", avgFrames                 > out2

        printf "\n   Summary data for hbond lifetime analysis:\n\n"             > out3
        printf "   Summed Avg Lifetime:    %10.4f\n", avgLT                     > out3
        printf "   Average Lifetime:       %10.4f\n", avgAvgLT                  > out3
        printf "      Summed Frames:  %10.0f\n", totFrames                      > out3
        printf "      Average Frames:      %10.4f\n", avgFrames                 > out3

        if (denom == 1) {
            x = "   Single HBOND event, no SD or SEM calculation possible!"
            print x     > out2
            print ""    > out3
            print x     > out3
        }
        else if (denom > 1) {
            sd_avgLT = sqrt(avgLTsq / (denom - 1))
            semAvgLT = (sd_avgLT / (sqrt(denom)))
            sd_totFrames = sqrt(avgFramessq / (denom - 1))
            semTotFrames = (sd_totFrames / (sqrt(denom)))

            printf "\n   SD lifetime:            %10.4f\n", sd_avgLT            > out2
            printf "   SEM lifetime:           %10.4f\n", semAvgLT              > out2
            printf "      SD Frames:           %10.4f\n", sd_totFrames          > out2
            printf "      SEM Frames:          %10.4f\n\n", semTotFrames        > out2

            printf "\n   SD lifetime:            %10.4f\n", sd_avgLT            > out3
            printf "   SEM lifetime:           %10.4f\n", semAvgLT              > out3
            printf "      SD Frames:           %10.4f\n", sd_totFrames          > out3
            printf "      SEM Frames:          %10.4f\n\n", semTotFrames        > out3

            if (denom == 2) {
                x = "   2 Hydrogen bond events found! No proper SD or SEM!"
                print ""        > out2
                print x         > out2
                print x         > out3
        }
    }

    print "\n\n----------------------------------------\n"                      > out3

    close(out1)
    close(out2)
    close(out3)
}

All of the above is untested of course since you didn't provide any sample input/output for us to test against but hopefully any bugs will be easy to spot and rectify.

7
  • Also note that GNU awk has an actual ENDFILE block.
    – Kusalananda
    Jun 21, 2019 at 8:26
  • 1
    @gusgrave you're welcome. Puting anything on some other site and expecting people trying to help you to be willing to go clicking around those other sites is not useful, you're expected to provide concise, testable sample input and expected output as text in your question when asking questions on this forum. When people put something in some external site it's invariably because they couldn't be bothered to come up with a minimal example to demonstrate their specific problem and want/expect everyone else to wade through some huge amount of data trying to understand it instead.
    – Ed Morton
    Jun 22, 2019 at 15:38
  • 1
    Wrt gender - I'm specifically not making the assumption that the person being referred to is male, that's all. When referring to a person whose gender you don't know "they" is contrived/awkward, "(s)he" is just as bad, and "it" is frowned upon ( :-) ) so that leaves us with "he" or "she" and too often people use "he". No-one has ever mentioned it when I use "he" but I frequently get complaints/questions/comments when I use "she" which I find interesting and encourages me to do it more to shake the stereotype of programmers being male.
    – Ed Morton
    Jun 22, 2019 at 15:48
  • 1
    Dear @EdMorton, thank you for clarifying the proper way to formulate a post. After consideration I incorrectly concluded that posting an input example to an external site was acceptable, avoiding adding to an already lengthy post. Now I know better and will make sure to include sample I/O in the post next time. Thank you
    – gusgrave
    Jun 23, 2019 at 18:44
  • 1
    Wrt gender - This is interesting indeed. After doing some research, the english language seems to be missing in the gender neutral/inclusive pronouns dept. so I do understand you position. In my native language there is a nationally adapted gender neutral/inclusive pronoun to use though I mostly just use usernames or "they/theirs". Thank you for a stimulating and developing exchange and your feedback, it is much appreciated.
    – gusgrave
    Jun 23, 2019 at 18:55

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .