Forum: Forums topas RSS
Sequential refinement results and gnuplot
rowlesmr #1
Member since Oct 2011 · 169 posts
Group memberships: Members
Show profile · Link to this post
Subject: Sequential refinement results and gnuplot
Hi all

I'd like to try and automate my plotting of refinement results, and due to John, gnuplot is the first thing to try.

Am I able to pivot data with gnuplot?

My data files look something like:

filename    seqno    phasename    a          b    c    scale     Rwp
blah_001.xye    1    corundum      3.00    3.00    12.0    0.001    3
blah_001.xye    1    silcon_NIST    5.40    5.40    5.4    0.002    3   
blah_002.xye    2    corundum      3.10    3.10    12.1    0.002    3.5
blah_002.xye    2    silcon_NIST    5.41    5.41    5.41    0.004    3.5   
blah_003.xye    3    corundum      3.20    3.20    12.2    0.001    3.1
blah_003.xye    3    silcon_NIST    5.42    5.42    5.42    0.002    3.1   
blah_004.xye    4    corundum      3.40    3.40    12.2    0.001    3.2
blah_004.xye    4    silcon_NIST    5.43    5.43    5.43    0.002    3.2   
blah_005.xye    5    corundum      3.50    3.50    12.3    0.001    3.9
blah_005.xye    5    silcon_NIST    5.44    5.44    5.44    0.002    3.9   
blah_006.xye    6    corundum      3.60    3.60    12.3    0.001    3.1
blah_006.xye    6    silcon_NIST    5.45    5.45    5.45    0.002    3.1

I would want to plot every data column (a,b,c,scale,Rwp) vs seqno for each phase (corundum & Si). The end result would be 10 plots -- Yes the Rwp plotted for corumdum and Si are redundent, but ease of gnuplot file wins out here.

I would like this to be as general as possible, so I wouldn't need to alter the gnuplot file everytime I make changes to my data file.

This is my first foray into gnuplot.

Thanks

Matthew
--
Matthew
johnsoevans (Administrator) #2
User title: John Evans
Member since Aug 2009 · 177 posts
Group memberships: Administrators, Members
Show profile · Link to this post
Matthew,

The script I give in Edinburgh is as below.  It assumes all the data is on one long line.  I'd therefore write your data out as just:

filename seqno a_cor b_cor c_cor scale_cor a_si b_si c_si scale_si Rwp

There is an "every" command that allows sample of points periodically.  You could do every second or third line with that format.  I don't tend to use it.

John

# Script to plot results file from multiple refinements
# Plot columns with no error first then columns with errors
# File has a dummy header line labelling the columns
# Temperature in col 1 Then data are in col 2 to 14; col 2 has no esd
# Change the column numbers to your example

reset

file = 'results.txt'
xcol     = 1       # x in column 1 is temperature
startcol = 2       # first column to plot is column 2
numnoesd = 1       # 1 column without an esd (Rwp in column 2)
numesd   = 27      # 27 columns with an esd

set style line 1 lt rgb "red"     linewidth 2 pt 7  

set key autotitle columnheader font "Arial,12" Left noenhanced
set term win font "Arial,24"
set ytics mirror
set xlabel 'Temperature/K'

#plot each column to the screen then to a gif file for columns with no esd
do for [col=startcol:startcol+(numnoesd-1):1] {
  plot file using xcol:col w linesp ls 1
  set term gif size 640,480 font "arial,12"; set output 'pic_col_'.col.'.gif'; replot; set term win
  pause -1 "<cr> to see the next plot"
}

#plot each column to screen then to a gif file for columns with associated esd
do for [col=startcol+numnoesd:startcol+numnoesd+2*(numesd-1):2] {
    plot file using xcol:col w linesp ls 1
  replot file using xcol:col:col+1 w errorbars ls 1 notitle
  set term gif ; set output 'pic_col_'.col.'.gif'; replot; set term win
  pause -1 "<cr> to see the next plot"
}
rowlesmr #3
Member since Oct 2011 · 169 posts
Group memberships: Members
Show profile · Link to this post
I've gone a cygwin solution.

