Skip to content

Instantly share code, notes, and snippets.

@jjhelmus
Created Apr 13, 2012
Embed
What would you like to do?
pysimpson
# modified example #1 which shows lost of precision in SIMPSON binary format.
spinsys {
channels 13C
nuclei 13C 13C
shift 1 0 6000 1 0 0 0
shift 2 0 6000 0 0 0 0
dipole 1 2 -1500 0 0 0
}
par {
method gcompute
spin_rate 2000
gamma_angles 20
sw spin_rate*gamma_angles
np 128
crystal_file rep320
start_operator Inx
detect_operator Inp
}
proc pulseq {} {
maxdt 5.0
delay 1e6
}
proc main {} {
global par
set iso 0
set f [fsimpson [list [list shift_2_iso $iso]]]
fsave $f reference_text.fid
fsave $f reference_rawbin.fid -raw_bin
fsave $f reference_binary.fid -binary
set f [fload reference_binary.fid]
fsave $f postread_text.fid
fsave $f postread_rawbin.fid -raw_bin
fsave $f postread_binary.fid -binary
# reference_text.fid and postread_text.fid will not match, but should
# reference_rawbin.fid and postread_rawbin.fid will not match, but should
}
"""
pysimpson: Python module for reading SIMPSON files
By: Jonathan J. Helmus (jjhelmus@gmail.com)
Version: 0.1 (2012-04-13)
License: GPL
"""
import math
import numpy as np
def read_raw_bin(filename):
"""
Read a SIMPSON file saved using: fsave $f filename -raw_bin
Returns a 1D complex array
"""
data = np.fromfile(filename, dtype='float32')
half = int(data.shape[0] / 2)
cdata = np.empty( (half, ), dtype='complex64')
cdata.real = data[:half]
cdata.imag = data[half:]
return cdata
def read_binary(filename):
"""
Reads a SIMPSON file saved using: fsave $f filename -binary
Note that the data read from a SIMPSON binary file will typically NOT match
the data in a raw_bin file as the character encoding process used by
SIMPSON causes a loss of one bit of accuracy.
Returns
-------
dic : dictionary
Dictionary of parameters contained in header
data : array
Array of data in file, complex and shaped appropiately.
"""
f = open(filename, 'r') # open the file
in_data_block = False # unset data block flag
chardata = "" # initalize data block string
dic = {} # initalize dictionary of parameters
# parse the file, pulling out the data characters and parameters
for line in f:
line = line.strip('\n')
if in_data_block:
if line == "END":
in_data_block = False
else:
chardata += line
elif line == "SIMP":
continue
elif line == "DATA":
in_data_block = True
elif "=" in line:
key, value = line.split("=")
dic[key] = value
else:
print "Warning, skipping line:",line
# DEBUGGING
#return dic, chardata
# convert every 4 character to 3 'bytes'
nquads, mod = divmod(len(chardata), 4)
assert mod == 0 # character should be in blocks of 4
bytes = []
for i in xrange(nquads):
bytes += chars2bytes(chardata[i * 4:(i + 1) * 4])
# DEBUGGING
#return dic, bytes
# convert every 4 'bytes' to a float
num_points, num_pad = divmod(len(bytes), 4)
data = np.empty( (num_points, ), dtype='float32')
for i in xrange(num_points):
data[i] = bytes2float(bytes[i * 4 : (i + 1) * 4])
if 'NI' in dic: # 2D data, reshape to NI, NP
return dic, data.view('complex64').reshape(int(dic['NI']), -1)
return dic, data.view('complex64')
BASE = 33
FIRST = lambda f,x: ((x) & ~(~0 << f))
LAST = lambda f,x: ((x) & (~0 << (8-f)))
def chars2bytes(chars):
""" Convert four characters from a data block into 3 'bytes' """
c0, c1, c2, c3 = [ord(c) - BASE for c in chars]
return [FIRST(6, c0) | LAST(2, c1 << 2),
FIRST(4, c1) | LAST(4, c2 << 2),
FIRST(2, c2) | LAST(6, c3 << 2)]
def bytes2float(bytes):
""" Convert four bytes to a string """
# the first 23 bits (0-22) define the mantissa
# the next 8 bits (23-30) the exponent,
# the last bit (31) defines the sign
b0, b1, b2, b3 = bytes
mantissa = ((b2 % 128) << 16) + (b1 << 8) + b0
exponent = (b3 % 128) * 2 + (b2 >= 128) * 1
negative = b3 >= 128
e = exponent - 0x7f
m = np.abs(mantissa) / np.float64(1 << 23)
if negative:
return -math.ldexp(m, e)
return math.ldexp(m, e)
def bytes2float_bitarray(bits):
""" byte2float function using BitArray for explanation """
from bitstring import BitArray
float_bytes = sum([BitArray(uint=b, length=8)[::-1] for b in bits])
mantissa = float_bytes[:23][::-1].uint # 23 bits
exponent = float_bytes[23:31][::-1].uint # 8 bits
negative = float_bytes[31] # 1 bit
e = exponent - 0x7f
m = np.abs(mantissa) / np.float64(1 << 23)
if negative:
return -math.ldexp(m, e)
return math.ldexp(m, e)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment