Skip to content
Snippets Groups Projects
Commit ec572915 authored by Germain Clavier's avatar Germain Clavier
Browse files

Added qp, GCcolors in plot and plot_table.py and lmp_ave_post.py to plot

and lammps directories.
parent e2151f06
Branches
No related tags found
No related merge requests found
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
'''
Takes the output file from some ave commands of LAMMPS
and returns average and sem.
'''
import argparse
import numpy as np
import sys
import logging
import utilsscript
def get_fields(infile, is_non_standard):
'''
Gets the fields of the file.
If non-standard, returns empty string.
'''
fields = []
if is_non_standard:
fields = ['']
else:
with open(infile, 'r') as f:
f.readline()
f.readline()
line = f.readline()
fields = line[1:].split()
logging.debug(fields)
return fields
def compute_mean(fields, infile):
'''
Computes data ave and sem.
'''
# Skip through header
nfields = len(fields)
with open(infile, 'r') as f:
line = f.readline()
while line[0] == '#':
line = f.readline()
# Initialisation
line = line.split()
step, nentries = int(line[0]), int(line[1])
nentries_ref = nentries
ave = np.zeros((nentries, nfields))
sem = np.zeros((nentries, nfields))
# Computing the averages
for i in range(nentries):
line = list(map(float, f.readline().split()))
for j, k in enumerate(fields):
ave[i, j] += line[j]
nrec = 1
while True:
try:
line = f.readline().split()
if not line:
break
step, nentries = int(line[0]), int(line[1])
if nrec % 1000 == 0:
logging.info(
"Reading record number {:<10d}, step {:<10.2f}."
.format(nrec, step)
)
if nentries != nentries_ref:
logging.error(
"Entries number changed between records. Check entry t={:10.2f}."
.format(step)
)
sys.exit(1)
for i in range(nentries):
line = list(map(float, f.readline().split()))
for j, k in enumerate(fields):
ave[i, j] += line[j]
nrec += 1
except EOFError:
break
logging.debug("Computes average over {:<10d} records.".format(nrec))
ave /= nrec
logging.info("Computing sem. Rereading file.")
# Computing sem
nrec = 0
f.seek(0)
line = f.readline()
while line[0] == '#':
line = f.readline()
while True:
try:
line = line.split()
step, nentries = int(line[0]), int(line[1])
if nrec % 1000 == 0:
logging.info(
"Reading record number {:<10d}, step {:<10.2f}."
.format(nrec, step)
)
if nentries != nentries_ref:
logging.error(
"Entries number changed between records. Check entry t={:10.2f}."
.format(step)
)
sys.exit(1)
for i in range(nentries):
line = list(map(float, f.readline().split()))
for j, k in enumerate(fields):
sem[i, j] += (line[j] - ave[i, j])**2
nrec += 1
line = f.readline()
if not line:
break
except EOFError:
break
sem = np.sqrt(sem/nrec)
return ave, sem
def write_output(outfile, fields, ave, sem):
'''
Self-explaining
'''
SPECIAL_CASES = [
'Chunk',
'OrigID',
'Coord1',
'Coord2',
'Coord3'
]
logging.info("About to write output file.")
with open(outfile, 'w') as f:
line = " ".join(('#', *["{:^17}"]*len(fields), '\n'))
f.write(line.format(*fields))
for i, (a, s) in enumerate(zip(ave, sem)):
# I found no nice way of doing this...
elems = []
for j, fi in enumerate(fields):
if fi in SPECIAL_CASES:
if fi == 'Chunk':
elems.append("{:^17d}".format(i+1))
else:
elems.append("{:>8.3f} {:>8.3f}".format(a[j], 0.))
else:
elems.append("{:>8.3f} {:>8.3f}".format(a[j], s[j]))
line = " ".join((*elems, '\n'))
f.write(line)
logging.info("DONE!")
return
def main():
parser = argparse.ArgumentParser(
description="Script to compute lifetime of tubes of symmetric star polymers."
)
parser.add_argument(
"-v",
"--verbose",
dest="verbosity",
default=0,
action="count",
help="display more information",
)
parser.add_argument(
"-f",
"--file",
dest="infile",
default="",
help="File to process",
)
parser.add_argument(
"-o",
"--output",
dest="outfile",
default="output.mean",
help="Output file name.",
)
parser.add_argument(
"-n",
"--non-standard",
dest="non_standard",
action="store_true",
help="Non standard LAMMPS headers.",
)
args = parser.parse_args()
##########
# Manage arguments
# -v/--verbose
utilsscript.init_logging(args.verbosity)
infile = args.infile
outfile = args.outfile
is_non_standard = args.non_standard
fields = get_fields(infile, is_non_standard)
ave, sem = compute_mean(fields, infile)
write_output(outfile, fields, ave, sem)
return
if __name__ == "__main__":
try:
main()
except KeyboardInterrupt:
raise SystemExit("User interruption.")
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
'''
Plot LAMMPS forces
'''
import argparse
import numpy as np
import os
import logging
import sys
import GCcolors
from matplotlib import pyplot as plt
def compute_energy(tp):
r = tp[0]
fo = tp[2]
e = np.zeros(r.shape)
for i, (ri, fi) in enumerate(zip(r, fo)):
if i == 0:
continue
dr = ri - r[i-1]
e[i] = e[i-1]-dr*fo[i-1]
e -= e[-1]
return e
def main():
parser = argparse.ArgumentParser(
description="Plots Lammps tabulated forces."
)
parser.add_argument(
"-v",
"--verbose",
dest="verbosity",
default=0,
action="count",
help="display more information",
)
parser.add_argument(
"-f",
"--file",
dest="infile",
default="",
help="File to read",
)
parser.add_argument(
"-x",
dest="xrange",
default="",
help="xrange : separated (for negative values: -x=-3:10)",
)
parser.add_argument(
"-y",
dest="yrange",
default="",
help="yrange : separated",
)
parser.add_argument(
"-e",
dest="extract",
action="store_true",
help="Extract the forces in separate files",
)
args = parser.parse_args()
##########
# Manage arguments
# -v/--verbose
try:
utilsscript.init_logging(args.verbosity)
except NameError:
logging.basicConfig(level=logging.WARNING)
logging.warning("utilsscript lib not found, using default logging at warning level.")
infile = args.infile
if not os.path.isfile(infile):
logging.error("Input file not found")
sys.exit(1)
toplot = []
with open(infile, 'r') as f:
lines = iter(f.readlines())
while True:
try:
r = []
force = []
ener = []
tok = []
while not tok:
tok = next(lines).partition('#')[0].rstrip()
logging.info("Found {} token".format(tok))
infos = next(lines).split()
npoints = int(infos[1])
dummy = next(lines)
for i in range(npoints):
line = next(lines).split()
r.append(float(line[1]))
ener.append(float(line[2]))
force.append(float(line[3]))
r = np.array(r)
ener = np.array(ener)
force = np.array(force)
toplot.append([r, ener, force, tok])
tok = []
dummy = next(lines)
except StopIteration:
break
for tp in toplot:
tp[1] = compute_energy(tp)
fig, axes = plt.subplots(1,2)
for tp in toplot:
axes[0].plot(tp[0], tp[1], label=tp[3], linewidth=3)
axes[1].plot(tp[0], tp[2], label=tp[3], linewidth=3)
if args.xrange:
xmin, xmax = list(map(float, args.xrange.split(":")))
axes[0].set_xlim(xmin, xmax)
axes[1].set_xlim(xmin, xmax)
if args.yrange:
ymin, ymax = list(map(float, args.yrange.split(":")))
axes[0].set_ylim(ymin, ymax)
axes[1].set_ylim(ymin, ymax)
# Setting axes 0
axes[0].set_title("Estimated energy", fontsize=30)
axes[0].legend(frameon=False, fontsize=30) # Fat font, no frame
axes[0].set_xlabel("r [A]", fontname="Hack", fontsize=30) # ylabel name, font is Hack, size 30pts
axes[0].set_ylabel("E [kcal/mol]", fontname="Hack", fontsize=30) # ylabel name, font is Hack, size 30pts
axes[0].tick_params(axis='both', which='major', labelsize=30) # Biggers ticks, bigger ticks label!
# Setting axes 1
axes[1].set_title("Tabulated force", fontsize=30)
axes[1].legend(frameon=False, fontsize=30) # Fat font, no frame
axes[1].set_xlabel("r [A]", fontname="Hack", fontsize=30) # ylabel name, font is Hack, size 30pts
axes[1].set_ylabel("F [kcal/mol/A]", fontname="Hack", fontsize=30) # ylabel name, font is Hack, size 30pts
axes[1].tick_params(axis='both', which='major', labelsize=30) # Biggers ticks, bigger ticks label!
plt.subplots_adjust(wspace=0.3)
plt.show()
if args.extract:
for tp in toplot:
outfile = "".join([tp[3], '.plot'])
logging.info("Writing file {}".format(outfile))
with open(outfile, 'w') as f:
f.write("# {} force extracted from {}\n".format(tp[3], infile))
f.write("# {:^20} {:^20} {:^20}\n".format('r', 'energy', 'force'))
for a, b, c in zip(tp[0], tp[1], tp[2]):
f.write("{:>18.16e} {:>18.16e} {:>18.16e}\n".format(a, b, c))
return
if __name__ == "__main__":
try:
main()
except KeyboardInterrupt:
raise SystemExit("User interruption.")
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
'''
Your docstring here
'''
import numpy as np
import re
from matplotlib import pyplot as plt
from matplotlib import cm
from matplotlib.colors import ListedColormap, LinearSegmentedColormap
# Found in a comment on Toward data science
# Tweaked my way
def colourGradient(fromRGB, toRGB, steps=50):
'''
colourGradient(fromRGB, toRGB, steps=50)
Returns a list of <steps> html-style colour codes forming a gradient between the two supplied colours
steps is an optional parameter with a default of 50
If fromRGB or toRGB is not a valid colour code (omitting the initial hash sign is permitted),
an exception is raised.
'''
# hexRgbRe = re.compile('#?[0–9A-F]{6}') # So we can check format of input html-style colour codes
# The code will handle upper and lower case hex characters, with or without a # at the front
# if not hexRgbRe.match(fromRGB) or not hexRgbRe.match(toRGB):
# raise Exception('Invalid parameter format')
# Tidy up the parameters
RGBA = np.ones((steps, 4))
rgbFrom = fromRGB.split('#')[-1]
rgbTo = toRGB.split('#')[-1]
# Extract the three RGB fields as integers from each (from and to) parameter
rFrom, gFrom, bFrom = [(int(rgbFrom[n:n+2], 16)) for n in range(0, len(rgbFrom), 2)]
rTo, gTo, bTo = [(int(rgbTo[n:n+2], 16)) for n in range(0, len(rgbTo), 2)]
# For each colour component, generate the intermediate steps
rSteps = ['{0:02x}'.format(round(rFrom + n * (rTo - rFrom) / (steps - 1))) for n in range(steps)]
gSteps = ['{0:02x}'.format(round(gFrom + n * (gTo - gFrom) / (steps - 1))) for n in range(steps)]
bSteps = ['{0:02x}'.format(round(bFrom + n * (bTo - bFrom) / (steps - 1))) for n in range(steps)]
# Values in 0, 1 format
Rcol = np.array([round(rFrom + n * (rTo - rFrom) / (steps - 1))/256 for n in range(steps)])
Gcol = np.array([round(gFrom + n * (gTo - gFrom) / (steps - 1))/256 for n in range(steps)])
Bcol = np.array([round(bFrom + n * (bTo - bFrom) / (steps - 1))/256 for n in range(steps)])
# Reassemble the components into a list of html-style #rrggbb codes
RGBA[:, 0] = Rcol
RGBA[:, 1] = Gcol
RGBA[:, 2] = Bcol
return RGBA, ["".join(['#',r,g,b]) for r, g, b in zip(rSteps, gSteps, bSteps)]
# My set of colors from colormind.io
GC_pink = "#FF01C0"
GC_dblue = "#122070"
# GC_blue = "#0013FF"
GC_blue = "#2CBDFE"
GC_green = "#558D3A"
GC_orange = "#FF8C00"
RGBA, GC_grad = colourGradient(GC_blue, GC_orange)
GCcmap = ListedColormap(RGBA)
color_list = [GC_pink, GC_dblue, GC_blue, GC_green, GC_orange]
plt.rcParams['axes.prop_cycle'] = plt.cycler(color=color_list)
plt.rc('axes', linewidth=2)
##### SOME EXAMPLES FOR TESTING PURPOSE ####
##### WITH SOME PERSONAL ADVICES
# 1) Better to set font and stuff for each axes
# 2) The bigger the better
##### Plotting linear data
# npoints = 10
# data = np.random.rand(npoints) # Random datas
# for i, _ in enumerate(color_list):
# plt.plot((i+1)*data, label=" ".join(['data', str(i+1)]))
# ax = plt.gca() # Get current axes to set properties
# ax.legend(frameon=False, fontsize=30) # Fat font, no frame
# ax.set_ylabel("Pouet", fontname="Hack", fontsize=30) # ylabel name, font is Hack, size 30pts
# ax.tick_params(axis='both', which='major', labelsize=30) # Biggers ticks, bigger ticks label!
##### Plotting heatmap
# data = np.random.rand(20, 20)
# plt.pcolormesh(data, cmap=GCcmap, rasterized=True) # Data is output from ListedColormap
# plt.colorbar(ax = plt.gca()) # Puts the colorbar on the side
##### Show result
# plt.show()
plot/qp 0 → 100755
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
'''
My module to plot a serie of files.
'''
import argparse as ap
import os
import sys
import re
import GCcolors
from itertools import cycle
import numpy as np
import scipy as sp
from scipy.signal import savgol_filter
from scipy.interpolate import interp1d
from matplotlib import pyplot as plt
# plt.rcParams.update({'lines.markeredgewidth': 4})
plt.rcParams.update({'font.size': 22, 'lines.linewidth': 3.})
def read_columns(fle, columns_list, fmt, skip):
'''
Reads input file and returns a list
containing columns' numbers followed
by array slices containing data.
'''
data = np.loadtxt(fle, skiprows=skip)
output = [fle]
for pair in columns_list:
partition = pair[0].split(':')
x = int(partition[0])
y = int(partition[1])
if fmt == 'xydy':
dy = int(partition[2])
elif fmt == 'xydx':
dx = int(partition[2])
elif fmt == 'xydxdy':
dx = int(partition[2])
dy = int(partition[3])
x_data = data[:, x-1]
y_data = data[:, y-1]
if fmt == 'xy':
output.append([x, y, x_data, y_data])
elif fmt == 'xydy':
dy_data = data[:, dy-1]
output.append([x, y, x_data, y_data, dy_data])
elif fmt == 'xydx':
dx_data = data[:, dx-1]
output.append([x, y, x_data, y_data, dx_data])
elif fmt == 'xydxdy':
dx_data = data[:, dx-1]
dy_data = data[:, dy-1]
output.append([x, y, x_data, y_data, dx_data, dy_data])
return output
def plot_sep_col(data, colors, fmt, is_log, legend):
'''
Makes separated ordered columns
wise: all 1:2 together, all 1:3 together
and so on...
'''
columns = []
for d in data:
for s in d[1:]:
columns.append([s[0], s[1]])
for i in range(len(columns)):
try:
x = columns[i][0]
y = columns[i][1]
j = i
while True:
try:
j += 1
pair = columns[j]
if pair[0] == x and pair[1] == y:
columns.pop(j)
except IndexError:
break
except IndexError:
break
plot_nb = len(columns)
if plot_nb % 3 == 0:
fig, axes = plt.subplots(plot_nb//3, 3, sharex=True, figsize=(20, 10))
elif plot_nb % 2 == 0:
fig, axes = plt.subplots(plot_nb//2, 2, sharex=True, figsize=(20, 10))
else:
fig, axes = plt.subplots(plot_nb, 1, sharex=True, figsize=(20, 10))
for pair, ax in zip(columns, axes.flatten()):
color_iter = cycle(colors)
suf = ":".join([str(pair[0]), str(pair[1])])
for d in data:
color = next(color_iter)
label = d[0]
for subdata in d[1:]:
x = subdata[0]
y = subdata[1]
if not(x == pair[0] and y == pair[1]):
continue
lab = label
x_data = subdata[2]
y_data = subdata[3]
if fmt == 'xy':
ax.plot(x_data, y_data,
label=lab, color=color
)
elif fmt == 'xydy':
dy = subdata[4]
plt.plot(x_data, y_data,
label=lab, color=color
)
up = y_data+dy
down = y_data-dy
plt.fill_between(
x_data, up, down,
alpha=0.2, edgecolor=color, facecolor=color,
linewidth=0
)
# ax.errorbar(x_data, y_data,
# yerr=dy,
# capsize=4, capthick=2,
# label=lab, color=color
# )
elif fmt == 'xydx':
dx = subdata[4]
ax.errorbar(x_data, y_data,
xerr=dx,
capsize=4, capthick=2,
label=lab, color=color
)
elif fmt == 'xydxdy':
dx = subdata[4]
dy = subdata[5]
ax.errorbar(x_data, y_data,
xerr=dx, yerr=dy,
capsize=4, capthick=2,
label=lab, color=color
)
if is_log:
ax.set_xscale('log')
ax.set_yscale('log')
color = next(color_iter)
if legend:
ax.legend(loc='upper left')
ax.set(title=suf)
return [fig, axes]
def plot_sep_fil(data, colors, fmt, is_log, legend):
plot_nb = len(data)
if plot_nb % 3 == 0:
fig, axes = plt.subplots(plot_nb//3, 3, sharex=True, figsize=(20, 10))
elif plot_nb % 2 == 0:
fig, axes = plt.subplots(plot_nb//2, 2, sharex=True, figsize=(20, 10))
else:
fig, axes = plt.subplots(plot_nb, 1, sharex=True, figsize=(20, 10))
for d, ax in zip(data, axes.flatten()):
color_iter = cycle(colors)
suf = d[0]
color = next(color_iter)
for subdata in d[1:]:
x = subdata[0]
y = subdata[1]
lab = ":".join([str(x), str(y)])
x_data = subdata[2]
y_data = subdata[3]
if fmt == 'xy':
ax.plot(x_data, y_data,
label=lab, color=color
)
elif fmt == 'xydy':
dy = subdata[4]
plt.plot(x_data, y_data,
label=lab, color=color
)
up = y_data+dy
down = y_data-dy
plt.fill_between(
x_data, up, down,
alpha=0.2, edgecolor=color, facecolor=color,
linewidth=0
)
# ax.errorbar(x_data, y_data,
# yerr=dy,
# capsize=4, capthick=2,
# label=lab, color=color
# )
elif fmt == 'xydx':
dx = subdata[4]
ax.errorbar(x_data, y_data,
xerr=dx,
capsize=4, capthick=2,
label=lab, color=color
)
elif fmt == 'xydxdy':
dx = subdata[4]
dy = subdata[5]
ax.errorbar(x_data, y_data,
xerr=dx, yerr=dy,
capsize=4, capthick=2,
label=lab, color=color
)
if is_log:
ax.set_xscale('log')
ax.set_yscale('log')
color = next(color_iter)
if legend:
ax.legend(loc='upper left')
ax.set(title=suf)
return [fig, axes]
def plot_data(data, colors, fmt, is_log, legend, names):
color_iter = cycle(colors)
plt.figure(figsize=(40, 30))
namelist=iter(names)
for d in data:
color = next(color_iter)
try:
label = next(namelist)
except StopIteration:
label = False
for subdata in d[1:]:
x = subdata[0]
y = subdata[1]
if not label:
label = d[0]
suf = ":".join([str(x), str(y)])
lab = label+suf
else:
lab = label
x_data = subdata[2]
y_data = subdata[3]
if fmt == 'xy':
plt.plot(x_data, y_data,
label=lab, color=color
)
elif fmt == 'xydy':
dy = subdata[4]
plt.plot(x_data, y_data,
label=lab, color=color
)
up = y_data+dy
down = y_data-dy
plt.fill_between(
x_data, up, down,
alpha=0.2, edgecolor=color, facecolor=color,
linewidth=0
)
# plt.errorbar(x_data, y_data,
# yerr=dy,
# capsize=4, capthick=2,
# label=lab, color=color
# )
elif fmt == 'xydx':
dx = subdata[4]
plt.errorbar(x_data, y_data,
xerr=dx,
capsize=4, capthick=2,
label=lab, color=color
)
elif fmt == 'xydxdy':
dx = subdata[4]
dy = subdata[5]
plt.errorbar(x_data, y_data,
xerr=dx, yerr=dy,
capsize=4, capthick=2,
label=lab, color=color
)
if is_log:
plt.xscale('log')
plt.yscale('log')
color = next(color_iter)
if legend:
# plt.legend(loc='upper left')
plt.legend()
return
def main():
'''
Main function.
'''
parser = ap.ArgumentParser(description='Script to plot a list of files.')
parser.add_argument('-v', '--verbose', dest='verbosity',
action='store_true', help='Verbosity.')
parser.add_argument('-sf', '--separated-files', dest='sf',
default=False, action='store_true',
help='Plot in separated windows based on filename. [default = False]')
parser.add_argument('-sc', '--separated-columns', dest='sc',
default=False, action='store_true',
help='Plot in separated windows based on columns. [default = False]')
parser.add_argument('-fmt', '--fmt', dest='fmt',
default='xy',
help='Format of the plot (xy, xydy, xydxdy, xydx). [default = xy]')
parser.add_argument('-i', '--input', dest='input',
nargs='+', action='append',
help='File name followed by columns and optionally lines to skip.')
parser.add_argument('-n', '--names', dest='names',
nargs='+', action='append',
help='Names for the legend.')
parser.add_argument('-l', '--loglog', dest='is_log',
action='store_true',
help='Plot in loglog scale.')
parser.add_argument('-t', '--title', dest='ttl', default='Your nice graph',
type=str, help='Title of the plot.')
parser.add_argument('-xl', '--xlabel', dest='xlab',
default=False, help='label for x')
parser.add_argument('-yl', '--ylabel', dest='ylab',
default=False, help='label for y')
parser.add_argument('-nl', '--no-legend', dest='no_legend', action='store_true',
help='removes the legend')
args = parser.parse_args()
verbose = args.verbosity
sep_files = args.sf
sep_col = args.sc
fmt = args.fmt
inpt = args.input
graph_title = args.ttl
xlab = args.xlab
ylab = args.ylab
is_log = args.is_log
names = args.names[0]
legend = not args.no_legend
if not inpt:
sys.exit('No file provided as input source.')
else:
for i in inpt:
if not os.path.isfile(i[0]):
sys.exit('Did not find file '+i[0]+'.')
if fmt == 'xy':
regex = re.compile(r'\d+:\d+')
elif fmt == 'xydy' or fmt == 'xydx':
regex = re.compile(r'\d+:\d+:\d+')
elif fmt == 'xydxdy':
regex = re.compile(r'\d+:\d+:\d+:\d+')
else:
sys.exit('Unknown fmt in -fmt option.')
regex_sk = re.compile(r'^\d+$')
colors = ['black', 'darkred', 'red',
'orangered', 'darkorange', 'orange',
'olive', 'darkolivegreen', 'chartreuse',
'forestgreen', 'lightseagreen', 'lightskyblue',
'steelblue', 'blue', 'darkblue', 'indigo']
colors = GCcolors.color_list
columns_list = []
data = []
for i in inpt:
fle = i[0]
columns_list = []
for column in i[1:]:
bingo = regex.findall(column)
if bingo:
columns_list.append(bingo)
temp = regex_sk.findall(i[-1])
if temp:
skip = int(temp[0])
else:
skip = 0
data.append(read_columns(fle, columns_list, fmt, skip))
if sep_col:
fig, axes = plot_sep_col(data, colors, fmt, is_log, legend)
elif sep_files:
fig, axes = plot_sep_fil(data, colors, fmt, is_log, legend)
else:
plot_data(data, colors, fmt, is_log, legend, names)
plt.xlabel(xlab)
plt.ylabel(ylab)
plt.title(graph_title)
plt.show()
return
if __name__ == "__main__":
try:
main()
except KeyboardInterrupt:
raise SystemExit("User interruption.")
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment