Bonsoir tout le monde,
Intéressé par un récent billet de blog sur B. Kernighan , j'essaie de scripter avec AWK une extraction de coordonnées cartésiennes de certains atomes dans un output de dynamique moléculaire.
Ce dernier est composée d'une ligne de titre, d'un header général, et d'un certain nombre de configurations (snapshots) du systeme pris a intervalles de temps réguliers (dénotés par un timestep). Plus précisément, chaque snapshot comporte 4 lignes de header (1 + 3 lignes pour les vecteurs de la cellule spatiale unitaire), et pour chaque atome : 1 ligne de header (élement, index, …) + 1, 2 ou 3 lignes comportant les positions/vitesse/forces.
Mon objectif est, pour l'instant de reconstruire un fichier de ce format, en ne gardant qu'une certaine molécule (dont je donne l'index du 1er atome, et le nombre d'atomes).
EDIT: finalement, le cœur de mon problème (voir 2e code) réside dans le fait que j'ai encore du mal a saisir comment et quand les conditions des actions sont vérifiées : je cherche a les "remettre sur le tapis" en changeant la valeur d'une variable dans la condition, en quelque sorte, émuler une (deux) boucle.
J'ai commencé naïvement par des boucles et des getline
, le plus embetant ayant été de faire passer des arguments (par le biai d'un surcouche bash):
extractHIS.sh:
#!/bin/bash
/home/francois/Sources/SCRIPTS/extractor.awk -v FirstAtomIndex=$1 NbAtomInMol=$2 HISTORY > redHISTORY
extractor.awk:
#! /bin/awk -f
#Compatible with DL_POLY4 HISTORY files
BEGIN{}
{
print "Extraction of ", NbAtomInMol, "atoms, starting from index #", FirstAtomIndex;
#TODO check for external parameters assignation
#Read the file specs and recreate it
getline;
CoordsWidth=$1+1;
NbAtomsInCell=$3;
NbStep=$4;
LastRecord=$5;
NewLastRecord=6+NbStep*(CoordsWidth*NewNbAtomsInCell+4);
print $1, $2, NbAtomInMol, $4, NewLastRecord;
#loop over every timesteps
for ( i=1;i<=NbStep;i=i++ ) {
getline;
print $1, $2, NbAtomInMol, $4, $5, $6, $7;
for ( j=0;j<3;j++ ) {
#cell specs
#loop over the number of cell vectors
getline;
print $0;
}
for ( j=1;j<=NbAtomsInCell;j++ ) {
getline;
if( j>=FirstAtomIndex && j<=(FirstAtomIndex+NbAtomInMol-1+CoordsWidth) ) {
print $1,(j-FirstAtomIndex+1),$3,$4,$5;
for ( k=0;k<CoordsWidth;k++ ) { getline; print $0; }
}
else {
for ( k=0;k<CoordsWidth;k++ ) { getline; }
}
}
}
}
END{}
2 Problemes:
* j'ai besoin de pouvoir echantillonner les timesteps, la prise en charge d'un nouvel argument d'intervalle est ajouté dans le script bash, mais les boucles utilisées sont mal adaptées
* les perfs sont sympas mais pas top. je mets ca sur le compte du fait que le fichier est parcouru avec des getlines
Puisque les numéros des lignes intéressantes peuvent être déduits du nombre d'atomes/de timesteps/… exploiter les NR
me semble etre la meilleure solution.
J'ai donc réecrit le script comme suit :
#! /bin/awk -f
#Compatible with DL_POLY4 HISTORY files
BEGIN{
print "Extraction of ", NbAtomInMol, "atoms, starting from index #", FirstAtomIndex; i1=0; i2=1;
}
NR == 2 {
CoordsWidth=$1+1;
NbAtomsInCell=$3;
NbStep=$4;
LastRecord=$5;
NewLastRecord=2+NbStep*(CoordsWidth*NbAtomsInMol+4);
print $1, $2, NbAtomInMol, $4, NewLastRecord;
getline;
}
#TODO check for external parameters assignation
#Read the file specs and recreate it
#timesteps
#{
#for ( i=1;i<=NbStep;i=i+Intrvl ) {
# if ( NR==( 2 + (i-1) * (4+NbAtomsInCell*(CoordsWidth+1)) +1 ) ) { print $1, $2, NbAtomInMol, $4, $5, $6, $7; }
# if ( NR==( 2 + (i-1) * (4+NbAtomsInCell*(CoordsWidth+1)) + (4 + (FirstAtomIndex-1)*(CoordsWidth+1) ) +1) ) {
# for (j=1;j<=NbAtomInMol;j++) {
# print $1,j,$3,$4,$5;
# for ( k=0;k<=CoordsWidth-1;k++ ) { getline; print $0; }
# getline;
# }
# }
#}
#
#}
NR == ( 2 + (i2-1) * (4+NbAtomsInCell*(CoordsWidth+1)) +1 ) { print $1, $2, NbAtomInMol, $4, $5, $6, $7; i1++; }
NR == ( 2 + (i1-1) * (4+NbAtomsInCell*(CoordsWidth+1)) + (4 + (FirstAtomIndex-1)*(CoordsWidth+1) ) +1 ) {
for (j=1;j<=NbAtomInMol;j++) {
print $1,j,$3,$4,$5;
for ( k=0;k<=CoordsWidth-1;k++ ) { getline; print $0; }
getline;
}
i2++;
print $1, $2, NbAtomInMol, $4, $5, $6, $7;
}
END{}
La partie commentée est fonctionnelle, mais très (très) lente! J'ai donc essayé de me servir des actions/patterns: cela s'avere beaucoup plus rapide, mais je ne parvient pas a transcrire a la fois les coordonnées des atomes, et les headers de chaque timesteps. Ici mon dernier essai d'incrémenter les indices i1 et i2 de façon croisée, pour essayer de forcer AWK a reconsidérer les pattern/actions mis a jour avec les nouvelles valeurs de variable (un seul index i partagé, ou deux croisés) mais ça a l'air d’être un peu hors de la philosophie AWK. Est-ce que quelqu'un aurait un conseil pour me remettre sur les rails?
Merci.
EDIT:
Bon, j'ai combiné un petit truc a base de modulos pour les boucles implicites, qui fait tout bien comme il faut.
Le script bash pour l'appel:
#!/bin/bash
./HISsampext.awk -v FirstAtomIndex=$1 -v NbAtomsInMol=$2 -v Intrvl=$3 HISTORY > redHISTORY
#TODO check that div done in header for Nb final records/Nb final steps corresponds to the modulo trick result in sampler
Le script AWK:
#! /bin/awk -f
#Compatible with DL_POLY4 HISTORY files
BEGIN{
print "Extraction of ", NbAtomsInMol, "atoms, starting from index #", FirstAtomIndex, "with interval timestep of ", Intrvl;
RS="\n";
a=1;
}
NR == 2 {
CoordsWidth=$1+1;
NbAtomsInCell=$3;
NbStep=$4;
LastRecord=$5;
NewLastRecord=2+NbStep*(CoordsWidth*NbAtomsInMol+4);
print $1, $2, NbAtomsInMol, $4, NewLastRecord;
a=(4 + NbAtomsInCell * (CoordsWidth+1) )*Intrvl;
b=1+(4 + (FirstAtomIndex-1) * (CoordsWidth+1) );
#c=NbAtomsInMol * (CoordsWidth+1);
#print a,b;
}
((NR-2) % a == 1) {
print $1, $2, NbAtomsInMol, $4, $5, $6, $7;
for ( k=1;k<=3;k++ ) { getline; print $0; }
}
((NR-2) % a == b) {
if(NR>2) {
for (j=1;j<=NbAtomsInMol;j++) {
print $1,j,$3,$4,$5;
for ( k=0;k<=CoordsWidth-1;k++ ) { getline; print $0; }
getline;
}
}
}
END{}
Ca marche pas trop mal, et je suis assez content de ma gueule d'avoir pensé (et réussi a faire marcher) le coup des modulo ^
Par contre, ca pourrait marcher un peu plus vite… il y a eu quelques idées/conseils dans les commentaires que je vais essayer de creuser.
# Python scientifique
Posté par Nonolapéro . Évalué à 3.
J'ai pas tout saisi à ton problème mais j'ai l'impression qu'en utilisant python et numpy c'est beaucoup plus simple.
Pour récupérer des colonnes depuis un fichier texte c'est :
Après si tes données sont un peu plus farfelues ou volumineuses ça doit valoir le coup de regarder du côté de pandas.
[^] # Re: Python scientifique
Posté par albahtaar . Évalué à 0.
Intéressant, bien que j'aimerais persévérer encore un peu avec AWK. Je ne connais pas encore python, mais il faudra bien que je m'y mette un jour ou l'autre.
Sinon j'avais pensé a utiliser R ou une base de données, au pire je peut toujours bricoler un truc en Fortran, mais je me méfie : j'ai un doute sur le comportement des OPEN (est-ce que tout le fichier est chargé en mémoire a ce moment la?).
[^] # Re: Python scientifique
Posté par Nonolapéro . Évalué à 2.
Il existe des solutions pour lire le fichier par morceaux et du côté de pandas ça semble ne plus poser de problème depuis la version 0.10. La dernière version est la 0.12.
http://stackoverflow.com/questions/11622652/large-persistent-dataframe-in-pandas/
http://wesmckinney.com/blog/?p=543
# Perl
Posté par Krunch (site web personnel) . Évalué à 3.
J'ai pas tout suivi non plus mais.
Bien souvent, quand awk est trop lent, Perl permet de résoudre le problème. Si tu veux juste tester en vitesse, il y a même a2p qui te permet de convertir un script awk en Perl. Après, si l'idée c'est d'apprendre awk ça ne t'avance pas beaucoup.
pertinent adj. Approprié : qui se rapporte exactement à ce dont il est question.
[^] # Re: Perl
Posté par albahtaar . Évalué à 0.
Je ne suis pas fan de Perl, j'ai un peu pratiqué mais ça ne me convient pas. Merci quand même.
# Awk
Posté par Sygne (site web personnel) . Évalué à 3.
Je ne peux pas t'aider beaucoup faute de pratiquer intensivement awk, et faute d'avoir sous les yeux les données à traiter, mais peut-être que les quelques remarques suivantes te seront utiles.
[^] # Re: Awk
Posté par albahtaar . Évalué à 1.
Merci pour ces infos, ça me donne en effet pas mal d'idée pour changer d'approche. Le top serait peut être alors de considérer chaque timestep comme un record, en utilisant un mot qui revient dans chacun ("timestep ) comme RS, et ensuite de gérer chaque atome comme un field…
Mais est-ce qu'un record de ~1000 lignes (comportant une 30aine de caracteres) serait viable? Mon gros probleme est finalement que je manque de perspective sur le fonctionnement interne de AWK.
Par exemple, dans la partie commentée du 2e code, je me demande si, a chaque boucle, AWK ne repart pas du début du fichier pour retrouver les nouvelles lignes demandées.
[^] # Re: Awk
Posté par Sygne (site web personnel) . Évalué à 2.
J'ai tendance à croire que Awk verrait les choses comme ça, ou peut-être autrement, mais pas à la façon d'un programme linéaire, si ce mot peut convenir. Après, je ne sais pas comment il va se comporter avec un enregistrement de 1000 lignes.
Je ne peux pas beaucoup t'aider. À ta place, je lirais quelques petits programmes awk pour commencer à penser en awk. Je crois que c'est un langage assez singulier, et c'est probablement pourquoi on lui préfère python et perl. Mais à mon avis, vu que chaque langage apprend à l'esprit de nouvelles façons de se structurer, il y a un intérêt à apprendre awk, justement parce qu'il est singulier.
Je n'ai pas pris le temps de le faire vraiment pour ma part, et c'est pourquoi je ne peux pas t'aider plus.
[^] # Re: Awk
Posté par albahtaar . Évalué à 1.
Pas de soucis, c'est déja bien urbain de ta part. Merci!
Suivre le flux des commentaires
Note : les commentaires appartiennent à celles et ceux qui les ont postés. Nous n’en sommes pas responsables.