The output from Topas goes through sort and awk and the transformed output then goes to gnuplot, with a little bit of grep thrown in. The current script is:
#!/bin/bash

file=$1
tmp="tmp.txt"
output="forGnuplot.txt"
header=$(head -n1 $file)

#sort the topas output file by the third column, and then sort
#  within each sorted section by the second column.
#  3rd = phase name, 2nd = seqno
# Unix stackexchange 52762
(head -n1 $file && tail -n +2 $file | sort -k 3,3 -k 2,2) > ${tmp}


#from that guy on the internet
#  get the name of the phase and prepend it to the data column header titles
awk -f script.awk ${tmp} > ${output}

rm ${tmp}

#get the number of columns in the data file
numcols=$(echo $header | awk '{print NF}')

#get the number of phases - count the number of empty lines and divide by two.
numphases=$(grep -cvP '\S' ${output})
numphases=$((($numphases/2)+1))

#if the 5th token ends in '_err', then skip
amIAnError=$(echo $header | awk '{print $5}')
if [[ "$amIAnError" == *\_err ]]; then
   skipError=2
else
   skipError=1
fi


rm figure.pdf

xcol=2

# numphases ; numcols ; every 1st/2nd col
gnuplot -c plot3.gp $numphases $numcols $skipError $xcol


where awk script is

BEGIN   { OFS = "\t" }

/^@/    {
    # save header fields into an array if the line starts with a "@"
    for (i = 1; i <= NF; ++i)
        header[i] = $i
    next
}

#If we're not on the first line, and the 3rd token of the current line
#  is not the same as the third token of the previous line
NR > 1 && $3!=p {
    # output two blank lines if needed
    if (print_blank) {
        print "\n"
    }
    print_blank = 1 #this gets set here to get the first line right

    # print first three headers as-is
    for (i = 1; i <= 3; ++i)
        printf("%s%s", header[i], OFS)

    # prepend column three to remaining headers
    for (i = 4; i < NF; ++i)
        printf("%s_%s%s", $3, header[i], OFS)
    printf("%s_%s%s", $3, header[NF], ORS)
}

# Otherwise, save value from column 3, and print everything
{ p=$3; print }

and the gnuplot script is

# ARG1 - number of phases
# ARG2 - total number of columns
# ARG3 - step interval in the column loop. 1= do every column, 2= skip the error column, if it's there.
# ARG4 - which column is my x-axis?

set print "-"  #set the print statment to go to STDOUT
set encoding iso_8859_1
set terminal postscript landscape noenhanced color 16
#set offset graph 0.0, 0.0, 0.05, 0.05
set output "figure.ps"
set key autotitle columnhead
set key outside
set key right above

do for [j=0:(ARG1-1)] {
  do for [i=4:ARG2:ARG3] {
    #the 3rd loop is only here to allow me to choose the column to plot as the x-ordinate
    #  It is the only way I could find to make the variable go into the plot command.
    do for [k=ARG4:ARG4] {
      # looking to see if the first datapoint is -1.
      #   if it is, then skip that column
      stats 'forGnuplot.txt' index j every 1:1:0:0:1:0 using i nooutput
      if ( STATS_min == -1 ) { continue; }

      plot 'forGnuplot.txt' index j using k:i with lp title columnhead
    }
  }
}

!ps2pdf figure.ps figure.pdf


I've also go a Windows batch file which will launch a cygwin process in the background and runs this script, so I don't need to by in cygwin for it to work.

@echo off

REM A constant to use later
SET cygwin=/cygdrive/

REM the folder where the batch file is being run from
SET winFolder=%cd%

REM replace the \ with / to go from Windows to Unix styles
SET unixFolder=%winFolder:\=/%

REM change the drive ":" to be a /
REM  and prepend the cygwin constant to get the
REM  cygwin directory
SET unixFolder=%unixFolder::=%
SET "unixFolder=%cygwin%%unixFolder%"

REM put the command together to be given to cygwin
REM  "replace" is the name of the script to run in cygwin
REM "%*" passess all the command line arguments
SET "command=cd %unixFolder%; ./replace %*"


REM actually run the script.
@echo on
c:\cygwin\bin\bash --login -c "%command%"


This turns data like:
@filename       seqno   phasename       a       b       c       scale   Rwp
blah_001.xye    1       corundum        3       3       12      0.001   3
blah_001.xye    1       silcon_NIST     5.4     5.4     5.4     0.002   3
blah_002.xye    2       corundum        3.1     3.1     12.1    0.002   3.5
blah_002.xye    2       silcon_NIST     5.41    5.41    5.41    0.004   3.5
blah_003.xye    3       corundum        3.2     3.2     12.2    0.001   3.1
blah_003.xye    3       silcon_NIST     5.42    5.42    5.42    0.002   3.1
blah_004.xye    4       corundum        3.4     3.4     12.2    0.001   3.2
blah_004.xye    4       silcon_NIST     5.43    5.43    5.43    0.002   3.2
blah_005.xye    5       corundum        3.5     3.5     12.3    0.001   3.9
blah_005.xye    5       silcon_NIST     5.44    5.44    5.44    0.002   3.9
blah_006.xye    6       corundum        3.6     3.6     12.3    0.001   3.1
blah_006.xye    6       silcon_NIST     5.45    5.45    5.45    0.002   3.1

into nice data for gnuplot:
@filename       seqno   phasename       corundum_a      corundum_b      corundum_c      corundum_scale  corundum_Rwp
blah_001.xye    1       corundum        3       3       12      0.001   3
blah_002.xye    2       corundum        3.1     3.1     12.1    0.002   3.5
blah_003.xye    3       corundum        3.2     3.2     12.2    0.001   3.1
blah_004.xye    4       corundum        3.4     3.4     12.2    0.001   3.2
blah_005.xye    5       corundum        3.5     3.5     12.3    0.001   3.9
blah_006.xye    6       corundum        3.6     3.6     12.3    0.001   3.1


@filename       seqno   phasename       silcon_NIST_a   silcon_NIST_b   silcon_NIST_c   silcon_NIST_scale       silcon_NIST_Rwp
blah_001.xye    1       silcon_NIST     5.4     5.4     5.4     0.002   3
blah_002.xye    2       silcon_NIST     5.41    5.41    5.41    0.004   3.5
blah_003.xye    3       silcon_NIST     5.42    5.42    5.42    0.002   3.1
blah_004.xye    4       silcon_NIST     5.43    5.43    5.43    0.002   3.2
blah_005.xye    5       silcon_NIST     5.44    5.44    5.44    0.002   3.9
blah_006.xye    6       silcon_NIST     5.45    5.45    5.45    0.002   3.1

Once I get my Sequential refinement macros working, I'll write it all up properly for the macros page.
--
Matthew
rowlesmr #4
Member since Oct 2011 · 169 posts
Group memberships: Members
Show profile · Link to this post
It's all working now, as far as I can tell.

http://topas.dur.ac.uk/topaswiki/doku.…?id=sequential_re…
--
Matthew
Close Smaller – Larger + Reply to this post:
Verification code: VeriCode Please enter the word from the image into the text field below. (Type the letters only, lower case is okay.)
Smileys: :-) ;-) :-D :-p :blush: :cool: :rolleyes: :huh: :-/ <_< :-( :'( :#: :scared: 8-( :nuts: :-O
Special characters:
Go to forum
Not logged in. · Lost password · Register
This board is powered by the Unclassified NewsBoard software, 20120620-dev, © 2003-2011 by Yves Goergen
Current time: 2018-10-22, 00:26:11 (UTC +00:00)