#!/bin/bash

ni=20; #number of images
nq=12; #number of quaternions
# list all log files of the zeroth replica
log="../output/job0.log.0 ../output/job1.log.0 ../output/job2.log.0 ../output/job3.log.0 ../output/job4.log.0"
#log="../sample-output/job4.log.0"

X=R; # replace with D or S as desired (R: reparametrized, D: drifted; S:smoothed)

# extracting the centers
cat $log | grep "TCL: Drifted " | awk -v ni=$ni '{for(i=4;i<=NF;i++)printf "%s ",$i;printf"\n";if($4==ni-1)print ""}' > cv-D.txt
cat $log | grep "TCL: Smoothed " | awk -v ni=$ni '{for(i=4;i<=NF;i++)printf "%s ",$i;printf"\n";if($4==ni-1)print ""}' > cv-S.txt
cat $log | grep "TCL: Reparametrized " | awk -v ni=$ni '{for(i=4;i<=NF;i++)printf "%s ",$i;printf"\n";if($4==ni-1)print ""}' > cv-R.txt

# calculating the string rmsd's
tail -$((ni+1)) cv-$X.txt > temp
cat temp cv-$X.txt | awk -v ni=$ni -v nq=$nq '
function acos(o){
    if (o>=1) {
	return 0
    } else if (o<=-1) {
	return 4*atan2(0.5,0.5)
    } else {
	return atan2(-o,sqrt(1-o*o))+2*atan2(0.5,0.5)
    }
}{
    if (NF==4*nq+1) {
	for (i=2;i<=NF;i++) c[$1,i-1]=$i
    }
    if ($1==ni-1) {
	if (!r) for (ij in c) c_last[ij]=c[ij]
	if (r==1) for (ij in c) c_first[ij]=c[ij]
	if (r>0) {
	    # with respect to first string
	    ds2=0; max=0
	    for (n=0;n<=ni-1;n++) {
		ds2i=0
		for (i=0;i<nq;i++) {
		    dot=0
		    for(l=1;l<=4;l++) dot+=c[n,4*i+l]*c_first[n,4*i+l]
		    acosdot=acos(dot)
		    ds2i+=acosdot*acosdot/nq
		}
		ds2+=ds2i;if (ds2i>max) max=ds2i
	    }
	    printf "%s %s %s",r,sqrt(ds2/ni),sqrt(max)
	    # with respect to last string
	    ds2=0; max=0
	    for (n=0;n<=ni-1;n++) {
		ds2i=0
		for (i=0;i<nq;i++) {
		    dot=0
		    for(l=1;l<=4;l++) dot+=c[n,4*i+l]*c_last[n,4*i+l]
		    ds2i+=acos(dot)*acos(dot)/nq
		}
		ds2+=ds2i;if (ds2i>max) max=ds2i
	    }
	    printf " %s %s",sqrt(ds2/ni),sqrt(max)
	    # with respect to previous string
	    if (r>1) {
    		ds2=0; max=0
		for (n=0;n<=ni-1;n++) {
		    ds2i=0
		    for (i=0;i<nq;i++) {
			dot=0
			for(l=1;l<=4;l++) dot+=c[n,4*i+l]*c_prev[n,4*i+l]
			ds2i+=acos(dot)*acos(dot)/nq
		    }
		    ds2+=ds2i;if (ds2i>max) max=ds2i
		}
		printf " %s %s",sqrt(ds2/ni),sqrt(max)
	    }
	    printf "\n"
	}
	for (ij in c) c_prev[ij]=c[ij]
	r++
    }
}' > s-rmsd-$X.txt

rm temp

if (( 1 ))
then
cat cv-R.txt | awk -v ni=$ni -v nq=$nq '
function acos(o){
    if (o>=1) {
	return 0
    } else if (o<=-1) {
	return 4*atan2(0.5,0.5)
    } else {
	return atan2(-o,sqrt(1-o*o))+2*atan2(0.5,0.5)
    }
}{
    for (i=2;i<=NF;i++) c[int((NR-1)/(ni+1)),$1,i-1]=$i
    if($1==0) n++
}END{
    for (j=0;j<n;j++) for (k=0;k<n;k++) {
	ds2=0
	for (m=0;m<ni;m++) for (i=0;i<nq;i++) {
	    dot=0
	    for (l=1;l<=4;l++) dot+=c[j,m,4*i+l]*c[k,m,4*i+l]
	    acosdot=acos(dot)
	    ds2+=acosdot*acosdot
	}
	print j,k,sqrt(ds2/(ni*nq))
    }
}' > pair-$X.txt
fi
