Scripts:get stats.py

From UC_Chemistry

---

Revision as of 16:32, 1 September 2009

Note: Make sure kT and h are correct for your particular system!

Script Source

#!/usr/bin/python

import sys
from ucgrad import *

kT = Boltzmann*298.15/4.184
h = 0.01

def main(argv):
        if len(argv) < 6:
                print "Usage: %s <p0.expot> <x0.expot> <neg dat> <pos dat> <out>"%(argv[0])
                return 1

        # lambda, beta mu, beta err
        p0 = read_matrix(argv[1])
        x0 = read_matrix(argv[2])

        # lambda, S, mu_DU, var_DU, bound, err, nb
        du0 = read_matrix(argv[3])
        du1 = read_matrix(argv[4])

        p0, x0, du0, du1 = match_lines(p0, x0, du0, du1)
        if len(p0) < 1:
                print "Unable to find any matching radii!"
                return -1

        # lambda, Upper bound+err, Upper Est.+err, Lower Est.+err, Lower Bound+err., LR Avg. Est.+err, HS+err, IS+err, total+err, err2
        M = zeros((len(p0),18))
        M[:,0] = p0[:,0]
        M[:,1:3] = du0[:,4:6]
        M[:,3] = du0[:,2]-0.5*du0[:,3]/kT
        M[:,4] = sqrt( (du0[:,3]+0.5*(du0[:,3]/kT)**2)/(du0[:,1]-1.0) )
        M[:,5] = du1[:,2]+0.5*du1[:,3]/kT
        M[:,6] = sqrt( (du1[:,3]+0.5*(du1[:,3]/kT)**2)/(du1[:,1]-1.0) )
        M[:,7:9] = du1[:,4:6]
        M[:,9] = 0.5*(du0[:,2]+du1[:,2])
        M[:,10] = 0.5*sqrt(du0[:,3]/(du0[:,1]-1.0)+du1[:,3]/(du1[:,1]-1.0))
        M[:,11:13] = p0[:,1:3]*kT
        M[:,11] *= -1.0 # should be pos., not neg.
        M[:,13:15] = x0[:,1:3]*kT
        M[:,15] = M[:,9]+M[:,11]+M[:,13]
        M[:,16] = sqrt(M[:,10]**2+M[:,12]**2+M[:,14]**2)
        M[:,17] = 0.5*abs(M[:,5]-M[:,3]) # Approx. to Gaussian approx. error
                        # being the difference between the 2-moment est.s.
        write_matrix(argv[5], M)

        return 0

def match_lines(*tbl):
        # round to int rep. using pre-set h
        rlist = []
        for t in tbl:
                rlist.append((t[:,0]/h+0.5).astype(int).tolist())
        inrad = set(rlist[0])
        for r in rlist[1:]:
                inrad = inrad & set(r)

        clip = []
        for rl, t in zip(rlist, tbl):
                t = t.tolist()
                for i in range(len(t)-1, -1, -1): # Check names with "in" list
                        if rl[i] not in inrad: # and be a good bouncer.
                                del t[i]
                sort(t, 0) # Numerically sort. 'em.
                clip.append(array(t))

        return tuple(clip)

def pfee(x):
        m, v, me, ve = stats(exp(x/kT))
        return kT*log(m), kT*me/m

def nfee(x):
        m, v, me, ve = stats(exp(-x/kT))
        return -kT*log(m), kT*me/m

if __name__=="__main__":
        main(sys.argv)