pysimpson

example.in

# 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.py

"""
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)